Seismogenic zonation and seismic hazard estimates in a Southern Italy area (Northern Apulia) characterised by moderate seismicity rates

The northernmost part of Apulia, in Southern Italy, is an emerged portion of the Adriatic plate, which in past centuries was hit by at least three disastrous earthquakes and at present is occasionally affected by seismic events of moderate energy. In the latest seismic hazard assessment carried out in Italy at national scale, the adopted seismogenic zonation (named ZS9) has defined for this area a single zone including parts of different structural units (chain, foredeep, foreland). However significant seismic behaviour differences were revealed among them by our recent studies and, therefore, we re-evaluated local seismic hazard by adopting a zonation, named ZNA, modifying the ZS9 to separate areas of Northern Apulia belonging to different structural domains. To overcome the problem of the limited datasets of historical events available for small zones having a relatively low rate of earthquake recurrence, an approach was adopted that integrates historical and instrumental event data. The latter were declustered with a procedure specifically devised to process datasets of low to moderate magnitude shocks. Seismicity rates were then calculated following alternative procedural choices, according to a “logic tree” approach, to explore the influence of epistemic uncertainties on the final results and to evaluate, among these, the importance of the uncertainty in seismogenic zonation. The comparison between the results obtained using zonations ZNA and ZS9 confirms the well known “spreading effect” that the use of larger seismogenic zones has on hazard estimates. This effect can locally determine underestimates or overestimates by amounts that make necessary a careful reconsideration of seismic classification and building code application. Correspondence to: V. Del Gaudio (delga@geo.uniba.it)


Introduction
Apulia region is the south-eastern end of the Italian peninsula and is constituted by an emerged portion of the Adriatic microplate, representing the foreland-foredeep area of the Southern Apennine chain.The northernmost part of Apulia, located between the Ofanto river and the Fortore river basin (Fig. 1), is occasionally affected by seismic events of moderate energy and has been historically hit by strong earthquakes which in some cases caused disastrous effects with thousands of fatalities.Recent estimates of seismic hazard conducted in Italy at national scale (Gruppo di Lavoro "Mappa della Pericolosità Sismica", 2004) were obtained through procedures based on Cornell (1968) approach, adopting a new subdivision of the Italian territory in seismogenic zones, the so called zonation ZS9.It was derived by modifying previous zonations to take into account advance in active tectonics knowledge in Italy and data derived from the most recent earthquakes, but also to solve the problem of the low number of events reported by seismic catalogues for several small seismogenic zones.The data shortage did not allow to well constrain the estimates of the seismicity rates, i.e. the number of events of different magnitude expected in a fixed time.Therefore in zonation ZS9, to enlarge the statistical bases of such estimates, areas previously belonging to distinct zones with relatively similar seismotectonic properties were joined together.
Some aspects of the ZS9 zonation are controversial.In particular, with reference to Northern Apulia, a recent study (Del Gaudio et al., 2007), on the basis of an integrated analysis of the characteristics of historical and instrumentally detected seismic events, raised doubts on the seismic homogeneity of the area included in the zone labelled as no.924: this zone extends for over 100 km crossing, from west to east, the Apennine chain front, the foredeep and the foreland (Fig. 1).In the aforementioned study its identification as a unique continuous strike-slip fault system responsible for major historical earthquakes was questioned on the basis of indications and constraints provided by data analysis and a subdivision of the area into zones with a differentiated seismic behaviour was consequently proposed.
In the present study such subdivision was assumed as basis for a locally more detailed seismogenic zonation to reevaluate the seismic hazard of Northern Apulia.The local adoption of smaller zones re-proposed the problem of the limited number of events reported by historical earthquake catalogue for each of such zones because, even though this area was historically affected by strong earthquakes, their temporal recurrence is not very frequent.To solve the problem of the weak constraints that historical seismicity provides to the assessment of seismicity rates, we adopted an approach integrating historical catalogue data with a database of low to moderate magnitude events, derived from the instrumental monitoring of seismicity during the last two decades.In particular the historical and instrumental data were integrated to obtain the coefficients of a unique Gutenberg-Richter (1944) relationship.The integration of these data required the application of declustering techniques to the instrumental data.
Seismic hazard estimates were then obtained by using the code SEISRISK III (Bender and Perkins, 1987).In seismic hazard assessment the outlining of the seismogenic zones has a considerable influence on the determination of seismicity rates.We were particularly interested in evaluating the in-fluence that the proposed local modification of seismogenic zonation would have on the hazard estimates for the study area and whether this influence is significant in comparison to the effect of other epistemic uncertainty factors affecting seismicity rate calculation.For this purpose estimates of seismicity rates based on integrated historical and instrumental datasets were carried out both for the Northern Apulia seismogenic zones defined in the new zonation and for the zone no.924 of the ZS9.The effect of other epistemic uncertainty factors was examined by parallely adopting alternative procedural choices, according to a "logic tree" approach (Kulkarni et al., 1984).Seismic hazard estimates were then obtained by using, for Northern Apulia zones, the seismicity rates calculated and, for the closest outer zones, those reported by the latest national scale hazard study (Gruppo di Lavoro "Mappa della Pericolosità Sismica", 2004).Finally, comparisons were carried out between the results obtained by using the ZS9 zonation and that locally modified.

