An Operational Earthquake Forecasting Experiment for Israel: Preliminary Results

Operational Earthquake Forecasting (OEF) aims to deliver timely and reliable forecasts that may help to mitigate seismic risk during earthquake sequences. In this paper, we build the first OEF system for the State of Israel, and we evaluate its reliability. This first version of the OEF system is composed of one forecasting model, which is based on a stochastic clustering Epidemic Type Earthquake Sequence (ETES) model. For every day of the forecasting time period, January 1, 2016 - November 15, 2020, the OEF-Israel system produces a weekly forecast for target earthquakes with local magnitudes greater than 4.0 and 5.5 in the entire State of Israel. Specifically, it provides space-time-dependent seismic maps of the weekly probabilities, obtained by using a fixed set of the model’s parameters, which are estimated through the maximum likelihood technique based on a learning period of about 32 years (1983–2015). According to the guidance proposed by the Collaboratory for the Study of Earthquake Predictability (CSEP), we also perform the N- and S-statistical tests to verify the reliability of the forecasts. Results show that the OEF system forecasts a number of events comparable to the observed one, and also captures quite well the spatial distribution of the real catalog with the exception of two target events that occurred in low seismicity regions.


INTRODUCTION
The frequent occurrence of deadly earthquake events highlights the importance of delivering reliable and skillful earthquake forecasts over different time windows (from days to decades) to support rational actions of risk reductions and enhancing preparedness and resilience. For this purpose, the International Commission for Earthquake Forecasting for Civil Protection, nominated by the Italian government after the M w 6.1 earthquake occurred in L'Aquila (Italy) on April 6, 2009 (Thomas et al., 2011(Thomas et al., , 2014, recommended the development of an Operational Earthquake Forecasting (OEF) system, which comprises procedures for gathering and disseminating authoritative information about the time dependence of seismic hazards, in order to help communities prepare for potentially destructive earthquakes. Specifically, OEF provides timely earthquake (probabilistic) forecasts over time-spaceintensity windows of interest for stakeholders (e.g., government agencies and departments).
OEF systems have already been implemented in several forecasting applications worldwide, such as in Italy, New Zealand and United States Marzocchi et al., 2014;Michael et al., 2020). In principle, OEF should deliver a continuous flow of information to avoid violating the so-called hazard/risk separation principle (Jordan et al., 2014). So far, this feature is in place only in the Italian OEF system, whereas in most cases OEF information is released upon specific subjective requests. In the wake of the Italian experiment, in this paper we describe the first version of the OEF system for the State of Israel, which is an active seismic region that experienced strong earthquakes in the past (e.g., M w 6.3 on 1927/07/11 in Jericho, M w 5.3 on 2004/02/11 in Israel, M w 5.1 on 2008/02/15 in south Lebanon). As part of the Middle East, Israel is in a thorny position, embedded between the four major tectonic plates: Nubia (Africa), Sinai, Arabia and Anatolia (see Figure 1). The collision between the Africa-Arabian plate with the European-Asian one resulted in faulting and folding of the sedimentary strata in Israel, intensified also by the fault zone developing along the Dead Sea-Red Sea Rift Valley.
Developing, implementing and delivering an OEF system for the State of Israel is indeed the final scope of the bilateral Italy-Israel project "Enhancing OPerational Earthquake foRecasting through innovative Approaches (OPERA)", on request of the Israeli steering committee for earthquake preparedness, in collaboration with the Geophysical Institute of Israel. The scope is to provide the very first real-time operating tool, which brings the State of Israel into alignment with several other countries' direction towards the release of reliable, live earthquake forecasts. Although the OEF system presented here is a prototype, it is an important first step which we trust can help the Israeli government agency and population to take pondered and calibrated actions aimed at reducing seismic hazard.
The first version of the OEF system presented here is based on the stochastic Epidemic Type Earthquake Sequence (ETES) model (Console and Murru, 2001;Giuseppe Falcone et al., 2010), which belongs to the general class of the self-exciting Epidemic Type Aftershock Sequence (ETAS) models firstly introduced by Ogata in 1988 (Ogata, 1988;Ogata, 1989;Ogata, 1998). The basic rational behind these models relies on the idea that seismic events cluster in space and time. In fact, this kind of spatiotemporal earthquake clustering models has been found to be the most reliable to track the probabilistic evolution of the seismic process (Marzocchi et al., 2017;Taroni et al., 2018). In this paper, we set up the ETES model specific for the Israeli seismicity, and we also check the consistency of the forecasts with the data contained in the Israeli seismic catalog (https:// earthquake.co.il/en/earthquake/searchEQS.php) carrying out two statistical tests proposed by the Collaboratory for the Study of Earthquake Predictability (CSEP), which is an international organization aimed at evaluating quantitatively earthquake predictions and forecasts at a global scale (https:// cseptesting.org/) .
The evaluation of the forecasts produced can help highlighting advantages and room for improvement of the OEF-Israel system, as well as giving possible hints to adjust the same experiment in other countries.

ISRAELI TECTONIC SETTING
The Dead Sea Transform (DST), known also as the Dead Sea Fault, is the largest fault system of the Eastern-Mediterranean area, place of the strongest earthquakes in the last centuries (Meghraoui et al., 2003). It has been formed during the Miocene due to the African-Arabian plate collision with the Gulf of Aden southward of the Red Sea, which caused the separation of the Sinai subplate. The Israeli region is located in the southern section of the DST fault system (see Figure 1), but it is affected also by the seismic activity occurring in the middle DST, that is a restraining bend characterized by transpression deformation (Quennell, 1959). The region consists of a topographic valley with raised flanks and normal fault boarders. Tectonic movements occur on oblique-and left strike-slip fault segments, delineating a string of rhombshaped, narrow and deep sectors releasing bends linked to orthogonal separation of the transform flanks on the surface, that could also extend well beneath the crust (Garfunkel and Ben-Avraham, 2001;Wetzler et al., 2014). This is the most active tectonic region in the Eastern-Mediterranean area, as witnessed by the several prehistorical, historical and more recent large Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 729282 earthquakes (Amit et al., 2002;Baer et al., 2008;Marco, 2008). A spatial distribution analysis of seismic events (Sharon et al., 2019) allowed to identify criteria for sorting the active faults in the region as "main strike-slip faults of the DST" and "faults with direct evidence of Quaternary activity," which are shown in the map available at the following link (see also Fig7 of Sharon et al. (2019)). As shown in the next subsection, and by comparing Figure 1 and the map of the url above, the complete dataset considered in this paper involves the main strike-slip faults of the DST and surrounding areas, as well as marginal faults, main branches and quaternary activity. The segments involved are Dakar fault, passing through Arava and Eastern/Western Dead Sea faults, until the further north Jordan Gorge and Roums faults.

The Philosophy of the OEF-Israel System
We adopt the same philosophy of the OEF-Italy system which is rooted in three main pillars. In particular, the OEF-Israel system 1) is fully transparent, such that the results are reproducible by anyone, expert or not in the field, who is interested in knowing and understanding the earthquake probabilistic forecast in the spatio-temporal window of interest; 2) is open to any modelers who want to contribute for future developments (the Italian system requires that any additional model has to be submitted to at least one CSEP experiment); 3) minimizes scientific controversies. To be more specific on this latter point, the idea is to use ensemble models built as the weighted combination of two or more models, where the weights are obtained from the objective statistical evaluation of each of the constitutive models' performance, thus significantly reducing the implicit subjectivity induced by selecting a single model. An ensemble model is not yet implemented in this pilot OEF-Israel system because we consider here only one model, the ETES model; however, this functionality can be activated when additional models will be added.

