Seismic hazard of the Iberian Peninsula: evaluation with kernel functions

The seismic hazard of the Iberian Peninsula is analysed using a nonparametric methodology based on statistical kernel functions; the activity rate is derived from the catalogue data, both its spatial dependence (without a seismogenic zonation) and its magnitude dependence (without using Gutenberg–Richter’s relationship). The catalogue is that of the Instituto Geográfico Nacional, supplemented with other catalogues around the periphery; the quantification of events has been homogenised and spatially or temporally interrelated events have been suppressed to assume a Poisson process. The activity rate is determined by the kernel function, the bandwidth and the effective periods. The resulting rate is compared with that produced using Gutenberg–Richter statistics and a zoned approach. Three attenuation relationships have been employed, one for deep sources and two for shallower events, depending on whether their magnitude was above or below 5. The results are presented as seismic hazard maps for different spectral frequencies and for return periods of 475 and 2475 yr, which allows constructing uniform hazard spectra.


Introduction
The seismic design of structures requires the definition of an adequate seismic action for the site.In Spain the official information for reference in this respect is the seismic hazard map included in the Spanish seismic code NCSE-02 (Ministerio de Fomento, 2003), already about a decade old and referring only to the Spanish territory.At about the same time, Peláez-Montilla and López-Casado (2002) presented a seismic hazard analysis (SHA) for the entire Iberian Peninsula based on the methodology by Frankel (1995).Since then, several regional studies within the Iberian Peninsula have been carried out: for the south-east of Spain, García-Mayordomo et al. (2007) and Gaspar-Escribano et al. (2008); for Andalusia, Benito et al. (2010); for the north of the Iberian Peninsula, Secanell et al. (2008) and Gaspar-Escribano et al. (2011); for Portugal, Vilanova and Fonseca (2007) and Sousa and Campos (2009); and finally Mezcua et al. (2011) have presented PGA (peak ground acceleration) results for the Spanish territory within the Iberian Peninsula.
Although national norms may lose jurisdiction at political boundaries, technical findings are not sensitive to them.It is therefore decided that the present study should cover the entire Iberian Peninsula.
In recent times important advances have taken place in several areas that strongly influence a SHA.
-New evaluation methods have been proposed, some of which forego the use of seismogenetic zones (Frankel, 1995;Woo, 1996a).
-The seismic catalogue includes new information: apart from recent events, the characterisation of older ones has been enriched with additional magnitude data or with the uncertainties associated with the different parameters.
-Research on attenuation models has led to better models and to guidance on their use.
At the same time, the goals pursued by seismic design have also evolved in a number of aspects.
-A design seismic input may now be required at unpopulated locations, such as maritime areas or coastal zones.
-The characterisation of the seismic input is expected to be more sophisticated; PGA or the felt intensity are no longer sufficient and uniform hazard spectra (UHS) are required.
-Increasingly, the design of structures and facilities is based on more than one level of probability, each level associated with a given set of performance requirements.
As a consequence, it is appropriate to reevaluate the seismic hazard for the whole Iberian Peninsula, incorporating all recent advances and satisfying the new engineering needs.
The present work was developed in the context of the Ph.D. thesis by Crespo (2011).Since then, the Spanish Instituto Geográfico Nacional has carried out a re-evaluation of the national seismic hazard map incorporating the present methodology (Ministerio de Fomento, 2013).

