International Journal of GeophysicsVolume 2012 (2012), Article ID 610712, 10 pagesdoi:10.1155/2012/610712

Research Article

Spatiotemporal Relationship between Geodetic and Seismic Quantities: A Possible Clue to Preparatory Processes of M ≥ 6.0 Inland Earthquakes in Japan

Masashi Kawamura,¹ Takeshi Kudo,² and Koshun Yamaoka³

¹Department of Earth Sciences and Graduate Institute of Geophysics, National Central University, Jhongda Road no. 300, Taoyuan, Jhongli 32001, Taiwan²General Education Division, College of Engineering, Chubu University, Matsumoto-cho 1200, Aichi, Kasugai 487-8501, Japan³Research Center for Seismology, Volcanology and Disaster Mitigation, Graduate School of Environmental Studies, Nagoya University, Chikusa-ku, Nagoya 464-8601, Japan

Received 23 June 2011; Revised 6 September 2011; Accepted 11 September 2011

Academic Editor: Rodol̀fo Console

Copyright © 2012 Masashi Kawamura et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Constructing a statistical model to characterize the physical conditions associated with large earthquake occurrence is crucial for disaster mitigation. With the aim to formulate such a model, we previously developed a statistical evaluation system which assesses the correlation of the spatiotemporal relationship between different kinds of physical quantities with the occurrence times of large inland earthquakes. In this study, we focused on assessing the relationship between geodetic and seismic quantities and attempted to find the pair of related quantities that most likely indicates preparatory processes of large earthquakes in Japan. We assessed the quantities prior to M ≥ 6.0 inland earthquakes for the period of 2001–2007 in terms of probability gains and error diagram. Our system revealed that the pair of absolute value of dilatation rate and seismic energy showed the highest statistical performance. Further validation of this result is required by updating the database of physical quantities.

1. Introduction

It is an urgent issue for disaster mitigation to develop a statistical model well reflecting crustal activities, which are expected to include the preparatory processes of large inland earthquakes. Resolution of this issue requires a comprehensive understanding and monitoring of crustal activities through more than one kind of physical observation. Therefore, keeping this in mind, in this study, we focused on examining the spatiotemporal relationships between different physical quantities, at least one of which is time-variable. By determining informative combinations of physical quantities, that is, those having some correlation, we aim to develop a monitoring index/indices in crustal activities which feature the preparatory processes of large inland earthquakes for Japan. Imoto [1, 2] insisted that their seismicity model shows higher performance by taking into consideration the correlation between different kinds of geophysical parameters. This implies that there are some correlations between different kinds of observations that may represent different aspects of crustal phenomena or show causality between apparently different crustal phenomena.

As the first step toward formulating a model to capture statistical and physical conditions for the occurrence of large inland earthquakes, we created a database of physical observations with spatially and temporally gridded formats (Figure 1). Furthermore, we constructed a statistical evaluation system to assess the relationship between the occurrence times of target events and interquantity relationship in the physical observation data, at least one of which is time variable [4]. This system comprises some calculation modules corresponding to respective steps in Section 2 and their subordinate functions (not shown in Section 2) to output several parameters which are required for next steps or to display the final result of statistical evaluation at (Table 1, Figure 2).

Table 1: Tabulated information on input parameters required at each step of statistical evaluation system. The related schematic explanation should be referred to Figure 2.

No.
Parameter

(1)
Range of region for analysis (inland Japan)

(2)
Interval between spatial grids (0.05°)

(3)
Interval between time grids (three month)

(4)
Time period referred to for calculating physical quantities such as dilatation rate for each time grid (two years)

(5)
Range of depth (0 ≤ depth ≤ 30 km)

(6)
Range of magnitude (M ≥ 1.0)

(7)
Smoothing (weight-adjusting) parameters for obtaining spatially gridded physical quantities (refer to Figure 2(a))

(8)
Search radius for applying smoothing parameters to area centered at a calculation grid (Figure 2(a))

(9)
Search radius for calculating statistical index for a calculation grid (60 km)

