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.

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 λ ( x , y , t , M | H t ) = f R ⋅ λ 0 ( x , y , M ) + ∑ i = 1 N H ( t − t i ) λ i ( x , y , t , M ) , (1) 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; λ ₀ (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) p ( M ) = 1 0 − b ( M − M c ) , M ≥ M c , (2) 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 time-independent background function λ ₀ (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 λ i ( x , y , t , M ) = Φ ( t − t i ) ⋅ Ψ ( r ) ⋅ p ( M ) , (3) where Φ ( t − t i ) = k ( t − t i + c ) p is the modified Omori law for the temporal decay, and Ψ ( r ) = d i 2 r 2 + d i 2 q , d i = d 0 ⋅ 1 0 α ( M i − M c ) , (4) 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 ₀, 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.

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 values for a selected cell or geographical area (Figure 2). All probability maps and the time evolution of the probability for each single cell are stored for future reference and quick access. In this way, users can query the database by name of location (inverse geocoding) or by direct selection of a point or area on the map, and plot the evolution of the probability in any time window of interest.

The instrumental earthquake data in the State of Israel have been collected since the beginning of the 20th century by seismic stations in Egypt, Lebanon, Israel [Geophysical Institute of Israel/GEOFOrschungsNetz Station Jerusalem (JER), GEOFOrschungsNetz Station Eilat (EIL)] and several tens of stations in Europe. More precisely, the complete earthquake catalog in the period 1900–1982 has been recorded by the International Seismological Summary (ISS), the International Seismological Centre (ISC), the National Earthquake Information Center (NEIC), and the compilation of Arieh et al. (1985). From 1982 on, seismic data have instead been recorded by the two national seismic networks of the State of Israel and Jordan. The first, named the Israel Seismic Network (ISN), was established in 1982 and is operated continuously by the Seismology Division of the Geophysical Institute of Israel (GII); the second, named the Jordanian Seismic Network (JSN), was established in 1983 and is operated continuously by the Jordan Seismo Obs and Geo Studies–Natural Resources Authority (NRA) of the Hashemite Kingdom of Jordan.

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.

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 log-likelihood 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 l o g L ( θ ) = ∑ i = 1 N l o g λ θ ( x , y , t , M | H t ) − ∫ 0 T ∬ S λ θ ( x , y , t , M | H t ) d x d y d t , 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 ₀, q). Since some physical investigations show that the static stress changes decrease with epicentral distance as r ⁻³ (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 ₀ also justifies this choice, in fact different pairs (q, d ₀) are shown to provide almost the same model’s likelihood (Kagan and Jackson, 2000).

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 7-days 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⁻³ 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.

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 f o r e , j } j ∈ J , the test is performed by computing the two quantile scores δ 1 = 1 − F N f o r e , j N o b s − 1 = P ( N f o r e , j ≥ N o b s ) , δ 2 = F N f o r e , j N o b s = P ( N f o r e , j ≤ N o b s )

The result of the test indicates an over–predicting or under–predicting forecast, respectively, when δ ₁ or δ ₂ 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 ζ = # S X | S X ≤ S , S X ∈ S S # S S , 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 N- and 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 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.