Earthquake forecasting using the rate-and-state friction model and a smoothing Kernel: application to Taiwan

In this work, two approaches were employed for estimating the spatiotemporal distribution of seismicity density in Taiwan. With the use of the rate-and-state friction model, a model for short-term forecasting according to the fault-interaction-based rate disturbance due to seismicity was considered. Another long-term forecasting model that involves a smoothing Kernel function is proposed. The application of the models to Taiwan led to good agreement between the model forecast and actual observations. Using an integration of the two approaches, the application was found to be capable of providing a seismicity forecast with a higher accuracy and reliability. To check the stability related to the regression the bandwidth function, the forecasted seismicity rates corresponding to the upper and lower bounds of the 95 % confidence intervals are compared. The result shows that deviations within the bandwidth functions had an insignificant impact on forecasting reliability. Besides, insignificant differences in the forecasted rate change were obtained whenAσ was assumed to be between 0.1 and 0.4 bars for the application of the rate-and-state friction model. By considering the maximum Coulomb stress change among the seismogenic depth, the model presents a better forecasting ability than that using any single fixed target depth. The proposed methodology, with verified applicability for seismicity forecasts, could be useful for seismic hazard analyses.


Introduction
In the past few decades, studies and interests in earthquake forecasting have increased.Among these studies, the rateand-state friction model was proposed by Dieterich (1994).
Based on this model, seismicity rate change is attributed to a change in the Coulomb stress caused by large earthquakes.Catalli et al. (2008) applied this model to the 1997 Umbria-Marche sequence, Italy, and found that it was able to illustrate the main features of the temporal evolution of seismicity.However, Chan et al. (2010) examined this model by applying it to the entire Italian region and concluded a marginal forecasting ability.Therefore, it is not clear if the results can be attributed to uncertainties from dispensable assumptions of the Coulomb stress calculation (Catalli and Chan, 2012) or to infrequent large earthquakes in Italy.In order to clarify these ambiguities, it was desirable to examine this model, once again, by applying it to a region with a higher seismicity rate.
Taiwan is a region that has a high amount of seismic activity and that also has a good earthquake catalogue.The operation of the modern seismic network, the Taiwan Telemetered Seismic Network (TTSN), began in the early 1970s (Tsai et al., 1981), and was operated with a total of 25 stations.During the operation period, approximately 4000 events were recorded each year in Taiwan.In the early 1990s, the TTSN stations were integrated into the Central Weather Bureau Seismic Network (CWBSN), with a total of 75 stations (Shin, 1992).Afterward, real-time digital recordings have been performed with the monitoring system.With the system, the arrival times of P and S waves are selected manually for the determination of earthquake location.The CWBSN records approximately 20 000 events each year in a region of roughly 400 × 550 km.Therefore, Taiwan is a good candidate region for evaluating earthquake-forecasting models.
In this study, we tested the feasibility of two forecasting models using a catalogue for the area surrounding Taiwan.First, to build the short-term seismicity rate change imparted by recent earthquakes, the rate-and-state friction model of Dieterich (1994) was introduced.Second, we built a longterm model using an epicentre-smoothing Kernel, as presented by Woo (1996), based on the distribution of past earthquakes.We tested their feasibilities by forecasting the distribution of earthquakes and evaluating their uncertainties.

Methodologies
In the following, we describe the Coulomb stress change, the rate-and-state friction model and the smoothing Kernel function that were employed for construction of the forecasting models.

