Abstract

The smoothed seismicity models plan to be introduced into the new Italian Probabilistic Seismic Hazard Maps 2017-2018. In this study we report progress on the use of smoothed seismicity models developed by using fixed and the adaptive smoothing algorithms and present an earthquake rate forecast developed from an ensemble smoothed seismicity model. Recent developments in adaptive smoothing methods and statistical tests for evaluating and comparing rate models prompt us to investigate the appropriateness of adaptive smoothing together with the fixed smoothing seismicity models for the new Italian seismic hazard maps. The gridded-seismicity models are based on both historical and instrumental earthquakes and assume that larger earthquakes occur at or near clusters of previous smaller earthquakes. In general, the approach of using spatially smoothed historical seismicity is different from the one used previously by Working group MSP04 (2004), and Slejko et al. (1998) for Italy, in which source zones were drawn around the seismicity and the tectonic provinces and is the first to be used for the new probabilistic seismic hazard maps for Italy. We develop two different smoothed seismicity models following the well-known and widely applied fixed (Frankel, 1995) and adaptive smoothing methods (Helmstetter et al., 2007) and compare the resulting models (Moschetti, 2015) by calculating and evaluating the likelihood test. In this framework, the smoothed seismicity models are constructed by using both the historical Catalogue Parametrico dei Terremoti Italiani, CPTI15, (1000-2014) (Rovida et al. 2015) and the instrumental Italian catalog (1981-20165) (Gasperini et al., 2013) and their associated completeness levels to produce a space-time forecast of the future Italian seismicity. We follow guidance from previous studies to optimize the smoothing seismicity parameters; the correlation smoothing distance (fixed smoothing) and the neighboring number (adaptive smoothing) by comparing model likelihood values, which estimate the likelihood that the observed earthquake epicenters from the recent catalog are derived from the smoothed rate models. We compare likelihood values from all rate models to rank the smoothing methods. We also compare two our models with the best two models of the several Italian Collaboratory for the Study of Earthquake Probability (CSEP) experiment models, to check their relative performances. Finally we create six an ensemble models combining using two different smoothing models (adaptive and fixed), which are weighted equally with different weights through a logic-tree approach to improve the forecast capability (Marzocchi et al., 2012; Taroni et al., 2014) and to estimate the uncertainty associated with the models.

Introduction

Probabilistic seismic hazard analysis (PSHA) quantifies the probability of ground shaking at a specified site that exceeds a specified intensity level (Cornell 1968, SSHAC 1997). It contains two basic ingredientsparts: (1) the earthquake rate and rupture models, for which specifyication of the statistical distribution of earthquakes in time, space and size (magnitude); and (2) the ground motion prediction equations, for which estimateion of the expected ground shaking level (and distribution) at a site for each earthquake rupture (Field et al. 2005). Recently developed eEarthquake rate models which are the inputs to PSHA, developed for PSHA, may include estimates of earthquake recurrence from historical, geological, and paleo-seismological observations, and inferencesas well as the fromestimates from the crustal deformations and seismicity rates (Stirling et al., 2002; Field et al., 2014; Petersen et al., 2014; Adams et al., 2015). In these models we mostly perform past earthquake rates estimating future rates from the spatial smoothing of earthquake locations from catalogs, especially in regions where we have little, limited or incomplete information about the active faults. Modern earthquake rupture forecasts commonly use past earthquake rates to estimate future rates. Particularly in regions where knowledge about active faults is limited and not complete, earthquake rate forecasts are mostly developed from the spatial smoothing of earthquake locations from catalogs (i.e., smoothed seismicity models). In fact the Ssmoothed seismicity models have been commonlybroadly performedapplied over the past two decades in developing earthquake rupture forecasts (Kagan and Jackson, 1994; Frankel, 1995; Wang et al., 2011; Akinci, 2010) and seismic-hazard models (Frankel et al., 1996; Petersen et al., 2014; Adams et al., 2015; Akinci et al., 2004; Akinci et al., 2009).

Early smoothed seismicity models employed spatially uniform smoothing parameters (i.e., fixed smoothing), such that the kernels used to smooth catalog-derived seismicity rates were invariant to spatial variations in seismicity rate. However, some recently developed methods (Stock and Smith, 2002; Helmstetter et al., 2007; Werner et al., 2011) spatially adapt the smoothing parameters to the earthquake rate (i.e., adaptive smoothing). Adaptive smoothing methods have the general effect of concentrating seismicity rates near the clusters of high activity and reducing them in areas of low-background seismicity compared to models derived from fixed smoothing.

