Short-term and long-term earthquake occurrence models for Italy : ETES , ERS and LTST

This study describes three earthquake occurrence models as applied to the whole Italian territory, to assess the occurrence probabilities of future (M ≥ 5.0) earthquakes: two as short-term (24 hour) models, and one as long-term (5 and 10 years). The first model for short-term forecasts is a purely stochastic epidemic type earthquake sequence (ETES) model. The second short-term model is an epidemic rate-state (ERS) forecast based on a model that is physically constrained by the application to the earthquake clustering of the Dieterich rate-state constitutive law. The third forecast is based on a long-term stress transfer (LTST) model that considers the perturbations of earthquake probability for interacting faults by static Coulomb stress changes. These models have been submitted to the Collaboratory for the Study of Earthquake Predictability (CSEP) for forecast testing for Italy (ETH Zurich), and they were locked down to test their validity on real data in a future setting starting from August 1, 2009.


Introduction
Despite the notable lack of success for reliable prediction of destructive earthquakes, there has been a recent resurgence of research into earthquake predictability that is now motivated by better monitoring networks and data of past events, new knowledge of the physics of earthquake ruptures, and a more comprehensive understanding of stress evolution and transfer.Unlike the old deterministic way of formulating earthquake predictions, which was based on various kinds of more or less popular precursors, much of the recent research into earthquake prediction is aimed at a quantitative specification of the uncertainty characterizing earthquake forecasts.These forecasts are formulated in statistical terms, and are currently applied to probabilistic hazard assessments of earthquake risk.Most importantly, the objective definition of the forecasting models and their capacity for computing the space-time density rate of future earthquakes allows the testing of such models in a prospective way, against observations of the real seismicity [Console 2001, Jordan 2006].
Stochastic short-term models that describe the phenomenon of earthquake clustering are achieving increasing success in the seismological community [e.g., Helmsetter et al. 2006].These models were proposed to answer the most common questions of the general public and the media that arise in particular after sizable events, such as, «What will happen next?» and, «What is the chance that another large earthquake will occur?».
Stochastic short-term models describe seismicity as a random point-process, for which a continuous space-time density distribution of the earthquake occurrence can be defined.A best-fit procedure based on the maximum likelihood criterion has been used for statistical analysis of random processes.In particular, Kagan andKnopoff [1976, 1987], Kagan [1991], Ogata [1999], Kagan and Jackson [2000], Console and Murru [2001], Imoto [2004], Rhoades and Evison [2006], and Helmstetter et al. [2006] have applied the likelihood method to earthquake occurrence studies.
A different approach is associated with forecasting the time of the next earthquake in the long term: that of the hypothesis of characteristic earthquakes.This probabilistic approach assumes that on a given seismogenic source, strong earthquakes will occur over the time of interest that have similar rupture areas, similar mechanisms, and similar magnitudes, with their time intervals characterized by a remarkable regularity.These earthquakes are often assumed to have similar hypocenters, similar displacement distributions within the rupture area, similar source-time functions (leading to similar seismograms), and a quasi-periodic recurrence.This approach can be applied assuming a renewal model with memory.In the renewal model, the elastic strain energy accumulates over a long period of time, after the occurrence of one earthquake and before the fault will release the energy as the next earthquake.This model of earthquake occurrence assumes that the probability of an earthquake is initially low following a segment-rupturing earthquake, and increases gradually as the tectonic processes reload the fault.
A direct implication of the characteristic earthquake hypothesis is that the occurrence of earthquakes on individual faults and fault segments does not follow a loglinear frequency-magnitude relationship of the form

Subject classification:
Earthquake interactions and probability, Statistical analysis, Stress, Historical seismology, Earthquake faults.