Methodology
The methodology employed is inspired in the non-parametric density estimation (NPDE).In NPDE, the objective is to find the density function from which a given sample derives, without specifying a priori a specific shape for the density function, such as a normal distribution or a Gamma one; instead, the shape of the distribution is expected to be provided by the sample itself.
The method consists in centring a density function on each element of the sample, adding all such functions and normalising their sum.It was initially proposed by Fix and Hodges (1951) but a more recent and very clear description is provided by Silverman (1986).
The specific application of NPDE to seismic data was proposed by Vere-Jones (1992).Later on, Woo (1996a, b) presented a way of using kernel functions for modelling seismic activity rates in SHA.
The mathematical definition of density estimated with kernel functions is as follows: where n is the number of elements in the sample, H is the bandwidth, which is a measure of the separation between sample elements, K is the kernel function, and x i is the position of event i.
For generating a seismic activity rate density λ k , two changes are introduced in the previous expression.
-The normalisation with respect to the number of events n is omitted, thus the result is expressed in terms of number of events.
-Each kernel function is divided by an effective period of detection T , so the density of events is per unit time.
With the above two changes, the expression becomes The effective period T is a measure of the detection probability of that event in past times.As noted in Eq. ( 2), each event can be assigned a different effective period T which makes the methodology very versatile.Typical event characteristics on which the effective period usually depends are the event magnitude, the time of occurrence and the characteristics of the epicentral location (onshore, offshore, unpopulated area at the time of occurrence, etc.), but other factors affecting the probability of the event being detected can also be incorporated through the effective period.
The resulting activity rate density λ k depends on location as well as magnitude through the bandwidth H.As can be seen, it is a summation of kernel functions K placed on each event of the catalogue with coordinates x i .Each function is weighted with an effective detection period T ; the normalisation is achieved through the bandwidth H which depends on the distance between events, as it will be seen later in Sect.5.1.The kernel function, the effective detection periods and the bandwidth are the three main parameters influencing the activity rate density.
Several kernel functions have been proposed for the cases where the sample is of seismic type, specifically the Gaussian kernel, the inverse bi-quadratic (IBQ) and a finite one that has zero values for distances greater than one bandwidth.The first two kernels were proposed by Vere-Jones (1992); later on, Woo (1996a) proposed employing the IBQ and alternatively a finite one.All three types of kernels can be seen in Figure 1.The kernels are axially symmetric.
The IBQ kernel originally proposed by Woo (1996b), apart from being the one in common to the proposals by Vere-Jones (1992) and Woo (1996a), has been employed here.Its mathematical expression is where λ IBQ is a parameter greater than 1 for which Vere-Jones (1992) suggests values between 1.5 and 2.0.As explained before, the effective period of detection must be defined for every earthquake in the catalogue.If the catalogue were complete, this effective period of detection would be the period of time covered by the events; since in most cases, including the present one, this is not so, effective periods of detection have to be established: the purpose is to capture correctly the temporal activity rate without eliminating any event, since this would affect the spatial distribution of the activity rate density.Details on the derivation of the effective periods in this case are presented in Sect.5.2.
The bandwidth is assumed to depend on magnitude through an exponential law.The relation proposed by Woo (1996a) was where c and d are parameters to be adjusted, and M is the magnitude, or the measure employed for the quantification of the event's size.This type of dependency follows the standard form of logarithmic correlation between the magnitude and other seismic parameters such as, for instance, the fault length.Again, specific details on the derivation of the parameters c and d in Eq. ( 4) are presented in Sect.5.2.
With the activity rate density already estimated and the attenuation model, the rate of occurrence of different levels of ground motion y j can be derived as follows: (5) The rate of occurrence is computed with Eq. ( 5), and follows exactly the summation specified in Eq. ( 2) for the construction of the seismic activity rate density; the procedure is quite simple though very intensive in calculations since for each magnitude it goes through the entire catalogue.The kernel function is given in Eq. ( 3) and the shape of the bandwidth in Eq. ( 4).These four equations give all the mathematical input needed to understand the construction of the rate of occurrence.
The seismic activity rate density has to be calculated for the entire area that is susceptible of generating earthquakes that are relevant for the location studied.The extension of this area depends on the minimum threshold of ground motion of interest and the maximum magnitude in the area; with these data and the attenuation relation, it is possible to estimate the extension of the area where the seismic activity rate density is needed.
Once the rate of occurrence of different levels of ground motion has been obtained, the seismic hazard curve at the site is derived assuming a Poisson process.
This methodology was originally implemented by Woo (1996b) in the Fortran program KERFRACT.In the context of the work by Crespo (2011), the code was modified in several aspects that are explained in detail in her work.However, the most important changes affecting the calculation algorithm are as follows.
-For each magnitude being integrated, the decision of whether an event is considered or discarded when computing the seismic activity is made dependent on the uncertainty associated to the event magnitude.This change is especially relevant for historical events that have large magnitude uncertainties.
-The computation of the activity rate may have the depth as an additional variable of integration.
This methodology has been already applied in the past in many parts of the world.Apart from the authors (Crespo and Martí, 2002;Crespo et al., 2003Crespo et al., , 2007)), Secanell et al. (2008) also employed it for studying the Pyrenean region, Menon et al. (2010) for South India, and Goda et al. (2013) for the UK.
A question may be posed as to the reliability of this description of the activity rate and its relative merits in comparison with traditional zoned procedures (Stein et al., 2011).Several considerations may be offered in this respect.A first advantage over zoned procedures comes from avoiding the imposition of constant seismic activity over zones with stepwise jumps across boundaries; this is particularly significant in areas of low to medium activity, with no clear association between seismicity and geological features, where the delineation of zones is a primarily subjective matter which has a strong influence on the results.In this respect Giner et al. (2002) presented five different zonifications, constructed by different authors for the same area of Spain, that produce very different hazard results.
Even for highly active regions like California, Jackson and Kagan (1999) presented their exercise in earthquake forecasting using a methodology based on a smoothing approach based on kernel functions, similar to the methodology proposed by Woo (1996a).These same authors participated in the RELM (Regional Earthquake Likelihood Models) test (Sachs et al., 2012) obtaining with their methodology the best mean forecast (Helmstetter et al., 2007).This further confirms the better performance of the zoneless methodology over traditional zoned approaches.
Further, as shown by the results presented here, there are areas like south-eastern Spain and the Pyrenean region where the zoneless approaches arrive at significantly higher acceleration values than the traditional methodologies.The recent activity felt in south-eastern Spain has been considerably higher than expected from available hazard maps, which again seems to confirm the limitations of the zoned approach.A similar situation has arisen in Italy regarding the recent activity around l'Aquila and Emilia Romagna (G.Woo, personal communication, 2013).