In this study, we apply both fixed and adaptive smoothed seismicity methods to develop a long-term earthquake rate model for independent events. The smoothed seismicity is an alternative to the approach that was used previously for the hazard calculation in Italy, in which area source zones were drawn around seismic or tectonic provinces. Special zones allow for local variability in seismicity characteristics within a zone (for example, hypocentral depth changes, changes in b-value, changes in maximum magnitude Mmax, and uniform seismicity characteristics). These models are combined to account for the suite of potential earthquakes that can affect a site. However, one advantage of the smoothed-seismicity methods is to avoid choosing zone boundaries that are sometimes poorly controlled by data and drawn by subjectively merging geological and seismological informationprobability. The delimitation and parametric characterization of small zones can be responsible for introducing uncertainties into the hazard evaluation. In its purest form, the smoothed-seismicity method simply assumes that patterns of historical earthquakes predict future activity, but it can easily be supplemented by tectonic – or geodetic – based zones or other model elements if there is reason to suspect that seismicity catalogs are insufficient. Even though in the approach by Frankel (1995) no seismicity source zones are needed, some model parameters can be taken as homogeneous throughout regional sub-zones. For example in this study, the study region is subdivided into seven broad zones mainly on the basis of the different catalog characteristics/completeness of the region.

In this study the smoothed seismicity models are constructed based on the spatial distribution of both the historical (Catalogue Parametrico dei Terremoti Italiani, CPTI15) and the instrumental Italian catalogues and their associated completeness levels produced for the Mappa di Pericolosita’ Sismica (MPS16) to attain a space-time forecast of the future Italian seismicity. The fixed smoothing model follows the method of Frankel (1995) and uses spatially uniform smoothing parameters (i.e., fixed smoothing), such that the kernels used to smooth catalog-derived seismicity rates are invariant to spatial variations in seismicity rate. However, the adaptive smoothing method (Helmstetter et al., 2007) defines a unique smoothing distance for each earthquake epicenter from the distance to the n-th nearest neighbor. We examined the ability of both models to predict the spatial distribution of seismicity from the recent part of the earthquake catalog using standard likelihood calculations. We followed the si

milar likelihood testing methods that have been spurred by the Regional Earthquake Likelihood Models (RELM) (Schorlemmer et al., 2007) and the Collaboratory for the Study of Earthquake Predictability (CSEP) (Zechar, 2010) testing centers for evaluation and comparison of earthquake forecast models (Schorlemmer et al., 2007; Zechar, et al., 2010). The aim of this effort is to develop and provide models for the seismicity rates using the most recently practiced methodologies along with updated databases (seismic catalogs) to be performed towards to new probabilistic seismic hazard maps for Italy.

Earthquake Catalogs and Completeness

The seismicity based source models benefit the smoothed seismicity rates calculated on a spatial grid platform from the two catalogs:

1) The parametric catalog of Italian earthquakes (Catalogue Parametrico dei Terremoti Italiani, CPTI15) that contains larger earthquakes up to magnitude 7.3 since 1000 up to year of 2014;

2) The instrumental catalog of Italian seismicity that contains the small earthquakes down to Mw 1.0, with a maximum of Mw 6.53, over the past 365 years (1981-201675).