(10)
Approval method for classifying the temporal change in statistical index R(t) into three categories (Figure 2(b))

(11)
Confidence level for statistical approval (95%)

(12)
Calculation period of total seismic energy released following calculation period of temporal change in statistical index R(t) (Figure 2(b))

(13)
Lower limit of total seismic energy released during period (12) which is regarded as seismically active (Figure 2(b))

(14)
Definition of large inland earthquakes (lower limit of magnitudes) (M ≥ 6.0)

Figure 1: Gridded maps for (a) dilatation rate, (b) maximum shear strain rate, (c) seismic energy [3] within 30 km search radius, and (d) the number of earthquakes within 30 km search radius, respectively, over the Japanese islands.

(a)

(b)

Figure 2: (a) Schematic explanation on the method of obtaining geodetic quantities (dilatation rate and maximum shear strain rate) and seismic quantities (seismic energy and number of earthquakes) for each grid, which were calculated by relating them to the displacement rates for GPS stations (open triangle) and the magnitudes for hypocenters (open circle) in circular area, respectively. We weighted the observation errors for GPS stations and seismic quantities for hypocenters depending on the distances between the calculation grid (black square) and the other grids in circular area, respectively. The vertical axis of each right figure in which weight curve is delineated shows the distances relative to the calculation grid. (b) Schematic explanation on the approach of obtaining Figures 3 and 4. See also text for detail.

Here, we actually operated the system to examine the relationship between two physical quantities and the occurrence times of mainshocks in Japan. We adopted the spatiotemporal relationship between geodetic quantities (dilatation rate and maximum shear strain rate) and seismic quantities (seismic energy and number of earthquakes) for representative physical observation. We look for the relationship for the earthquake from 2001 through 2007 for inland Japan. In evaluating this relationship, the system computes the statistical performance of a proposed monitoring index which would reflect an aspect of the preparatory process of large earthquake occurrence. We searched for the most informative combination between geodetic and seismic quantities based on evaluating the statistical performance of the proposed monitoring index for inland mainshocks. Evaluation of the statistical performance revealed that the time variation of the monitoring index created based on the relationship between the absolute value of dilatation rate and seismic energy was most closely associated with the occurrence times of inland mainshocks. Our result revealed the relationship between the absolute value of dilatation rate and seismic energy would best reflect the crustal activities involving preparatory processes of the occurrence of M ≥ 6.0 inland mainshocks. An important future issue is to further validate the highest performance of the proposed monitoring index by updating the database to increase the number of samples, or large inland earthquakes.

2. Statistical Evaluation of the Relationship between Geodetic and Seismic Quantities

The statistical evaluation system developed by Kawamura et al. [4] statistically assesses the relationship between the spatiotemporal connection in any two kinds of physical quantities and the occurrence times of target physical events of interest We actually operated the system as follows to find a monitoring index, which is created based on a pair of items (a) to (d) in Figure 1, or a pair of geodetic and seismic quantities, such that the index can feature crustal activities which are expected to reflect the preparatory processes of large inland earthquakes.

2.1. Step 1: Making a Database with Spatially and Temporally Gridded Formats

We created a database comprising items (a) to (d) in Figure 1. It shows spatially and temporally gridded maps of four physical quantities: (a) dilatation rate, (b) maximum shear strain rate, (c) seismic energy [3] within 30 km search radius, and (d) the number of earthquakes within 30 km search radius, respectively, over inland Japan. Items (a) and (b) were computed from GEONET (GPS Earth Observation NETwork) data (Figure 2(a)). Items (c) and (d) were calculated from unified hypocenter data operated by Japan Meteorological Agency (Figure 2(a)).