ETES Model
The earthquake forecasting model adopted in the OEF system for Israel is the clustering, self-exciting Epidemic Type Earthquake Sequence (ETES) model proposed in Giuseppe Falcone et al. (2010) and previously introduced in Console and Murru (2001) and Console et al. (2007). According to this purely stochastic model, the conditional probability density function quantifying the occurrence rate of a seismic event (x, y, t, M) is given by where: H t is the past history of the current event; f R is the "failure rate" function, i.e., the ratio between the expected number of independent events and the total number of shocks; λ 0 (x, y, M) is the time-invariant "spontaneous" background seismicity, obtained by analyzing a real seismic catalog using a smoothing algorithm; H (·) is a step function and λ i (x, y, t, M) is the punctual contribution of the past ith event (x i , y i , t i , M i ). The summation performed over the contributions of all the previous shocks, that is the second term in the right hand side of the ETES rate (1), provides the time-varying "triggered" aftershock component of the seismic sequence. According to the ETES model, the magnitude distribution of both the spontaneous and triggered components of rate (1) is separable with respect to time and space, and it is given by the well-known decreasing exponential Gutenberg-Richter law (Gutenberg and Richter, 1944) where M c is the completeness magnitude, that is the value such that all the events with a higher magnitude are surely recorded in the earthquake catalog (Mignan and Woessner, 2012). The parameter b is the so-called b-value, which controls the slope of the magnitude-frequency distribution in a semi-log plot: the higher is this parameter, the lower is the fraction of larger shocks with respect to the smaller ones; some authors claim that a decrease of the b-value can be interpreted as a precursor for forthcoming large shocks, but this is actually an open debate in the literature due to the potential sources of bias that can affect this parameter's estimation (Marzocchi et al., 2020). The timeindependent background function λ 0 (x, y, M) is obtained through the method introduced by Frankel (1995), with an exponential kernel distribution adopted for the smoothing algorithm. The triggered component λ i (x, y, t, M) of rate (1) is given by where p is the modified Omori law for the temporal decay, and is the isotropic spatial function of the distance r between the location of the current event (x, y) and the location of the previous event (x i , y i ). Following Giuseppe Falcone et al. (2010), in the function d i we set α 1 2 , so that the average triggering distance of the aftershock zone is proportional to the square root of the mainshock rupture area, as observed for real data (Kagan, 1991). We note that when α is not fixed a priori, but it is estimated, and when the spatial function is Gaussian, the model is switched to the classical ETAS (Ogata, 1998).
The ETES parameters (p, c, k, d 0 , q, b) are assumed to be positive and typically estimated through the maximum likelihood function technique (Aki, 1965;Bender, 1983;Ogata, 1988;Ogata, 1998;Marzocchi and Sandri, 2003). These estimates are performed by considering a learning period of the seismic catalog, and they are then used in the model to forecast future earthquake events.

