Aftershock Hazard Magnitude , Time , and Location Probability Forecasting

This study combines branching aftershock sequence (BASS) and modified Omori’s law to develop a predictive model for forecasting the magnitude, time, and location of aftershocks of magnitude Mw ≥ 5.00 in large earthquakes. The developed model is presented and applied to the 17:47 20 September 1999 Mw 7.45 Chi-Chi earthquake Taiwan, 09:32 5 November 2009 (UTC) Nantou Mw 6.19, 00:18 4 March 2010 (UTC) Jiashian Mw 6.49 earthquake sequences, Taiwan, and 05:46 11 March 2011 (UTC) Tohoku Mw 9.00 earthquake, Japan. The estimated peak ground acceleration (PGA) results are remarkably similar to calculations from the recorded magnitudes in both trend and level. This study proposes an empirical equation to improve the aftershock occurrence forecast time. The forecast time results were greatly improved. The magnitude of aftershocks generally decreases with time. It was found that the aftershock forecast probability of Mw ≥ 5.00 is high in the first six days after the main shock. The results will be of interest to seismic mitigation specialists. Spatial and temporal seismicity parameters to the aftershock sequence investigation into the 17:47 20 September 1999 (UTC) Mw 7.45 Chi-Chi earthquake, Taiwan found that immediately after the earthquake the area closest to the epicenter had a lower b value. This pattern suggests that at the time of the Chi-Chi earthquake, the area closest to the epicenter remained prone to large magnitude aftershocks and strong shaking. With time, however, the b value increased, indicating a reduced likelihood for large magnitude aftershocks.


1.INTRODUCTION
Aftershocks following large to moderate earthquakes are potentially hazardous.They are also rich in seismic data.The 17:47 20 September 1999 (UTC) M w 7.45 Chi-Chi earthquake struck Central Taiwan, causing 2489 confirmed deaths with 50 people still missing, and damaging a total of 106159 buildings (Tsai et al. 2001).This study provides spatial and temporal analysis of the p (the power law decay of modified Omori law) and b (a constant of the Gutenberg-Richter relation) values of the aftershock sequence to the Chi-Chi earthquake and uses these parameters to estimate the aftershocks associated with the 09:32 5 November 2009 (UTC) Nantou M w 6.19, and 00:18 4 March 2010 (UTC) Jiashian M w 6.49 earthquake sequences.
The primary goal of most earthquake research is hazard mitigation (Holliday et al. 2008).While the main shock is generally the most destructive phase of an earthquake sequence, aftershocks can compound the initial destruction, interfere with rescue efforts, and even create further destruction.Large aftershocks can be particularly hazardous if they occur in populated areas such as the M w 5.80 aftershock to the 1999 M w 7.40 Izmit, Turkey, earthquake.This aftershock killed 7 and injured a further 420 people (Wiemer et al. 2002).Toda et al. (1998), and Wiemer and Katsumata (1999) systematically analyzed the spatial distribution of seismicity within aftershock zones.Their results demonstrate that seismic activity, a, magnitude-frequency distribution of events, b, and the decay rate, p, all show significant spatial variation.These results are consistent with our own.Forecasting the probability of aftershocks to large earthquakes has progressed considerably (Ebel 2009).Earthquake aftershock sequences have been found to approximately satisfy three empirical scaling relations: (1) the Gutenberg-Richter frequency-magnitude scaling (Gutenberg and Richter 1954); (2) Bath's law for the difference in main shock magnitude and its largest aftershock (Båth 1965); and (3) modified Omori's law for temporal decay in aftershock rates (Shcherbakov et al. 2004).In this paper, we replace the frequency-magnitude relation with the branching aftershock sequence (BASS) model to predict aftershock locations and magnitude.The BASS model recognizes that each earthquake has an associated aftershock sequence with the parent shock creating generation upon generation of aftershocks.The model follows Bath's law and is fully selfsimilar.The most important parameter is the magnitude difference.According to Bath's law the magnitude difference is the difference between the parent shock and the largest expected aftershock.If this value is positive the number of aftershocks will eventually die out.It is theoretically possible to determine the total number of aftershocks.In this study, we calculated magnitude difference as follows: the largest aftershock to the 17:47 20 September 1999 M w 7.45 Chi-Chi earthquake was M w 6.90.The largest aftershock to the 09:32 5 November 2009 Nantou M w 6.19 is M w 5.74; and the largest aftershock to the 00:18 4 March 2010 Jiashian M w 6.49 was M w 5.72.Therefore, the average difference in magnitude between these main shocks and their largest aftershocks is 0.59 ± 0.16.The value 0.59 ± 0.16 is almost the same as 0.55 (the magnitude difference between the main Chi-Chi earthquake shock and its largest aftershock), therefore, we set the magnitude difference at 0.50.
Besides predicting the magnitude and location using the BASS model, we apply Omori's law to illustrate the temporal decay of aftershocks.We focus primarily on M w ≥ 5.00 aftershocks of the Chi-Chi, Nantou, and Jiashian earthquakes and analyze the difference between the aftershock predictions and occurrence observations.The predicted versus actual peak ground acceleration trends are examined for the Nantou and Jiashian areas.A successful predictive model could be of value to rescue and reconstruction crews in earthquake zones.