Geological and seismotectonic setting
From a structural geological point of view, Northern Apulia consists of three different zones (Fig. 1): a) A foreland area to the NE constituted by the Gargano promontory, a horst elongated towards the Adriatic sea, generated by the uplift of a carbonate plateau and delimited by steep scarps.This is the most elevated part of a carbonate platform that towards the inland sinks under the front of the Apennine chain.
b) A large central alluvial plain named "Tavoliere", constituting a local enlarged section of the Apennine chain foredeep, characterised by soft sediment deposits overlying the carbonate platform.
c) The marginal external front of the Southern Apennine chain, named Dauno Sub-Apennine, consisting of northwest-southeast oriented thrusts of tectonically deformed turbiditic formations.The thrust belt front rests on a clastic wedge which, in turn, overlays the marginal part of the Apulian foreland carbonate platform.
Despite the structural similarity to the rest of the Apulia region, mostly constituted by the same carbonate platform outcropping in the Gargano promontory, seismicity of its northern part, documented by historical data and instrumental monitoring, appears much more active.It includes at least three events (1361,1627,1731) that caused effects up to X degree on the Mercalli-Cancani-Sieberg (MCS) scale and casualties in the order of thousands, but severe damages and an uncertain number of victims were reported also for other events (e.g. in 1223, 1414, 1646).In recent years seismic shocks that caused slight damages were recorded in the Gargano area, like the event of moment magnitude M w =5.2 that hit the north-eastern part of promontory on 30 September 1995.On the contrary, south of the Ofanto river no event has caused damaging effects comparable with those of major Northern Apulia earthquakes and also instrumental seismicity shows much lower rates of seismic event occurrence.Figure 2 shows the spatial distribution of historical and instrumentally detected earthquakes, analysed in recent studies (Del Gaudio et al., 2005, 2007) that provide more details on the characteristics of regional seismicity.
In the regional differences of seismic behaviour an important role is played by a spatial variation in structure and thickness of the lithosphere between northern and southern part of the Adriatic microplate (Venisti et al., 2004(Venisti et al., , 2005)), in correspondence of a belt crossing the central Adriatic sea from the Gargano promontory to the opposite Croatian coast.Along this belt a concentration of intra-plate seismic activity has been observed, possibly as consequence of the structural weakness represented by the mentioned structural heterogeneity, which can determine a focusing of seismic energy release.
The integrated analysis of historical earthquake characteristics, instrumental seismic data and geological structural elements, conducted in a previous study (Del Gaudio et al., 2007), pointed out significant differences also within Northern Apulia.According to this study, the foreland area is characterized by a major concentration of events compatible with a transpressive stress field having an approximately NW-SE compression axis (P) and NE-SW extension axis (T).The spatial distribution and energetic characteristics of event foci suggested a possible differentiation between the southeastern part of the foreland, constituted by the core of the Gargano plateau, and the northernmost part, partially submerged by the sea, extending through the area of the Tremiti Islands, the coastal Lesina lake and the mouth-lower course of the Fortore river.In comparison to the Gargano plateau, the latter area is characterized by relatively shallower events and a 30 • -40 • counter-clockwise relative rotation of regional T and P axis.This area was hit by the strongest earthquake historically documented, i.e. the magnitude 6.7 event of 30 July 1627, which affected the Fortore-Lesina area and generated a tsunami which submerged the town of Lesina.
If compared to the foreland, the foredeep area appears relatively less seismically active and the regional stress field shows a transitional character towards the prevailingly NE-SW extensional regime of the Appennine chain domain.Indeed the focal mechanisms of events in south-central Tavoliere plain, even having still a prevailingly strike-slip nature, show an accentuation of the relative weight of NE extension with respect to NW compression, probably as effect of a reduced efficiency in the transmission of axial compression along the less rigid border of the Adriatic microplate (Del Gaudio et al., 2007).Bath and Duda (1964) formula for earthquakes extracted from the CPTI 2004 catalogue (CPTI Working Group, 2004).The circles corresponding to the three major historical earthquakes are identified by the year of their occurrence.

