A double-branching model applied to long-term forecasting of Italian seismicity ( M L ≥ 5 . 0 ) within the CSEP project

The purpose of this study was to apply a double-branching model to forecasting of moderate-to-large Italian seismicity within the Collaboratory Study for Earthquake Predictability project. This project is designed for statistical evaluations and comparisons of various forecasting models, on both the global and regional scales. This proposed doublebranching model is time-dependent, as it assumes that each earthquake can generate, or is correlated with, other earthquakes through physical mechanisms that act on different spatio-temporal scales. Specifically, it consists of an application of two branching processes, in which any earthquake can trigger a family of later events on different space-time scales. In our recent study [Lombardi and Marzocchi 2009], we applied this model to a declustered historical database that included the strong Italian seismicity from over the last few centuries. This catalog only allowed us to describe the long-term time evolution of moderate-to-strong seismicity. Here, we have applied this double-branching model to a new database that has allowed us to describe both short-term clustering and long-term features at the same time. As the model can produce forecasting calculations of future seismicity, we provide some probability maps of occurrence of predicted events over different temporal windows.


Introduction
The modeling of the spatio-temporal distributions of moderate-to-large earthquakes is a major goal of geophysical investigations.The huge social impact of risk mitigation has promoted the formulation of both physically based and stochastic models, which are all directed towards providing better understanding of the driving features of earthquake occurrence and/or towards forecasting of future events [Dieterich 1994, Ogata 1998, Stein 1999, Kagan and Jackson 2000, Gerstenberg et al. 2005].Most of the proposed models hypothesize a substantially stable temporal behavior of the seismic rate around a nearly constant mean, with short-term random fluctuations arising from aftershocks.This assumption explains the current use of the stationary Poisson paradigm in many practical applications that are related to long-term behavior of earthquake occurrence, such as in the formulation of probabilistic seismic-hazardassessment methodologies based on the Cornell method [Cornell 1968] and in the evaluation of earthquake prediction/forecasting models [e.g., Kagan and Jackson 1994, Frankel 1995, Varotsos et al. 1996, Gross and Rundle 1998, Kossobokov et al. 1999, Marzocchi et al. 2003].
Over the last two decades, considerable attention has been focused on long-term trends of seismic rates and on the multiyear stability of earthquake occurrence that has been called into question [Kagan and Jackson 1991a, Kagan and Jackson 1991b, Rhoades and Evison 2004, Lombardi and Marzocchi 2007, Marzocchi and Lombardi 2008].The longterm memory of the seismogenetic potential of a region can be ascribed to different physical causes, such as possible stress perturbations on spatio-temporal scales a lot larger than the duration of the sequences [Stein 1999], or combined effects of fault recurrence and fault interactions [Marzocchi et al. 2009].
In Marzocchi and Lombardi [2008], we formulated a new model, known as a double-branching model (DBM).The DBM takes into account long-term modulation of earthquake occurrence, as well as short-term clustering of earthquakes.This model has shown better earthquake forecasting performances for large earthquakes on both worldwide [Marzocchi and Lombardi 2008] and regional [Lombardi and Marzocchi 2009] scales, with respect to models with time-independent background rates.
Following on from this, we applied the DBM to Italian seismicity using a declustered historical catalog that allowed recovery of about 400 years of seismic history above a moment magnitude (M w ) of 5.5 [Lombardi and Marzocchi 2009].Here, we now apply this model to events above a local magnitude (M L ) of 4.5 that have occurred in Italy over about the last century.This last dataset also includes aftershock sequences, and it allows us to model the short-term features of the seismicity, which was ruled out in our previous study [Lombardi and Marzocchi 2009].
The DBM can provide earthquake forecasts for participation in the initiative known as the Collaboratory Study for Earthquake Predictability (CSEP; www.cseptesting.org)[Jordan 2006, Zechar et al. 2009].Here, we present probability maps of future earthquakes with M L ≥ 5.0 in Italy over two different temporal windows: 5 years and 10 years.These forecasts have been submitted to the CSEP EU-Italy Testing Center, to be evaluated by suitable statistical tools.