Sources
The main catalogue used in the study is the seismic database maintained by the Spanish Instituto Geográfico Nacional (IGN, 2010): more than 95 % of the final catalogue employed here is constituted by events from the IGN database.However, some areas that are relevant for the seismic hazard of the Iberian Peninsula are not fully covered by the IGN catalogue.This is especially true for the south-east of France; also some large events in north-western Italy may influence the seismic hazard of the north-east of Spain.The IGN catalogue has been completed with data from two international organisations (USGS, 2011;ISC, 2010), data from the catalogues of neighbouring countries (BRGM, 2010;Gruppo di lavoro CPTI, 2004) and published works that describe the seismicity of surrounding regions (Vilanova and Fonseca, 2007;Peláez et al., 2007).
Nearly 60 000 earthquakes have been considered, which include all magnitudes/intensities registered and dependent events.Most of them, more than 97 %, come from the IGN catalogue; some 1500 were added from the ISC, most of them located in France and in the north of Africa and belonging to the instrumental period; finally, just 45 earthquakes were taken from the USGS, all of them catalogued as "significant events" and located outside the Spanish territory in areas not covered by the IGN catalogue.
The most ancient event with an assigned intensity dates back to 1048 and belongs to the IGN catalogue, while the most recent events considered date from November 2010.

Magnitude unification
The quantification of events has to be homogeneous throughout the catalogue.The homogenisation has been conducted in terms of magnitude, specifically the moment magnitude M w , which is becoming the standard measure for seismic events and is the one used in the majority of modern attenuation relations.Some 10 % of the events in the IGN catalogue are only quantified with epicentral intensity, as is frequently the case for events from the historical period.For these historical earthquakes, if they had an assigned magnitude either by the IGN (2010) or by specific studies (Martínez-Solares and López-Arroyo, 2004;Vilanova and Fonseca, 2007) that value has been considered; in the other cases the correlation by L. Cabañas (personal communication, 2010) specifically developed for the Iberian Peninsula has been employed: where M w is the moment magnitude, I 0 is the epicentral intensity, and σ is the standard deviation.
During the instrumental period most events have an assigned magnitude; generally this magnitude is either m b or m bLg , and since 2002, some of them have also a moment magnitude.The type of magnitude is specified in the catalogue provided by the IGN (2010).The conversion of both epicentral intensity and magnitude m b or m bLg into moment magnitude has been carried out using correlations specifically developed for the Iberian Peninsula by the IGN (L.Cabañas, personal communication, 2010).Specifically, three subsets of earthquakes are considered, each one having a different correlation.
-For earthquakes characterised with magnitude m b : -For earthquakes before March 2002 and characterised with magnitude m bLg : -For earthquakes after March 2002 and characterised with magnitude m bLg : Moment magnitudes already assigned by the IGN ( 2010) have been maintained, and some other moment magnitudes from the literature have also been incorporated (Stich et al., 2003(Stich et al., , 2010)).As it will be explained in Sect.6, the magnitude integration range starts at M w = 3.5.Magnitude uncertainties are considered on an event-by-event basis assuming a Gaussian distribution with the mean being the magnitude in the catalogue.Hence it is necessary to consider lower magnitude events in the process of integration.In this case the minimum magnitude threshold in the catalogue is M w = 3.0.

Dependent events
It is assumed that earthquake events are Poisson distributed, which requires removing the dependent events.The methodology followed for this purpose is the traditional one of placing a time-space window around the main events to identify those that need to be discarded.The methodology followed is inspired in the one proposed by Gardner and Knopoff (1974), but modified in two aspects.
-The window size depends on the event magnitude according to the indications given by Peláez et al. (2007); namely, for M w = 3.0 the spatio-temporal window was of 20 km and 10 days, and 100 km and 900 days for M w = 8.0.
-The catalogue is scanned in decreasing magnitude order, so that in the potential series the first event identified is always the main shock (Crespo, 2011).
As a result of the pruning process conducted, 36 % of the events were considered dependent events (either foreshocks or aftershocks) and were discarded from the catalogue.