DATA COLLECTION AND ANALYSIS
This study uses the data from Chen and Tsai's research (2008).We converted the original Taiwan' earthquake catalog 1900 -2006 magnitudes into homogenized M w magnitudes.The catalog was updated in this study updated to include earthquakes up to 2010.The aftershock distribution for the 17:47 20 September 1999 M w 7.45 Chi-Chi earthquake Taiwan is shown in Fig. 1.The figure shows aftershocks for the period 17:47 20 September 1999 to 17:47 20 September 2000 (Tajima and Kanamori 1985).In order to under-stand the spatial variations in b values, we divide the aftershock zone (Fig. 1) into three equal sections (based on equal range): (1) north of the epicenter; (2) at the epicenter; and (3) south of the epicenter (Fig. 2).We used the total recorded aftershock data to estimate a, b value variations spatially and temporally (Fig. 3).The data were then sorted according to the divided area or time range.We then used the Gutenberg-Richter relation to estimate the a and b values.Data from magnitudes ≥ 3.00 were used to understand spatial variations in the p value (Reasenberg and Jones 1989).We then sorted the data according to the divided area and use the modified Omori's law to estimate the p value (Fig. 4).
The next step in this analysis involved determining if it is possible to predict the observed patterns using the BASS model and modified Omori's law.That is, we are attempting to forecast the time, magnitude, and location of aftershocks.

BRANCHING AFTERSHOCK SEQUENCE MODEL AND MODIFIED OMORI'S LAW
The BASS model recognizes that aftershock sequences are self-similar multi-generational series.The governing parameter in the BASS model is the magnitude difference.If the magnitude difference is positive the earthquake sequence will then diminish with time.If we specify a minimum aftershock magnitude, we can theoretically determine the total number of aftershocks following the main shock.
The BASS formulation requires that the frequencymagnitude distribution of each order of aftershocks (Holliday et al. 2008) be expressed as: where m d is the magnitude of an aftershock; N d (≥ m d ) is the number of aftershocks of magnitude equal to or larger than m d ; and a d and b d are the a and b values of the distribution.Shcherbakov and Turcotte (2004) introduced a definition of magnitude difference (Δm * ) whereby Δm * is the difference in magnitude between the main shock and its inferred largest aftershock.It is less than the magnitude of the main shock, m p , such that: Substituting this condition to satisfy Eq. ( 1) whereby ( ) This relation fully specifies the frequency-magnitude distribution of each family of aftershocks.The formula implies an infinite number of small earthquakes.To eliminate this singularity, we must prescribe a minimum magnitude earthquake m min to allow us to obtain the total number of aftershocks; i.e., set m d = m min in Eq. ( 3) : This relation is the essential feature of the BASS model.
If we use the simple aftershock model of Reasenberg and Jones (1989, 1990, 1994), we can obtain an expression for the aftershocks in terms of time: .From which, we can get the a value.Therefore, this equation describes the parameters p, a, and b.Using Eq. ( 4) to replace the numerator in Eq. ( 5), we can get an expression for aftershocks greater than m min in terms of time: The above combined expression has never been previously used.From Eq. ( 6), we are able to obtain the theoretical number of aftershocks of magnitude equal to or larger than m min by expected time, marking the importance and usefulness of this expression.By consequently using Eqs.( 3) and (4), we obtain the cumulative distribution function P cm for aftershock magnitude: where m d is the aftershock magnitude.In order to simulate the aftershock magnitude, we randomly choose a value for P cm in the range 0 < P cm < 1.By rearranging Eq. ( 7) we obtain the aftershock magnitude m d : Next, we determined occurrence time.Because aftershock sequences satisfy the general form of Omori's law (Shcherbakov et al. 2004), then: where R(t d ) is the rate of aftershock occurrence, and x , c, and p are parameters.The number of aftershocks that occur after a time t d is given by (Holliday et al. 2008): We can obtain the total number of aftershocks by setting t d = 0, and ( 1) Using Eqs. ( 10) and ( 11), we obtain the cumulative distribution function P ct giving the time of occurrence for aftershocks: In order to simulate the aftershock occurrence time, we randomly choose a value for P ct in the range 0 < P ct < 1.The aftershock occurrence time is then given by: Finally, we require the radial distance of aftershock occurrence relative to the main shock.According to Felzer and Brodsky (2006), the cumulative distribution function P cr for the radial distance r d of each aftershock is: In order to simulate the aftershock occurrence locations, we randomly choose P cr in the range 0 < P cr < 1.The radial aftershock occurrence distance relative to the main shock is determined as: In order to fully determine the aftershock locations their direction relative to the main shock θ d must be specified.Therefore, the aftershock direction is chosen at random from 0 < θ d < 2π.Equations ( 6), ( 8), (13), and (15) allow us to forecast the magnitude, time, and location of each aftershock of magnitude ≥ M.