The shallow background seismicity rates are calculated on a spatial grid using the historical CPTI15 catalog (Rovida et al. 2015), consisting of 442718 records (M≥>4.0) in the time window from 1000 AD to 2014. Since PSHA assumes a Poissonian process of earthquakes, where seismic events are considered temporally independent, the CPTI15 catalog is declustered using Gardner and Knopoff (1974) method to remove large fluctuations of seismicity rates in space and time due to aftershock sequences, and we select total of 28971594 main shocks for the final computation with depth ≤30 kKm. The completeness periods of the CPTI15 catalog were identiﬁed by the Working Group 2015 (http;//emidius.mi.ingv.it/CPTI15-DBMI15/data/CPTI15_descrizione.pdf), based on both historical and the statistical analyses. The completeness magnitude threshold was deﬁned over 157 magnitude bins, and represents the centers of each magnitude class with a width of 0.232, starting from a minimum magnitude of Mcw 3.734.0 (as the central value of the ﬁrst magnitude class) to a maximum magnitude of Mcw 7.41410 for different periods of time using both the historical and statistical approach (Albarello et al., 2001). For example, the lower magnitude limit 4.650 in the complete catalog represents the center of the 4.5353.96 to 4.76518 class. Table 1a and b shows the completeness of the catalog in the ﬁve six zones in terms of the variability at the beginning of this completeness time over seventeentwelve magnitude ranges from the historical and the statistical approaches, respectively. Figure 1a shows these seven six zones, which are indicated by #zone numberdifferent colors for different historical completeness time intervals and magnitude thresholds, as given in Table 1a. These zones have been used also to estimate b values through the Gutenberg-Richter (GR) relation (Gutenberg and Richter, 1949), for the fixed smoothed seismicity model (see following section). Only earthquakes with focal depths less than 30 km are considered in this study.

The instrumental catalog 1981-201567 (actually it starts on 1/1/1981 and ends on 30/4/2017) is based on the Gasperini et al (2013) catalog, updated for the new seismic Italian hazard map. After declustering (method by Gardner and Knopoff 1974), from this catalog we select 2560412 shallow earthquakes (depth < 30 Km) inside the CSEP Italy zone, with magnitude M≥3.0 (Figure 1b2). Completeness magnitude Mcw 3.0 is the same suggested by the authors of the catalog (Gasperini et al, 2013).

Constructing Smoothed Seismicity Models

The gridded seismicity sources assume earthquake rate models with generating cumulative earthquake counts by assuming the magnitude-frequency distribution with the truncated form of the GR relation (Gutenberg and Richter, 1949),

Log10 N (MMc) = a – bMcM, (2)

where N(Mc) is the number of events with magnitude equal to or greater than McM, and the b-value is the slope of the distribution that describes the relative frequency of small and large earthquakes within the completeness range MminMcM≤Mmax. The Weichert maximum-likelihood method (Weichert, 1980) was used to obtain the 10a values, with the completeness magnitude thresholds over different periods of time as given by the MSP16 Working Group (2016) (see Table 1a and b), so that the rare large and the smaller earthquakes with short recording times were obtained more accurately. Based on the assumption that the seismicity rate is constant with time, the method included each earthquake that was counted as of equal weight in the rate calculations. Thus, the seismicity rates during time periods that have more countable earthquakes (for example; time periods with lower completeness thresholds) will be affected by the completeness period much more strongly than the seismicity rate at other times (Felzer 2008). Cumulative numbers of events are sampled on a 0.1-by-0.1 degree grid by counting the number of earthquakes with Mmin≥4.5Mc. , However which is the minimum magnitude threshold stabilized by in the MPS16 Project, is M4.5 in each cell of the grid across the Italian territory and surrounding region for the final earthquake rate model to be used in the probabilistic seismic hazard calculations. In order to calculate the probabilistic ground motions from our smoothed seismicity models, however, we convert from cumulative recurrence to incremental recurrence using established relations (Herrmann 1977) which also requires a b-value. Therefore, the resulting “agrid” gives the annual rate of earthquakes in each grid cell as an incremental 10a in the Gutenberg-Richter notation, between plus or minus 0.05 magnitude units or a 0.1 bin-width centered on M= 0 in each grid cell.

We use two different approaches to estimate the b-value for the two catalogs. For the historical catalog CPTI15 the b-value is assumed to be variable throughout the sixeven regions (Figure 1a,b2) and calculated using the Weichert (1980) formulation for for both catalog from different magnitude-time completeness interval (Table 2). We start our computation from the third bin of magnitude (M 4.535) because the completeness of the first two bin is not entirely reliable (as suggested by an internal report of the project).

To estimate the maximum magnitude for the historical catalog, we use a classical approach, where maximum possible magnitude is equal to the maximum observed magnitude plus a constant, that can be computed using statistical consideration, or using the standard error associated with the maximum observed magnitude (Kijko 2004). In our case, the maximum observed magnitude is M7.3 (Val di Noto earthquake 1693) and the associated standard error is 0.2; so we chose a maximum magnitude equal to 7.5. We use the same maximum magnitude for the whole Italian territory.

For the instrumental catalog, we use the tapered version of the Gutenberg-Richter law (Kagan and Jackson 2000). We chose to use a different distribution for this catalog to explore the epistemic uncertainty relative to the choice of the type of the Gutenberg-Richter law that must be used (truncated or tapered). For this distribution we need to estimate the b-value and also the corner magnitude (that is a “soft” version of the maximum magnitude): due to the difficulty of estimation of the corner magnitude with a catalog of few years, like the Italian instrumental one, we also use some information of the historical catalog. Then we merge this information to estimate the parameters (Kagan and Schoenberg 2001). In this estimation, the total rate is the observed number of events with magnitude M≥ 4.5. In Appendix A we give more detail of this computation.

The rate i in the each cell ith is calculated following equation;

(3)

where , , and are the smoothing kernel, smoothing distance

, Mc value and event rate respectively, for the jth earthquake from a total of N events, and rij is the distance between the jth event and ith cell.

We then smooth seismicity rates using an isotropic, 2-D Gaussian function as the smoothing kernel:

(4)

The kernel is normalized by the coefficient, C (σj) that is proportional to σj -2. In equation (4), σj is the smoothing distance associated with an earthquake j in the spatial domain as the horizontal distance between event j and its n-th closest neighbor. So that the epicentral distance as “neighbor distance” corresponds to a value N as the “neighbor number”. Adaptive kernels assume that σj vary as a function of the density at location of each earthquake. So that σj smoothing distance decreases when the density of the seismicity is high at the location ri of the earthquake j, providing higher resolution where the density of seismicity is higher. The adaptively smoothed seismicity models make use of spatially varying smoothing distances. In fact the fixed and adaptive smoothed seismicity models differ only in the application of the smoothing distance applied to each earthquake (Moschetti, 2015). Fixed smoothed seismicity models use a single smoothing distance for all earthquakes while the adaptive smoothed seismicity models use unique smoothing distances to each earthquake. Since the adaptive smoothed models vary smoothing distance with the local seismicity density, the resulting smoothing distances are relatively smaller in regions of high seismicity than regions with sparse seismicity. In this study we implemented a maximum smoothing distance of 200 km following Moschetti, (2015) where we have to impose a physical constraint on the smoothing in regions of sparse seismicity in parts of Italy.

Smoothing parameters are based on earthquake location of earthquakes uncertainties and spatial trends seen both in historical and instrumental seismicity. One problem with the smoothing method is apparent in Italy where seismicity that occurs in narrow linear zones is over smoothed into nearby aseismic regions and become important to select optimal and proper smoothing parameters to be used for the earthquake rate forecasting in the study region.

In this study we use both the fixed and adaptive smoothing seismicity parameters in the application of the smoothing distance. For the fixed smoothed seismicity we tested the set of smoothing distances (5,10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95,.. and 100 km) for producing optimized fixed smoothed models, using historical seismicity from the CPTI15 catalog for Mcw≥4.5 events with the completeness magnitude thresholds over different periods of time (see Table 1a and b), so that the rare large and the smaller earthquakes with short recording times were obtained more accurately. However we selected neighbor numbers 1-10 and computed intermediate neighbor distances within this neighbor-number range for identifying optimal adaptive smoothed models. Then we follow the methodology by implementing an adaptive smoothing scheme that provides some smoothing but generally keeps the modeled hazard seismicity rates closer to the original seismicity using respective correlation smoothing distances for directions parallel and normal to dominant seismicity trends.

To construct the adaptive seismicity models, we employ the method of Helmstetter et al. (2007), which determines the smoothing distance for each earthquake from the distance to neighboring earthquake epicenters. We refer to the smoothing models that result from the use of smoothing distances to the n-th nearest neighbor as N-neighbor-number models. We perform compare the use of the instrumental catalog, in which completeness magnitude is smaller Mw 3.0, compare to that of the historical catalog (Mw 4.5), which so that contains more events in the catalog (N=1430 M≥3.0 between 1981-2000 for the instrumental compare to N=675 M≥4.5 between 1895-2000 for the historical catalog) and may be better identifyies the active fault structures and the spatial distribution of the local seismicity that may be locations of the larger events especially in the case of the adaptive smoothing approach. In the following we give optimized smoothing parameters (correlation smoothing distance and the N-neighbor number) obtained through the L-test likelihood testing follows the general CSEP testing methodology (Werner et al., 2011). Finally, for the recurrence time of the seismic activity, we assumed a time-independent (Poisson) model in which the occurrence of an earthquake does not change the probability of occurrence for following events. For a Poisson process, the probability P of occurrence of one or more events, in a time period T of interest, is given by Equation (5):

(5)

where ni is the rate of earthquakes and is the inverse of the average recurrence time.

More to add into this section??? MORGAN

Likelihood calculations for model optimization

We test compare the ability of the a smoothed seismicity models generated from the fixed and adaptive methods—and for different correlation lengths and neighbor-numbers, respectively— to predict forecast the location of earthquakes in the recent part of the earthquake catalog in order to optimize the smoothed seismicity models using through a likelihood parameter. In particularFor the likelihood parameter, we make use of a variant of the L-test, which is commonly used by several studies for the comparison of seismicity rate forecast models (Field et al., 2014) and the Collaboratory for the Study of Earthquake Predictability (CSEP) projects (Schorlemmer, 2010) and in similar studies for optimizing smoothed seismicity models (e.g., Werner et al., 2011). The L-test likelihood calculation is performed by separating the earthquake catalog into two parts: (1) An early part which we refer to as the forecasting catalog; and (2) A more recent part of the catalog, which we refer to as the testing catalog. From the forecasting catalog, wWe generate adaptive smoothed seismicity models from the forecasting catalog using neighbor-numbers N=1 to N=10 and fixed smoothed seismicity models using correlation smoothing distances 10 5 km to 100 km. We and calculate the number of events with magnitudes Mw ≥4.5 that occur in 5, 10- and 15-year window within each grid cell. Forecast earthquake counts are made for the four different completeness levels for the historical and instrumental catalogs. To calculate the likelihood that the spatial distribution of earthquakes in the testing catalog derives from an underlying distribution given by the smoothed seismicity model, we assume the Poisson distribution for earthquake occurrence:

(5)

where i and ni are the predicted and observed number of earthquakes in grid cell, i, respectively. Following others (e.g., Werner et al., 2011), wWe then normalize the total number of earthquakes in each smoothed seismicity model to the number of observed earthquakes in the testing catalogs, Ntot:

(6)

The joint log-likelihood value associated with a smoothed seismicity model, m, is computed from the sum of log-likelihood values computed for all grid cells:

(7)

and comparison of smoothed seismicity models employs the per-earthquake information probability gain parameter:

(8)

where Gm, is the information probability gain for the smoothed seismicity model, m, and LLref is the log-likelihood computed for a seismicity model with uniform rates in all cells.

Forecast earthquake counts are made for the four selected completeness levels and magnitudes for the historical catalogs, considering two historical and statistical completeness’s (see Table 1), and for the instrumental catalog.

We also varied the magnitude threshold of the forecasting catalog to test whether including small earthquakes results in greater predictability of future Mw ≥ 4.5 earthqua

kes. and magnitudesThen the correlation smoothing distance and the neighbor-number are optimized using four different completeness levels, and threshold magnitudes which are defined separately for the instrumental and the historical catalogs. We also changed the testing periods to investigate the robustness of the results and constructed from the most recent 5, 10 and 15 years of the catalogs (Table 3a and b). Although the procedure is followed over the two catalogs with different magnitude thresholds, the same testing periods for five-year of 2011-2015, ten-year of 2006-2015 and fifteen-year of 2001-2015 are performed over the two catalogs. We counted 18, 345 and 59 earthquakes Mw ≥4.5 occurred in Italy every 5, 10 and 15 years used for the testing catalogs. We determined the optimal smoothing parameters by choosing the value that maximized the likelihood that the forecast rate models gave rise to the spatial distribution of the of the testing earthquakes in the testing catalog. given the smoothed spatial density estimated from forecasting catalog through the information gain parameter. Events in the forecasting catalogs are used to compute the smoothing parameters and the predicted forecast earthquake density in all bins. The number of earthquakes occurs within the Italian territory examined and respective smoothing distances for these completeness levels, identified from the smoothed seismicity model corresponding to the highest information probability gain value (Eq. 6) for observed earthquakes with magnitudes of Mw ≥ 4.5 and applied to all forecast models. We varied only the minimum magnitude for events in the forecasting catalog fixing the minimum magnitude Mw4.5 for events in the testing catalog. Variation in information probability gain is presented with optimized correlation smoothing distance and neighbor number for the fixed and the adaptive smoothing, respectively. The results are summarized in Table 3a and b, which presents the information probability gain values from the Mw ≥ 4.5 testing catalogs as a function of testing period, forecast catalog, and of the type of the smoothing parameter for the adaptive and fixed smoothing models (Fixedhist and Fixedstat).

Results

Optimized smoothed models and parameters and models

The highest information probability gain value calculated for the adaptively smoothed seismicity models arises from the use of the most recent completeness level; Mw ≥ 3.0, from 1981 in the case of the instrumental catalog. However, it is different in the case of the historical CPTI15 catalog. Note that the information probability gain values are sensitive only to the spatial distribution of seismicity rates relative to the spatial distribution of earthquakes observed in the testing catalog since the seismicity rates are normalized in the log-likelihood calculations (Eq. 4). It is observed that the information probability gain values over the four completeness levels are effected affected by either its regency or by the magnitudes composing the catalogs.

The earlier part of the instrumental catalog with the larger-magnitude completeness levels (Mw ≥ 4.5, 1871-2000) results in a relatively smoother seismicity models because there are fewer earthquakes at this completeness level. For example, neighbor number is larger (N ~5-6) with respect to those obtained from the recent part of the instrumental catalog with relatively smaller threshold magnitude (Mw ≥ 3.0, 1981-2000) result with the smaller neighbor numbers (N ~2-3) (Table 3). The information probability gains are maximized for the Mw ≥ 4.5 testing catalog, using the second nearest neighbor (N=2) with smaller earthquakes from the longer completeness level of the instrumental catalog (Figure 3a). Information Probability gain is uniformly higher for the 10 years period Mw ≥ 4.5 target catalog for all completeness levels than for the 5 and 15 years period of target catalog (Table 3a). However in general including smaller earthquakes increase the probability gain for each three period of target catalogs compare the larger events.

In the case of the historical catalogs the situation gets contrary; the model obtained from the larger-magnitude completeness levels (M ≥5.2, 1787-2000) relatively better predicts the observed seismicity rates for M ≥4.5 events in a 5 year window from 2011 to 2015 both for the historical and the statistical completeness catalog. However the information probability gain decreases with increasing correlation smoothing distance for the 10 and 15 years target catalog. In general results show that more recent and short period of earthquake catalogs produce higher informationprobability gain values. The largest Gm is obtained as 4.4923 with 30 km correlation smoothing distance for the forecasting catalog from Mw ≥ 5.2 threshold magnitudes. In general fixed smoothed seismicity seems to provide more stable smoothing parameters with correlatiosmoothing distances around 30- 35 km than the adaptive smoothing neighboring numbers in which result with some fluctuations between 2 and 8. Adaptive smoothing distances are smaller than the fixed smoothing distances with little smoothing where density of seismicity is high while adaptive smoothing distances are larger than the fixed smoothing distances with large degree of smoothing where density of seismicity is low. Moreover, others (e.g., Moschetti, (2015) have highlighted demonstrated that the rate comparisons in the context of the information probability gain reveal the inherent trade-offs that occur with spatial smoothing methods; the broader kernels may better forecast a region encompassing more future catch the future earthquakes at the expense of the accuracy achieved by narrow smoothing kernels. better but reduce the accuracy of the spatial distribution of rates whereas the narrow kernels miss some of the future earthquakes but predict the average rates better.

Therefore in order to better predict forecast better the spatial distribution of earthquakes and earthquake rates we chose both the adaptive smoothed seismicity model that results from 1) the larger number of small magnitude earthquakes (Mw ≥ 3.0) from the instrumental catalog and 2) the fixed smoothed seismicity models that takes into account longer catalogs with larger magnitudes. This latter is constructed not only from the catalog M>5.2 but through the Weichert (1980) approach considering the magnitude-time completeness, including each earthquake in seven six macro zones, and and counted as of equal weight in the rate calculations. Thus, the seismicity rates during time periods that have more countable earthquakes (for example; time periods with lower completeness thresholds) will be affected much more strongly than the seismicity rate at other times (Felzer and Cao, 2007; Felzer, 2008).