Uncertainties
The methodology used requires the catalogue to incorporate uncertainties with respect to magnitude as well as to location (epicentre and depth).
Regarding magnitude, the uncertainties quoted in the various catalogues have been taken into account.When using a correlation, its uncertainties have been combined with those in the catalogue.It should be pointed out that the correlations provided by the IGN (L.Cabañas, personal communication, 2010) have been derived with an RMA (reduced major axis) methodology which takes into account the variability of the two variables being adjusted.These correlations have been published by the IGN (Ministerio de Fomento, 2013).
As for the location uncertainty, the catalogues provide this type of information for many earthquakes, especially in the instrumental period.For events lacking this information in the catalogue, an uncertainty has been estimated.It seems logical that this uncertainty be higher for events that occurred at older times, and that it be higher for epicentres at sea.The uncertainties assigned to the epicentral location follow the recommendations by G. Woo (personal communication, 2000) and are also in agreement with the suggestions by other authors (Giner, 1996;Molina, 1998;García-Mayordomo, 2005).The uncertainties finally assigned to sea and land epicentres appear in Table 1.

Attenuation model
Major research on the attenuation of ground motion has been taking place for at least half a century, but important advances have been produced recently, both in relation with the number of models available and with their functional forms.
Very relevant work has also been conducted on the criteria for selecting attenuation models when conducting a SHA.Cotton et al. (2006) proposed a series of seven criteria that in practice were rather difficult to satisfy.Later on, Bommer et al. (2010) reformulated the criteria into a set of ten recommendations.
Attenuation studies at present tend to focus on medium to high magnitudes, a range that is representative of only a small fraction of the seismicity in low to medium seismicity areas such as the Iberian Peninsula.Bommer et al. (2007) alerted to the problems that can arise when an attenuation model is employed out of the range for which it was developed.They also presented some preliminary results obtained with an attenuation model that employs a database spanning magnitudes 3.0-7.5;their results in the lower range are consistent with those of Bragato and Slejko (2005).
In the present work, three different attenuation models have been employed.
- Bragato and Slejko (2005), for shallow earthquakes with moment magnitudes below or equal to 5.0.
The boundary between shallow and deep earthquakes has been fixed at 35 km, all deep earthquakes occurring between 35 and 200 km, except for a few around Granada which are around a 600 km depth.These very deep earthquakes are not included in the calculation.In the same way, deep earthquakes are accounted for just for magnitudes above M w = 5.0, which is the range of applicability of the attenuation relation and a sensible lower bound for deep earthquakes; hence, deep earthquakes around the Pyrenean region and many around the south of the peninsula also fall outside the group of interest.
Most modern attenuation models are well suited for dealing with shallow earthquakes with magnitude above 5.0.The reasons for choosing that by Ambraseys et al. (2005) are that it was developed from a wide database of European records, including some from the Iberian Peninsula; that it has been used in a number of SHA, including recent ones in the Iberian Peninsula (Garcia-Mayordomo, 2005;Mezcua et al., 2011); and that it satisfies 9 of the 10 points proposed by Bommer et al. (2010).The model by Akkar and Bommer (2010), based on the same database as the one by Ambraseys et al. (2005), was also initially considered but eventually discarded because its sampling of spectral periods is equally spaced throughout the range covered instead of gradually increasing with the period as the amplification becomes smaller.In an area of medium-low seismicity like the Iberian Peninsula, the maximum spectral amplification in terms of acceleration is expected to appear at low periods.Since one of the objectives of the work was to capture this spectral amplification, the sampling provided by Ambraseys et al. (2005), smaller than other attenuation models and consistent with that of Bragato and Slejko (2005), was considered more appropriate.Additionally it is consistent with the model by Bragato and Slejko (2005), with which it will be combined.In relation with the magnitude, It is not important that the model by Ambraseys et al. (2005) does not include the non-linear magnitude dependence proposed by Bommer et al. (2010), because it will only be used in part of the magnitude range.
The number of attenuation models generated with databases that cover the low magnitude range, specifically below 5.0, is considerably smaller.That by Bragato and Slejko (2005) is one of them and, as already mentioned, is consistent with the preliminary results presented by Bommer et al. (2007).Additionally, it satisfies 8 out of the 10 criteria proposed by Bommer et al. (2010), the remaining two being of lesser importance since the model will be employed in only part of the magnitude range.The minimum magnitude at which the integration of the seismic activity rate density begins is M w = 3.5.This value was finally adopted after repeated calculations showed that further lowering of this threshold only affects the results in areas with lower hazard and for return periods below 475 yr.
There are studies based on macroseismic intensity data that conclude different attenuations for different regions within the Peninsula (Sousa and Oliveira, 1996;López Casado et al., 2000).Here, only one attenuation model for shallow seismicity is employed for the whole Iberian region; this decision follows the recommendations by Bommer et al. (2010), who claim that there is no strong evidence for persistent regional differences in ground motions among tectonically comparable areas.Their recommendations are based on the works by Douglas (2007) and Stafford et al. (2008) which rely on recorded strong ground motions rather than macroseismic data.
Attenuation models for deep earthquakes are very scarce and those available refer to subduction zones like Chile, west coast of Mexico, Peru or Japan.To account for deep seismicity in the calculations, a model must be used that includes focal depth as an independent variable; there is then little choice other than a subduction model, even if the tectonic patterns represented are not the ideal ones.Among the available subduction models, that by Youngs et al. (1997) was selected because it scores higher in the list by Bommer et al. (2010).
The attenuation models have to be harmonised for consistency in their input parameters and in the results produced.In relation with magnitude, the catalogue had been homogenised in terms of M w , which is already the measure employed by Ambraseys et al. (2005) and Youngs et al. (1997).The model by Bragato and Slejko (2005) was originally developed in local Richter magnitude M L but, for the range in which this model is applied here, both magnitudes can be considered equivalent, as also assumed by Bommer et al. (2007).
The distance considered is the epicentral distance one.The model by Bragato and Slejko (2005) is based on that distance.That by Ambraseys et al. (2005) employs the Joyner-Boore (JB) distance but, for moderate seismicity, the JB distance can be assumed to be equal to the epicentral one.Youngs et al. (1997) make use of a similar equivalence between the epicentral distance, adequately combined with depth, and the distance to the rupture surface originally employed in the model.
The type of soil for which the results are derived is "stiff soil" with a shear wave velocity between 360 and 750 m s −1 : this is the range for which the Bragato and Slejko (2005) model was investigated in detail and is also one of the soil types considered by Ambraseys et al. (2005).It corresponds also to one of the two types of ground studied by Youngs et al. (1997).
Finally, the acceleration measure used in the regression, which is also the output acceleration, may vary; indeed all three models use different measures of acceleration.The present results employ the same measure as Ambraseys et al. (2005), which uses the maximum horizontal component.For Bragato and Slejko (2005) a conversion has been carried out using their own model guidelines, and in the case of Youngs et al. (1997), the correlations by Beyer and Bommer (2006) have been incorporated.
Figures 2 and 3 present part of the attenuation model constructed for shallow seismicity, which includes the correlations by Ambraseys et al. (2005) and Bragato and Slejko (2005), with the aforementioned modifications.As can be seen, the combination of these two models, with the appropriate homogenisation, yields good continuity at the M w = 5.0 boundary.As a qualitative verification, the maximum Spectral Acceleration (SA) can be seen to move towards the lower frequencies with increasing magnitude, as could be expected.