OEF System Interface
The forecasts are available to users of the system by means of a comprehensive and interactive dashboard composed of an embedded Google Maps showing 1) the current weekly probability map, 2) a timeline graph showing the probability history for a single cell or area, and 3) the current probability The complete dataset used for magnitude analysis counts 31,322 events in the time span from 1900 to 2020 covering a large area extending also beyond the Israeli boundaries, as shown in Figure 1. Since we focus on the area of the State of Israel, we restrict the analysis to the region included in the latitudes 29.4°N -34°N and longitudes 33.9°E-36.3°E (red box in Figure 1). The number of events in this smaller area is 10,340. We also integrated the catalog with 5 events with M w in the range [4.8, 5.3] reported by the "Global Centroid Moment Tensor" (Global-CMT) database in the selected area. The analysis of the events' magnitudes in the considered catalog has shown a lot of heterogeneity. Duration, body and moment magnitudes have been adopted singularly or jointly to describe the size of the events, depending on the geographical location or on the type of energy released. To eliminate this source of uncertainty that would invalidate the reliability of the forecasts, we have developed a new catalog where the magnitudes have been homogenized through the General Orthogonal Regression (GOR) technique (Lolli and Gasperini, 2012). For the details, see the Supplementary Appendix.

Learning Phase
We selected the period from January 1, 1983 to December 31, 2015 for the learning period, as it was found to reproduce the best fit of the model parameters. It contains 1,690 events within 30 km depth and of magnitude ≥2.2, the latter being the completeness threshold value computed through the maximum curvature method (Wiemer and Wyss, 2000). The correlation distance of 9 km was adopted in the exponential kernel distribution of the smoothing algorithm, as determined by maximizing the loglikelihood function of the seismicity contained in half of the catalog, under the time-independent model obtained from the other half. We recall that the log-likelihood function of the ETES model in the spatiotemporal domain [0, T] × S is given by where θ is the set of parameters. Since a future earthquake could occur also in an area of low or null observed seismicity, we labeled the cells of a square lattice covering the whole region with the 1% of the total rate, divided by the total number of cells (surprise coefficient), as in Kagan and Jackson (2000). The seismicity rate of the learning catalog is shown in Figure 3, together with the estimates for the model parameters (p, c, k, d 0 , q). Since some physical investigations show that the static stress changes decrease with epicentral distance as r −3 (Hill et al., 1993;Antonioli et al., 2004), for the sake of simplicity and following the line of several papers in the literature (Lombardi and Marzocchi, 2010b;Giuseppe Falcone et al., 2010) here we impose that q 1.5. This is a value close to that obtainable from the maximum likelihood best fit, and which follows the theory of elasticity when the spatial distance from the triggering event tends to infinite; the recognized trade-off between q and the other spatial parameter d 0 also justifies this choice, in fact different pairs (q, d 0 ) are shown to provide almost the same model's likelihood (Kagan and Jackson, 2000).