FORECASTED VERSUS RECORDED DATA
We used the total data from the Chi-Chi earthquake (as shown in Fig. 1) to estimate the a and b values in (Fig. 5a), and M w ≥ 3.00 to estimate the p and c values (Reasenberg and Jones 1989) in Fig. 5b in order to increase the confidence in our forecast data.We used these p, c, and b parameters to forecast the magnitude, time, and location of each aftershock of M w ≥ 5.00.Selection of the magnitude of predicted aftershocks is based on smaller magnitude events being less likely to present seismic hazard.
Destructive earthquakes such as the 17:47 20 September 1999 M w 7.45 Chi-Chi earthquake are typically followed by large aftershocks.These aftershocks can impede rescue efforts and endanger the lives of rescue teams and survivors.Naturally, rescuers would greatly benefit from a system that could predict aftershocks and give information on peak ground acceleration (PGA).Therefore, we combined Eqs. ( 6), ( 8), (13), and (15) to forecast the magnitude, time and location of larger aftershocks in Fig. 6.From this figure, we can see that most large aftershocks occur within the first hours after the main shock.Next we applied the acceleration attenuation relationship for Taiwan (Liu and Tsai 2005) to calculate the peak ground acceleration for each of the forecasted and recorded magnitudes in each grid, and then sorted the maximum PGA for each grid (Fig. 7).From Figs. 7a and b through to Figs. 7k and l, we can see that the pattern and level of forecasted results are similar to the results calculated from recorded magnitudes, standard deviation, and correlation coefficients.These are: Day 1: 30.14 gal, 0.936; Day 2: 36.25 gal, 0.903; Day 3: 8.94 gal, 0.880; Day 4: 17.14 gal, 0.689; Day 5: 13.38 gal, 0.848; and Day 6: 54.04 gal, 0.895.The correlation coefficients are high.The accuracy of the modeled data compared to the actual data given in Figs. 6 and 7 shows that this model could potentially benefit rescue teams by predicting the location of aftershocks and the corresponding degree of ground shaking.

IMPROVED FORECAST TIMES FOR AFTERSHOCKS
To improve the accuracy of forecasting aftershock occurrence times, we used the first day of aftershocks of M w ≥ 5.00 to calculate the empirical relation as follows: . . .log T M M 0 471 0 527 7 12 Where T is the f aftershock occurrence time in minutes; M min is magnitude of the smallest aftershock considered; and M p is the magnitude of the main shock or preceding aftershocks.M p is given by Eq. ( 8).The results comparison is a random process (Figs.8a, c, e, g, i, and k) with the empirical relation results (Figs. 8b, d, f, h, j, and l) given in Fig. 8. From Fig. 8, it is evident that the spread range of the empirical relation is better than that of the random process.To understand the probability of aftershock occurrence, we use the following equation: Where m d is the magnitude of the forecast aftershocks; and t is the improved aftershock occurrence forecast time.The result is shown as Fig. 9. From Fig. 9, we can see the probability of forecasting aftershocks of magnitude equal to or larger than 5.0 is high in the first six days after the main shock.This result provides important information for rescue workers.