The double branching model
In the present study, we apply the stochastic model proposed by Marzocchi and Lombardi [2008], which consists of the sequential application of two branching processes in which any earthquake can trigger a family of later events.The main goal of our model is to highlight the interactions between events that are due to different physical processes and that involve largely different spatio-temporal domains.
In the first step of our modeling, we applied a version of the well-known epidemic type-aftershocks sequences (ETAS) model [Ogata 1998] to describe the short-term clustering of seismic events in space and time.The intensity function of this first branching model is given by: (1) where H t = {(t i , x i , y i , M i ); t i < t} is the observation history up to time t, M c is the minimum magnitude of the catalog, K 1 , c and p are parameters of the modified Omori Law that define the decay in time of short-term triggering effects, and a 1 defines the dependence of the triggering capability on the magnitude of an earthquake.The parameters d 1 and q 1 define the spatial distributions of triggered events as functions of the distance r i between a general location (x, y) and the epicenter of the i-th earthquake (x i , y i ) (C d1,q1 is a normalization constant).The function u 1 (x, y) is the probability density function of the locations of spontaneous events [Ogata 1998].Finally, b = b•ln( 10) is a parameter of the well-known Gutenberg-Richer Law [Gutenberg and Richter 1954], which is assumed as the distribution for the magnitude of all events.Using ETAS modeling, it is possible to compute for each earthquake a probability of being a triggered event, and consequently to formulate a stochastic declustering algorithm that can remove the short-term triggered events from the original dataset [see Zhuang et al. 2002, for details].
The second step of our procedure consists of reapplying a branching process to the filtered database that is obtained using the ETAS-derived declustering procedure proposed by Marzocchi and Lombardi [2008].This step is based on the stochastic technique of Zhuang et al. [2002], but unlike this stochastic technique, we apply an algorithm that produces a single declustered catalog.We have verified in the present and in previous analyses that the significance of the results does not depend on the adopted declustering procedure [see Marzocchi and Lombardi 2008, for details].This second branching model works on larger space-time scales compared to the smaller domains involved in shortterm clustering, which are removed after the first step.The time-dependent conditional rate of earthquake occurrences ( , , , / ) u ( , ) / for the second step of our procedure is given by: (2) The residual seismicity that is obtained by filtering the original database of these short-term triggering effects is therefore ascribed to the superposition of two physical processes: the time-independent and spatially variable tectonic loading n 2 = o 2 •u 2 (x, y), and the long-term coupling between events.This latter is described in time by an inverse exponential function with a characteristic time x.Similar to the first branching step, the dependence of a hazard function with the magnitude of an exciting event is assumed to be of an exponential type, i.e. proportional to e a2M .Finally, the spatial decay of long-term stress variation is described by an inverse power law probability density function, the function of distance r i from the epicenter of i-th perturbing event, and the parameters d 2 and q 2 (C d2,q2 is a normalization constant).To estimate the parameters of both of these branching processes, we used an iteration algorithm developed by Zhuang et al. [2002], which is based on the maximum likelihood method and on a kernel estimation of the total seismic rate.Further details on this model and on the estimation of its parameters can be found in Marzocchi and Lombardi [2008].
From a technical point of view, we chose to follow a step-wise procedure rather than to apply a single complex model, for two reasons: (1) the simultaneous estimation of a large number of parameters involves problems that do not have easy solutions; and (2) it has been shown that the description of additive complex processes through a step-wise procedure can lead to significantly better modeling [e.g. the boosting approach of Bühlmann 2003].