We here explain the method of making physical quantities (a) to (d) (Figure 2(a)). To obtain dilatation rate (a) and maximum shear strain rate (b) with grid formats, we first removed abnormal values from the GEONET coordinate data, where the coordinate data with a standard deviation >3 over half a month before and after each time was regarded as abnormal value. We next calculated the displacement rate at each observation point. We assumed that the displacement rate on an arbitrary coordinate grid can be represented by using a linear function of the displacement rate at its nearby observation point. Then, the displacement rate at observation point is related to the position of grid point as follows: where and are east-west and north-south components of displacement rate at each observation point; and are east-west and north-south components of the position of observation point relative to the position of grid point. and are the error term of respective components. Then, dilatation rate and maximum shear strain rate ∑ can be obtained as follows: Furthermore, the error term , was weighted depending on the distance between observation and grid points as follows: where is each component of the estimation error of the displacement rate at observation point. R is the distance between observation and grid points. D is weight-adjusting parameter. In the following sections, we show the system operation result obtained by setting and .

To obtain seismic energy (c) and the number of earthquakes (d) for each grid point, the earthquakes within a radius of 60 km of grid point with depths shallower than 30 km and magnitudes ≥1.0 were collected (Figure 3(b)). Then, the total seismic energy and the total number of earthquakes were calculated, where the seismic energy and the number of each earthquake were weighted depending on the distance between observation and grid points as follows: where is physical quantity related to the th earthquake, to which 1 is assigned in the case physical quantity is defined as the number of earthquakes, n is the total collected number of earthquakes for grid point, and R is the distance between observation and grid points. is weight-adjusting parameter. In this way, by taking the distance-depending weight into consideration, the essential event collection radius becomes 30 km. In the following sections, we show the system operation result obtained by setting = 10 km and 0 ≤ R ≤ 60 km.

(a)

(b)

(c)

(d)

Figure 3: Spatial distributions of linear trends of statistical index R(t) (pure red: significant decrease, pure blue: significant increase, pure green: insignificant change), which is obtained by comparison between (a) the absolute value of dilatation rate and seismic energy per two years, (b) maximum shear strain rate and seismic energy per two years, (c) absolute value of dilatation rate and number of earthquakes per two years, (d) maximum shear strain rate and number of earthquakes per two years. Time period for calculating linear trends of R(t) is from 1 April 1998 through 30 September 2004, which is shown in panel (c) and (d). Superimposed in dark colors are the same distributions of areas in which seismic energy more than or equal to M4.0 earthquake (6 × 10¹⁰ Nm) was radiated during one year after the time period shown in panel (c) and (d). Linear trends of R(t), which were obtained by least squares fitting, were classified using T approval with a confidence level of 95%. The numbers of dark-red, dark-blue, and dark-green grids, which are normalized by those of their corresponding pure-color ones, are here defined as no. PTN1, no. PTN2, and no. PTN3, respectively, for convenience in Figures 4 and 5.

2.2. : Comparison between Geodetic and Seismic Quantities

After the selection of one of the four possible pairs (here, the absolute value of dilatation rate and seismic energy (per two years), maximum shear strain rate and seismic energy, absolute value of dilatation rate and the number of earthquakes (per two years), and maximum shear strain rate and the number of earthquakes), its relationship was examined using scatter diagrams for circular regions, each having 60-km radius centered in 0.05°-by-0.05° spatial grids, over inland Japan; each pair was compared for two-year time windows centered in a three-month-interval time grid.

2.3. : Quantification of the Relationship Between Geodetic and Seismic Quantities

Quantification of the relationship between geodetic and seismic quantities requires us to define an index well characterizing the relationship. We here define the index as a statistical index R(t), where represents time. One of the simplest indices is correlation coefficient. Scatter diagrams for respective time windows, however, implied that the relationship between geodetic and seismic quantities showed relatively complicated patterns, which forced us to define another index as follows.(1)For a particular circular region centered in a spatial grid and for a two-year time window centered in a time grid, the mean of the absolute value of dilatation rate or the maximum shear strain rate was calculated from the grids ranking in the top 20% in terms of seismic energy or number of earthquakes. We defined this parameter as statistical index R(t), where represents the midpoint of each two-year time window.(2)Process (1) was carried out for each time window which is moved by three months. The linear trend of R(t) was then calculated for each time period from 1 January 2001 through 1 October 2007 with reference to the fixed time grid of 1 April 1998. The linear trend of R(t) thus obtained was classified into three categories on the basis of trend significance with a 95% confidence level: significant decrease, significant increase, and insignificant change.(3)Processes (1) and (2) are carried out for every spatial grid with an interval of 0.05°.