Comparisons with the CSEP Italian models and Ensemble models

Figure 4 and 5 shows the seismic earthquake rate model for Mw ≥ 4.5, which was calculated using the spatially smoothed locations of Mw ≥ 4.5 earthquakes with a fixed 30 km smoothingcorrelation distance. The smoothing parameter maximized the information probability gain calculated for the Mw ≥ 4.5 testing catalogs. The optimized smoothed seismicity model, respectively the one based on historical and statisical completeness, forecasts an annual rate of 5.0 and 5.44.8 for Mw ≥ 4.5, 1.6 and 1.81.59 for Mw ≥ 5.0, and 0.17 and 0.220.19 for Mw ≥ 6.0. Figure 6 shows the seismic earthquake potential rate model for Mw ≥ 4.5, which was calculated using the spatially smoothed locations of Mw ≥ 4.53.0 earthquakes computed from the adaptive smoothing model with the second nearest neighbor number. The expected total annual number of earthquakes for Mw ≥ 4.5,and ≥ 5.0 and ≥ 6.0 events are 4.37 , 1.3 and 0.131.78. The general spatial patterns of seismicity are similar to those from fixed smoothed seismicity especially where the seismicity distributed along the belt of the central Apennines. Although the total seismicity rates are quite similar

both for the adaptive and smoothed seismicity models, seismicity rates near high density of seismicity are higher in the adaptive smoothed model than fixed smoothed seismicity model and then greatly reduced through sparse seismicity regions.