Seismic activity rate
The seismic activity rate is calculated according to Eq. ( 2).The kernel function employed is the IBQ, as stated in Sect. 2.
In this section particular choices made for the determination of the bandwidth and the effective detection periods are described first, which are the two fundamental elements on which the seismic activity rate relies on.
The seismic activity rate presented in this section represents an intermediate result of the overall calculation in contrast with the seismic activity represented by the Gutenberg-Richter (1944) relationship, which is an input in the traditional zoned approach.There is no need to explicitly obtain these results for deriving the final hazard output, in fact the original implementation by Woo (1996b) did not offer the possibility of retrieving them; however it has been considered a useful exercise in order to help understanding the differences between the traditional and the kernel approach.

Bandwidth
The bandwidth H(M w ) depends on magnitude as indicated in Eq. ( 4).This type of logarithmic relation between moment magnitude and bandwidth is commonly assumed to be applicable between moment magnitude and other seismic parameters.
The derivation of c and d is performed via a least-square fit.
-Events are classified in groups according to their magnitude.
-For each event, the distance to the nearest epicentre within the same magnitude range is determined.
-All minimum distances calculated for each magnitude range are averaged.
-A least-square fit is conducted in order to obtain the two parameters, c and d, that appear in Eq. ( 4).
The derivation of the parameters c and d is done independently for shallow and deep earthquakes.If Eq. ( 4) is adjusted using the complete catalogue, the resulting correlation is that presented in Figure 4.However, since the fit is conducted at all points where the hazard is computed, the spatial distribution of c can also be produced, as shown in Figure 5.As reflected in the figure , c shows stable values between 0.25 and 0.50 for the areas of higher activity.In the central part of the Peninsula and in the Balearic Islands, the values increase to around 3.0, which is consistent with the fact that a lower seismicity implies earthquakes more distant from each other.

Effective detection periods
The effective detection periods, parameter T (x i ) in Eq. ( 2), have been calculated considering whether an earthquake is shallow or deep, and for the shallow ones, whether the epicentre is on land or at sea.
The magnitude and year of occurrence of each earthquake are taken into account for determining the effective period of the earthquake as follows.
-History is divided into time intervals as a function of the means of detection available at each time.
-For each interval with duration D i and each level of magnitude a probability of detection p im is estimated.
-A reference year A m is established for each magnitude level: where A 0 is the most recent year with records.Probabilities of earthquake detection over time need to be assigned.This can be done by direct estimation based on the possibilities offered by the available technology.However, in this case, the estimation of the probability has been made by comparison of activity rates with the one observed during the completeness period: the probability is considered to be proportional to the observed activity rate in the catalogue, being 1.0 in the completeness period and lower in the previous times.
Table 2 shows the reference years obtained by the previous procedure for the calculation of the detection periods.