Seismogenic zonation
The aforementioned seismotectonic data suggested the possible identification in Northern Apulia of four separate seismogenic zones (Fig. 1): two foreland areas, corresponding to the Fortore lower course -Lesina lake -Tremiti Islands (Zone 1) and to the Gargano promontory (Zone 2), a foredeep area corresponding to the Tavoliere plain (Zone 3) and the external front of the Apennine chain corresponding to the Dauno Sub-Apennine (Zone 4).
To support this hypothesis of zonation (henceforth named ZNA) under the aspect of the temporal pattern of seismic energy release, this was comparatively analysed for all the four zones by examining a catalogue of events recorded for 20 years from 1985, generated during a previous study (Del Gaudio et al., 2007).To avoid a possible bias due to temporal change in event list completeness, the minimum magnitude was evaluated for which the catalogue can be assumed complete: this threshold was identified as the minimum magnitude for which the diagram of the logarithm of the cumulative number N of events exceeding a local magnitude M L deviates from the linear decrease expected according to the Gutenberg-Richter relationship (Gutenberg and Richter, 1944) (Fig. 3).On this basis a magnitude value of 1.9 was assumed as completeness threshold.(Castello et al., 2006): deviation from linear decrease at low magnitudes was used to recognize completeness threshold.
The monthly cumulative energy released per unit area was calculated for each of the four zones and also for the entire study area through a relation proposed by Gutenberg and Richter (1942): where E is the energy released in joules and M L is the local magnitude of events.
On the whole, during the considered 20 years, the seismic energy released through the overall study area was equal to 8.4•10 12 joules, however its spatial and temporal distribution was extremely irregular and variable for the different zones.About 60% of this energy was released by a seismic sequence including two major shocks of M L =5.4 and 5.3, which in 2002 hit the left side of the Fortore river middle course, at the border between Molise and Apulia regions.It was located in Zone 4, which, on the other hand, before that event had released only 0.2•10 10 joules through an average of 3-4 events per year, mostly of magnitude 2.1-2.4,with an energy release rate of 0.4•10 5 joules/km 2 •year −1 .Thus this area was practically quiescent before to be hit by the 2002 earthquake.The seismic record of this zone (both historical and instrumental) reports only another event of similar energetic characteristics, i.e. an earthquake of intensity IX MCS occurred in 1125.These observations suggest that Zone 4 is characterised by a very weak seismic activity punctuated by rare moderately energetic events having rather long mean return time.
A different behaviour is shown by the other 3 zones which present a more frequent activity, even though with different rates: in the examined 20 years, these zones have released seismic energy with rates per unit area and per year respectively of 0.15 (Zone 1), 4.48 (Zone 2) and 0.08•10 7 (Zone 3) joules/km 2 •year −1 .Examining the time variation (Fig. 4) one can see that in Zone 3 the energy release per unit area has had a relatively constant rate and is systematically lower than the regional average (black line in Fig. 4).In the other two zones the energy release appears to proceed through major bursts occurring episodically, which in some periods confer to these zones energy release higher than the regional average.For Zone 1, two thirds of the total energy were released from 1986 to 1990 through events of magnitude around 4.0 isolated or included in a sequence, occurred in the area between the Lesina lake and the Tremiti Islands.For Zone 2, 95% of energy was released through the seismic sequence of 1995 (with a main shock of M L =5.4).As a consequence of this temporal pattern, these two zones have exchanged the role of the area releasing most of the seismic energy (Fig. 4).
On the whole also this analysis supports the idea that the four outlined zones may have significant differences in the pattern of seismic energy release, which justifies their separation in the seismogenic zonation used for the following tests.