We compare these selected smoothed models with the two best performingprevious CSEP models in Italy. To compare models, we use the official CSEP cataloginstrumental catalog with magnitude Mw ≥ 4.95, depth < 30 km, as a testing catalogue from 2010 to 20145, 56 years window length using log-likelihood and spatial log-likelihood methods (Table 4); we also perform the same computation by including the last central Italy seismic sequence, i.e. with a catalog from 2010 to April 2017.. The log-likelihood is computed by using equation (7); the spatial log-likelihood (LL???spatial/N) uses the same equation (7) and by summing the values along magnitude bins and finally normalized it to match the number of observed earthquakes following the procedure described by Zechar et al., (2010). We calculate the spatial distribution of the earthquakes up to year 2009 for our two models then compare those with the best time-dependent, HAZ-BPT (Akinci at al., 2009) and the best time-independent, MPS04a (Meletti et al., 2014) CSEP models recalibrated for a 65 and 7.33 years comparison. Since the CSEP model was build using the Ml magnitude instead the Mw one, we compare results only for the spatial-loglikelihood, not for the loglikelihood, to have a fair comparison.

In table 4 the likelihood scoretests indicate that both models constructed in this study perform similarly. Hence, we conclude that, at the present state of knowledge, both earthquake rate models may be considered equally skilled. Therefore we define the best model as the ensemble model that represents the average of the two models. Moreover, we know that the CPTI15 catalog is longer than the instrumental one, but some recent experiment, e.g. the RELM California experiment, demonstrate that model that use smaller earthquake (Helmstetter et al, 2007) have better performance to forecast larger earthquakes (Zechar et al., 2013). We then merge the fixed seismicity models (Fixedhist and Fixedstat) with the adaptive seismicity model by using equal equal weights (0.5 for each model, linear combination). In order to check whether the ensemble model is more skilled than any of the single models, we applied the same likelihood score tests described in the previous section. The results of the likelihood test are reported in Table 4, showing the higher skill leading a better performance of the ensemble model respect to the single ones in longer experiments (Taroni et al. 2014). Also in this case, by retrospective comparing our ensemble to the single models with data from 2010 to 2015, we find an improving of the performance (Table 4).

