A Forecasting Procedure for Plate Boundary Earthquakes Based on Sequential Data Assimilation

A forecasting procedure is proposed for plate boundary earthquakes in subduction zones. It is based on spatio-temporal variation in slip velocity on the plate interface, which causes interplate earthquakes. Model outputs are not only information about the occurrence of great earthquakes but also information about the physical state evolution that causes earthquakes. To overcome the difficulty in forecasting earthquake generation resulting from uncertainty both in the physical model and in the observation data, we introduce a type of sequential assimilation for crustal deformation data. We are currently constructing a prototype, applying to the plate boundary along the Nankai Trough.

atmosphere is monitored continuously so that its evolution can be forecast.Both the atmosphere's monitored (current) and forecast states are important outputs in weather forecasting.Our ultimate goal in earthquake forecasting is to provide outputs similar to those in weather forecasting-monitored and forecast physical states that cause great interplate earthquakes.In this paper, we describe a procedure for forecasting earthquake generation based on incorporation of real-time observations into a physical model.
Mathematical models for earthquake generation that correspond to atmospheric models exist only for interplate earthquakes (e.g., Rice, 1993;Lapusta et al., 2000;Hashimoto et al., 2013).These simple models can, at least qualitatively, reproduce a range of slip behaviors on the plate interface, including those that generate great interplate earthquakes (e.g., Kodaira et al., 2006;Kato, 2008;Nakata et al., 2012;Noda et al., 2013).In these models, earthquakes are represented as unstable frictional slip on a plate interface.Strain energy that will be released during an earthquake is stored around the coupled area of the interface, where the slip velocity of the interface is much slower than the velocity of plate motion.In the model, we use a simple law of friction derived from rock experiments because we cannot measure the friction at the plate interface.Because of such limitations in the model and also limitations in observed data, we cannot forecast the slip behavior in a deterministic way.
We can, however, use physical models and run as many different scenarios as possible to avoid overlooking any risks.
Despite these limitations, our model is worth pursuing because improvements can be introduced easily, and soon we will be able to include more real-time observations, especially offshore, where great interplate earthquakes nucleate.
In this article, we show that physical state evolution that causes earthquakes on a plate interface can be modeled as spatio-temporal variation in slip velocity along the interface.We then introduce a simple mathematical model for plate-interface earthquakes and explain our forecasting procedure.Finally, we introduce a model prototype and discuss remaining issues.

EARTHQUAKE GENER ATION MODEL FOR PL ATE BOUNDARIES
Earthquake generation can be represented in a model as a repetitive sequence of stick-slip events on a plate interface.
The coupled zone is loaded by the velocity difference between the area that is stuck and surrounding areas where slow sliding occurs."Slow" means that the sliding velocity is an order of magnitude smaller than that of seismic slip, which causes elastic wave radiation.Our model

INTRODUCTION
Earthquake forecasting or prediction attempts to estimate when, where, and how large earthquakes will occur.
Long-term earthquake forecasts provide probabilities of earthquake occurrence within a few decades.This information is used, for example, in seismic hazard assessment for building codes, insurance rate calculations, and policymaking.However, short-term forecasts-whether an earthquake will occur within several years or less-are still at the trial stage (e.g., Jordan, 2006).Short-term deterministic predictions have been tried, but they have not been successful (e.g., Geller, 1997;Geller et al., 1997;Kagan 1997).Note that the outputs of most earthquake predictions and forecasts, whether based on statistics or physics, are information on earthquake occurrence, such as occurrence time, magnitude (M), hypocenter location, frequency magnitude distribution, and recurrence time probability density functions (e.g., Tullis et al., 2012).
In contrast, weather forecasts provide not only the probability of rainfall but also the state of the atmosphere that causes rainfall.The state of the ABSTR ACT.A forecasting procedure is proposed for plate boundary earthquakes in subduction zones.It is based on spatio-temporal variation in slip velocity on the plate interface, which causes interplate earthquakes.Model outputs are not only information about the occurrence of great earthquakes (time, place, and magnitude) but also information about the physical state evolution that causes earthquakes.
To overcome the difficulty in forecasting earthquake generation resulting from uncertainty both in the physical model and in the observation data, we introduce a type of sequential data assimilation.In this method, we compare observed crustal deformation data to simulations of several great interplate earthquake generation cycles.We are currently constructing a prototype, applying this forecasting procedure to the Nankai Trough, Southwest Japan, where great interplate earthquakes have occurred and are anticipated.Geodetic data can confirm that coupled areas on plate interfaces correspond to source areas where M>7 earthquakes have generated tsunamis (Hashimoto et al., 2009(Hashimoto et al., , 2012)) Peaks in slip deficit that coincide with the source regions of large, historical interplate earthquakes are observed (Figure 1a; Hashimoto et al., 2012).More importantly, from these data it is possible to determine the extent of the coupled area that is consistent with the large slip area of the 2011 Tōhoku-Oki earthquake (Figure 1b; Hashimoto et al., 2012), suggesting that the source areas of future M9 great earthquakes could be estimated using observational data.For example, Figure 1a indicates that a wide coupled area exists along the Kuril Trench off northernmost Japan, which corresponds to a past tsunami source area as large as that of the 2011 Tōhoku-Oki earthquake, whose recurrence is anticipated based on geological data analyses (Nanayama et al., 2003;Sawai et al., 2009).
Uncertainties remain about the size of the coupled area, especially offshore near the trench if only on-land data are available (e.g., Nishimura et al., 2004).
The relationship between the coupled area and potential earthquakes is also uncertain because large amounts of slip can occur in slow sliding areas during interseismic periods (Noda and Lapusta, 2013).Note that the uncertain coupled area corresponds to the relatively large fault slip area during the 2011 Tōhoku-Oki earthquake (e.g., Fujiwara et al., 2011).
During interseismic periods, afterslip along faults triggered by the main earthquakes (Heki et al., 1997) and slow slip events (SSEs; Hirose et al., 1999) have been observed.An SSE denotes spontaneous slow sliding on the plate interface.
A few years before the 2011 Tōhoku-Oki earthquake, SSEs were observed near the Tōhoku-Oki hypocenter (Ito et al., 2012), and M7 earthquakes with large afterslips were observed within the source area (Suito et al., 2011).SSEs and large afterslips may be preseismic signals that could have indicated that the 2011 Tōhoku-Oki earthquake was about to occur, although it is also possible that such phenomena occur without being followed by a larger earthquake (Ohtani et al., 2014).The SSEs occurred within the subseafloor source area of the 2011 Tōhoku-Oki earthquake so that GNSS observations couldn't detect them.They Trench (Kaneda, 2014, in