The Coulomb stress change
In order to evaluate the short-term seismicity rate using the rate-and-state friction model, we first computed the Coulomb stress change CFS caused by source events.Based on the constant apparent friction law (Harris, 1998;Cocco and Rice, 2002), the general expression of the Coulomb stress change CFS can be represented as follows: where τ is the shear stress change along the slip direction, µ is the apparent friction coefficient and σ n is the normal stress change on the assumed plane (note that a positive σ n represents unclamping).According to CFS, a positive stress change encourages the occurrence of conse-quent events; otherwise a negative stress change inhibits future seismicity.
Application of the Coulomb model requires knowledge of the rupture parameters for source events, such as the geometry of the rupturing fault and the size of the slip.For this requirement, we adopted a homogenous slip model using the dimensions and average slips derived from scaling laws, for the general focal mechanism as proposed by Wells and Coppersmith (1994), as follows: log (L) = −2.44 + 0.59M W ; (2) where L is the rupture length in kilometers, M W is the moment magnitude, W is the rupture width in kilometres, and AD is the average slip in metres.
Another key parameter for the CFScalculation is the receiver fault mechanism.In general, receiver faults can be represented in the following forms: (1) optimally oriented fault planes according to the combination of the regional stress field and a stress change caused by the source event (King et al., 1994); (2) geometries and slip rakes of active faults (Toda et al., 1998); and (3) fixed focal mechanisms of earthquakes in sub-regions (Chan and Stein, 2009).Hainzl et al. (2010) investigated the effect of considering more realistic fault systems in CFS and found that considering earthquake nucleation on multiple receiver fault orientations significantly changed the predicted spatial stress change pattern and the Nat.Hazards Earth Syst.Sci., 12, 3045-3057, 2012 www.nat-hazards-earth-syst-sci.net/12/3045/2012/ total number of triggered events.In this study, we followed the suggestions of Hainzl et al. (2010) and assumed a spatially variable receiver fault plane for each calculation grid node.We introduced focal mechanisms within the time window from 1991 to 2007, collected by Wu et al. (2010), as a reference for the focal mechanisms (Fig. 1a).We assumed that a receiver fault plane for each calculation grid node (Fig. 1b) consisted of the reference focal mechanism with the shortest epicentral distance.In other words, we assumed a temporally immovable fault orientation within the study region.The procedure is in accordance with applications obtained from previous studies (Toda et al., 2008;Hainzl et al., 2010;Chan et al., 2010;Catalli and Chan, 2012).For each grid node, we evaluated CFS based on both nodal planes and reported the higher one.
In order to minimise the depth uncertainty for the CFS calculation, Catalli and Chan (2012) evaluated CFS among the seismogenic depth and reported the maximum one for each calculation grid.In this study, we first follow the suggestions of Catalli and Chan (2012) and then discuss the influence of forecasting ability using depth uncertainties to point out the importance of maximum CFS for earthquake forecasting.We estimated the CFS within a homogeneous half-space by applying the COULOMB 3.2 programme (Toda and Stein, 2002).

The rate-and-state friction model
To quantify the impact of CFS on the seismicity rate, we used the rate-and-state friction model (Dieterich, 1994).We followed Chan et al. (2010) and represented the evolution of the seismicity rate R (M, x, t) using CFS with the n'th source event CFS n (x) at the site of interest x as a function of magnitude, M, and time, t, as follows: where λ (M, x) is the long-term seismicity rate, and R n−1 (M, x) is the short-term seismicity rate promptly before the occurrence of the n'th source event (i.e., R 0 = λ (M, x)), Aσ is a constitutive parameter of the model as described by Dietrich (1994), t n is the occurrence time of the n'th source event and t na is the aftershock duration.The relationship describes the short-term seismicity rate change by considering a series of source events.

The smoothing Kernel function
We estimated the long-term seismicity rate λ (M, x) at the site of interest, x, as a function of the magnitude, M, as described by Woo (1996), as follows: where K (M, x − x i ) is the Kernel function as a function of the magnitude and the distance between the site of interest, x, and the epicentre of the i'th earthquake, x i , T M represents the time within the complete catalogue for the magnitude, and N M represents the total number of earthquakes with a magnitude within the earthquake catalogue.In order to acquire the mean annual seismicity rate, the Kernel functions for all of the events in the earthquake catalogue K (M, x − x i ) were summed and divided by the duration of the complete catalogue, T M .In this study, we followed Woo (1996) and described the Kernel function, K (M, x − x i ), as follows: where PL denotes the power law index.The bandwidth function H (M) is defined as the mean of the distances between each event of magnitude M and its nearest neighbour, and can be represented as follows: where c and d are constants that can be obtained by regression from the earthquake catalogue.Based on the magnitudedependent bandwidth function, forecasting models can illustrate variations in seismic densities for different magnitudes.