Note that the ensemble slightly outperforms, in terms of log-likelihood, the models that compose it. In general, all our ensemble three modelss (Fixed, Adaptive and Ensemble) outperform the best time independent model of the CSEP Italian experiment (MPS04a), but the time-dependent HAZ-BPT model is even better.

Finally we produce six ensemble models merging the adaptive and fixed smoothing models using different percentage combinations as given in Table 5. We also associated each of those with different weights in order to take into account the uncertainties in the model constructions. The final total seismicity rates presented by six ensemble models match quite well to the historic rate of earthquakes from the earthquake catalog of the CPTI15 and present the model rate uncertainties (Figure 7ab).

The one of the final ensemble earthquake rate models for Italy from smoothed seismicty, ensemble model-3 is presented in Figure 8 and as the gridded, annual rate of Mw ≥ 4.5 earthquakes, available in the electronic supplement together with other ensemble models (Supplemented Materials). This ensemble The final smoothed seismicity model forecasts an annual rate of 1.471.68???? for M ≥ 5.0 earthquakes; this rate corresponds to one Mw≥4.5.0 earthquake approximately every 8 monthsone and half years. The higher rates in Figure 8 mainly concentrate along the Apennines chain, eastern Alps, Calabria and eastern Sicily, Messina areas.

The mean occurrence probabilities to have at least one earthquake with Mw≥ 4.5 in the next fifty years range from 20% to 30% in those regions under the Poisson model (Figure 9). Eastern Alps show clearly high earthquake occurrence probability with respect to the rest of the Alpine belt. This is consistent with the recent tectonic activity of the central Mediterranean region, where the Adriatic plate moves northward, creating compression along the Alps. All models consistently define the axial belt of the central Apennines as one of the most hazardous regions in Italy. This derives from both the frequent M ~ 6 earthquakes that have occurred in the last few decades and from some large (6.5-7.0) infrequent events (2016, 1703, and 1915). It is worth noting that many large earthquakes (ML 6.5–7) of past centuries have occurred in the southern section of the Apennines and in Calabria (Akinci et al., 2016). Eastern Sicily shows medium probability values in all three models. It is worth remarking that at least two large historical earthquakes did occur in this area. Central and western Sicily, where some moderate events occurred in the Tyrrhenian offshore and in the Belice area (as in 1968) show low to moderate probabilities for the next decade.