Methodology
The influence of the proposed zonation modification on seismic hazard estimates was evaluated by calculating the spatial distribution of a parameter commonly used in hazard assessment (see Giardini, 1999), i.e. the peak ground acceleration having an exceedence probability of 10% in 50 years (PGA 0.10, 50 ), which is a usual reference for seismic building codes.PGA 0.10, 50 values were calculated through a standard program (SEISRISK III, Bender and Perkins, 1987) on a grid of points spaced by an angular distance of 0.05 degrees both in N-S and in E-W direction.The calculation procedure requires the definition of the seismicity rates for each seismogenic zone relevant for the assessment of hazard in the study area.Furthermore, an attenuation relationship is used to provide probabilistic estimates of PGA expected at a given distance for an earthquake of given magnitude: at this regard we adopted the relationship obtained by Sabetta and Pugliese (1996), calibrated on an Italian accelerometric database.
The calculation of seismicity rates poses some problems.Hazard estimates conducted in Italy to provide reference for building code have been based on a quite rich historical documentation available in form of earthquake catalogues spanning approximately the last 1000 years.However the reduced size and the relatively moderate rate of earthquake occurrence characterising the seismogenic zones considered in this study would lead to derive seismicity rates from small numbers of historical events reported, which might provide unreliable results.
To overcome this problem we adopted a procedure which integrates historical and instrumental data to constrain the coefficient a and b of the Gutenberg-Richter law (1944): where N (M) is the number of events of magnitude equal to or larger than M occurring in a given area during a fixed time interval.The instrumental data recorded in the last two decades provide constraints for this relationship at low magnitudes (down to 2-3), whereas historical data do the same at the highest magnitudes observed in each zone (∼6 or more).
Integrating both kinds of data one can obtain better constrained values of a and b, which can be used to estimate the number of events expected also at intermediate magnitude (4-5) for which both historical and instrumental catalogue might not provide reliable estimates of seismicity rates (the historical catalogues for incompleteness, the instrumental ones for temporal shortness).Since a basic assumption in the application of SEIS-RISK III code is a Poissonian model of seismic event time recurrence, earthquake historical catalogues developed in Italy for hazard assessment were declustered by removing aftershocks and foreshocks and including only main shocks (Slejko et al., 1998).
Therefore, to integrate instrumental data with historical ones, a declustering procedure needs to be applied also to the instrumental catalogue.A standard procedure used for the identification of cluster of interdependent events is that proposed by Reasenberg (1985).The procedure is based on the identification of an "interaction zone", i.e. a space and a time span around the location and time of each event, calculated as function of event magnitude, such that following events occurring within these space-time limits are considered as events "stimulated" by the previous one, that is as "genetically dependent" events belonging to a same cluster.Generally declustering procedures proposed in literature have been developed for catalogues in which main shocks are expected to have moderate to high magnitude (e.g.not less than 4.0), so that their effectiveness, when applied to low magnitude events, cannot be taken for granted.Indeed, preliminary tests carried out on a catalogue of Northern Apulia recent events showed that the results of declustering with the Reasenberg procedure have still an excess of short inter-event times (and a deficit of longer ones) in comparison to what expected for a Poissonian distribution having the same mean inter-event time (Fig. 5).Therefore as alternative procedure we tested a new purposely developed technique.It is based on the assumption that the sequence of main shocks must have the properties of a Poissonian process, which implies that the probability P that one or more independent events occur during a time interval t is: (3) where T represents the mean interval between successive events (mean inter-event time).
The proposed method adopts an iterative procedure implemented in a code named DECLPOI, consisting of the following steps: 1. the cumulative frequency distribution F ( t) of each the inter-event times t i observed in the catalogue is compared with the value P ( t) expected, according to the Eq. ( 3), for a Poissonian distribution having the same mean inter-event time of the catalogue and the t i ' value is found for which the difference F ( t i )−P ( t i ) is maximum; 2. for all the couples of events whose t i is smaller than or equal to t i ', a space-time distance d ST is calculated according to the criterion proposed by Davis and V. Del Gaudio et al.: Seismogenic zonation and seismic hazard in Northern Apulia Frohlich (1991), which associates time separation τ and space distance d through the expression where C is a transformation coefficient of time separation into space distance and is equal to 1 km/day; if in previous iterations the two events have been identified as belonging to clusters, then d ST is calculated considering the minimum time separation and space distance between the corresponding clusters, rather than the distance between the single events; 3. the couple of events with d ST minimum is identified as belonging to a common cluster and the event of smaller magnitude is marked to be excluded in the final catalogue; then the frequency distributions of t i values (both for the catalogue and for the comparative Poisson distribution) are recalculated excluding the marked events; 4. for the t i values of the modified catalogue the coefficient of variation C V (i.e. the ratio between the standard deviation and the mean of t i values) is calculated and if larger than 1 (i.e. the value expected for a Poissonian distribution) the steps from 1 to 4 are repeated, otherwise the procedure stops.
At the end the declustered catalogue for which the quantity |C V −1| is minimum is adopted as final catalogue.Figure 5 shows that, applying such a procedure, the final declustered catalogue shows inter-event time frequency in very close agreement with what expected for a Poissonian distribution.Seismic events located in each seismogenic zone of Northern Apulia were then extracted both from the declustered instrumental datasets and from the historical catalogue and the number of the events of different magnitudes were grouped into fixed intervals.The number of events in each interval was normalised multiplying it by a factor 100/ T , where T is the temporal extension in years of the used part of catalogue, which is defined according to the result of completeness analyses carried out both on historical and instrumental catalogues.These data were used to determine the coefficients a and b of the Eq. ( 2) for each zone both by using a simple linear regression (LS) and by applying the "maximum likelihood" method (MLM) (Aki, 1965;Bender, 1983).
A recent study pointed out that the MLM technique should be preferred (Sandri and Marzocchi, 2007) because LS estimates are affected by a bias causing an underestimate of the coefficient b, particularly in case of small datasets.According to these authors, this bias is mainly caused by the logarithm transformation of the number of events used in regression and by the underrepresentation of negative fluctuation at higher magnitudes, due to the exclusion of magnitude classes having zero events.
However, considering the peculiar characteristics of the datasets used in this study, a potential source of bias might affect also MLM estimates in consequence of the underrepresentation of intermediate magnitude classes at the passage from the magnitude range covered by instrumental data to that covered by historical ones.Indeed, it is to take into account that the MLM estimates are based on the calculation of the average differences (M −M min ) between the dataset event magnitudes and their minimum value (see Aki, 1965) and provide increasing b values as such average decreases.Since the number of events diminishes exponentially with magnitude increase, possible lack of data at intermediate magnitudes, due to the temporal shortness of instrumental catalogue and to the incompleteness of the historical ones, would cause an underestimate of the average (M −M min ) and a consequent overestimate of b.
Considering that the two methods may produce results affected by opposite sign errors, they were both used to derive alternative estimate of the Gutenberg-Richter relationship coefficients.Then this relationship was used to calculate the number of events of different magnitude intervals expected in 100 years, to be provided as input to SEISRISK III.

Data processing
The adoption of different zonations has an obvious influence in the calculation of the seismicity rates of the seismogenic zones, which can considerably modify the results of hazard assessment.Therefore, one of the main purposes of the test carried out in this study was to evaluate whether the modifications of seismogenic zonation ZS9 suggested according to the indications of our previous work (Del Gaudio et al., 2007) imply significant variations of seismic hazard estimates for Northern Apulia, i.e. variations which would modify seismic classification and building code application in some part of this territory.For this purpose, the previously described methodology was used to calculate PGA 0.10, 50 values both by adopting the zonation ZS9 and, in alternative, by replacing in it the Northern Apulia zones (Zones 924 and 925) with the four zones outlined for the zonation ZNA.For the Zone 924, which largely overlaps the study area, the seismicity rates were also recalculated following the same procedure adopted for ZNA, in order to avoid the introduction of differences related to the peculiarity of the calculation methods followed (in particular the integration of historical and low magnitude instrumental data).For the Zone 925, which covers only marginally the southern part of the study area, and for external zones the seismicity rates used were those reported by the published national scale estimates (Gruppo di Lavoro "Mappa della Pericolosità Sismica", 2004).
Since at some stages of the procedure of seismicity rate calculation alternative choices could be made in data selection or methodology, epistemic uncertainties associated to these choices were explored through a "logic tree" approach (Kulkarni et al., 1984): at each stage proposing alternative choices, the calculation procedure branches to follow all the possibilities, so that at the end PGA values are obtained with any combination of procedural choices.The use of this approach allowed to compare the influence of the seismogenic zonation choice with that of other factors of seismicity rate uncertainty.