2.4. : Statistical Evaluation of the Relationship between Geodetic-Seismic Quantity Connection and Large Event Occurrence Times

The relationship between the occurrence times of large inland mainshocks and different types of statistical index R(t) was evaluated by calculating the probability gains of alarm rate (AR) and success rate (SR) for M ≥ 6.0 inland mainshocks and drawing an error diagram based on a contingency table. To calculate AR and SR for each pair of physical quantities, we need to tabulate the relationship between the proposed potential condition, which is proposed for the possible occurrence of physical events of interest, and the actual occurrence of the event (Table 2). The items in the row and column of contingency table correspond to the potential condition for event occurrence and the actual occurrence, respectively. Defining the time period for positive potential condition in which the events will be predicted to occur as the period of A and the period of negative potential condition as the period of B, each element of (a, b, c, d) in Table 2 can be explained as follows:(a)the number of cases in which physical events of interest occur during the period of A,(b)the number of cases in which physical events did not occur during the period of A,(c)the number of cases in which physical events of interest occur during the period of B,(d)the number of cases in which physical events did not occur during the period of B.

Table 2: Model contingency table for explaining the relationship between potential conditions for/against the occurrence of physical events such as large inland earthquakes and their actual occurrence. The relation can be classified into the following four elements: a: the number of “yes” warnings accompanied by the occurrence of events, b: the number of “yes” warnings not accompanied by the occurrence of events, c: the number of “no” warnings accompanied by the occurrence of events, d: the number of “no” warnings not accompanied by the occurrence of events.

Potential condition for event occurrence
Event occurrence

Yes
No

Positive
a
b

Negative
c
d

By this definition, a+b and a+c can be interpreted as follows.(a+b):Total number of the periods of A.(a+c):Total number of periods in which physical events of interest occurred.

If an event occurs during the time period of A, one count is added to element . Using elements of (a, b, c, d), alarm rate (AR), and success rate (SR) can be calculated as follows [5, 6]: The role of this step is to find a useful monitoring index/indices which feature the preparatory crustal activity of a large inland earthquake through a statistical evaluation process. Such a monitoring index is defined based on the information on grid counts showing a specific (ex. negative significant) trend of R(t) for a time window which are accompanied by high seismic activities. It should be noted that the performance of the monitoring index must be verified and validated in prospective way by updating the database.

We here focused on time variations of spatial grid counts which are accompanied by one-year release of total seismic energy larger than or equal to M4.0 earthquake (6.0 × 10¹⁰ Nm); the seismically active grid counts are specifically defined as no. PTN1 for those of significant decrease in R(t), no. PTN2 for those of significant increase in R(t), and no. PTN3 for those of insignificant change in R(t). In addition, the three kinds of grid counts were normalized by their respective total grid counts (pure-color plus dark-color grids in Figure 3) for each time period of calculation. Step tabulates thus calculated spatiotemporal distributions of no. PTNs1–3 and outputs the following figures Figure 4(a) spatial distributions of no. PTNs1–3 for each time period (Figure 3), Figure 4(b) temporal changes in no. PTN1, no. PTN2, and no. PTN3 (Figure 4), and Figure 4(c) temporal changes in no. PTN1/no. PTN2 (Figure 5); this figure represents superiority or inferiority in the normalized area of PTN1 and PTN2. Assuming the temporal changes in no. PTN1/no. PTN2 reflect the preparatory processes of large inland earthquakes, we attempted to formulate a positive potential condition for the occurrence of large inland earthquakes in terms of a way of temporal changes in no. PTN1/no. PTN2 as shown in the first row of Table 3.