MATHEMATICAL MODEL OF SLIP VARIATION AT THE PL ATE INTERFACE
We need to model the geometry of the plate interface in order to model stick-slip behavior there.The interface is divided into small subfaults.The stress transferred to each small subfault from slip on other subfaults is calculated using a Green's function for an elastic half-space medium (Okada, 1992).
The interface is assumed to obey a simple law of friction derived from rock experiments.This law of friction describes the relationships between slip velocity, shear stress in the slip direction, and fault strength (Nakatani, 2001).If the shear stress is less than the fault strength, the slip velocity is much smaller than the relative plate motion.
However, if the stress is higher than the fault strength, the velocity can reach seismic slip rates as high as 1 m s -1 .
Fault strength is not constant but rather evolves, depending on slip history.Fault strength declines when slip occurs (slip weakening) and recovers when the slip rate is low (logarithmic healing).Laws for fault strength evolution have been proposed (e.g., Dieterich, 1979;Ruina, 1983;Kato and Tullis, 2001).Stick-slip and slow sliding areas can be identified from the heterogeneous distribution of friction law parameters.For the stick-slip area, the rate of fault weakening is high and fault strength recovery time is short.
Slow sliding areas have the reverse properties.Recent experiments indicate that fault strength becomes very low when the sliding velocity reaches the order of 1 m s -1 (e.g., Wibberley et al., 2008).
When the mathematical model qualitatively (Kodaira et al., 2006;Kato, 2008).The model of an earthquake generation cycle for an M9 earthquake, with preceding earthquakes of M7~M8 within the source area, is qualitatively consistent with the characteristics of the 2011 Tōhoku-Oki earthquake (Hori and Miyazaki, 2011;Ohtani et al., 2014).In addition to earthquake generation, slow sliding behavior on the plate interface has also been reproduced qualitatively using this mathematical model.For example, the coexistence of afterslip and slow slip events for an M7 earthquake was reproduced for the Hyuga-nada region in Southwest Japan (Nakata et al., 2012).The above simulations are based on quasi-dynamic models in which elastic wave propagation is not included.
Although fully dynamic models are more realistic (e.g., Lapusta et al., 2000;Noda and Lapusta, 2013;Noda et al., 2013), especially for rupture propagation, fully dynamic simulations with realistic plate geometry require very large computational resources, so their introduction remains to be tackled in the future.

A NEW EARTHQUAKE GENER ATION FORECASTING SYSTEM CONCEPT
As described above, we can quan-

PROTOT YPE FORECAST PROCEDURE FOR NANK AI TROUGH EARTHQUAKES
We are now constructing observation and simulation databases for the Nankai Trough, where a great earthquake is anticipated (Hori et al., 2013).For   Figure 4 shows how we changed the model parameter distributions for the various scenarios.We changed (1) the deeper limit of the seismogenic slip area (seismogenic zone) and its strength for 56 scenarios (Figure 4a), (2) the upper limit of the seismogenic zone for 30 scenarios (Figure 4b), and (3) the strength of the barrier for slip propagation between eastern and western areas of the Kii Peninsula for 30 scenarios (Figure 4c).An SSE patch is included for two scenarios (Figure 4d).
The resultant recurrence time intervals of great earthquakes for the different scenarios were from 50 to 210 years, with an average of 115 years, which is comparable to historical data on recurrence of M8 class earthquakes along the Nankai Trough (Ishibashi, 2004).
From the scenarios of spatio-temporal