Discussion and Conclusions

We estimate the probabilities of future earthquakes of magnitudes M≥4.5 in Italy using an a set of ensemble smoothing seismicity models in Italy. We consideredestimated the spatial distribution of future seismicity by smoothing the locations of past earthquakes listed in two Italian catalogs: short instrumental catalog, M≥3.0 and longer historical catalog, M≥4.5. Therefore we used also small earthquake to identify the active fault structures that might be locations of the larger events. For the adaptive model based on shorter instrumental catalog with small magnitude of events, the average number of events Mw4.5 is 4.34.71 per year, whereas 5.013 is the yearly rate for the longer historical catalog with larger magnitude of events, calculated from fixed smoothed seismicity model. Results indicate some probable potential zones for future earthquakes of higher magnitudes in Italy, like the central Apennines, Friuli, the Irpinia region, and eastern Sicily with about 20-30% probability of occurrence within 100 km2 areal zone in the next fifty-year.

The bandwidth of the fixed and adaptive spatial kernel is estimated by optimizing the predictive power of kernel estimate of the spatial earthquake density in pseudo-retrospective forecasts. We used also small earthquake to identify the active fault structures that might be locations of the larger events. Three different catalog durations (5, 10 and 15 years) are used to calibrate these models on two catalogs. We observed that the optimal smoothing parameter is not varied for different target periods but it varied substantially for different forecasting periods with different threshold magnitude especially for the adaptive spatial kernel. Our results show that more recent earthquake catalogs with lower Mc≥3.0 values produce higher probability gain values (4.46) for the forecasting catalog between 1981-2000 while it is slightly lower, 4.19 for the same forecasting catalog with events Mc≥3.5. Meanwhile,time other research (Werner et al., 2010) has also highlighted that the demonstrated that the insufficient information overvarying Mc values among the different completeness levels precludes resolving the effects of minimum magnitude and the recency of the catalogmost recent seismicity: it is not possible to and threshold magnitudes in the forecasting catalogs which make hard to determine whether the improved forecasts result because they use recent seismicity better foreca