Table 3: Contingency table showing the relationship between proposed potential conditions associated with the occurrence of M ≥ 6.0 inland mainshocks for four pairs of geodetic and seismic quantities and their actual occurrence times. The four figures for each element correspond to the following four pairs of geodetic and seismic quantities from left: the absolute value of dilatation rate and seismic energy, maximum shear strain rate and seismic energy, absolute value of dilatation rate and the number of earthquakes, and maximum shear strain rate and the number of earthquakes. A proposed potential condition occurring prior to the occurrence of inland mainshocks, Δ(no. PTN1/no. PTN2)<0, implies that the trend of no. PTN1/no. PTN2 is negative (Figure 5).

Occurrence of inland earthquakes (M ≥ 6)
Nonoccurrence
Alarm rate (AR)

(no. PTN 1/no. PTN2) < 0
5/5/4/3
7/12/7/10
0.417/0.294/0.364/0.231

(no. PTN 1/no. PTN2) > 0
0/0/1/2
13/8/13/10

Success rate (SR)
1.00/1.00/0.80/0.60

Probability gains of AR and SR: 3.56/2.49/3.00/1.97 and 3.58/3.49/2.73/2.13.

(a)

(b)

(c)

(d)

Figure 4: Temporal variation in no. PTN1, no. PTN2, and no. PTN3 (ref. Figure 3). Black arrows correspond to the occurrence date of M ≥ 6.0 inland mainshocks. Also shown are the error bars which were obtained by 100 times of shuffling of no. PTN1, no. PTN2, and no. PTN3 in time, respectively. The horizontal axis denotes the time corresponding to the end of the time period referred to for calculation of linear trend of R(t). The beginning of each time period is fixed to April 1, 1998. Panels (a) to (d) correspond to different four pairs of geodetic and seismic quantities, which are the same as those in Figure 3.

(a)

(b)

(c)

(d)

Figure 5: Temporal variation in the ratio of no. PTN1 to no. PTN2 (ref. Figure 3). Black arrows and the horizontal axis denote the same parameters as in Figure 4. Also shown are the error bars which were obtained by calculating the ratio between no. PTN1 and no. PTN2 shuffled in Figure 4. Panels (a) to (d) correspond to different four pairs of geodetic and seismic quantities, which are the same as those in Figures 3 and 4.

Statistical evaluation process at requires the probability gains [7, 8] of alarm rate (AR) and success rate (SR). To calculate the probability gains of AR and SR, we prepared reference AR and SR for each pair of geodetic and seismic quantities in the following way.(1)We applied the bootstrap method to shuffle the spatial distributions of no. PTNs1–3 grids and their corresponding seismically inactive grids for a particular time window.(2)Process (1) was carried out by moving the time window with its beginning time grid fixed. By repeating this, we obtained the temporal variation of the shuffled spatial distributions.(3)We applied the bootstrap method to shuffle the temporal variation obtained in process (2).(4)Processes (1) to (3) were repeated 100 times.(5)The 100 ARs and SRs obtained were averaged.

Thus, calculated references AR and SR were considered as being due to a random phenomenon.

3. Most Informative Pair of Geodetic and Seismic Quantities

Among the four pairs of geodetic and seismic quantities examined, the pair comprising the absolute value of dilatation rate and seismic energy was found to be the most closely related in terms of probability gains of both alarm rate (AR) and success rate (SR) for M ≥ 6.0 inland mainshocks (3.56 and 3.58, resp.) (Table 3). This pair showed the best performance even from the viewpoint of the error diagram (Figure 6). This is because the actual plot was located closest to position of the origin in Figure 6, which corresponds to the best performance of a proposed potential condition for event occurrence. The condition is represented in the first row of Table 3 or as the horizontal axis of Figure 6. This result was the same regardless of the approval methods and confidence levels. Figures 4(a) and 5(a) and the leftmost figures of respective elements in Table 3, which correspond to the pair comprising the absolute value of dilatation rate and seismic energy, indicate that M ≥ 6.0 inland mainshocks are most likely to occur during the time period of the decrease in no. PTN1 relative to no. PTN2 or the increase in no. PTN2 relative to no. PTN1.