The earthquake catalogue
The completeness magnitude (M c ) is a key factor, within the earthquake catalogue, needed for earthquake forecasting studies.In order to check the quality of the catalogue in Taiwan, the maximum curvature approach (Wiemer and Wyss, 2000) was employed for calculating the spatiotemporal evolution of M c .For the analysis, events with a focal depth less than 40 km from 1973 to 2009 were applied.
Since the CWBSN became a denser seismic network and switched in operation from a trigger to a continuous recording mode, increasing the detection of small events (Wu et al., 2008), we found that M c decreased significantly after 1994.Thus, we separated the catalogue into two periods (one from 1973 to 1993 and the other from 1994 to 2009) for the spatial M c analysis (Fig. 2).Events within the time window from 1973 to 1993 were major as recorded in the TTSN, so we named the catalogue TTSN.Using the same logic, we named events from 1994 to 2009 as the CWBSN catalogue.
The M c for the CWBSN catalogue (Fig. 2a) was lower than the TTSN (Fig. 2b).Regions with a M c ≤ 4.0 for the TTSN and a M c ≤ 3.0 for the CWBSN nearly occupied the same area.For our study area (Fig. 2c), we used the intersection of the two catalogues, regions with M c ≤ 4.0 for TTSN and M c ≤ 3.0 for CWBSN.The time window from 1973 to 2007 was considered to be the "learning period" for establishing the long-term forecasting model using the smoothing Kernel function to here, as the "testing period" for the retrospective forecast and for establishing a short-term forecasting model using the rate-and-state friction model.