Computation of the seismic activity rate
Having defined the kernel function K(x), the bandwidth H(M w ), and the effective detection periods T (x i ), the activity rate can be calculated using Eq. ( 2).The resulting activity depends on location and magnitude.
For a given magnitude level, the seismic activity rate can be contoured; similarly, for a given location, the activity rate can be plotted.This second result is equivalent, except for a factor of area, to the traditional activity rate derived with the Gutenberg and Richter (1944) relationship.Figure 6 shows, for the south-west of the Iberian Peninsula, a contour map of the activity rate for earthquakes with magnitude above 3.5, which is consistent with the seismicity of this area.This plot constitutes the first type of intermediate results that can be produced.ALthough this type of results have been obtained for the whole Peninsula, we include here the ones for the South-West, which is one of the most illustrative areas.
Figure 7 shows plots of the activity rate for the four locations marked in the previous map.These curves represent the number of earthquakes per unit area and per year for a given location, similar to the concept used in the Gutenberg-Richter relationship.However, here this is just a byproduct of the calculations, while in the traditional procedure it is a necessary step on which some assumptions, like the analytical shape of the correlation, will be based.In the present case, the shape, not necessarily linear, arises directly from the catalogue data.
The graphs in Figure 7 include also a linear fit of the activity rate derived from the calculations.The slope of this line, referred to here as b k , may be compared with the b parameter of the Gutenberg-Richter relationship.For each point of the map, similar fits could be performed leading to a value of the b k parameter; the results are plotted in Figure 8.As can be seen, the values found for b k are consistent with the typical ones already predicted by Gutenberg and Richter (1944) for the b parameter.

General considerations
Once the activity rate λ k has been calculated (Sect.5) and an attenuation model has been adopted (Sect.4) the seismic hazard can be calculated.
The magnitude integration range for shallow seismicity starts at M w = 3.5 and for deep earthquakes at M w = 5.0.The lower limit for integration depends on the return period of interest and on the minimum acceleration levels that need to be correctly captured.In the present case it was verified that, for the 475 yr return period, the integration should start at M w = 3.5 if accelerations as low as 0.04 g must be determined: using a lower limit for integration does not change the results significantly, while higher values would affect the results.This observation confirms the stability and consistency of the attenuation model constructed.The question of whether or not the minimum accelerations derived are of engineering interest should be decided by the engineer who will be making use of them and will depend on the type and purpose of the calculation being performed.What has to be guaranteed at this stage is that the minimum accelerations being derived are reliable, which requires an adequate selection of the magnitude integration threshold.
The upper limit of integration arises from the maximum magnitude recorded in the catalogue and its uncertainty, specifically, M wi + 2σ Mwi , σ Mwi being the uncertainty associated with event i.Again, it was verified that the results are not sensitive to increases in this upper limit, since each event contribution is weighted with the probability that its magnitude falls within the integration range, assuming a lognormal magnitude distribution.
For deep earthquakes, the range of integration used is 35 to 200 km.

Contour maps
The PGA has been calculated for a rectangular area that spans the Iberian Peninsula for return periods of 475 and 2475 yr.The resulting contours appear in Figs. 9 and 10.
The reason for choosing the 2475 yr period is that it is becoming an important reference in seismic design (ASCE, 2010;NRC-IRC, 2010).
It is also around this return period that the contribution to the seismic activity rate arising from geological considerations may start being significant.Given the convergence rate between the Iberian and African plates it is expected that faults in the Iberian Peninsula are also slow, so morphogenic earthquakes (hence with magnitude above 6.0-6.5)have recurrence periods of the order of few thousand years, which exceeds the period covered by the seismic catalogue: note that the first event in the catalogue dates from the fourth century BC, and the first quantified event from the fifth century AD.The explicit inclusion of geological data would be carried out adding earthquakes that are representative of the characteristic magnitude of the fault and with an effective period consistent with the fault recurrence period.This inclusion affects two aspects: the location where the activity rate is modelled, which concentrates around specific geometrical features or faults; and the activity rate itself, which is enriched with geological information that the seismic catalogue does not reflect.
The effect on the distribution of the activity rate is partly incorporated using the kernel methodology: the activity rate is not uniformly spread over zones, but this is only guided by the epicentral locations.All this should be taken into account when considering the results.
Looking at the 475 yr map, the areas with higher hazard levels are Granada (south of Spain) and the French slope of the Pyrenees, where the PGA values reach 0.30 g, followed by the south-east of the Peninsula, around the city of Alicante, with a value of 0.27 g.Outside the Iberian Peninsula, the area around Algiers has 0.25 g.Finally, Lisbon and the north of Catalonia attain values of 0.20 g and 0.15 g, respectively.
The results for a 2475 yr return period are presented in Figure 10.It can be seen that the distribution of the hazard is similar to that observed for 475 yr, with approximately double the acceleration values.This ratio is consistent with the dependence on the return period proposed in the Spanish seismic code (Ministerio de Fomento, 2003).
Two additional hazard maps are shown in Figs.11 and 12: they correspond to the same return periods mentioned earlier but refer to a spectral period of T = 0.1 s.As clarified in the next section, this spectral ordinate is the more representative one for deriving the maximum spectral accelerations, which are those of the plateau.This spectral ordinate peaks in the same places as the PGA: around Granada it reaches 0.90 g for 475 yr and 1.90 g for 2475 yr.
The ratio of the plateau acceleration to the PGA for 475 yr turns out to be between 2.5 and 2.8 in most of the Iberian Peninsula (Figure 13), which is in agreement with the value of 2.5 frequently found in seismic codes.