Earthquake catalogue selection
The first stage of data processing consisted in the extraction of events located in the examined seismogenic zones from historical and instrumental earthquake catalogues.With regard to historical seismicity, the catalogue CPTI04 (CPTI Working Group, 2004) was used adopting as magnitude values those reported as moment magnitude.CPTI04 is a recent version of an Italian seismic catalogue specifically designed for hazard assessment, i.e. satisfying two main requirements: i) the inclusion of only the main shocks of seismic sequences with damaging effects and ii) the estimation of time interval of completeness at different magnitude levels.About the latter aspect, the completeness intervals for different magnitudes were assumed according to the results of the estimates carried out prevailingly on the basis of historical analysis, indicated in the published national scale estimates (Gruppo di Lavoro "Mappa della Pericolosità Sismica", 2004) as the preferred criterion.
Concerning the instrumental catalogue, a dataset of 1789 events occurred from 1985, relocated with a local velocity model in a previous work (Del Gaudio et al., 2007) was taken into consideration.Local magnitudes M L of these events were derived from three sources: the catalogue CSTI -"Catalogo Strumentale dei Terremoti Italiani dal 1981 al 1996" (Gruppo di Lavoro CSTI, 2001), reporting events recorded up to 1996; the catalogue CSI -"Catalogo della Sismicità Italiana 1981-2002" (Castello et al., 2006) covering a time span up to 2002; the seismic bulletins published by the Italian "Istituto Nazionale di Geofisica e Vulcanologia -Centro Nazionale dei Terremoti" (INGV -CNT), available online at: ftp://ftp.ingv.it/pro/bollet/.
The catalogue CSI was the result of an extension and revision of the previous catalogue CSTI: the revision affected particularly the magnitude attribution to seismic events occurred until 1996, which, according to the results obtained by Gasperini (2002), appears to have been overestimated for magnitude smaller than 2.5 and underestimated for those larger than 5.0.Even though the most recent catalogue should be considered more reliable, we assumed the differences between the two catalogues as representative of uncertainties affecting magnitude estimate methods, thus we carried out parallel seismicity rate calculations using both sources, in order to evaluate the influence of such uncertainties on the final results.Therefore two separate datasets were prepared, one based on CSTI and the other on CSI, both integrated by the INGV -CNT data for the period until 2004.The choice of using only events from 1985 was motivated by the outcome of the cited study (Del Gaudio et al., 2007) which showed as the poor seismic network coverage in Southern Italy existing before that year makes the catalogue rather incomplete and affected by large location uncertainties.
A completeness analysis of both datasets was carried out using two methods, i.e.: 1. the analysis of the deviation from linearity expected for log N (M) according to the Eq. ( 2) (as made before in Sect.3, preliminarily to the energy release estimates): such deviation at low magnitudes is considered to reflect dataset incompleteness; 2. the examination of the slope change in the cumulative number of events as function of time, for different magnitude thresholds, according to the method proposed by Tinti and Mulargia (1985).
Both methods converge to indicate that the seismic datasets of Northern Apulia are complete from 1985 for magnitude greater than or equal to 2.0 or 1.9 (depending on whether CSTI or CSI magnitudes are adopted).