Short-term and long-term earthquake occurrence models for Italy: ETES, ERS and LTST
Giuseppe Falcone 1,* , Rodolfo Console 1 , Maura Murru 1 described by Gutenberg and Richter [1944] (logN = a − bM).The characteristic earthquakes are assumed to be large enough to dominate the seismic moment release, and to substantially reduce the average stress.This approach to earthquake forecasting has been widely applied as the basis for long-term forecasts of future seismic activity, particularly in Japan and the United States.
Progress is also being made with physically constrained models that link stress changes to seismicity rate changes using the Dieterich rate-and-state model [Ruina 1983, Dieterich 1986, Dieterich 1992, Dieterich 1994, Stein et al. 1997, Parsons 2004, Parsons 2005].
Validation of earthquake forecasting/prediction models is the main rationale behind some recent international efforts, including Regional Earthquake Likelihood Models (RELM) and Collaboratory for the Study of Earthquake Predictability (CSEP).The validation process consists of two steps: 1) running of all of the codes simultaneously to forecast future seismicity in well-defined testing regions; and 2) comparison of the forecasts through a suite of tests.These tests are mostly based on the likelihood score, and they evaluate both the time and space performances of the forecasting models.All of these tests rely on some basic assumptions that have never been deeply discussed and analyzed.In particular, it is assumed that in each spatiotemporal bin, any forecasting model is expected to produce a number of earthquakes according to a Poisson distribution (with the rate estimated by the forecasting model), and independent of adjacent bins.
Here, we give a brief description of three different forecasting models that we have submitted on the Test Center for earthquake forecasting/prediction experiments of the CSEP initiatives at ETH in Zurich.All of these models cover the whole of Italy as our testing region.The data for the learning phase have come from the seismic catalog collected by the Istituto Nazionale di Geofisica e Vulcanologia (INGV), which cover from July 1, 1987, to January 1, 2009, for a minimum magnitude of 2.1.We present two 24-hour forecasts that are based on the epidemic type aftershock sequence (ETAS) models [e.g., Ogata 1998, Console andMurru 2001], with the starting time for these forecasts of August 2009.The first model for short-term forecasts is a purely stochastic epidemic type earthquake sequence (ETES) model, where the temporal aftershock decay rate is governed by the modified Omori Law [Ogata 1983], and the distance decay follows a power law.The second short-term epidemic rate-state (ERS) forecast uses the Dieterich [1994] rate and state model to predict the spatial and temporal decay of the aftershocks.This more physical approach requires knowledge of the stress changes caused by each event.Rather than computing this directly for each triggering event (e.g., using the focal mechanism with Coulomb stress calculations), we assume that it follows a radially symmetric power-law decay with distance.For further information on the ETES and ERS models adopted, reference can be made to the studies of Console and Murru [2001], Console at al. [2003Console at al. [ , 2006aConsole at al. [ , 2006bConsole at al. [ , 2007Console at al. [ , 2010aConsole at al. [ , 2010c] ] and Murru et al. [2009].Both of these epidemic models assume that each event can spawn its own sequence of events (including aftershocks of aftershocks) and a Gutenberg-Richter distribution of events.The ETES and ERS models have been tested on real seismicity in Italy, the United States, Greece and Japan, through comparisons with a plain timeindependent Poisson models based on likelihood methods; these have shown the validity of these models in a retrospective way [Console et al. 2003[Console et al. , 2006a[Console et al. , 2006b[Console et al. , 2007]].
The ETES model has been applied since 2006 in an automatic way to the real-time data of the Italian Earthquake Data Center that is operated by the INGV [Murru et al. 2009].Moreover, a real application was started for the first time soon after the strong earthquake that struck the city of L'Aquila (central Italy) on April 6, 2009, at 01:33 UTC (M w ≥ 6.3).For a duration of a number of months, 24-hour earthquake forecasts were produced in near real time through an algorithm based on the ETES model, which were provided every morning for the Italian Agency of Civil Protection, for their use in the planning of rescue activities.At the time of the preparation of this article, this whole real-time process based on both the ETES and the ERS models was being developed at the INGV, as a real-time automatic information system that will be accessible as an experimental and confidential system through a specific website.This system is designed to become a public information service in the future.
We also present here our 5-year and 10-year earthquake forecasts that are based on a long-term model that considers the perturbations of earthquake probability for interacting faults by static Coulomb stress changes, which we refer to as long-term stress transfer (LTST).A similar approach has already been applied in Italy by Console et al. [2010b].

A short-term forecasting model based on earthquake clustering probability: the ETES model
Here we consider the short-term clustering properties of earthquakes, and we provide a brief outline of a statistical method for modeling the interrelationships between any earthquake and any other, as applied to the real-time Italian data at the ETH Testing Center; for further details, see the citations in the Introduction.This method is based on algorithms pertaining to the ETAS model that we have reported previously.Here, we have named this as the ETES model, to distinguish it from other ETAS models of the ETH Testing Center.
This stochastic process is characterized by a limited number of free parameters, as f r , k, c, p and d 0 in Equation ( 1), which allow the computation of the expected occurrencerate density as a continuous function in space and time, according to the definition of Ogata [1998], for the assessment of their maximum likelihood values [Kagan 1991].The expected occurrence-rate density is modeled as the sum of the independent, or time-invariant "spontaneous", activity f r m 0 (x, y, m), where f r is a factor known as the «failure rate» (i.e., the ratio between the expected number of independent events and the total number of events), and the contribution of every previous event using a kernel function that correctly takes into account: (i) the magnitude of the j-th triggering earthquake; (ii) the spatial distance r from the triggering event; and (iii) the time interval (t-t j ) between the triggering event and the instant considered for the computation, such that: (1) Although the coefficient f r is constant in space and time, the ratio between the background seismicity and the total seismic rate is not the same at any given time.During a sequence, the occurrence rate is dominated by the second term on the right-hand side of Equation ( 1).The failure rate f r is just the average proportion of independent events to the total number of events, such that Equation (1) holds when it is integrated over a large space-time volume.The smooth total time-independent rate-density function m 0 (x, y, m) is computed using the method introduced by Frankel [1995], which was described in detail by Console and Murru [2001] and then successively modified by Console et al. [2010cConsole et al. [ , 2010d]].The correlation distance used in the exponential kernel distribution of the smoothing algorithm, which was found to be 29 km, was determined by maximizing the likelihood of the seismicity contained in half of the catalog, under the time-independent model obtained from the other half of the catalog.
In our smoothing algorithm, the events received a weight proportional to the probability of being independent, as in the method introduced by Zhuang et al. [2002].These weights were adjusted following an iterative procedure that was similar to that adopted by Marsan and Longliné [2008].The rate distribution m 0 (x, y, m) was obtained from the distribution based on the weighted catalog, by normalizing this to the actual total number of events in the catalog.
During a sequence, the occurrence rate is dominated by the second term on the right-hand side of Equation (1).This does not mean that the f r needs to change in time and space.The coefficient f r is just the average proportion of independent events, such that Equation (1) holds when it is integrated over a large space-time interval.
In our algorithm, parameter f r is not determined from the best fit in the learning phase, because it is related to all of the other parameters of the model by the condition that Equation (1) must give the total number of observed events when it is integrated over the whole space-time volume spanned by the catalog.
Figure 1 shows the smoothed seismicity in the Italian region for period from July 1, 1987, to January 1, 2009.We fixed the value of q (the exponent of the spatial kernel of triggered events) at 1.5, which is a value close to that obtainable from the maximum likelihood best fit, such that this limits the number of free parameters in the model and makes the inversion procedure more robust.The b value of the Gutenberg-Richter magnitude distribution is assumed to be constant over the geographical area spanned by the catalog, and it has been obtained independently of the other free parameters, the meaning of which is reported in Table 1.The b value was determined from the maximum likelihood method of Aki [1965] (b = 0.9): (2) where M mean is the mean magnitude, and M 0 = M c − 0.05 (0.05 is half the magnitude bin unit of 0.1).
The average triggering distance of the aftershock zone, d j , is related to the magnitude of the main shock, m j , such that: (3) According to Equation (2), in our ETES model the average density of the triggered events does not depend on the magnitude of the triggering event.This allows a reduction by one free parameter with respect to other versions of the ETAS model [see for example, Ogata 1998,  Ogata and Zhuang 2006].We fixed the value of a at 1.0, to reduce the number of free parameters to be determined in the inversion, which increases the stability of the data.Indeed, parameter a is strongly correlated to parameter k.
The learning dataset for the best fit of the model parameters was the Italian seismic catalog from July 1, 1987, to January 1, 2009, over a rectangular area of 1,000 km × 1,200 km centered on the point with the geographical coordinates (42.0N, 13.0E).This provided 14,083 events of magnitude ≥2.1 (assumed this as the lower magnitude threshold of both the triggering and target events).A comprehensive analysis of the completeness for this catalog and its variations in time justified the choice of magnitude ≥2.1, as carried out by Console et al. [2010a].Here, we ignored the effects of increases in completeness magnitude soon after large mainshocks.
The maximum-likelihood best-fit parameters of the ETES model provided the results shown in Table 1.It can be noted here that the log-likelihood values are positive, while the values from most other studies are negative.This comes from the units chosen for m.As in all of our previous studies [Console and Murru 2001;Console et al. 2003Console et al. , 2006aConsole et al. , 2006bConsole et al. , 2007Console et al. , 2010aConsole et al. , 2010c;;Murru et al. 2009], we have multiplied m by a large factor V 0 , the dimensions of which are the inverse of those of m.This is equivalent to adding a large positive constant to the LogL.
Even if the data shown in Table 1 do not represent a test, as this is not the purpose of this study, it can be noted that, as expected, the likelihood of the catalog under the time-dependent clustering hypothesis is greater by far than that of the simple time-independent, spatially variable Poisson hypothesis.

A physically constrained epidemic model: the ERS model. Application of the rate-state model to earthquake clustering
We have also submitted a slightly different epidemic model to the ETH Testing Center: the ERS model, for short-term forecasting of moderate and large earthquakes in Italy.This model merges into a single algorithm the classical concept of ETAS (the purely stochastic model) and the rate-state constitutive law for the seismicity rate, as introduced by Ruina [1983] and Dieterich [1986Dieterich [ , 1992Dieterich [ , 1994]].In practice, our ERS model derives from the more typical ETES model by the substitution of the Omori law with Equation (3) for the temporal decay of the triggered events, and from some physical constraints put on the spatial distribution of the triggered events.The final model is stochastic, which allows computation of the likelihood of a seismic catalog, but which also reflects to at least some extent the physics of the earthquake processes [Console et al. 2006a, Console et al. 2007].According to Dieterich [1994] the rate m(t) of earthquakes after a Dx stress change at time t = 0 is given by: (4) where m 0 is the previous reference rate density, and A, v and t a are the physical parameters of the constitutive law.For all of our applications, A and v always appear as a multiplicative product, so we consider this product as only one independent parameter.It is also possible to show that t a is equal to Av divided by the stressing rate of the seismic area, , and the latter is itself related to the reference rate m 0 [Dieterich 1994, Console et al. 2006a].Rather than computing the stress change Dx caused at any point by any earthquake, for which it would be necessary to know the source parameters of all of the earthquakes, we have introduced a shortcut that allows use of the most common of the catalog information: the origin time, the epicenter coordinates, and the magnitude.Empirically, we hypothesize that the stress change produced by an earthquake is given by: (5) where Dx 0 is a free parameter that represents the maximum shear stress produced by the fault at its epicenter, r and q have the same meaning as in Equation ( 1), and d j is as in Equation (3).
Equation (4) implies that the stress increases everywhere, which may appear not very physical.However, experience has shown that frequently it is difficult to detect the effects of negative DCFF in terms of seismicity shadows.Several explanations have been proposed to justify this circumstance, among which there are the heterogeneity of slip on the causative fault and the effects of dynamic stress.In synthesis, seismicity increase is much more commonly observed than seismicity decrease.Moreover, in our case, there is also spreading of the triggered events due to location errors, which at least for weak triggering events, might dominate the geometric distribution.
It must be noted also that while the rate-state model specifies that the influence of successive stress steps multiply each other [see e.g.Dieterich 1994], our ERS model simply adds the contributions from each earthquake, as for all of the models of the ETAS group.
For numerical applications, it is necessary to define the values of the various parameters.In this study we aimed to reduce the number of free parameters as much as possible, so we arbitrarily fix the values of some of the parameters that can be reasonably estimated based on prior experience.These included: • d 0 is the radius of a circular fault of magnitude m 0 , and it was derived from the equivalent seismic moment M 0 through expression of the seismic moment for a circular fault in terms of stress drop and source radius [Keilis-Borok 1959], and of the relationship between the seismic moment M 0 and the magnitude m of an earthquake.For a magnitude threshold m 0 = 2.1, the computation gave d 0 = 64 m.
• q was fixed at 1.5, as in the previous section, and for consistency with the theory of elasticity (Dx decays as r −3 when r → ∞, and the total number of triggered events is proportional to Dx, according to Equation 3).
• b was estimated from the catalog analyzed.
• the geographical distribution of m 0 was as adopted for the ETES model.This left us with only Dx 0 and the product Av as free parameters in this new model, to be determined by a maximum-likelihood best fit.For the ERS model, the best fit provided the following values: Dx 0 = 0.08 MPa; Av = 0.012 MPa.
The maximum log-likelihood under the ERS and Poisson models were 152,092.69(lnL 1 ) and 144,926.00(lnL 0 ), respectively, with f r and dlog (= ln(L 1 /L 0 )) of 0.84 and 7,166.69,respectively.The comparison between the maximum log-likelihood of the previous ETES model and that of this ERS model gives a dlogL of 18,287.51.
As mentioned for the data with the ETES model, these data cannot be considered as a test of the models (which is the purpose of the CSEP Testing Center).However, a comparison of the data for the ETES and ERS models in the learning phase shows a clear difference in favor of the ETES model.This was not the case in a previous study that was carried out on Japanese seismicity [Console et al. 2006a].As far as we understand, the explanation for this lies in the smaller value that was adopted in the ERS model for the spatial parameter d 0 , which was not a free parameter: d 0 = 64 m, compared to d 0 = 134 m obtained from the best fit of the ETES model.The latter value accounts for the geometrical spreading due to uncertainties in epicenter locations.This circumstance was not so important for the larger magnitudes of the Japanese earthquakes analyzed by Console et al. [2006a].
To submit our models to the Testing Center of the CSEP at ETH in Zurich, two codes were written to supply an output file in the format specified by the Testing Center: for each magnitude bin of 0.1 unit in the range from 4.0 to 9.0, the codes supply a Table containing the expected number of events in adjacent cells covering the test area of 0.1˚× 0.1˚.
Figure 2 shows examples of how the ETES and ERS models work with real data that shows the spatial changes of the expected occurrence rate density over 24 hours starting from March 21, 2009 (00:00 UTC).The parameter values used in these cases are reported in Table 1 for the ETES model and in this section for the ERS model.In Figure 2b (ERS), a background seismicity that is higher than that in Figure 2a (ETES) stands out, which is consistent with a larger number of independent events with respect to the total number of events.Conversely, the ERS model shows a SHORT AND LONG-TERM EARTHQUAKE MODELS seismicity that is characterized by greater peaks, but is spatially concentrated.A small zone can be seen in central Italy where the expected rate of M ≥ 5.0 events in these next 24 hours was larger than 0.1 (events/day/cell).This zone shows the effects of the L'Aquila seismic sequence that started in December 2008, and that strongly hit the city of L'Aquila in the Abruzzo region, central Italy, on April 6, 2009, at 01:33 UTC (M w ≥ 6.3).

Long-term stress transfer model
The methodology adopted with the LTST algorithm was based on fusion of a statistical renewal model with a physical model that considered fault interactions that in real circumstances can either increase or decrease the future earthquake probability, with respect to what would be expected by a simple renewal model.We considered the fault interactions by computation of the co-seismic static permanent Coulomb stress change (DCFF) caused by all of the earthquakes that occurred after the latest characteristic earthquake on a given investigated fault segment and in the surrounding sources.
According to the methodology developed over the last decade [Stein et al. 1997, Toda et al. 1998, Parsons 2004], and also applied by Console et al. [2010b], the probability of the next characteristic earthquake on a known seismogenic structure (for which the mechanism and size of the characteristic event is given) in a future time interval (in this case, 5 and 10 years) starts from the estimate of its occurrence rate, which is conditioned by the time (t) that has elapsed since the previous characteristic event.To do this, two parameters are necessary: the expected mean recurrence time, T r , and the aperiodicity, a (also known as the coefficient of variation) of the renewal process [Mc Cann et al. 1979, Shimazaki andNakata 1980].In the present study, from among various statistical renewal models, we adopted the Brownian passage-time distribution (BPT) [Matthews et al. 2002] to represent the inter-event time probability distribution, f (t), for earthquakes on individual sources.
This distribution function that provides the instantaneous number of events is expressed as: (6) The influence of the stress change on the probability of an impending characteristic event was computed by introduction of a permanent shift (clock advance or delay) in the time, t, that had elapsed since the previous earthquake, obtaining a modified time t' according to: (7) where DCFF is the static stress change, and is the tectonic stressing rate (assumed to be unchanged by the stress step).
As DCFF is not uniform on the receiving fault, we adopted its average value, which was given by the sum of the values computed on small elements of the triggered fault, divided by the number of elements.
We now need to consider how the stress changes caused by earthquakes that occurred on neighboring faults might affect the probability of occurrence of the next possible event  Basili et al. 2008, DISS Working Group 2009].Some sources in the investigated area were discarded, as they lacked the date of the last event.This reduced the number of sources to 104.Among the causative events that potentially changed the stress conditions on the 104 faults included, we considered the following: 1.The characteristic events associated with the seismogenic sources themselves (as reported in DISS 3.1.0);2. The events from 1895 to 2004 with M w ≥ 5.0 reported in the parametric catalog of the historical Italian earthquakes CPTI04 [CPTI Working Group 2004] associated with the DISS 3.1.0.seismogenic areas (192 events with M w ≥ 5.0).To this database, we added the recent largest earthquakes from May 2004 up to the main shock of the L'Aquila 2009 sequence (April 9, M w 6.4) (reported in http://emidius.mi.ingv.it).
Obviously, all of these events were considered only once if they were reported by more than one information source, with preference given to the information coming from the source in the order as they are listed.
The stress change DCFF on an individual fault was computed by adding the contributions from all of the other sources that ruptured after the latest known earthquake on the fault considered.The computation was carried out at the hypocenter depth of this latest earthquake.
The implementation of the method outlined in the previous section requires some information that is available in DISS, including the hypocenter coordinates, the expected magnitude, the focal mechanism, the fault size, the average slip, the mean recurrence time, and the date of the latest event.Following Ellsworth et al. [1999] and Matthews et al. [2002], in the present study we adopted the value a = 0.5 for the coefficient of variation.
The fault parameters of strike, dip, rake, dimensions and average slip are needed for all of the triggering sources.The fault mechanism is also needed for the triggered source (the receiving fault), to resolve the stress tensor on it.As we were dealing mainly with pre-instrumental events, for which details such as the fault shape and slip heterogeneity were not known, we assumed rectangular faults with uniform stress-drop distributions.We computed the clock change, Dt,  according to the ratio between DCFF and (tectonic stresschange rate).The clock change was positive (the fault comes closer to failure) if the Coulomb stress change was positive.In the opposite case (negative change), the fault moves farther from failure, where it is possible that in extreme situations the elapsed time can be reset to zero.The values of have been computed for each source according to the strain rates obtained by S. Barba for the INGV-Dipartimento della Protezione Civile (DPC) S2 project (2005)(2006)(2007).These values of the strain rate were obtained by different numerical models of deformation for Italy: using a finite-element method and calculating the nodal velocity by the weighted residual method, through the use of the SHELLS software [Bird 1999], which was suitably modified to include the representation of the seismogenic Italian faults.
The strain tensor was resolved on the specific source, taking into account the mechanisms of its characteristic earthquakes.The values of the shear-strain component so obtained allowed the computation of through multiplying it by the shear modulus n = 3 × 10 4 MPa.
SHORT AND LONG-TERM EARTHQUAKE MODELS  Once the distribution function f (t) was estimated, the expected number of events N (rate) over a given time interval (t, t + Dt) was computed on the cells defined by the CSEP as 0.1˚× 0.1˚ by integration: (8) With the hypothesis of a generalized Poisson process, we can estimate the probability of earthquake occurrence in the given time interval as: (9) Figure 3 shows the expected probabilities plotted on a map of seismogenic sources for the 5 years to 10 years from August 2009, as the probabilities obtained from the permanent effects only.This forecast relates to only the individual seismic sources, and so it is not testable according to the CSEP methodology, which is based on a grid approach.Considering the forecast over the next 10 years (Figure 3b), the highest probabilities are related to the following sources: Selci_Lama (Umbria region, central Italy) (7.37%), Aspromonte Northwest (Calabria region, southern Italy) (6.45%), Camerino (Marche region, central Italy) (5.16%) and Pontremoli (Liguria region, northwest Italy) (4.61%).The last characteristic events occurred on these four sources, on: September 30, 1789; February 6, 1783; July 28, 1799; and February 14, 1834; respectively.
To provide the CSEP Testing Center with a grid-based forecast, the expected number of events obtained for each of the sources was uniformly subdivided on the basis of the cell numbers included inside the sources.The center of these cells must be inside the rectangle that delimits the vertical projection of the sources.Moreover, two maps of the expected rates for cells of 2 km × 2 km were generated using the historical catalog (CPTI04, January 11, 1895, to May 2004; M w ≥ 5.0) and the Italian Seismological Instrumental and Parametric database (ISIDE, http://iside.rm.ingv.it;July 1, 1987, to June 2005, M 1 ≥ 2.1).These rates derived from the historical and instrumental catalogs were computed by the method introduced by Frankel [1995].We preferred to use a dense grid of 2 km × 2 km for the rate calculations, instead of the CSEP 0.1˚× 0.1˚ grid, for better resolution of the fault areas.Three rate values were obtained for each cell (one provided by the sources, the second from historical data, and the third from the instrumental catalog).We used the rate values corresponding to the sources when the cell was inside it.Otherwise, we used the average rate values obtained from two catalogs if the cells were outside of the sources.To submit the LTST algorithm to the CSEP Testing Center at ETH in Zurich, code was written to supply an output file in the format specified by the Testing Center.Thus, for each magnitude bin of 0.1 unit in the range from 5.0 to 9.0, the code supplies a Table containing the expected numbers of events in adjacent cells of 0.1˚× 0.1˚ covering the Test area.Figure 4 shows the expected event rates for M ≥ 5.0 in cells of 0.1˚× 0.1˚ for the next 5 years and 10 years.
As indicated for the short-term forecasts above, the testing of these models is the task of the CSEP Testing Center, to which our models have been submitted.Therefore, we have not discussed this issue further in this study, instead limiting our attention to the descriptions of  these models.Moreover, testing long-term earthquake forecasting models on new and independent data, as is necessary for any objective validation purposes, poses great difficulties due to the long times needed to collect the observational data.A description of a retrospective test of this LTST model can be found in Console et al. [2010b].

Conclusions
We have submitted three grid-based models to the CSEP ETH Testing Center that were applied to the whole Italian territory for the forecasting of impending shocks, starting from August 1, 2009.The first two 24-hour short-term forecasting models are based on two epidemic-type earthquake sequences.The ETES model is only a purely stochastic model, while the ERS model is also based on physical constraints through the application of the Dieterich rate-and-state model for earthquake nucleation.The LTST model is a 5-10-year, long-term forecasting physical model that considers the interactions among the seismic sources through the computation of the Coulomb failure functions (DCFF).While the first two models only use the information contained in a seismic catalog (time, latitude, longitude, depth, magnitude), the LTST model also uses geological and geodetic information.These models were produced to test the validities of their assumptions in a truly prospective test against the observed seismicity, and to explore which models among these submitted will be preferable for use in seismic hazard and risk assessment.

Figure 1 .
Figure 1.Smoothed distribution of the Italian seismicity (from July 1, 1987, to January 1, 2009) through the smoothing algorithm using 29 km as the correlation distance.The color scale represents the average number of earthquakes (M ≥ 2.1 and h ≤ 70 km) in an area of 2 km × 2 km over this time period.

Figure 3 .
Figure 3. Conditional occurrence probabilities obtained by the LTST algorithm for the next characteristic earthquakes in a future timeframe starting on August 1, 2009.Only the permanent effects have been considered.In all, 104 seismogenic sources were considered.(a) Probabilities over the next 5 years.(b) Probabilities over the next 10 years.

Figure 4 .
Figure 4. Forecast maps of modeled expected occurrence-rate density (events M ≥ 5.0, per 5 years per cell of 0.1˚× 0.1˚) for the whole Italian Test area, starting on August 1, 2009, provided by the LTST algorithm.(a) Probabilities over the next 5 years.(b) Probabilities over the next 10 years.