Uniform hazard spectra
Results were also produced in terms of uniform hazard spectra for the eight locations marked with red crosses in Figure 9.They are shown in Figure 14 for selected sites from the south of the Iberian Peninsula and in Figure 15 for those in the north.
The maximum amplification appears between periods T = 0.1 s and T = 0.15 s, which is consistent with spectra corresponding to low magnitudes derived directly from the attenuation model (Figure 3).
In the figures, the dashed lines represent the results of the seismic hazard analysis, while the solid lines correspond to an analytical approximation based on constant acceleration, velocity and displacement branches.The more adequate periods for anchoring those branches were found to be T = 0.1, 0.4 and 2.0 s, respectively.However, the latter corresponds, in fact, to the highest period for which results are available; hence the possibility exists that a longer period might have been preferable.

Comparison with previous studies and seismic codes
It is interesting to compare the results shown above with those of other recent regional studies.Benito et al. (2010) presented a study for Andalusia.The geometry of their contour lines is consistent with that shown in Figure 9 and their numerical values are also consistent, possibly slightly lower, after correcting for the fact that their results correspond to rock conditions.For Granada the present study calculates 0.30 g (for stiff soil) while they find 0.22 g for rock; this difference may be partly explained from the consideration of different soil types.Secanell et al. (2008) studied the area around the Pyrenees.They used the kernel methodology proposed by Woo (1996a) in one of the branches of a logic tree.As in the previous case, after correcting for the soil type, there is good consistency with the esults presented here, both qualitative and quantitative.
For the area around Alicante, we calculate 0.28 g, a value that is significantly higher than the one reported using zoned methods (García-Mayordomo, 2005;García-Mayordomo et al., 2007;Benito et al., 2006), but that is in agreement with Peláez-Montilla and López-Casado (2002), who employed a variant of the zoneless approach by Frankel (1995) working in terms of felt intensity.
For these two latter areas (Pyrenees and Alicante) good consistency is found with previous studies based on zone-less approaches but not otherwise.This suggests that the spatial gradient of the activity rate in these areas is high, so the zonation needed in order to represent this seismic activity rate would require small zone dimensions, which might be difficult to characterise, in terms of magnitude-frequency relationships, with the seismic information available.
For the area of Portugal the present study reaches 0.20 g around Lisbon and 0.13 g in the southern coast of Portugal (Algarve region).The work by Vilanova et al. (2007) finds for rock 0.19 g for Lisbon and 0.20 g for the Algarve region; after accounting for the different soil type, the value around Lisbon would be slightly lower in our study, with a somewhat larger difference for the south of Portugal.The attenuation relationships employed by Vilanova et al. (2007) are probably responsible for this discrepancy.Vilanova et al. (2007) used three different attenuation models in the branches of their logic tree and the lowest values were derived using Ambraseys et al. (1996), not very different from Ambraseys et al. (2005) employed here.Considerably larger are the values found by Peláez-Montilla and López-Casado (2002), but this is likely a consequence of the fact that they work in terms of intensity and translate the final results into accelerations.
Table 3 compares the calculated PGA values for eight cities in the Iberian Peninsula and compares them with those reported elsewhere.The locations are indicated with red crosses in Figure 9.
Regarding the return period of 2475 yr, although it is becoming a more frequent reference (OPCM, 2008;NRC-IRC, 2010;ASCE, 2010), no results have been published recently for it so comparisons cannot be offered.
Regarding the UHS results, the range of spectral periods with maximum amplification appears to be narrower than indicated by the Spanish seismic code (Ministerio de Fomento, 2003) for soil type II, which is equivalent to the stiff soil considered here; however it is in agreement with the periods used to limit the plateau in Eurocode 8 (CEN, 2004) for earthquakes with magnitude below M s = 5.5.
It should also be noted that a number of seismic codes construct the design spectra by using several spectral ordinates.ASCE (2010) takes T = 0.2 s as the reference period for the plateau and T = 1.0 s that for the constant velocity branch.These values are probably adequate for regions more active than the Iberian Peninsula (see again Figure 3): Figures 14  and 15 suggest that the period T = 0.2 s should more likely belong to the constant acceleration branch.
Using several spectral ordinates to construct design spectra is probably the best way of taking advantage of all the information produced with the new attenuation models.The Italian seismic code (OPCM, 2008) is a good example in Europe.