Estimates of seismicity rates
To estimate the seismicity rates of each zone, data from the historical catalogue and from the instrumental datasets were integrated after having declustered the latter.Two declustering procedures were tested as alternatives, i.e. the technique reported by Reasenberg (1985) and the new method described in Sect. 4 (code DECLPOI).Declustering was preliminarily carried out on the whole instrumental dataset of Northern Apulia and then separate datasets for each zone were obtained cutting away the events with magnitude below the previously found completeness thresholds (M L <2.0 and 1.9, respectively for datasets with CSTI-and CSI-based magnitudes).
Table 1 shows the results of the declustering procedure in terms of number of cluster identified, events left and events removed from the catalogue for the two datasets and the two techniques adopted.One can notice that the DECLPOI procedure removes 2-3 times more events and identifies 5-6 times more clusters than the Reasenberg method.Most of these clusters consists in couple of events occurred at close location and within a short interval of time (from few hours to few days), so that the number of events per cluster turns out much lower for DECLPOI (approximately one half in comparison to Reasenberg method).This result suggests that DECLPOI is more efficient in recognising inter-dependence between couple of events of low magnitude, which could belong to clusters whose other events, having a still lower magnitude, may have not been recorded or located.Furthermore DECLPOI is less sensitive to difference of magnitude estimates: for the two magnitude sources used (CSTI and CSI) the number of events left in the declustered datasets differs www.nat-hazards-earth-syst-sci.net/9/161/2009/Nat.Hazards Earth Syst.Sci., 9, 161-174, 2009 only by 4 units on a total of ∼600 events, since magnitude affects the removals only in case that difference of estimate causes an exchange between the removed event and the left one in a couple of inter-dependent shocks of similar magnitude.Comparatively Reasenberg technique shows a larger differences in the number of events removed (about 5-6% of the events left in the declustered dataset) because magnitude affects space-time extension of the "interaction" zone.Declustered and historical datasets were then used together to calculate the coefficients a and b of the Gutenberg-Richter relation (Eq.2) both with least squares (LS) and "maximum likelihood" method (MLM).For this purpose events were grouped into magnitude classes at intervals of 0.2.The number of events for each class, derived from historical data at magnitudes larger than 4.5 and from instrumental data at lower magnitudes, was normalised on a time interval of 100 years.As previously specified, estimated moment magnitude values were considered for historical data, whereas local magnitudes M L were used for instrumental events.The merging of these two kinds of magnitude is justified by the observation that at magnitude lower than 4.5 the M L values are comparable with the M W ones within the few tenths of unit uncertainty commonly affecting magnitude estimates (see Utsu, 2002).Figure 6 shows an example of the event number distribution with magnitude.These diagrams suggest that the completeness thresholds of instrumental catalogues might be slightly higher than those (M L =1.9 or 2.0) derived by analysing cumulative frequencies before declustering.However it should be taken into account that grouped magnitude frequency distribution shows a larger scattering (see Fig. 6) related to the superimposition, on the linear trend expected according to the Gutenberg-Richter law, of statistical fluctuations that cumulative frequencies tend to smooth.Such fluctuations become more prominent examining the event number included within narrow magnitude intervals, particularly when, as in the studied case, datasets contain a limited number of events.Thus the consequent data scattering can determine a misidentification of the "slope change point" defining the completeness threshold.Nevertheless, these sources of uncertainties have a limited influence on the determination of seismicity rates, considering that the Gutenberg-Richter coefficients can be quite well constrained by using frequency estimates spanning a large magnitude interval.
Table 2 shows the results obtained calculating Gutenberg-Richter a and b with both LS and MLM techniques, together with 95% confidence interval derived for least squares from standard deviation and for "maximum likelihood" by the method of Tinti and Mulargia (1987).Data declustered with the Reasenberg technique provided systematically higher values of both a and b in comparison with those processed with DECLPOI: this is clearly an effect of the larger number of low magnitude events left in the datasets by the first procedure.Furthermore, LS and MLM estimates are in good agreement for Zones 1 and 2, whereas in the other two zones LS gives significantly lower values (beyond the 95% confidence intervals) than MLM and quite similar to those of the previous zones.Considering that the less seismically active Zones 3 and 4 provided poorer datasets (particularly for the historical part covering the higher magnitude range), the discrepancies observed between LS and MLM results can be explained as a consequence of the maximisation of the diverging bias effects discussed above (see Sect. 4).
Finally, the a and b coefficients were used to calculate the seismicity rates of each zone in terms of number of events expected in 100 years at different magnitude intervals from 4.6 to a value corresponding to the maximum magnitude historically documented for each zone.The minimum magnitude considered in seismicity rate definition was chosen for homogeneity with the value adopted in national scale estimates (Gruppo di Lavoro "Mappa della Pericolosità Sismica", 2004), in view of the planned comparative analysis.However a test carried out extending the magnitude intervals down to 4.0 showed that, in the studied area, the contribution to hazard estimate from events of magnitude lower than 4.6 is negligible, causing localised increases of PGA 0.10,50 at most by 0.02 g.
The seismicity rates obtained are shown in Table 3.In Zones 1 and 2 the values obtained from Gutenberg-Richter coefficients calculated with LS are generally lower than those derived using MLM, whereas the opposite occurs for Zones 3 and 4: these differences reflect the differences in b values (influenced by the aforementioned biases), the seismicity rates resulting higher as b values are lower.Furthermore, except for the case of Zone 4 datasets processed with MLM, seismicity rates derived from datasets declustered by DE-CLPOI are smaller at lower magnitudes in comparison with data declustered with Reasenberg technique, whereas the situation tends to be inverted at the highest magnitudes: this reflects the more efficient capacity of small event removal characterizing DECLPOI.Finally the use of datasets with magnitude estimates based on the catalogue CSTI appears to produce higher rates, especially for lower magnitudes, in comparison to those based on CSI data, possibly as effect of the systematic differences of magnitude estimates between the two catalogues (see Sect. 5.1).

