Interactive comment on “ Importance of three-dimensional grids and time-dependent factors for the applications of earthquake forecasting models to subduction environments

The paper evaluates future earthquake occurrence in two regions, one in the northeast of Taiwan (Ryukyu) and the other off Japan (Kanto) with two different models: one timeindependent (Woo’s method) and one time-dependent (rate-and-state friction model), in 2D and 3D. I found the paper basically clear, but there are several elements that need clarification, as I point out in what follows. Formula (1) gives a rate lambda that has the following dimensions: number of earth-


Introduction
Earthquake forecasting models generally provide essential knowledge for seismic hazards mitigation; i.e. they point out the regions with high seismicity activity and provide fundamental information regarding seismic hazard assessment (Marzocchi et al., 2003;Lombardi and Marzocchi, 2009).Therefore, studies and general interest in this issue have significantly increased and many forecasting models have been proposed.
However, most forecasting studies focus on the crustal earthquakes; i.e. the credibility of these models remains con-troversial for the application to subduction environments.Such regions include non-vertical seismogenic structures; depth-independent grid cells thus become crucial for forecasting models.Besides spatial distribution, temporal evolution of seismicity is another factor that dominates forecasting precision.Wiens et al. (1997) concluded that in comparison to the sequence of crustal earthquakes, a smaller number of aftershocks follows the occurrence of a subduction event.Such differences in temporal behaviour might result in a forecasting bias.
Therefore, this study applies several forecasting models and discusses their feasibility for applications to subduction regions.To model the behaviours of non-vertical seismogenic structures precisely, approaches with threedimensional grids cells are developed.By comparing these approaches with those with two-dimensional cells, I show the importance of the depth component.To reveal the role of the temporal factor for forecasting, I evaluate the forecasting ability of time-independent and renewal models.I apply the models to two subduction regions, the southwestern portion of the Ryukyu and Kanto.
C.-H. Chan: Earthquake forecasting models for the subduction zone 2.1 The smoothing kernel function Woo (1996) proposed a forecasting model, which described time-independent seismicity rate density λ (M, x) at the site of interest, x, as a function of magnitude, M, as follows: where K M, x − x i represents a smoothing kernel as a function of magnitude and distance between the site of interest, x, and the location of the ith earthquake, x; T M represents the period of a complete catalogue with a magnitude threshold; and N M represents the total number of earthquakes with magnitudes larger than the threshold.The unit for λ (M, x) can be per area or volume, depending on the dimensions of forecasting models.Although λ (M, x) is a function of magnitude, for simple representation, forecasting seismicity rates in different magnitude ranges will be shown in this study.This study follows the procedure of Woo (1996) and describes the kernel function K M, x − x i as follows: where PL represents the power law index and H (M) represents the bandwidth function, defined as the nearest distance to other events for each magnitude bin, M. The function can be represented as follows: where c and d are constants and c is a length, obtained from regression analysis of earthquake spatial distribution.Molina et al. (2001) recommend a PL of between 1.5 and 2.0, and Chan et al. ( 2010) concluded insignificant differences between the results when the PL is assumed to be between the recommended values.This study thus assumed an intermediate value of 1.75.The kernel function represents seismicity rate as a function of magnitude, and its feasibility has been proven through implementation in various regions (e.g.Molina et al., 2001;Beauval et al., 2006;Chan et al., 2010Chan et al., , 2012)).
The smoothing kernel function forecasts seismicity rate based on the seismic activity during an observation period; i.e. this model minimizes the factor of temporal evolution and provides a time-independent model.