Forecasting Phase
The forecasting test we carried out in the OEF experiment for Israel has been implemented from January 01, 2016 till November 15, 2020, in the spatial extent [29.4°N, 34°N] latitude × [33.9°E, 36.3°E] longitude, gridded in a 0.1°× 0.1°square lattice covering the entire region. As preliminary OEF experiment presented here, for the sake of simplicity we use in the whole forecasting period the same parameters set estimated in the learning phase. This is also justified by the fact that the ETES model has been shown to be able to track the seismic variability even with constant parameters, as shown for example in the case of Tohoku M w nine earthquake, when ETES performed better than models with varying parameters (Nanjo et al., 2012). In any case, our future scope (and work in progress) is to improve the reliability of the forecast by following the Bayesian approach, which consists in daily or weekly updating the parameter estimations (Omi et al., 2013;Omi et al., 2015).
In this OEF experiment, each day of the temporal window specified above we produce weekly forecasts of the expected rates (or probabilities) of earthquakes with magnitudes ≥ M w , where M w 4.0 and 5.5. These two threshold can be re-set in the future for any magnitude or macroseismic intensity threshold, that are of interest for the stakeholders. We also stress that the choice of 7days forecasts has been inherited from OEF-Italy, as an explicit request from Italian Civil Protection to have short-term forecasts within a temporal window reasonable to activate state of alert procedures; anyway, the forecasts are produced at the midnight of ever day, and after the occurence of any M w ≥ 3.5 event recorded by the seismic network, therefore the government agency can modify the forecating time window according their specifici needs. The output produced by the system is the same as illustrated for Italy in Marzocchi et al. (2014): time-dependent maps showing the weekly probability for M w ≥ M w earthquakes occurring in each cell of the spatial grid covering the analyzed region. As an example, we show in Figure 4 the forecast of M w ≥ 4.0 events produced for the week starting from August 31, 2018, 00:00 UTC; see also Figure 2. The map shows that on August 31 the OEF forecast have conferred a maximum of probability in the cells of the grid covering the Sea of Galilee, on the city of Tiberias, each cell having a probability of about 10 −3 for M w ≥ 4.0 earthquakes. This maximum is caused by the occurrence, 1 month before, by a M w 4.7 event which induced an increase of the moderate earthquake rate. Besides this maximum in probability, the map shows also a peristent higher probability along the Death Sea valley as expected by the higher rate of historical seismicity.

Testing Phase
In this section we discuss in detail the testing analysis of the produced forecasts' reliability.
The testing phase of the model is carried out following the guidelines of the Collaboratory for the Studies of Earthquake Predictability (CSEP). Specifically, the N-test and S-test are performed to assess the reliability of the delivered forecast for Israel, i.e., to test the consistency of the OEF forecasts with the earthquakes occurred.
The N-test evaluates the consistency between the number N fore of forecasted earthquakes in all space-time-magnitude bins, and the number N obs of events observed over the entire testing region within any forecasting time window, and over any magnitude larger than the selected threshold. The two-tailed p-value is obtained by assuming that the target earthquakes follow the Poisson distribution with N fore mean. More precisely, given the collection J of synthetic catalogs representing the forecast, and the relative empirical cumulative distribution functions {F N fore,j } j∈J , the test is performed by computing the two quantile scores The result of the test indicates an over-predicting or under-predicting forecast, respectively, when δ 1 or δ 2 goes below a critical threshold value, which in our case corresponds to the 0.01 significant level (Taroni et al., 2018).
The S-test evaluates consistency of spatial occurrence of target events with respect to the model's normalized spatial forecast. The spatial component of the forecast is isolated, and the forecast is normalized so that its sum matches the observation. The test is then summarized by the quantile score that is the fraction of simulated synthetic catalogs of target earthquakes having spatial log likelihoods S X { } smaller than the observed spatial log likelihood S, calculated with the observed target earthquakes ( S S { } is the set of all the simulated spatial likelihoods). This quantile score represents also the p-value of the S-test. An inconsistent forecast is then obtained from the latter when ζ falls below the significance level fixed for this test, which we set here at 0.01. More precisely, a small p-value indicates that the fit between forecasts and data is worse than the one expected if the model was generating the data (Taroni et al., 2018). We stress that we use the 1% significance level since the Nand S-tests are susceptible to the use of a Poisson assumption; in fact, target earthquakes are expected to follow overdispersed distributions but, from its perspective, the Poisson assumption could induce the rejection of models that capture such overdispersion (Lombardi and Marzocchi, 2010a), like it happens for the ETES model adopted in the OEF-system for Israel.
We apply the N-and S-tests for all the 255 weekly windows entirely included in the forecasting phase, that is, 255 weeks starting from the first Sunday after January 1, 2016, occurred on January 3. The target events are those with magnitude M w ≥ 3: we use the rescaling technique to obtain the relative probabilities from the forecasts produced for the events ≥4.0. In the 255 weekly testing windows we found 48 target events, distributed as in Figure 5A. The consistency of the ETES model is found positive for N-and the p-value obtained is shown in Figure 5B. The S-Test results are presented in Figure 5C, and indicate that the cumulative p-value trend is not always positive and that in some cases it falls below the significance level 0.01. This could be due to the occurrence of small earthquakes with M w 3 in regions with expected low seismicity rate. In any case, we also calculate the S-test by considering as target the few events with M w ≥ 4, that is the same threshold adopted in the Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 729282 7 learning period. In this case, the resulting p-values fall completely above the significance level 0.01 for the pair of events occurring within the 255th week ( Figure 5D). The corresponding N-test is shown in Figure 5E.