"Logic tree" development
The alternative choices considered with regard to instrumental dataset magnitude source, declustering procedure and Gutenberg-Richter coefficient determination were adopted through a "logic tree" approach which led to define eight possible combinations for the calculation of seismicity rates, according to the scheme represented in Fig. 7.This scheme was followed both for the seismogenic zonations ZS9 and ZNA.Thus, eight different input datasets were prepared for each zonation to calculate, through the code SEISRISK III, the PGA 0.10, 50 values in a grid of nodes covering the study area.Figures 8-9 show the map of minimum and maximum PGA values among those obtained with the eight calculation schemes both for ZNA and ZS9, respectively.
Following the approach adopted in Italy for recent assessment of seismic hazard (Gruppo di Lavoro "Mappa della Pericolosità Sismica", 2004), weighted medians were also calculated for PGA at each node of the grid, by attributing relative weights to the alternative choices at each stage of the seismicity rate calculation, and thus obtaining by multiplication a final relative weight for each combination of procedural choices.Different weighing schemes were tested ranging from the attribution of equal weights to all the choices, to weights unbalanced in favour of choices deemed preferable (i.e.CSI vs. CSTI, DECLPOI vs. Reasenberg, MLM vs. LS).The results obtained turned out quite similar: Fig. 7 shows one of the weighing scheme adopted and Fig. 10 shows the maps derived from this scheme for the two zonations.Finally, Fig. 11 shows the spatial distribution of differences Fig. 7. Logic tree scheme adopted to calculate seismicity rates with any combination of procedural alternative choices at the stages of the selection of the instrumental data magnitude source (catalogue CSTI or CSI), of the instrumental data declustering (DECLP=by DECLPOI; REAS=by Reasenberg, 1985) and of the estimate of the Gutenberg-Richter coefficients a and b (LS=by least squares; MLM=by maximum likelihood method).The relative weights associated to each procedural choice and the weight relative to each branch are those attributed to the results of each procedural sequence for the calculation of the weighted medians mapped in the following Fig. 10. between the values obtained for the two zonations, deriving them from the maps of Fig. 10.

Discussion and conclusions
The results obtained with the currently national reference zonation (ZS9) and with the one hypothesised in this work by including local modifications for Northern Apulia (ZNA) were comparatively examined.This comparison shows that the variations of PGA 0.10,50 estimates using different pro-cedural sequences are locally at most of 0.10 g.Major effects on these variations appear to derive from the choice of the method for Gutenberg-Richter coefficient determination, with LS causing, in comparison to MLM, decreases in foreland zones and increases in chain-foredeep zones by up to 0.06 g.The adoption of magnitudes from CSI catalogue, instead that from CSTI, tends to cause a generalised decrease also up to 0.06 g throughout the study area if Gutenberg-Richter coefficients are obtained with MLM, but produces minor differences (negative in foreland and posi- tive in foredeep-chain zones) if LS is the technique adopted for the calculation of a and b.The choice of the declustering method seems to have less influence on the results, causing at most differences by ±0.03 g, particularly in the two foreland zones.
In comparison to the effects of the procedural choices adopted in hazard estimates, differences related to the choice of seismogenic zonation appear however more prominent, being evident the zonation influence on the "geometry" of the spatial distribution of PGA 0.10,50 values and on the maxima reached by them.The zonation ZNA provides more pronounced maxima: in some areas (northern coast of Gargano promontory, Tremiti Islands) even the minimum PGA 0.10,50 estimate obtained with ZNA is not less than the maximum estimate obtained with ZS9 (see Figs. 8 and 9).Elsewhere in the study area (e.g. in some parts of the Dauno Sub-Apennine) the opposite can also occur, with maximum of ZNA estimates lower than ZS9 minimum: this concerns areas where seismic hazard obtained with ZNA zonation is below the regional average.
For a synthetic comparison, median values mapped in Fig. 10 can be considered.The results mapped in Fig. 10b appear in quite good agreement with those reported for the study area by the reference national hazard estimates (Gruppo di Lavoro "Mappa della Pericolosità Sismica", 2004), which confirms the compatibility of the adopted approach with the standard procedures used in that study.
The comparison between the results obtained with ZS9 and ZNA was focused on the area to the north of the Ofanto river because to the south of this limit major discrepancies are caused by the exclusion in zonation ZNA of Zone 925, marginally overlapping the southern part of Zones 3 and 4 (see Fig. 1): this determines an underestimate of hazard for the area of Zone 925 extending outside the limits of the zonation ZNA.Actually there are controversial aspects also in the definition of this zone and of its seismic characteristics: however a more specific study of the hazard for this area is beyond the scope of this paper.
Overall the PGA 0.10, 50 values obtained with the proposed zonation ZNA are not excessively dissimilar from those derived from ZS9, discrepancies being comprised in a ±0.05 g interval for most of the study area (see Fig. 11), however local significant differences can be found.These are the consequence of the well known "hazard spreading effect" of the Cornell approach, which is enhanced by the use of more ex- tended seismogenic zones: the assumption of homogenous seismicity rates in a seismogenic zone tends to "distribute", throughout the entire zone extension, the hazard associate to seismic activity observed in a delimited portion, providing more uniform hazard estimates and smoothing local differences.Therefore, the adoption of smaller zones, thanks to the integrated use of data derived from historical catalogue and instrumental datasets, can determine more territorially differentiated hazard estimates.Areas where seismic activity is more frequently observed (like the zone of Lesina lake and Tremiti Islands) show an increase of PGA 0.10, 50 , whereas a decrease is observed where seismic activity is more rare (like the Dauno Sub-Apennine zone).
Indeed, in comparison to the results based on ZS9, those derived with zonation ZNA provide increases and decreases that locally can be in the order of 0.1 or even more: the maximum difference is found for the Tremiti Islands because they are outside any seismogenic zone of ZS9, whereas are well inside an active zone, according to ZNA.Considering that the seismic classification of Italian territory is based on PGA 0.10, 50 intervals of 0.1 g, such differences, if confirmed, would imply a modification of classification for some localities of the study area.
The latest approach adopted by Italian regulations for seismic building code reduced the importance of seismic classification which, previously, fixed the scale factor of the design spectrum on the basis of a single PGA 0.10, 50 value adopted throughout the extension of a municipality.The new regulations, at present, prescribe criteria for the calculation of design spectra, based on a dense grid of calculated PGA 0.10, 50 values covering the entire national territory.This approach reduces possible abrupt variations of building criteria across the boundaries of administrative territorial units and should www.nat-hazards-earth-syst-sci.net/9/161/2009/Nat.Hazards Earth Syst.Sci., 9, 161-174, 2009 allow to better fit the spatial variation of hazard factors.However a full exploitation of this approach requires a more detailed recognition of such variations.In this regard, our results demonstrate the importance of a better comprehension of seismic behaviour of seismogenic structures and of the recognition of inhomogeneities in this behaviour.For this purpose an important role can be played by the integration of historical datasets with low energy event instrumental datasets properly processed (e.g. through declustering procedures devised specifically for this kind of data).
It is opportune to stress that caution should be adopted in integrating historical data with instrumental ones, particularly if temporal coverage of the instrumental catalogue is limited.Actually, low energy seismicity might show significant variations in time and it is not easy to recognise whether the time interval of the available datasets is sufficiently representative of long term behaviour.However a reasonable consistency in seismicity rates observed for strong historical events and for small instrumental shocks (e.g. in terms of alignment of magnitude frequency distribution according to the Gutenberg-Richter relationship: see Fig. 6) can make one more confident on the obtained results, even though a reliable confirmations will be possible only with a long term continuation of seismicity monitoring.