The rate-and-state friction model
Another implemented forecasting approach is the rate-andstate friction model (Dieterich, 1994), which evaluates seismicity rate evolution based on earthquake Coulomb stress changes.According to the constant apparent friction law (Harris, 1998;Cocco and Rice, 2002), Coulomb stress change by a source event, n, CFS n (x) at the site of interest x as a function of magnitude, M, and time, t, as below: where τ n (x) represents the shear stress change along the slip direction; µ represents the apparent friction coefficient; and σ n (x) represents the normal stress change on the assumed plane.The law suggests that a positive Coulomb stress change enhances the occurrence of subsequent events, while a negative one inhibits future seismic activity.According to the law, however, the Coulomb stress change model does not quantify seismicity rate changes.
To quantify seismicity rate evolution, Dieterich (1994) proposed the rate-and-state friction model.This model presents the evolution of the seismicity rate R n (M, x, t) by considering the Coulomb stress change by the nth source event CFS n (x) as below: , where λ (M, x) represents the time-independent seismicity rate shown in Eq. ( 1); R n−1 (M, x) represents the seismicity rate change just before the occurrence of the nth source event (i.e., R 0 = λ (M, x)); Aσ represents a constitutive parameter of the model with the dimension of a stress; t n represents the occurrence time of the nth source event; and t a represents the aftershock duration.The rate-and-state friction model forecasts the temporal evolution of seismicity rate after occurrence of large earthquakes.
3 Forecasting application to the Ryukyu region

Tectonic setting and earthquake catalogue
The southwestern portion of the Ryukyu trench near Taiwan is seismically active as the Philippine Sea Plate subducts under the Eurasian Plate from the south (Fig. 1).In addition to high seismic activity, this region contains a good-quality earthquake catalogue.The Taiwan Telemetered Seismic Network (TTSN), the modern seismic network in the Taiwan region, was initiated in the early 1970s (Tsai et al., 1981).Since its operation, approximately 4000 earthquake events have been recorded each year.After the early 1990s, TTSN stations were integrated into the Central Weather Bureau Seismic Network (CWBSN), which records approximately 20 000 events each year (Shin, 1992).With a large amount of seismic activity and high-quality earthquake catalogues, the region is an ideal site for earthquake forecasting tests.Implementing a complete portion of an earthquake catalogue is a key factor for precise forecasts.I checked the magnitude of completeness, M c , for the catalogues using the maximum curvature approach (Wiemer and Wyss, 2000).Due to station coverage, both the TTSN and CWBSN catalogues (represented in Fig. 1a and b, respectively) obtain better observation quality inland than in the offshore region.The M c for the CWBSN catalogue (Fig. 1b) was lower than that for the TTSN (Fig. 1a), and the regions with a M c ≤ 4.0 for the TTSN and a M c ≤ 3.0 for the CWBSN are nearly the same.Thus, the intersection of the two catalogues, regions with M c ≤ 4.0 for TTSN and M c ≤ 3.0 for the CWBSN (Fig. 1c), determines the study region and the magnitude thresholds.I applied the earthquakes before 2009 for model construction and referred those in 2010 and 2011 as forecasting events for retrospective tests.Thus, T M is 37 and 14 years, while N M is the total number of earthquakes with magnitudes larger than 4.0 and 3.0 for TTSN and CWBSN cases, respectively.
Based on the parameters of the earthquakes before 2009, the linear regression determined that the c and d values of the bandwidth function in Eq. (3) were 0.0174 km and 1.1209, respectively.