(a)

(b)

(c)

(d)

Figure 6: Error diagram (Molchan diagram) for evaluating the significance of the performance of a proposed potential condition for the occurrence of large inland earthquakes.

On the assumption that this tendency holds for even a small region influenced by a change in Coulomb failure stress (CFS), a possible physical occurrence mechanism of a M ≥ 6.0 mainshock would be related to its preparatory process over a fault system. A mainshock would be triggered at the stage after microseismicity becomes active in the fault system, where stresses are considered to be accumulated on its asperities. Such a region may be accompanied by smaller strain rates. Because only a few samples of M ≥ 6.0 inland mainshocks were available to us, we need to continue to collect observations and further evaluate and validate the condition for event occurrence and the closest relation between geodetic and seismic quantities in the same manner as the one described in this study. We emphasize that this manner must not be changed in validating the condition for event occurrence. Furthermore, it is important to operate the statistical evaluation system to search for other informative pairs of physical quantities to gain a physical understanding of crustal activities as preparatory processes of large inland earthquakes.

4. Summary

Based on the notion that there are some correlations between different kinds of observations that may represent different aspects of crustal phenomena, we began with creating a database of various physical quantities and further constructed a basic system for statistically evaluating and validating a monitoring index which would reflect crustal activities associated with the occurrence of physical events such as large inland earthquakes. We have indicated that the most informative pair of geodetic and seismic quantities comprises the absolute value of dilatation rate and seismic energy in terms of probability gains of alarm rate (AR) and success rate (SR) and error diagram for M ≥ 6.0 inland mainshocks. The probability gains were calculated on the basis of contingency tables which show the relationship between the temporal changes in linear trends of statistical index R(t) in seismically active areas and the occurrence times of M ≥ 6.0 inland mainshocks. Further evaluation and validation of the best pairs and the search for other informative combinations of physical quantities are necessary to gain a physical understanding of crustal activities, particularly, the preparatory processes of large inland earthquakes.

Acknowledgments

The authors deeply acknowledge valuable comments provided by two anonymous reviewers, which highly contributed to the improvement of the manuscript. They are also grateful to Geospatial Information Authority of Japan, and Japan Meteorological Agency for letting them utilize GEONET (GPS Earth Observation NETwork) data and unified hypocenter catalog, respectively.

References

M. Imoto, “Earthquake probability based on multidisciplinary observations with correlations,” Earth, Planets and Space, vol. 58, no. 11, pp. 1447–1454, 2006.

M. Imoto, “Information gain of a model based on multidisciplinary observations with correlations,” Journal of Geophysical Research B, vol. 112, no. 5, Article ID B05306, 2007.

B. Gutenberg and C. F. Richter, “Magnitude and energy of earthquakes,” Annali di Geofisica, vol. 9, pp. 1–15, 1956.

M. Kawamura, T. Kudo, K. Yamaoka, and M. Furumoto, “Understanding of crustal activity with statistical approaches: I. Construction of spatiotemporal database with a unified data format,,” Journal of the Seismological Society of Japan (Second Series), vol. 63, pp. 21–33, 2010 (Japanese).

I. T. Jolliffe and D. B. Stephenson, Forecast Verification, John Wiley & Sons, New York, NY, USA, 2003.

S. Uyeda and A. Kumamoto, “Evaluation of the Kushida method of short-term earthquake prediction,” Proceedings of the Japan Academy Series B, vol. 80, no. 3, pp. 140–147, 2004.

K. Aki, “A probabilistic synthesis of precursory phenomena,” in Earthquake Prediction: An International Review, D. W. Simpson and P. G. Richards, Eds., Maurice Ewing Seires 4, pp. 566–574, AGU, Washington, DC, USA, 1981.

A. Helmstetter, Y. Y. Kagan, and D. D. Jackson, “High-resolution time-independent grid-based forecast for M ≥ 5 earthquakes in California,” Seismological Research Letters, vol. 78, no. 1, pp. 78–86, 2007.