VFig. 1 .
Fig.1.Location of the study area and delimitation of the seismogenic zones according to the zonation proposed for Northern Apulia (solid lines) or the ZS9 zonation (dashed lines): 1=Fortore lower course -Lesina lake -Tremiti Islands; 2=Gargano promontory; 3=Tavoliere plain; 4=Dauno Sub-Apennine.The two ZS9 zones outlined are the 924 (to the north) and the 925 zone (to the south).

VFig. 3 .
Fig. 3.Diagram of cumulative number of events located in Northern Apulia as function of magnitude provided by catalogue CSI(Castello et al., 2006): deviation from linear decrease at low magnitudes was used to recognize completeness threshold.

Fig. 4 .
Fig. 4. Cumulative seismic energy release per unit area as function of time in the four seismogenic zones outlined and, globally, in the whole study area.

Fig. 5 .
Fig. 5.Cumulative frequency distribution of inter-event times of the original dataset and of those declustered with the code DECLPOI (DECLP) and with theReasenberg (1985) technique (REAS).Each frequency distribution is compared with that expected for a Poissonian distribution having the same mean inter-event time T (reported in brackets).

ZONEFig. 6 .
Fig. 6.Diagrams of grouped magnitude frequency distribution normalised on 100 years for Zone 1: (a) and (b) are referred to datasets whose magnitudes were derived from the catalogue CSTI, (c) and (d) to datasets with magnitudes derived from CSI; (a) and (c) are referred to datasets declustered with the code DECLPOI, (b) and (d) to datasets declustered with the Reasenberg technique.Full circles represent frequencies derived from historical data, open circles are referred to instrumental data.Linear trends according to theGutenberg-Richter (1944) relationship and derived using the least squares (LS) and the "maximum likelihood" method (MLM) are shown with solid and dashed lines, respectively.

Fig. 8 .
Fig. 8. Maps of the minimum (a) and maximum values (b) of PGA with a 10% exceedence probability in 50 years (PGA 0.10,50 ), obtained using the zonation proposed in this study for Northern Apulia (ZNA).

Fig. 10 .
Fig. 10.Maps of the medians of PGA 0.10,50 weighted according to the scheme shown in Fig. 7 for the seismogenic zonation ZNA (a) and ZS9 (b).

Table 1 .
Results of declustering with the new procedure described in the text (DECLPOI) and with Reasenberg technique: CAT=catalogue used as source for event magnitude; N. ev=total number of events in the initial dataset; Rem.=number of removed events; Res.=number of residual events left in the dataset; Cluster=number of clusters identified; N./clust=mean number of events per cluster.
V. Del Gaudio et al.: Seismogenic zonation and seismic hazard in Northern Apulia

Table 2 .
Estimates of a and b coefficients of the Gutenberg -Richter relationship, obtained with least square (LS) and with "maximum likelihood" approach (MLM) from integrated historical and instrumental datasets, for the four seismogenic zones of the new zonation (ZNA): M SOURCE=source of magnitude estimates for instrumental data (catalogues CSTI or CSI); DECLUST.PROCED.=procedureusedfordeclustering(DECLP=DECLPOI, REAS=Reasenberg, 1985).For each coefficient 95% confidence intervals are reported (conf.int.a, conf.int.b) and for LS the R 2 determination coefficients of regressions are shown.

Table 3 .
Seismicity rates obtained for the four zones (Z) of the zonation ZNA, expressed as mean number of events expected in 100 years for magnitude classes defined at intervals of 0.2 and identified by the central value of each class.Seismicity rates are reported for each combination of choices between instrumental dataset magnitude sources (M.S.=CSTI or CSI), declustering procedure (DECLUST.PROCED.=DECLP or REAS, for DECLPOI or Reasenberg, respectively) and Gutenberg-Richter coefficient estimate (G-R=LS or MLM for least squares and "maximum likelihood" method, respectively).