The rate-and-state friction model
To calculate Coulomb stress change ( CFS), rupture behaviours of source earthquakes and source parameters of receiver fault planes are two important factors.For the source earthquake parameters, the hypocentre location, the moment magnitude and the focal mechanisms were obtained from the Broadband Array in Taiwan for Seismology (BATS) website (http://bats.earth.sinica.edu.tw/), and fault dimension and magnitude of slip were determined through the scaling laws of Yen and Ma (2011).For source parameters of receiver fault planes, both location and focal mechanism should be defined and utilized for the calculation.The loca- tions are assumed to be the coordinates and the target depth of each calculation grid.For receiver fault mechanisms, I followed the procedure of Catalli and Chan (2012) and assumed a spatially variable receiver fault plane for each calculation grid.A receiver fault plane for each grid is consistent with the closest reference focal mechanism determined by Wu et al. (2010).For each grid node, I evaluated the CFS on both nodal planes and reported the higher value.To minimize depth uncertainty, this study followed the procedure of Catalli and Chan (2012) that evaluated the CFS for seismogenic depth and reported the maximum value for each calculation grid.Since earthquakes with small magnitudes do not significantly influence the current seismicity rate within the model (Catalli et al., 2008), I only analysed the CFS for the M ≥ 4.5 events (Table 1).The CFS evaluation is based on the assumption of a homogeneous half space with an intermediate value of µ = 0.4, the Poisson's ratio of 0.24 and Young's modulus of 8 × 10 5 by applying the COULOMB 3.3 program (Toda and Stein, 2002).Application of the rate-andstate friction model (Eq.5) requires parameters of Aσ and t a .Previous studies (e.g.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.I assumed a fixed Aσ of 0.2 bars, corresponding to the assumption of previous studies (e.g.Chan et al., 2012Chan et al., , 2013)).t a was assumed to be a function of the moment magnitude (M w ), as proposed by Burkhard and Grünthal (2009), described as follows: The unit of t a is days.t a is determined based on the magnitude of each source event (Table 1).

The two-dimensional models
I first represent the forecasting models based on twodimensional calculation cells with a 0.1 • × 0.1 • size (i.e. the depth-independent model).For application of the smoothing kernel function, x − x i in Eq. ( 1) was the epicentre distance between the site of interest and the epicentre of earthquakes (i.e.depth-independent).The model forecasted higher rates along the coastline of Taiwan and for the area east of latitude 122.5 • , which correspond to the distribution of the forecasting events during 2010 and 2011 (Fig. 2a).
For the CFS calculation on the two-dimensional grids, the target depth corresponds to the hypocentral depth of each source event (Table 1).Through the rate-and-state friction model, I calculated the time-dependent rate evolutions for different moments (Fig. 3).In comparison with the spatialtemporal pattern of the forecasting events (open circles in Fig. 3), many of the consequent earthquakes are in the region with rate decrease (green stars in Fig. 3); i.e. it is difficult to confirm the feasibility of this model.

The three-dimensional models
I then propose the forecasting models based on threedimensional cells with 0.1 • × 0.1 • × 10 km sizes (i.e. the depth-dependent model).For the smoothing kernel function application, x − x i in Eq. ( 1) was the distance between the site of interest and the hypocentre of earthquakes (i.e.depthdependent).Two profiles along the longitudes of 122.0 and 122.5 • (Fig. 2b and c, respectively) presented higher forecasted rates above the depth of 30 km and along the subduction slab dipping to the north, which fit the distribution of the forecasting earthquakes (the open circles in Fig. 2b and c) well.
For the rate-and-state friction model application, I evaluated the maximum CFS along the seismogenic depth for each cell and modelled the corresponding seismicity rate evolution (Fig. 4).Departing from the outcomes of the twodimensional model (Fig. 3), a significant rate increase near the epicentre of each source event corresponds to the distribution of forecasting events (Fig. 4).
4 Forecasting application to the Kanto region

Tectonic setting and earthquake catalogue
The Kanto region, Japan, is an area with complex tectonic settings.Most parts of this region sit on the Eurasian Plate, under which the Philippine Sea Plate subducts from the south.At greater depth, the Pacific Plate subducts from the east (Toda et al., 2008).The complex plate interactions in this region result in seismic activity.Fortunately, similarly to the Ryukyu region, the activity has been well recorded by a high-quality seismic network.The modern seismic network maintained by the Japan Meteorological Agency (JMA) was initiated in 1923 and has been modernized over time (Nanjo et al., 2010).In addition, significant change in the seismicity behaviour in the Kanto region followed the 2011 M9.0 Tohoku earthquake (Ishibe et al., 2011;Toda et al., 2011).Such spatial-temporal conditions provide an ideal environment for testing the credibility of forecasting models in respect to both depth and time factors.

Procedure of application
Due to the seismic network modernization, M c of the JMA catalogue decreased dramatically after 1980 and 1990 (Nanjo et al., 2010).Thus, I analysed the catalogue with various magnitude thresholds in the three periods: magnitudes 4.5, 3.5 and 2.5 since 1923, 1980 and 1990, respectively.The thresholds correspond to the M c determined by Nanjo et al. (2010) through the entire magnitude range method (Woessner and Wiemer, 2005).I input the complete part of sion and magnitude of slip through the scaling laws proposed by Wells and Coppersmith (1994).A receiver fault plane for each grid is consistent with the closest reference focal mechanism from the F-net catalogue.

The two-dimensional model
This study first forecasts in two-dimensional calculation cells with 0.2 • × 0.2 • sizes defined by the Collaboratory for the Study of Earthquake Predictability (CSEP) Japan Testing Center for the Kanto region (Tsuruoka et al., 2012).The target depth for the CFS calculation is 47.5 km, which corresponds to the averaged hypocentral depth of the earthquakes in the region.Combining the smoothing kernel function and the rate-and-state friction model, the models represent higher seismicity rates for smaller magnitude ranges (e.g.Fig. 6a) than for larger ones (e.g.Fig. 6d), consistent with the Gutenberg-Richter law (Gutenberg and Richter, 1954).In addition, due to being stress-enhanced by the Tohoku sequence (Ishibe et al., 2011;Toda et al., 2011), a high seismicity rate is forecasted along the coast and offshore in the Pacific Ocean 40 km northeast of Tokyo.Note that the 5.0≤M≤5.96.0≤M≤6.97.0≤M≤7.94.0≤M≤4.95.0≤M≤5.96.0≤M≤6.97.0≤M≤7.9stress shadow zone at the target depth by the source events (including the 2011 Tohoku mainshock) resulted in some low forecasted rate zones offshore in the Pacific Ocean.

The three-dimensional model
I then proposed a depth-dependent model based on threedimensional cells with 0.2 • × 0.2 • × 5 km sizes.Compared with the forecast with two-dimensional cells, the threedimensional model illustrated more detailed patterns along the depth (Fig. 7b-d); e.g. the high seismicity rate 40 km northeast of Tokyo is identified at depths between 30 and 70 km (Fig. 7b), and the high rate along the coast and offshore in the Pacific Ocean is located at depths of 25-75 km (Fig. 7d), consistent with the boundary between the Pacific and Eurasia plates (Toda et al., 2008).

Forecasting ability of each forecasting model
To validate the forecasting ability statistically, I compared models with the distribution of forecasting earthquakes using the Molchan diagram (Molchan, 1990, and   therein).The diagram was designed for evaluating forecasting ability through presenting the fraction of alarm-occupied space vs. the fraction of failure in forecasting by considering the locations of the forecasting earthquakes with respect to the distribution of forecasting seismicity density rate.The "fraction of alarm-occupied space" indicates the percentage of events within the study region with a forecasting level equal to or higher than "alarm".The "fraction of failure in forecasting" represents the percentage of consequent earthquakes with a lower forecasting level than the alarm, corresponding to "miss rate" defined by some previous studies (e.g.Zechar and Jordan, 2008).For each event, the area with a forecasting rate equal to or smaller than that at the location of the forecasting earthquake is extracted and represented as a percentage of the entire study area.The events are then sorted according to percentage of area and plotted against event count, represented as the percentage of the total number of forecasting events.In the diagram, when data points are distributed along a diagonal line, the distribution of target earthquakes is independent of forecasting; convex distribution suggests that the majority of consequent earthquakes occur within regions with a lower forecasted rate, whereas concavity suggests that the majority of consequent earthquakes are within a high forecasted rate area.An optimistic forecasting is represented in the Molchan diagram by the condition of having the lowest fraction of alarm-occupied space, and the lowest fraction of failure in forecasting.
Through Molchan diagrams, the forecasted seismicity rate obtained using different models was compared with the locations of earthquakes for the Ryukyu and Kanto cases (shown in Figs. 8 and 9, respectively).Most of the models show concavity distribution, suggesting good forecasting ability, except the case of the combining two models in   [2010][2011] for the Ryukyu case.Blue and yellow dots represent the results for the models using the smoothing kernel models in two-and three-dimensional grids, respectively; red and green dots represent the results for the models using the rate-and-state friction models in two-and three-dimensional grids, respectively; grey dots represent the 99 % significance level determined by 1640 forecasting events.two-dimensional grids in Kanto (the yellow dots in Fig. 9).Such an exception can be attributed to slip model misfit and hypocentral depth uncertainties of forecasting earthquakes (Catalli and Chan, 2012).To further confirm the significance of the forecasting ability for the rest of models, the null hypothesis (Zechar and Jordan, 2008) is implemented.The 99 % significance level, i.e. α = 1 % in Eq. (3) of Zechar and Jordan (2008) for both Ryukyu and Kanto cases (the grey dots in Figs. 8 and 9, respectively), is plotted based on the number of the forecasted earthquakes (1640 and 703 earthquakes, respectively).Most of the models cannot be rejected by the 99 % confidence level of the null hypothesis (the data below the grey dots in Figs. 8 and 9), confirming their robustness.

Importance of the temporal factor
Since the smoothing kernel function averages the seismic activity during the observation period, it can be regarded as a time-independent model.In contrast, the rate-and-state friction model forecasts temporal evolution of seismicity rate disturbed by a series of source events and can be renewed with time.The comparison between the two models may indicate the importance of the temporal factor for forecasting.
The Molchan diagram in the Ryukyu case shows a lower fraction of failure in forecasting for the smoothing kernel model using the two-dimensional grids (the blue dots in Fig. 8) than that for the rate-and-state friction model (the red dots in Fig. 8) when fraction of alarm-occupied space is fixed.Such a result confirms a better forecasting ability for the smoothing kernel model.A similar conclusion can be obtained for the three-dimensional models (the yellow and green dotes in Fig. 8).This finding corresponds to the conclusion of Chan et al. (2012), obtained from forecasting experience in entire Taiwan.For the Kanto case (Fig. 9), by contrast, departing from conclusion of the Ryukyu case, the rate-and-state friction model provides a better performance.
The discrepant conclusions between the two cases might be attributed to the effect of recent earthquakes.For the cases of Ryukyu and Chan et al. (2012), there was no significant short-term rate perturbation during the forecasting periods.For the Kanto case, on the contrary, the time dependency becomes a crucial factor in forecasting the consequences of the 2011 Tohoku earthquake (Ishibe et al., 2011;Toda et al., 2011).

Importance of the depth factor
Both of the Ryukyu and Kanto cases have qualitatively shown that three-dimensional models have a better performance in Figs. 2 and 7, respectively.The smoothing kernel function using three-dimensional grids forecasted distribution along depth in detail.For the rate-and-state friction model, the three-dimensional applications presented significant rate increase for most regions near the epicentre of each source event, corresponding to the locations of forecasting events.
For the Ryukyu application, comparing between forecasted rates obtained using the smoothing kernel function and the locations of target earthquakes using the Molchan diagram (blue and yellow dots in Fig. 8), the three-dimensional forecasting model had a better forecasting ability, i.e. a smaller fraction of failure to predict.For the rate-and-state friction model, the three-dimensional applications grids also provide a better forecasting ability (green dots in Fig. 8) in comparison to the two-dimensional ones (red dots in Fig. 8).
The application to the Kanto region also confirmed the better performance of the three-dimensional models (yellow and green dots in Fig. 9).In addition, in the Kanto case, the Molchan diagram raised the disadvantage for the rate-andstate model using two-dimensional calculation grids.When the fraction of alarm-occupied space is large, convex distributions are presented (red and yellow dots in Fig. 9); i.e. some earthquakes took place in the region with low/no forecasted rates (region in white in Fig. 6), suggesting forecasting failure.In contrast to the two-dimensional models, the threedimensional ones have proved their forecasting ability (green dots in Fig. 9).The conclusion is consistent with the findings of Catalli and Chan (2012), and confirms that the depth fac-tor is one of the upmost important parameters for Coulomb stress calculation.
Through the applications to different forecasting approaches, I have confirmed that models with threedimensional grids always obtain better forecasting ability.In addition, the importance of depth dependency for forecasting models is determined, especially for the application to a subduction environment or to a region with non-vertical seismogenic structures.

Data availability
In the Ryukyu case, the earthquake parameters were obtained from the catalogue of the Central Weather Bureau Seismic Network (last accessed December 2013); the focal mechanisms parameters were obtained from Broadband Array in Taiwan for Seismology (BATS) (last accessed December 2013).In the Kanto case, the earthquake parameters were obtained from the catalogue of the Japan Meteorological Agency (JMA) network (last accessed April 2012); the focal mechanisms parameters were obtained from the F-net catalogue (last accessed April 2012).The detailed slip dislocation model of the 2011 M9.0 Tohoku earthquake was derived from Fujii et al. (2011) (http://iisee.kenken.go.jp/staff/fujii/OffTohokuPacific2011/tsunami_inv.html).The Coulomb stress change is calculated from COULOMB 3.3 (http://usgsprojects.org/coulomb/).

Figure 1 .
Figure 1.The magnitude of completeness (M c ) for (a) the TTSN and (b) the CWBSN catalogues.(c) The study area is shown in grey, with the intersection of regions with M c ≤ 4.0 for the TTSN shown with dashed lines and with M c ≤ 3.0 for the CWBSN shown with solid lines.The M c for each grid is determined according to the events that occurred within 30 km of the centre of each grid.
Figure 2. (a) A map view and (b-c) profiles of the forecasted seismicity rate for M ≥ 3.0 modelled by the epicentre-smoothing kernel function.White dots denote the earthquakes from 2010 to 2011.Earthquakes within 25 km of each side of the profiles are presented.

Figure 3 .
Figure 3.The seismicity rate evolution at different time moments.The target depth for the CFS calculation corresponds to the hypocentral depth of each source event (Table 1).Source events from 2010 to 2011 are shown as open black stars.Open circles denote earthquakes during each time sequence.

TheFigure 5 .
Figure 5. Distribution of the M ≥ 6.0 earthquakes, which took place in or near the Kanto region during January of 2010 and August of 2011.The coseismic dislocation model of the 2011 M9.0 Tohoku earthquake is obtained by tsunami waveform inversion (Fujii et al., 2011), whereas the focal mechanisms of the others are obtained from the F-net catalogue.

Figure 6 .
Figure 6.The two-dimensional forecasting models for the magnitudes between (a) 4.0 and 4.9; (b) 5.0 and 5.9; (c) 6.0 and 6.9; and (d) 7.0 and 7.9, respectively, at the end of August 2011.Black dots denote the M ≥ 4.0 earthquakes during January 2010 and August 2011.
Figure 7. (a) Map view and (b)-(d) profiles of the three-dimensional forecasting models at the end of August 2011 and the distribution of the target earthquakes during 2010-2011.Black dots represent the M ≥ 4.0 earthquakes during January 2010 and August 2011.Earthquakes within 5 km of each side of the profiles are presented.

Figure 8 .
Figure8.The Molchan diagram used for investigating the correlation between different forecasting models and earthquakes during the forecasting period (2010-2011) for the Ryukyu case.Blue and yellow dots represent the results for the models using the smoothing kernel models in two-and three-dimensional grids, respectively; red and green dots represent the results for the models using the rate-and-state friction models in two-and three-dimensional grids, respectively; grey dots represent the 99 % significance level determined by 1640 forecasting events.

Figure 9 .
Figure9.The Molchan diagram used for investigating the correlation between different forecasting models and earthquakes during the forecasting period (January 2010-August 2011) for the Kanto case.Blue and red dots denote the results for the models using solely the smoothing kernel function and the rate-and-state friction model, respectively; yellow and green dots denote the results for the combination of the two models in two-and three-dimensional grids, respectively; grey dots represent the 99 % significance level determined by 703 forecasting events.

Table 1 .
Source parameters for the source events used for the inputs of the rate-and-state friction model.Earthquakes with M w ≥ 4.5 that occurred in 2010 and 2011 were considered.Parameters were determined based on the Broadband Array in Taiwan for Seismology (BATS).