MODEL APPLICATION TO THE NANTOU AND JIASHIAN EARTHQUAKES
We applied the parameters calculated from 17:47 20 September 1999 M w 7.45 Chi-Chi earthquake, Taiwan to the 09:32 5 November 2009 M w 6.19 Nantou earthquake, Taiwan and 00:18 4 March 2010 M w 6.49 Jiashian earthquake, Taiwan (Fig. 10).From Fig. 10, we can see that the forecasting of aftershocks in terms of time is still not perfect for the Nantou and Jiashian earthquakes.However, the aftershock location and magnitude forecast approximates the recorded data.Figure 11 gives the forecasted PGA for aftershocks to the Nantou and Jiashian earthquakes.Figures 11a and b the show standard deviation and correlation coefficient to be 5.20 gal and 0.976 for the Nantou earthquake, respectively.Figures 11c and d give the standard deviation and correlation coefficient as 13.42 gal and 0.977 for the Jiashian earthquake, respectively.Note the correlation coefficients are high.The forecasted PGA matches the trend well in calculated PGA from the actual data for each of the grids.This is an important facet of replication as PGA describes the degree of ground shaking rescuers may experience.

MODEL APPLICATION TO THE 11 MARCH 2011 JAPAN EARTHQUAKE
An earthquake of magnitude M w 9.00 struck Tokohu, Japan at 05:46 on 11 March 2011 (UTC).The massive offshore earthquake triggered a devastating tsunami that destroyed a good deal of Fukushima prefecture's foreshore and inland farming and industrial regions.To help understand how this model may be of use in mitigating seismic hazards and reducing the loss of life and property from such extreme events, we apply the model [Eqs.( 6), ( 8), ( 15)] and empirical relation [Eq.( 16)] using parameters calculated from Chi-Chi aftershocks, Taiwan, to forecast the magnitude, time, and location of aftershocks of magnitude M ≥ 6.00 to the 11 March 2011 Tokohu earthquake, Japan.The result is shown in Fig. 12. From Fig. 12, we can see the result is remarkably good.The estimated peak ground acceleration is shown in Fig. 13.Although, the distribution range of aftershocks is broad (1200 × 1200 km), after comparing the PGA contour between the forecasted PGA (from forecasting magnitude) and calculated PGA (from recorded magnitudes) (USGS Web site) in Figs.13a and b through to Figs. 13i and j, we can see that the estimated peak ground acceleration results are remarkably similar to calculations from recorded magnitudes (USGS Web site) in both trend and level.The respective standard deviation and correlation coefficients are: Day 1: 73.08 gal, 0.678; Day 2: 4.01 gal, 0.952; Day 4 : 4.61 gal, 0.981; Day 5: 1.32 gal, 0.998; and Day 6: 1.72 gal, 0.913.(Three days after the main shock, there were no aftershock events with magnitude ≥ 6.00.)The correlation coefficients are high.Such good results should be of interest to seismic mitigation specialists and rescue crews.
Fig. 6.The comparison of magnitude, time, and location of forecasted aftershocks with recorded aftershock data respectively to the 21 September 1999 as follows: 1 day after main shock (a) and (b); days after main shock (c) and (d); days after main shock (e) and (f); days after main shock (g) and (h); days after main shock (i) and (j); and days after main shock (k) and (l).(a) (e) (i)

DISCUSSION AND CONCLUSION
One area that has seen progress in terms of earthquake prediction is the prediction of aftershocks to strong earthquakes as evidenced by the work of Ebel (2009).However, Ebel's valuable work only describes the range of a, b, and p parameters, and does not forecast the magnitude, time, and location of aftershocks.This study, first established the parameters a, b, and p using the total recorded data from the 17:47 20 September 1999 Chi-Chi earthquake Taiwan.Figure 2, shows the b value for aftershocks is lowest at the epicenter of the Chi-Chi earthquake.The low b value suggests large aftershocks and strong shaking were to be expected close by the epicenter of the Chi-Chi earthquake (Wyss and Stefansson 2006).This result is consistent with the distribution of aftershocks shown in Fig. 1. Figure 3 shows the temporal variations in b values.It is evident from the figure that b values were lowest in the first 7 days after the earthquake and then increased with time.This indicates event occurrence potentially being more likely at lower magnitudes with time.Figure 4 shows spatial variations in p values for aftershocks to the Chi-Chi earthquake.The results of this analysis show that p values are larger south of the epicenter than at the epicenter or to its north.This evidence suggests the possibility of aftershock decay being faster to the south of the epicenter.Figures 2 to 4 show that the area closest to the epicenter of the Chi-Chi earthquake was at the highest risk for aftershocks.In fact, many aftershocks of M w ≥ 5.00 did occur within the central epicenter area.We used these parameters to predict aftershocks to the Chi-Chi earthquake as well as the 09:32 5 November 2009 M w 6.19 Nantou earthquake, Taiwan, 00:18 4 March 2010 M w 6.49 Jiashian earthquake, Taiwan, and 05:46 11 March 2011 M w 9.00 Tokohu earthquake, Japan.The predictive model used the BASS model in place of the Gutenberg-Richter frequency-magnitude relationship.This allows for the prediction of magnitude and location of aftershocks.Modified Omori's law gives the temporal decay.Though mathematical manipulation, we achieved Eqs. ( 6), ( 8), (13), and (15) above.These equations provided the model predictive skill for the magnitude, time decay, and location of each aftershock of magnitude ≥ M. We made use of the empirical relation [Eq.( 16)] to improve aftershock time forecasting.The empirical relation results are better than those derived from random processes.Finally, we predicted the PGA for the Nantou and Jiashian aftershock sequences.The resulting PGA trends match the actual PGA trends well.We also applied these parameters to forecast the magnitude, time, and location of aftershocks of magnitude M ≥ 6.00 to the 05:46 11 March 2011 Tokohu earthquake, Japan.Although the distribution range of aftershocks is broad (1200 × 1200 km), after comparing PGA contour maps between forecasted PGA and calculated actual PGA, we think the result is remarkably good, and the estimated peak ground acceleration results are remarkably similar to the actual calculations from recorded magnitudes in both trend and level.The results of this study will be of interest to hazard mitigation specialists, especially those concerned with rescue efforts following large destructive earthquakes.