Summary and conclusions
A seismic hazard assessment of the Iberian Peninsula has been carried out using the current state-of-the-art in seismic engineering and taking into account the present engineering needs.
From the point of view of the methodology, the seismic activity rate is constructed via a non-parametric estimation based on kernel functions.The resulting rate presents a continuous variation in space and magnitude and its shape is not assumed but is derived in the process.No seismogenic zones are needed and the uncertainties in magnitude and position are incorporated in the methodology.
The main seismic catalogue employed in the calculation is the IGN database.However, the construction of the activity rate in certain areas requires information from a wider area than the one covered by the IGN, particularly the south of France.As a consequence, the IGN catalogue has been supplemented with information from other catalogues or published studies.
For the homogenisation of the catalogue, made in terms of moment magnitude M w , it has been of great help the existence of a growing number of events in the IGN catalogue that have more than one type of magnitude assigned, information that permits to establish suitable correlations between different types of magnitude.Once homogenised, dependent events have been identified and discarded.
The seismic activity rate compiled has a continuous variation with respect to location and magnitude: its characteristics are derived solely from the information contained in the catalogue.
The attenuation model has been constructed combining relations from three different authors.It has proven to be stable with respect to the variations of the minimum and maximum magnitudes between which the integration of the activity rate density is performed.
Results have been obtained for the Iberian Peninsula for four spectral frequencies, two return periods, and a soil with a shear wave velocity between 360 and 750 m s −1 .
The results for PGA and 475 yr are generally consistent with different recent studies performed at a regional scale.There are two regions, specifically the south-east of Spain and the Pyrenean region, where consistency is achieved only with studies that followed a zoneless approach, being the accelerations obtained here significantly higher than the ones obtained with traditional procedures; this fact suggests that the spatial gradient of the activity rate in these areas is high so the zonation needed to represent this seismic activity rate would require small zone dimensions, which might be difficult to characterise with the available seismic information.
The results for 2475 yr are consistent with those for 475 yr, with the acceleration ratios between the two return periods being in the usual range.
Uniform hazard spectra have been derived for eight locations.The spectral shapes present the maximum amplifica-tion of acceleration around 0.1 s, which should be expected for a medium-low seismicity area like the Iberian Peninsula.The best points for anchoring branches of constant acceleration, velocity and displacement branches appear to be 0.1, 0.4 and 2.0 s, respectively; they are lower than usually found in seismic codes, but are consistent with the low seismicity in the area.

Fig. 1 .
Fig.1.Probability densities of various kernel functions: the Gaussian kernel (GAUSS), the inverse bi-quadratic kernel (IBQ), and a finite kernel that vanishes in one bandwidth (FIN).

Fig. 4 .
Fig. 4. Derivation of the bandwidth parameters using the catalogue data, separately grouped into shallow and deep events.

Fig. 5 .Fig. 6 .
Fig. 5. Distribution of the bandwidth parameter c; the higher values (red and purple) that appear in low seismicity areas (Central Peninsula and Baleares Islands) indicate a greater distance between events of similar size.

Fig. 7 .Fig. 8 .
Fig. 7. Activity for each of the four locations identified in Figure 6: (a) point A; (b) point B; (c) point C; (d) point D. The relation is approximately straight.

Fig. 9 .
Fig. 9. Distribution of the accelerations obtained for T = 0 s and 475 yr return period.

Fig. 10 .
Fig. 10.Distribution of the accelerations obtained for T = 0 s and 2475 yr return period.

Fig. 11 .
Fig. 11.Distribution of the accelerations obtained for T = 0.1 s and 475 yr return period.

Fig. 12 .
Fig. 12. Distribution of the accelerations obtained for T = 0.1 s and 2475 yr return period.

Fig. 13 .
Fig. 13.Maximum amplification of the accelerations for 475 yr return period, measured by the ratio of the plateau to the PGA.

Fig. 14 .
Fig. 14.Comparison of uniform hazard spectra and fitted code spectra for cities in the south of the Iberian Peninsula: (a) Huelva; (b) Granada; (c) Alicante; (d) Lisbon.

Fig. 15 .
Fig. 15.Comparison of uniform hazard spectra and fitted code spectra for cities in the north of the Iberian Peninsula: (a) Orense; (b) Pamplona; (c) Gerona; (d) Oporto.

Table 2 .
Reference years for different magnitudes and types of events.

Table 3 .
Results of PGA (g) for selected locations.Results of the present study are compared with values given by the official Spanish seismic code as well as other published studies.