CONCLUSIONS
The forecasting procedure presented here is a prototype, but we consider it important to demonstrate a test forecast using it.Although earthquake predictability is low, monitoring and continuous communication with the public about slip variations on plate interfaces and their possible results may improve public scientific literacy with regard to earthquake generation.Further, testing our procedure will identify issues that need to be solved for forecasting and public communication of earthquake information.We need to improve our S P E C I A L I S S U E O N U N D E R S E A N AT U R A L H A Z A R D SA Forecasting Procedure for Plate Boundary Earthquakes Based on Sequential Data AssimilationB Y TA K A N E H O R I , M A M O R U H Y O D O , R Y O K O N A K ATA , S H I N ' I C H I M I YA Z A K I , A N D Y O S H I Y U K I K A N E D A Oceanography | Vol. 27, No. 2 94A snapshot of slip velocity (V) distribution for M8 earthquake occurrence in an earthquake generation cycle simulation along the Nankai Trough, off Southwest Japan.Red, yellow, and blue colors indicate high-speed slip (seismic), moderate-speed slip (comparable to the plate motion = V pl ), and slow speed slip (much slower than the plate motion), respectively.

Figure 1 .
Figure 1.(a) Slip-deficit zones estimated from interseismic global navigation satellite systems (GNSS) data (from Hashimoto et al., 2012).The blue and red contours represent, respectively, interseismic slip-deficit and slip-excess rates at intervals of 3 cm yr -1 .The green circles and yellow stars indicate tsunami source regions and epicenters, respectively, of large interplate earthquakes that occurred in the last century.(b) Coseismic slip distribution estimated from inversion of GEONET (GNSS Earth Observation Network System) data (from Hashimoto et al., 2012).The green contours represent the inverted coseismic slip distribution on the plate interface at intervals of 4 m.For comparison, the blue contours show the interseismic slip-deficit rate distribution at intervals of 3 cm yr -1 .
were revealed during processing of ocean bottom pressure gauge data recovered after the earthquake.Ocean bottom cabled observatories that are now operating off Kumano in the eastern part of the Nankai Trough (Kaneda, 2010) and being constructed in the western part of the Nankai Trough and in the Japan described above is applied to slip behavior at subduction zone plate interfaces, it successfully reproduces various earthquake occurrence patterns and slow sliding behaviors.For example, variations in the size of and the time intervals between significant earthquakes along the Nankai Trough, offshore Southwest Japan, and along the northern part of the Japan Trench have been reproduced

Figure 2 .
Figure 2. Schematic diagram of the earthquake forecasting procedure proposed in this study.
database, we show the preliminary results of 118 scenarios whose parameter variations do not yet cover an appropriate range of uncertainty.Each scenario corresponds to one set of model parameters.To calculate one scenario, it takes about 12 hours using 128 nodes of the K computer.

Figure
Figure 3a shows the standard case: V pl is the plate motion rate, and L, A and B are frictional parameters.L determines the slip amount necessary for fault weakening.A-B control the stability of steady sliding.Figure 3b shows some examples of earthquake slip distribution.Note that the moment magnitude can reach 9.0 in the simulated scenarios, although interpretations of historical data indicate only M8 earthquakes occurred (Ishibashi, 2004).See Hyodo and Hori (2013) for a more detailed description of model parameter distribution and simulated results of earthquakes along the Nankai Trough.
coseismic stress drop and lower limit (blue arrows) of seismogenic zone Strength of the barrier (blue ellipse) assumed at Kii peninsula Upper limit of the seismogenic zone Inclusion of SSE patch (red ellipse) in the Bungo channel

Figure 4 .
Figure 4. Classification of scenarios based on the heterogeneity of frictional properties.

Figure 6
Figure 6 shows temporal variations in vertical displacement at Kushimoto and Muroto for GEONET data and simulated scenarios in which the initial values are optimally fit to the data.Each of 10 scenarios with similar likelihood

Figure
Figure 5. (a) Vertical (above) and horizontal (below) velocity distribution of observation data (black), the best scenario (red), and the worst scenario (blue).(b) Slip velocity distribution on the plate interface for the best case scenario (above) and the worst case scenario (below).Blue indicates a coupled area.A large difference can be seen beneath the Kii Peninsula, indicated by a red ellipse.

Figure
Figure6.Time variation of observation data (dots) and simulated scenarios with their likelihoods indicated by color at the Kushimoto and Muroto stations (located in Figure5).One of ten scenarios was picked from among those with similar likelihoods.Abrupt deformation indicates the timing of a great earthquake in each scenario.