DISCUSSION AND CONCLUSION
Since deterministic earthquake prediction is still an elusive enterprise, growing efforts are dedicated to deliver reliable earthquake forecasts which may help communities to be prepared and possibly to reduce the number of casualties. In line with this goal, in this paper we presented the development of the first Operational Earthquake Forecasting system for the State of Israel. In fact, the governmental Israeli steering committee for earthquake preparedness expressed interest in developing and implementing a real-time OEF system to reduce seismic hazard in Israel, on the basis of the prototipe developed for Italy. This resulted in the bilater Italy-Israel OPERA project, to whose final scope this paper responds.
This pilot OEF-Israel system is based on the stochastic clustering ETES model, and it is implemented similar to the Italian system . We have considered the time window 1900-2020 and the spatial extent [29.4°N, 34°N] latitude × [33.9°E, 36.3°E] longitude, which covers the land of Israel.
The set of the model parameters and the background seismicity have been estimated through a maximum likelihood approach using a learning period of about 32 years. Then we have produced earthquake forecasts for Israel. In each day between January 3, 2016 and November 15, 2020 we have delivered weekly forecasts for M w ≥ 4.0 and M w ≥ 5 events in the area [29.4°N, 34°N] latitude × [33.9°E, 36.3°E] longitude, providing timedependent seismic maps of the relative expected rates. An example is illustrated in Figure 4 for the M w ≥ 4.0 events produced in the week starting from August 31, 2018, 00:00 UTC. This map shows an increase in probability around the Sea of Galilee that was caused by the M w 4.7 event which occurred in this region a few days before.
The reliability of the OEF forecasts have been checked by adopting the CSEP guidelines. In particular, we have performed the N-test and the S-test to analyze the outcomes relative to M w ≥ 3 target events occurred during the 255 weeks entirely contained in the forecasting temporal window. We have obtained satisfactory results for both tests, with only a 2% overestimation of the predicted events with respect to the observed ones. The spatial distribution is also well captured by the model most of the time. However, we have found a few cases in which the p-value is particularly low for the S-test; this is due to the occurrence of a moderate seismicity in the testing area in areas where a low activity is expected. This calls for the need of considering improved alternative models (stand-alone and/or in an ensemble approach), for the aim of obtaining more reliable spatial distributions. Some new models have already been proposed in the literature by some Israeli colleagues (Zhang et al., 2021), and this will be the object of future further studies. As a preliminary experiment, we are confident that the OEF system we developed in this paper for the State of Israel could be an additional tile that, reinforced in its weaknesses, will allow pursuing the aim of delivering real-time probabilistic forecasts worldwide.

AUTHOR CONTRIBUTIONS
All authors contributed to the conception and design of the study, which has been developed in the framework of the Italy-Israel bilateral project OPERA coordinated by SH and WM. GF organized the database and implemented the software codes; GF, IS, and WM performed the statistical analysis. GF and IS wrote most of the paper with contributions from all the other authors. All authors reviewed and approved the submitted version.

FUNDING
This work has been supported by the bilateral Italy-Israel project "Enhancing OPerational Earthquake foRecasting through innovative Approaches (OPERA)", which has been funded by the Ministry of foreign affairs and international cooperation of the Italian republic and the Ministry of science, technology and space of the State of Israel, and by the Real-time Earthquake Risk Reduction for a Resilient Europe (RISE) project, funded by the European Unions Horizon 2020 research and innovation program (Grant Agreement Number 821115).