DATA AND RESOURCES
All data sources were taken from published works listed in the References.Some plots were made using the Generic Mapping Tools version 4. 3.1 (www.soest.hawaii.edu/gmt;Wessel andSmith 1998, last accessed August 2006.).

Fig. 1 .Fig. 2 .
Fig. 1.The epicenter distributions of Chi-Chi earthquake aftershocks studied from 21 September 1999 to 21 September 2000 (local time).White star is the Chi-Chi main shock M w 7.45, and the red lines are the fault lines.
Fig. 3.The b value varies temporally.The vertical dashed line is the range of regression, and the thicker dashed line is the extrapolation part.(a) Aftershock 0 -7 days, (b) aftershock 7 -30 days, (c) aftershock 30 -365 days.
of aftershocks with magnitude equal to or larger than the magnitude threshold, M c ; M m is the magnitude of the main shock, occurring at time t; b is derived from the Gutenberg-Richter relation; and 10 ( ) a b M M m + -l is equal to the value of k in the modified Omori's law:

Fig. 4 .
Fig. 4. The p value spatial variation.The range is the same as Fig. 2's respective areas.
Fig. 5. (a) The a and b values calculated from total data (Fig. 1), (b) p, and c parameter calculated from total data (Fig. 1), and selected with magnitude ≥ 3.0.The parameters are applied to Nantou M w 6.19, Jiashian M w 6.49 earthquakes, Taiwan and 05:46 11 March 2011 (UTC) Tohoku M w 9.0 earthquake, Japan.

Fig. 7 .
Fig. 7.The comparison of maximum PGA of each grid calculated from forecasted magnitudes and from recorded magnitude with respect to Fig.6as follows: 1 day after main shock (a) and (b); 2 days after main shock (c) and (d); 3 days after main shock (e) and (f); 4 days after main shock (g) and (h); 5 days after main shock (i) and (j); and 6 days after main shock (k) and (l).

Fig. 8 .Fig. 9 .Fig. 10 .
Fig.8.The comparison of forecasting time of occurrence of aftershocks between random processes and empirical relation as follows: 1 day after main shock (a) and (b); 2 days after main shock (c) and (d); 3 days after main shock (e) and (f); 4 days after main shock (g) and (h); 5 days after main shock (i) and (j); and 6 days after main shock (k) and (l).

Fig. 11 .Fig. 12 .
Fig. 11.The comparison of maximum PGA of each grid calculated from the respective forecasted magnitudes and from recorded magnitudes with respect to Fig. 10: (a) 1 day after main shock; (b) 1 day after main shock; (c) 1 day after main shock; (d) 1 day after main shock.

Fig. 13 .
Fig. 13.The comparison of maximum PGA of each grid calculated from the respective forecasted magnitudes and recorded magnitudes with respect to Fig.12: 1 day after main shock (a) and (b); 2 days after main shock (c) and (d); 4 days after main shock (e) and (f); 5 days after main shock (g) and (h); and 6 days after main shock (i) and (j).(Because 3 days after main shock, there are no events of aftershock with magnitude ≥ 6.00, and we do not compare it.)