The datasets: the CPTI04 and CPTI08 earthquake catalogs
Italy is characterized by well tested experience in compiling historical databases [Camassi and Stucchi 1997, Boschi et al. 1995, Boschi et al. 1997, Boschi et al. 2000].The most recently compiled historical catalog of Italian seismicity is the parametric catalog of Italian earthquakes (Catalogo Parametrico dei Terremoti Italiani; hereinafter referred to as CPTI04); this catalog was published at 1999 [CPTI Working Group 1999] and revised in 2004 [CPTI Working Group 2004].The CPTI04 catalog was obtained by merging some previously collected datasets and reports of earthquakes with intensities I 0 ≥ V/VI on the Mercalli-Cancani-Sieberg (MCS) scale that occurred in Italy (or in neighboring zones) from 217 B.C. to 2002.Given the strong hazard-oriented nature of this catalog, foreshocks and aftershocks were removed from the original datasets that were merged into the CPTI04 catalog.Specifically, the declustering procedure consisted of removing all of the events that occurred in a 30-km and 90-day spatio-temporal window with respect to a shock with a larger magnitude.Historical investigations relating to the completeness of the CPTI04 catalog have identify a substantially complete recording for events with magnitude M w 5.5 from 1600 (217 events) [Stucchi et al. 2004].
The CPTI04 catalog is suitable for hazard analysis and for studies of the long-term behavior of seismicity, although it cannot be used to model short-term clustering of events.Recently, a new catalog was proposed in its first version, known as CPTI08 [Rovida and the CPTI Working Group 2008].This CPTI08 catalog includes aftershocks for which A DOUBLE BRANCHING MODEL FOR ITALY  Given the strong uncertainty associated with magnitude values, we analyzed the completeness magnitude of the CPTI08 catalog using the simple maximum curvature (MAXC) method [Wiemer and Wyss 2000].The overall estimate predicted a completeness magnitude M w 4.6, although the analysis in time shows a decreasing completeness magnitude from M w 4.8 to M w 4.6.So we selected the events above a magnitude M w 4.8 so as to have a homogeneous and complete catalog (see Figure 1).
The major goal of the present study is to present a model that can forecast future seismicity in Italy, in the framework of the CSEP project.Thus, we converted the moment magnitudes reported in the CPTI08 catalog into local magnitudes (M L ).In this way, we conform to the CSEP in their tests of proposed models of the seismic bulletin of the INGV, Istituto Nazionale di Geofisica e Vulcanologia (Bollettino Sismico Italiano; http://bollettinosismico.rm.ingv.it/and http://iside.rm.ingv.it),where they have adopted M L as their official magnitude scale.The magnitude conversion was obtained through the relationship proposed by the MPS Working Group [2004]: (3) By using this relationship, the completeness moment magnitude M w 4.8 is transformed into the completeness local magnitude M L 4.5.In Figure 2, we show the locations of the CPTI04 catalog events (M w ≥ 5.5) and the CPTI08 catalog events (M L ≥ 4.5).

Application of the double-branching model to the CPTI08 catalog.
We applied the ETAS model, the first step of the DBM, to the Italian seismicity recorded in the CPTI08 database.Following the procedure proposed by Zhuang et al. [2002], we estimated the model parameters, together with the spatial distribution of the seismicity that was not coseismically triggered (u 1 (x, y); see Equation 1).Table 1 lists the inferred values of the model parameters, together with their errors and the associated log-likelihood values.Figure 3 shows the ratios between the short-term triggering and the total rates for the whole period recovered from the CPTI08 database.These rates were obtained by integrating in time the intensity function m 1 (t, x, y, m/H t ) of the first branching model (Equation 1) and the contribution provided only by the short-term triggering effects (sum in Equation 1).As can be seen, in most of the zones, the short-term triggered rate is below 30% of the total seismic rate.In some areas however (Friuli, central Apennines, Irpinia, Messina Straits, western Sicily), there were well-identified sequences that reach 90% of the total rate.
The application of the second-order branching model (DBM) to the CPTI08 dataset provided the parameters reported in Table 2, together with their relative errors and the maximum log-likelihood values.Given the short duration of the CPTI08 dataset, we set the x parameter to a value (30 years) estimated from the longer CPTI04 catalog [see Lombardi and Marzocchi 2009].Moreover we fixed a 2 = 0, as seen in previous studies on both global and local scales [see Marzocchi andLombardi 2008, Lombardi andMarzocchi 2009].As discussed in Lombardi and Marzocchi [2009], this a 2 value might be due to the low number of events.In any case, we have verified that no different combinations of these parameters provides better log-likelihoods.To determine the reliability of the DBM, we also report here the mean rates and the maximum log-likelihoods obtained by Poisson modeling (see Table 2).The time-dependent model is significantly better than the Poisson process considering the Akaike information criterion [Akaike 1974].This result was   = confirmed by applying residual analysis [Ogata 1998, Lombardi andMarzocchi 2009] to the CPTI08 dataset.This diagnostic technique confirms that if a point process model with intensity m(t) describes the temporal evolution of the seismicity well, the transformed data are expected to behave like a stationary Poisson process with a unit rate [Ogata 1998].We tested this hypothesis by the Runs and the Kolmogorov-Smirnov tests [Lombardi and Marzocchi 2007], for independence and exponential distribution of residuals, respectively.The P value of the Kolmogorov-Smirnov test leads to the rejection of the ETAS model (P < 0.01), but not the rejection of the DBM (P = 0.2); the Runs test does not reject the independence hypothesis of residuals for either ETAS or the DBM.
The most reliable cause of this finding is a probable long-term time-dependence of the occurrence of damaging shocks in Italy.Indeed, the temporal scale recovered by the second model (with a characteristic time x of 30 years) suggested that the interactions revealed by the second branching model (DBM) are not likely to be due to poor ETAS model formulation.
To represent the significance of the long-term interactions revealed by the DBM, we show in Figure 4 the map of the proportion of seismicity due to long-term triggering effects, together with the map of the spatial distribution of tectonic-driven seismicity (u 2 (x, y)).The model defines the maximum effects of long-term interactions in the Friuli region (see Figure 4b), and a larger tectonic loading along the northern-central Apennines range and in the Friuli region (see Figure 4a).In any case, as the CPTI08 dataset recovers a time period well below the reloading time of the faults in Italy, it will not be representative of the real seismogenetic potential of these different zones.By considering the model estimated by Lombardi and Marzocchi [2009] (the first application of DBM to Italian seismicity) with the CPTI04 catalog, we obtained the map of the probability distribution u 2 (x, y)  To provide forecasts of future seismicity, we decided to use the parameters estimated in the present study and the spatial distribution of tectonic loading (u 2 (x, y)) estimated by Lombardi and Marzocchi [2009] on the CPTI04 catalog (Figure 5), which we judged to be the more representative.This choice implies that the distribution of M w ≥ 5.5 earthquakes is the same as the M L ≥ 4.5 events.For a very preliminary check of the forecasting capability of our model, we show in Figure 6 the map of the predicted number of events for the period January 1, 2007, to July 31, 2009.We also plot the locations of the events with M L ≥ 4.5 and greater than 30 km in depth that occurred over the same period in the CSEP testing region that were collected by the official bulletin of the INGV (http://iside.rm.ingv.it)(13 events).The basic parameters of these events are reported in Table 3.All of these events occurred in cells with relatively high forecast rates.For 10 of the 13 events, the proportion of cells with a higher forecast rate (with respect to the cell where the earthquake occurred) is lower than 5% (Table 3, last column).We note here that the forecasted rates were computed for the whole period, without taking into account the triggering effects of the real events that occurred within this period.In this way, we intended to reproduce the forecasts made for the CSEP, for which we have no information on future earthquakes and we can only simulate the triggering effects of these events.By including the real triggering effects for forecasting calculations, the results in Table 3 are further improved (at least for events 2 and 9).

Forecasting maps
The model formulated and tested in the previous sections allowed us to compute the forecasts for future seismicity in Italy.The predictions are in a form of the probabilities of the occurrence in Italy of at least one earthquake with a M L ≥ 5.0 within a cell of 0.1˚× 0.1˚ over the next 5 years (August 2009 to July 2014) and 10 years (August 2009 to July 2018).These are obtained by integrating in time and space a combination of the intensity functions of two branching models (Equations 1 and 2).Specifically for each cell C i and for each forecasting   period T j , we computed the relative rates of occurrence R ij according to the following integration: (4) where w i is the probability that the i-th event is not coseismically triggered, computed by the first branching ETAS model [see Marzocchi and Lombardi 2008].The forecasted rates above M L 5.0 were obtained by rescaling the rate of earthquakes above M L 4.5, in agreement with the Gutenberg-Richter relation.Equation ( 4) shows that time-dependent modeling such as the DBM can also take into account the triggering effects of the seismicity that occurred before and that are expected during the forecast interval.So, in the past history, we included all of the seismicity with magnitudes greater than M L 4.5 and depths greater than 30 km that occurred up to July 31, 2009, and were collected in the CPTI08 catalog ( January 1, 1901to December 31, 2006) and INGV catalog ( January 1, 2007to July 31, 2009).To generate the forecasting maps, we simulated 10,000 different stochastic realizations for each of two time windows, using the thinning method proposed by Ogata [1998] and the intensity function formulated in Equations ( 1) and (2).Here, no significant variations were found using a much lower number of simulations (up to 1,000).Then we averaged the predictions coming from each of these synthetic catalogs.The results relative to both of these time periods are shown in Figure 7.The expected numbers of events with M L ≥ 5.0 over the next 5 and 10 years are 9.4 and 18.4, respectively.The probabilities of one or more events over next 5 years with M L ≥ 6.0 and ≥7.0 are 40% and 4%, respectively.The DBM identifies the most dangerous zones as the regions that have been affected by the most recent earthquakes.Above all, there are the Abruzzo region, which was hit by the recent M w 6.3 L'Aquila earthquake (April 6, 2009; see Figure 5) and the Irpinia region that have the highest probabilities of experiencing a shock in the next 5 years and 10 years, respectively (see Figure 7).Smaller maxima of probability are seen in the northerncentral Apennines, Friuli, Calabria and the eastern part of Sicily.

Conclusions
The main goal of the present study was to describe the DBM as applied to forecasting Italian seismicity for the CSEP project.This model represents an implementation of the model described in Lombardi and Marzocchi [2009]; the A DOUBLE BRANCHING MODEL FOR ITALY main difference is in the seismic catalogs that were used to set-up the models.From a seismological point of view, the results obtained in the present study confirm the main findings of Lombardi and Marzocchi [2009].First, large earthquakes in Italy tend to cluster in time and space, and on a temporal scale (from decades to a few centuries) that is a lot longer than what would be expected for a typical aftershock sequence [see also Faenza et al. 2003, Cinti et al. 2004].Secondly, the time-dependent model proposed here is a significantly better model, in which the only time-dependent feature is short-term clustering.
The main improvements to this model, as compared to that described in Lombardi and Marzocchi [2009], is linked to the use of the first version of the CPTI08 catalog.This has allowed us to apply the whole DBM to Italy, including also the modeling of short-term behavior (Equation 1).In Lombardi and Marzocchi [2009], we applied only the second step of the DBM (see Equation 2) as the CPTI04 catalog, that we used in that study, is a declustered dataset.
The comparison of the forecasts proposed in the present study with the maps presented by Lombardi and Marzocchi [2009] highlights some differences.Specifically, the map of the probability distribution of occurrence over next 10 years (Figure 6b) proposed in the present study shows some hazardous regions that were not identified in our previous study: northern Apennines, the Gargano region and Eastern Sicily [see Figure 5b in Lombardi and Marzocchi 2009].These findings are due to the larger proportions of moderate events that have occurred in these zones over last decades, with respect to the stronger seismicity of the last four centuries that was collected in the CPTI04 catalog.We stress here that the spatial distribution of the tectonic driven seismicity (u 2 (x,y)) used in both of these studies is the same, and that the differences are due to long-term features that are modeled by the second step of our model (see Equation 2).
Finally we note that an intrinsic limitation of the DBM is that it is based on information coming from only a limited historical catalog.In this case, we argue that a few centuries could not be enough to get all of the time features to accurately describe the long-term variations; moreover, the inclusion of the relevant tectonic/geological information might represent a future direction for improvements to these models.

Figure 1 .
Figure 1.Completeness magnitude of the CPTI08 catalog (January 1, 1901, to December 31, 2006; 1591 events) by the maximum curvature method (MAXC).(a) Frequency magnitude distribution for the whole dataset: the MAXC Method provides M c = 4.6 (moment magnitude); (b) M c as a function of time.
known.The first published version of CPTI08 collected the shocks that occurred from 1901 to 2006 (1591 events).

Figure 3 .
Figure3.Map of ratios between short-term triggered and total seismic rates for the CPTI08 dataset.

Figure 4 .
Figure 4. Maps of (a) spatial distribution of tectonic-driven seismicity u 2 (x, y); and (b) ratio between long-term triggered rate and total seismic rate for the CPTI08 declustered catalog.

Figure 5 .
Figure 5. Map of spatial distribution of tectonic-driven seismicity u 2 (x, y) for the CPTI04 catalog.

Figure 7 .
Figure 7. Maps of probability of occurrence of one or more events with magnitude above M L 4.5 per cell of 0.1˚× 0.1˚ over the next 5 (a) and 10 (b) years.Green circles: spatial bin with highest probability.

Figure 6 .
Figure 6.Map of seismic rates (number of events per cell of 0.1˚× 0.1˚) predicted by the DBM for the period January 1, 2007, to July 31, 2009, inside the CSEP testing region (dashed line).Blue circles: locations of 13 events that occurred over the same period (from the INGV dataset).