The forecasting model using the rate-and-state friction model
In order to forecast the distribution of aftershock sequences or triggered earthquakes, we proposed a short-term forecasting model using the CFS incorporated into the rate-and-state friction model.Based on this model, earthquakes with small magnitudes or those that had occurred far in the past do not have a significant influence on the current seismicity rate (Catalli et al., 2008;Chan et al., 2010Chan et al., , 2012)), we considered M W ≥ 4.5 earthquakes that occurred during the testing period (Table 1) as source events for the rate change calculation.The focal mechanism, the depth and the M W of each event were determined using the moment tensor inversion provided on the website for the Broadband Array in Taiwan Source events have occurred Target events in each time span Before Eq.1 Before Eq.2 Before Eq.3 Before Eq.8 Before Eq.9 Before Eq.10 Before Eq.4 Before Eq.5 Before Eq.7 Cal. time Time span Beginning of 2008 to Eq.1 Eq.1 to Eq.2 Eq.2 to Eq.3 Eq.7 to Eq.8 Eq.8 to Eq.9 Eq.9 to Eq.10 Eq.3 to Eq.4 Eq.4 to Eq.5 Eq.5 to Eq.7 for Seismology (BATS, http://bats.earth.sinica.edu.tw/).The magnitude completeness of this catalogue is 3.9 since 2001.
For evaluating CFS, we considered an intermediate value of µ = 0.4.The value is in good agreement with a range of µ between 0.2 and 0.5, being referenced from the study of earthquake focal mechanisms in Taiwan (Hsu et al., 2010).Applying the rate-and-state friction model, previous studies (Toda and Stein, 2003;Toda et al., 2005;Catalli et al., 2008) have suggested that the physically reasonable range for Aσ is between 0.1 and 0.4 bars.We first assume a fixed Aσ of 0.2 bars and then discuss its influence on the calculations.The t a was assumed to be a function of the magnitude as proposed by Burkhard and Grünthal (2009) and Grünthal et al. (2009).Here, we evaluate the rate change solely imparted by source events without considering a longterm seismicity rate λ (M, x).In other words, we assume Before Eq.11 Before Eq.12 Before Eq.13 Before Eq.20 Before Eq.21 Before the end of 2009 Before Eq.14 Before Eq.15 Before Eq.18 Cal. time Time span Eq.10 to Eq.11 Eq.11 to Eq.12 Eq.12 to Eq.13 Eq.18 to Eq.20 Eq.20 to Eq.21 Eq.21 to the end of 2009 Eq.13 to Eq.14 Eq.14 to Eq.15 Eq.15 to Eq.18 λ (M, x) = 1 and report the seismicity rate change at each time point.The importance of a long-term seismicity rate will be discussed in the later part of this study.By considering the corresponding completeness in time and magnitude (Fig. 2), the seismicity rate change was calculated on a 0.2 • × 0.2 • grid.The calculated seismicity rate changes for different moments (Fig. 3), indicating that the seismicity rate increases near source events that have just occurred and then retreats with time.Based on this model, the distribution of earthquakes during the testing period (referred to as target earthquakes) can be associated with the triggered aftershock sequence.
In order to validate our results, we compared forecast seismicity rates with the distribution of target earthquakes using the Molchan diagram (Molchan, 1990(Molchan, , 1991)).The diagram was designed for evaluating earthquake forecasting ability and is presented as the fraction of failure to predict versus the fraction of space occupied by the expectation.Here, we present the "fraction of space occupied by alarm" as the proportion of the study area having a forecast seismicity rate equal to or higher than the threshold, defined as "alarm".The  "fraction of failure to predict" indicates the proportion of target earthquakes that had a lower forecast seismicity rate than the alarm.In other words, when data points are distributed along the diagonal line, the distribution of target earthquakes is uniform or independent of the forecast rate.When a convexity is present it suggests that the majority of the target earthquakes occurred within regions with a lower forecast rate as compared to the entire area.When a concavity is present it suggests that the majority of the target earthquakes occurred in an area with a higher forecast rate.An optimistic forecast result is represented by a condition of having the least space occupied by alarms, and the lowest percentage of target earthquakes with a failure to predict.We compared the forecasted seismicity rate obtained using the rate-andstate friction model with the locations of target earthquakes in the Molchan diagram (the black dots in Fig. 4).We evaluated seismicity rate changes in the study region at the occurrence time for each target earthquake and compared them with the corresponding grid cell of the target earthquakes.The Molchan diagram confirms the forecast ability, and only 28 % of target earthquakes were located within the study area having a low (<50 percentile) forecast rate change.

The forecasting model using the smoothing Kernel function
We forecast the long-term seismicity rate as a function of magnitude using the smoothing Kernel function by analysing earthquakes that occurred during the learning period.We determined c and d values of the bandwidth function using a linear regression of ln(H (M)) (Fig. 5).The analysis indicated that c and d values were 0.053 and 0.8653, respectively.With respect to uncertainty for regression, we obtained a 95 % confidence interval of 0.19 (dashed lines in Fig. 5).Below, we first construct a forecasting model using the regression curve of the c and d values and then discuss the influence of this uncertainty.Recommended values for the power law index in Eq. ( 7) are between 1.5 and 2.0, corresponding to the cubic or quadratic decay of the seismic activity within a hypocentral distance (Molina et al., 2001).Chan et al. (2010) suggested that differences between the results are insignificant when the power law index is assumed to be between the recommended values.Therefore, in this work, we assumed an intermediate value of 1.75.
The calculation grids are consistent with those of the forecasting model using the rate-and-state friction model (0.2 • × 0.2 • ).The highest seismicity rate (Fig. 6) was determined along the eastern coastline and offshore in the northeast, corresponding to a high crustal deformation rate along the boundaries between the Eurasia and Phillipine Sea Plates (Yu et al., 1997).Higher seismicity rates are determined for smaller magnitude ranges (e.g., Fig. 6a) rather than for larger ones (e.g., Fig. 6d), as the Gutenberg-Richter law (Gutenberg and Richter, 1954).The locations of high rate regions were slightly different in each magnitude range and could be associated with the magnitude-dependent bandwidth function (Fig. 5).
In order to validate forecasting ability, we compared the estimated seismicity rate change obtained by the smoothing Kernel with the distribution of target earthquakes using the Molchan diagram (the gray dots in Fig. 4).We found a positive correlation between the two.Furthermore, only 18 % of target earthquakes were located within the study area that has a low forecast seismicity rate.Generally, the forecasting ability of this approach was better than the rate-and-state friction model (the black dots in Fig. 4).

The forecasting model using a combination of the smoothing Kernel function and the rate-and-state friction model
In the sections above, we proposed forecasting models using the rate-and-state friction model and the smoothing Kernel function, and demonstrated their abilities.Here, we present another forecasting model using a combination of these approaches.In Eq. ( 5), we assume the forecasting model from the smoothing Kernel function as the input parameter, λ (M, x), for the rate-and-state friction model.For validating forecasting ability, we compared the calculated seismicity rate with the distribution of target earthquakes using the Molchan diagram (the light gray dots in Fig. 4).We evaluated the seismicity rate in the study region at the time of the occurrence of each target earthquake and compared the corresponding grid cell of the target earthquakes.Since only 16 % of target earthquakes were located within the study area having a low forecast rate, we found good agreement using this approach.Additionally, in a comparison of forecasting models using either the rate-and-state friction model (the black dots in Fig. 4) or the smoothing Kernel (the gray dots in Fig. 4), the integrated approach displayed the best forecasting ability.

Discussion
The long-term seismicity rate was evaluated based on the distribution of earthquakes during the learning period by the smoothing Kernel function.The results indicate that in the region without neighbouring earthquakes quiet seismic activity is expected in the future.Therefore, this approach is applicable to Italy (Chan et al., 2010) and Taiwan (this study), since they have a long catalogue period and a high seismicity rate, respectively.However, when the approach was applied for forecasting large events with long return period or/and to the region with short observation period, the limitation of this approach were exposed that seismicity rate may be underestimated.
Using the smoothing Kernel function, the model worked in combination with the bandwidth function by considering earthquakes in adjacent regions.The bandwidth function for the study region was obtained from the earthquake catalogue using linear regression (Fig. 5).In order to check the uncertainty related to this regression, we compared the forecasted seismicity rates corresponding to the upper and lower bounds of the 95 % confidence intervals of the bandwidth function (Fig. 7).The results indicated that, at most, 25 % of the studied regions had more than a 5 % difference for the forecast rate within the confidence interval.By associating this distribution with forecasting events, only a few forecasting events occurred in the region that had significantly different forecasting rates (e.g., the northern tip of Taiwan).In contrast, the difference was insignificant in regions where the majority of target earthquakes occurred.We suggest that de-viations within the bandwidth functions had an insignificant impact on forecasting reliability.
In order to obtain short-term forecasting ability and to calculate the fault-interaction-based disturbance on seismicity, we introduced the CFS that was incorporated into the rate-and-state friction model.Here, the rate-and-state friction model was applied by considering a fixed Aσ of 0.2 bars.In order to check the forecast uncertainty related to Aσ , we compared the rate changes by assuming a Aσ of 0.1 and 0.4 bars, corresponding to the upper and lower bounds of the physically reasonable ranges, respectively (Toda and Stein, 2003;Toda et al., 2005;Catalli et al., 2008).At the end of 2009, differences in the forecasted rate change were found to be insignificant when Aσ was assumed to be between 0.1 and 0.4 bars (Fig. 8).At least 73 % of the studied region had a less than 0.1 % difference within the physically reasonable range (the white regions in Fig. 8).In contrast, only 2 % of the studied regions had more than a 5 % difference (the red regions in Fig. 8).Chan et al. (2010) considered a fixed target depth for the Coulomb stress calculation and concluded a marginal forecasting ability for the rate-and-state friction model.Above, we considered the maximum CFS among the seismogenic depth instead and presented better forecasting ability than that using different target depth.Here, we prove the importance of target depth for CFS calculations.We evaluated the seismicity rate change imparted by the source events (Table 1) and compared them with the distribution of target  events using the Molchan diagram.The procedure shows the variance of forecasting qualities by assuming different target depths (Fig. 9).For any single fixed target depth (smaller dots in Fig. 9), the forecasting ability was marginal, and corresponded to the conclusion of Chan et al. (2010).On the other hand, better forecasting ability was acquired by considering the maximum CFS among seismogenic depths (larger dots in Fig. 9).The result can be attributed to depth uncertainties from rupture geometries, especially when homogenous slip models are applied (Catalli and Chan, 2012).

Conclusions
We applied the rate-and-state friction model and the smoothing Kernel function for forecasting to the region surrounding Taiwan.The results indicated good agreement between the forecasting models and observations.The smoothing Kernel function approach displayed a better forecasting ability than the rate-and-state friction model approach.However, a combination of these two methods provided the best result.Through this approach, we obtained not only a long-term background rate, but also the short-term rate evolution.The application of this approach to the 2010 Darfield, New Zealand, earthquake sequences has presented its utility (Chan et al., 2012).The Darfield sequence was initiated at the occurrence of the 4 September 2010 M = 7.1 Darfield earthquake.On 21 February 2011, the M = 6.3 Christchurch earthquake took place 40 km east of the epicenter of the Darfield earthquake.During the 2011 Christchurch earthquake, a large peak ground acceleration (PGA) was recorded in downtown Christchurch and resulted in severe damage and fatalities.If only a long-term rate model is considered, Christchurch is in a low seismic hazard region, despite the occurrence of the 2010 Darfield earthquake.Chan et al. (2012) evaluated the seismic rate evolution through the rate-and-state friction model.They obtained a significantly higher seismic rate expected after the 2010 Darfield earthquake.The application could provide a warning before the occurrence of consequent earthquakes and would be valuable for seismic hazard mitigation.

Fig. 1 .
Fig. 1.(a) The focal mechanisms from 1991 until 2007 collected by Wu et al. (2010).(b) The focal mechanisms as receiver faults for the CFS calculation.As shown in Fig. 2c, the dashed line illustrates the study area.

Fig. 2 .
Fig. 2. The magnitude of completeness (M c ) for (a) the TTSN and (b) the CWBSN catalogues within the time window from 1973 to 1993 and from 1994 to 2009, respectively, for shallow earthquakes (with a focal depth ≤ 40 km).(c) The study area, which is the intersection of regions with M c ≤ 4.0 for the TTSN (the dashed lines) and M c ≤ 3.0 for the CWBSN (solid lines), shown in gray.Note that a lower M c between Green Island and Orchid Island for the TTSN catalogue can be attributed to an additional station at Green Island.During the CWBSN period, the CWB did not install an instrument at Green Island and caused the effect.

Fig. 3 .
Fig. 3.The short-term seismicity rate change acquired using the rate-and-state friction model for various periods.Source events during 2008-2009 are shown as open green stars.Target earthquakes during each time span are shown as gray dots.Source parameters of the source events necessary for calculating the seismicity rate change are shown in Table1.

Fig. 4 .
Fig. 4. The Molchan diagram for investigating the correlation between different forecasting models and target earthquakes.The dashed line denotes half of the fraction of space occupied by the alarm, and the corresponding fraction by the failure to predict using each model presented.

Fig. 5 .
Fig.5.Bandwidth functions within the 95 % confidence interval for the study area using the linear regression of earthquakes that occurred following catalogue completeness by considering the time and magnitude.Note that the mean distance of the nearest event is scaled in ln(km).

Figure 6 Fig. 6 .
Figure 6 Fig. 6.The distribution of the reference seismicity rate acquired for different magnitude ranges.Forecasting earthquakes are shown as white open circles.

Fig. 7 .
Fig. 7.The difference in the seismicity density rate by considering the upper and lower bounds of the 95 % confidence interval for the bandwidth function, as shown in Fig. 5. Target earthquakes are shown as open circles.

Fig. 8 .
Fig. 8.The difference of the seismicity rate change for different constitutive parameters (Aσ ) within the rate-and-state friction model.Source events are shown as open stars.Target events are shown as open circles.The source events for calculating the seismicity rate change are shown as stars.As shown in Fig. 2c, the dashed line illustrates the study area.

Table 1 .
Source parameters of the source events for calculating seismicity rate change using the rate-and-state friction model.No.Year Month Day Longitude (• ) Latitude ( • ) Magnitude Depth (km) Strike ( • ) Dip ( • ) Rake ( • )