sts the location of earthquakes or whether theor lower threshold minimum magnitudes of completeness. from more recent completeness level and result reduction in smoothing distances with higher information gain. There exists trade-offs between two methods the use of broad smoothing kernels produces forecasts that catch more future earthquakes but reduce the accuracy of the spatial distribution rates, whereas the use of narrow smoothing kernels produces forecasts that miss some future earthquakes but better predict average rates. In this study we propose a new approach demonstrating that the ensemble models that combine different seismicity models may perform better than single model in longer periods.

In the present study, the probabilities of future large earthquake occurrence of future large earthquake in the next fifty years is are calculated based on the assumption that earthquake processes have no memory, i.e., the occurrence of a future earthquake is independent of the occurrence of previous earthquakes from the same source. Although this is the most widely used hypothesis in probabilistic seismic hazard analysis, its assumption is not physically valid for individual fault sources, given that the process of stress builds up and its release is inherently time dependent. Moreover, other factors such as clustering, static-elastic fault interactions, dynamic-stress changes and viscoelastic stress transferred from earthquakes on nearby faults can also influence the short-term and long term probabilities for earthquake recurrence. Important questions remain for the development of future seismicity rate models that the appropriateness of the Poisson probability distribution for characterizing earthquake occurrence. Even though it is still a challenge to answer what probability distribution and occurrence model best describes earthquake behavior it is important to discriminate the available models that differ significantly from one another, and to evaluate differences in a rigorous and quantitative way.

Acknowledgments

This study is under the framework of the Mappa di Pericolosita Sismica, MPS16 Project supported by Centro di Pericolosita Sismica (CPS), Instituto Nazionale di Geofisica e Vulcanologia. We thank Warner Marzocchi, Carlo Meletti, Bruno Pace and Francesco Visini for their contribution and collaboration under the framework of the MPS16 project.

**...(download the rest of the essay above)**