Abstract
We present a time-independent gridded earthquake rate forecast for the European region including Turkey. The spatial component of our model is based on kernel density estimation techniques, which we applied to both past earthquake locations and fault moment release on mapped crustal faults and subduction zone interfaces with assigned slip rates. Our forecast relies on the assumption that the locations of past seismicity is a good guide to future seismicity, and that future large-magnitude events occur more likely in the vicinity of known faults. We show that the optimal weighted sum of the corresponding two spatial densities depends on the magnitude range considered. The kernel bandwidths and density weighting function are optimized using retrospective likelihood-based forecast experiments. We computed earthquake activity rates (a- and b-value) of the truncated Gutenberg–Richter distribution separately for crustal and subduction seismicity based on a maximum likelihood approach that considers the spatial and temporal completeness history of the catalogue. The final annual rate of our forecast is purely driven by the maximum likelihood fit of activity rates to the catalogue data, whereas its spatial component incorporates contributions from both earthquake and fault moment-rate densities. Our model constitutes one branch of the earthquake source model logic tree of the 2013 European seismic hazard model released by the EU-FP7 project ‘Seismic HAzard haRmonization in Europe’ (SHARE) and contributes to the assessment of epistemic uncertainties in earthquake activity rates. We performed retrospective and pseudo-prospective likelihood consistency tests to underline the reliability of our model and SHARE's area source model (ASM) using the testing algorithms applied in the collaboratory for the study of earthquake predictability (CSEP). We comparatively tested our model's forecasting skill against the ASM and find a statistically significant better performance for testing periods of 10–20 yr. The testing results suggest that our model is a viable candidate model to serve for long-term forecasting on timescales of years to decades for the European region.
1 INTRODUCTION
A basic component of any probabilistic seismic hazard assessment (PSHA) is an earthquake source model that defines earthquake activity rates, that is the occurrence rates of events as a function of space, time and magnitude. Earthquake source models commonly use catalogued earthquakes as input, combined with knowledge of faults, geology, tectonics and/or other indicators of the seismogenic potential of a region. Such models are essentially earthquake rate forecasts that can be compared with seismicity following the period of data these models have been developed from. For a PSHA, they are combined with ground motion prediction equations, resulting in exceedance probabilities for a specified ground motion intensity measure in a given time period.
There is no unique solution to building a rate forecast model, and multiple approaches to build an earthquake source model for a PSHA exist, with the most often used approaches being (1) the regionalization of the study region following the traditional area source model (ASM) concept (including the consideration of fault sources; Cornell 1968; Wiemer et al.2009a) and (2) the description of earthquake activity for geologically assessed faults in combination with a background seismicity described by some smoothing algorithm (Petersen et al.2008; Field et al.2009; Field & Page 2011; Stirling et al.2012). In addition to these approaches that involve subjective expert judgments, efforts have been made to introduce more objective, data and algorithm driven models. In particular, kernel-smoothed seismicity approaches, starting on a prominent yet experimental attempt to map the seismic hazard within the central and eastern United States (Frankel 1995), have been suggested as an alternative model branch to address epistemic uncertainties of earthquake occurrence within a PSHA.
Following the pioneer models of Kagan & Jackson (1994) and Frankel (1995), several other kernel-smoothed seismicity approaches have been proposed (e.g. Cao et al.1996; Lapajne et al.2003; Woo 1996; Jackson & Kagan 1999; Stock & Smith 2002; Helmstetter et al.2007; Zechar et al.2010a; Werner et al.2011) which all rely on using only seismicity as input to estimate future earthquake locations and rates, with differences in the shape and functional form of the smoothing kernel. Helmstetter & Werner (2012) used space–time kernels to obtain spatial densities of earthquakes, thereby circumventing the relatively subjective choice of a declustering algorithm. Some efforts were also made to include tectonic knowledge in terms of regionalization schemes (e.g. Burkhard & Grünthal 2009). Pseudo-prospective forecast experiments have been used to optimize the kernel width for a given time period by separating the earthquake catalogue into learning and target periods, finding the optimal kernel width based on likelihood evaluation procedures (e.g. Zechar et al.2010a). Such data-driven smoothed seismicity models perform well in prospective testing of earthquake forecasts within the framework of the Collaboratory Study for Earthquake Predictability (CSEP, Jordan 2006) for the intermediate-term forecast period of 5 yr (e.g. the Regional Earthquake Likelihood Models experiment, RELM, Field 2007; Schorlemmer et al.2010; Zechar et al.2013). This is a result of one of the basic assumption of these models: earthquakes occur at or very close to locations of previous seismicity.
Despite the relative success to forecast seismicity rates, skepticism still exists within the hazard community towards the use of this method in PSHA due to its basic assumptions: (1) future seismicity occurs in areas close to where past seismicity occurred, (2) seismic catalogues limited to the past few decades to centuries are sufficient to model future seismicity on time spans of centuries to millennia, (3) tectonic structures are not considered, assuming that active features are mapped via the catalogued seismicity alone and (4) the occurrence of the rare largest events can be forecast using the numerous smaller ones, hence earthquake size scaling holds up to the maximum event size considered.
The first assumption has proven to be very useful along well-defined plate boundaries. For intraplate seismicity evidence has been presented supporting the idea that events delineate wider zones with a higher likelihood for events to occur in the proximity of former events (Kafka 2007); however, the assumption has been challenged not to be a very good proxy for future seismicity in intraplate regions as seismicity may jump between fault structures that have been inactive for long time spans (Swafford & Stein 2007; Stein et al.2009; Liu et al.2011). The additional assumptions arise from the successful application of the models in CSEP-testing experiments and from evidence for no breaks in scaling over a large magnitude range (e.g. Main 2000; Kwiatek et al.2010). In fact, the seismicity record in many regions is not long enough to support these assumptions albeit the efforts of prolonging this record with historical, archeological and palaeoseismological studies (Meghraoui et al.2001; Fäh et al.2006; Camelbeeck et al.2007; Hinzen & Reamer 2007) which are used when available.
Hiemer et al. (2013) introduced a stochastic earthquake source model that addresses assumption (3) of earlier gridded smoothed seismicity models. By combining information on active faulting with a smoothed seismicity approach, they build an alternative source model for California. The model applies in essence the kernel density estimation technique to both, past seismicity and fault moment release, with the latter being estimated from slip rates on mapped active fault structures. We use the term SEIFA for such a combined smoothed SEIsmicity and FAults model. The resulting forecast relies on data-driven likelihood optimization techniques and is thus less dependent on subjective expert judgments compared to other source model types used in PSHA—though some cannot be avoided. A similar model has been presented for New Zealand (Rhoades & Stirling 2012) pointing out shortcomings within the New Zealand national seismic hazard assessment (Stirling et al.2012).
In this paper, we adapt and improve the SEIFA-approach presenting its application to Europe, a much larger region, more diverse in terms of tectonic regimes and data availability. We use the European Database of Seismogenic Faults (EDSF; Basili et al.2013), which includes both crustal faults and the subduction zones of the Calabrian, Hellenic and Cyprus arcs. The implementation at the European scale is possible thanks to the large-scale community effort made in the European Union project SHARE (Seismic HAzard haRmonization in Europe) to homogeneously parameterize the fault sources and to provide the necessary input data for PSHA. The EDSF comprises not only structures in active tectonic regions similar to California (e.g. the North Anatolian fault) or moderately active tectonic regions (e.g. the Apennines), but also in regions of low activity such as the Lower Rhine embayment (Camelbeeck et al.2007; Hinzen & Reamer 2007) or regions like the Po Plain typically characterized by blind faulting (Burrato et al.2012). In addition, the SEIFA-approach can treat crustal faults and subduction sources (i.e. subduction interfaces) in the same way as both are characterized by geometry and slip rate data that can be used to estimate geological moment release.
In the following sections, we outline the methodological concept of the stochastic earthquake source model and detail choices for the construction of the model. We focus on the pathway to construct a gridded earthquake rate forecast that is readily employed for earthquake rate forecasting experiments and PSHA purposes. The main differences with respect to the original model by Hiemer et al. (2013) are: (1) we use an optimized variable-sized kernel instead of a fixed-sized kernel to estimate the spatial seismicity density following Helmstetter et al. (2007); (2) we introduce pseudo-prospective likelihood-based forecast experiments to constrain the weighting function between the seismicity- and fault-based spatial model components; (3) we estimated the a- and b-values of a truncated Gutenberg–Richter magnitude distribution based on a maximum likelihood approach that considers the spatial and temporal completeness history of the entire seismic catalogue; and (4) we adapt the model to treat subduction zones. Finally, we present properties and features of the resulting long-term earthquake rate forecast and its sensitivity to some model choices. We perform retrospective and pseudo-prospectively likelihood tests à la CSEP to evaluate and compare the consistency of all model's dimensions against independent catalogues providing moment magnitudes and compare the performance to the ASM rate forecast of the SHARE project (Giardini et al.2014).
2 DATA
2.1 Earthquake catalogue and declustering
We used the earthquake catalogue that was compiled in the framework of project SHARE and publicly released as the SHARE European Earthquake Catalogue (SHEEC; Grünthal & Wahlström 2012; Stucchi et al.2012; Giardini et al.2013; Grünthal et al.2013). By combining macroseismic data and instrumental seismological data in the period 1000–2006, the catalogue provides location and moment magnitude estimates m along with uncertainties for both (Fig. 1a). SHARE provided a completeness assessment of the catalogue within large completeness zones (Fig. 2) that were drawn based on historical constraints, that is considering the documentation history throughout Europe, not tectonic regionalization or seismic network constraints. These large zones consequently do not reflect many local variations in the completeness history, yet are a good approximation for the entire region. The SHEEC catalogue was declustered with a windowing technique (Grünthal 1985; Burkhard & Grünthal 2009).
(a) Epicentres of earthquakes within the study region from 1000 to 2006 (declustered SHARE European Earthquake Catalogue). Black polygons outline the 23 superzones for completeness assessment (cf. Fig. 2). (b) Traces of the upper edge of fault sources, projected to ground surface, for the entire Euro-Mediterranean area (colour-coded by maximum slip rate). Inset shows fault discretization example: grey lines correspond to depth contours. The estimated corresponding moment rate point sources are labelled by their focal mechanisms (A, downdip extent; sr, slip rate; ra, rake).
Magnitude of completeness histories within 23 superzones. Each horizontal bar is annotated by its corresponding magnitude of completeness value, and its background colour denotes the involved number of events. For a detailed description see SHARE deliverable D3.2 (http://www.share-eu.org).
We used the complete part of the declustered SHEEC catalogue as input to estimate the spatial and magnitude probability density of seismicity. Due to the limited quality of the depth information in the SHEEC catalogue, earthquakes in subduction zones and crustal seismicity were separated at 40 km depth. Accordingly, we applied the density estimation separately to crustal seismicity and subduction zone seismicity.
2.2 Database of seismogenic faults
The fault sources of the European Database of Seismogenic Faults (EDSF; Basili et al.2013) are the second input to our model. Building on previous experiences in this field (Basili et al.2008, 2009; Haller & Basili 2011), new standards have been adopted for the definition and characterization of active faults, including subduction zones, to ensure a homogenous data set for use in hazard assessment. This data set is homogeneous in the sense that all records are characterized by a common mapping strategy and the same set of parameters. For each fault parameter, the compilers documented the origin of the information being provided, along with the pertinent scientific reference to a publication, where applicable. Data uncertainties are handled by assigning a range of values to each parameter of a seismogenic source to capture its aleatoric variability (Basili et al.2008). The current version encompasses 1128 crustal faults with a slip rate of at least 0.01 mm yr−1 (Fig. 1b) and three subduction zones.
2.2.1 Crustal faults
The EDSF adopt the concept of the composite seismogenic source for crustal faults, which is a generalized, 3-D representation of a dipping surface in the Earth's crust on which fault slip occurs and where most of the seismic energy is released during an earthquake. Such fault sources are idealized as a uniformly dipping surface constrained between two horizontal lines that define the top and bottom edge of the source itself. Pairs of latitude–longitude geographic coordinates in decimal degrees with positive values for north/east define the locations of seismogenic sources. Depth values are positive downwards and are given in km. We applied a tessellation of the fault planes with elements of 5 × 5 km. Each fault element is associated with a set of geometrical and kinematic parameters, such as average strike, dip, rake and minimum and maximum annual slip rate; the latter are always positive and measured in millimetres per year along the direction of movement.
2.2.2 Subduction zones
Subduction zone seismicity is known to generate earthquake ruptures of several types with varying size and faulting styles at the slab interface or within the slab. These earthquakes follow different scaling laws (Strasser et al.2010) and rheology (Bilek & Lay 1999) from those of crustal earthquakes.
In the EDSF, a subducting plate is mapped as a collection of free-form subelements with a variable number of nodes. Each single subelement represents a portion of the entire surface of the subducting plate and is bounded by lines of constant depth except for the uppermost line when it coincides with the seafloor. Subduction sources are characterized by geometric (strike, dip and depth) and kinematic (rake, and slip rate taken equal to convergence rate) parameters. These parameters are given for all subelements together with their range of variability for the entire subduction zone (SHARE deliverable D3.4, www.share-eu.org). Thus the homogeneity and consistency of the EDSF data format facilitates the same treatment of both subduction and crustal fault data.
3 METHOD
3.1 Overview
3.2 Kernel density estimation
Kernel smoothing is a density estimation technique when no parametric density function is known. Thus it is a suitable technique to capture spatial clustering properties of observed seismicity. The locations of earthquakes are redistributed in space, where the kernel function and its bandwidth govern the shape and the amount of the redistribution (Stock & Smith 2002). Different kernels (e.g. Gaussian, power law) generate qualitatively similar densities, whereas the choice of an appropriate bandwidth has a crucial impact on the results (e.g. Silverman 1986; Wiemer et al.2009b). A spatially varying bandwidth accounts for the clustered nature of earthquake occurrence and depicts a better representation of the spatial distribution of seismicity than estimates from spatially invariant bandwidths (Stock & Smith 2002). The optimal local or regional kernel width is found through retrospective forecast experiments by dividing the data into a learning and target period (e.g. Zechar et al.2010a).
Location probability density for the declustered SHEEC catalogue based on two different kernel density estimation techniques: (a) constant kernel size of 10 km and (b) variable kernel size using the two spatially nearest events. The insets show the results of retrospective forecast experiments for optimizing the kernel size and the number of nearest neighbours, respectively. For two distinct 5-yr target periods, we used all m ≥ 5.0 earthquakes as target events. The probability gain per earthquake was estimated with respect to a spatially uniform earthquake density model as reference.
Location probability density for SHARE's database of crustal faults based on moment-weighted kernel density estimation. The upper inset shows the 9.2 per cent of all cells that contain 97.5 per cent of the density (as indicated by the dashed line in lower left-hand panel).
3.3 Spatial density weighting
The SEIFA model requires the combination of two distinct location probability densities μseis and μfault to a single spatial density μ. Hiemer et al. (2013) introduced a simple magnitude-dependent linear weighting function, such that the contribution from the fault-based density linearly increases from 0 to 1 with increasing magnitude (in the magnitude range of 6.5 ≤ m ≤ 8.0). Such a weighting function resembles the assumption that larger future earthquakes are more likely to occur in the vicinity of a fault. A magnitude-dependent location probability density leads to an apparent variability in the local frequency–magnitude distribution. Previously, the weighting function's linear slope and magnitude range was adjusted to minimize strong variations of local apparent b-values. Note that such an apparent b-value is not a b-value estimated from the catalogue, but it follows from the weighting of the scaled spatial probability density functions, with the scaling factors being normalized to match the global productivity, or a-value, of the study region.
(a) Seismicity-fault spatial density weighting as a function of target magnitude for the target period 1997–2007. In each magnitude bin the log-likelihood per target event is normalized with respect to its corresponding maximum. (b) Log-likelihood per target event for two different selected target magnitudes. Black symbols highlight maximum for different 10-yr target periods (including results shown in a). (c) Overview of results for all target periods. Grey line denotes weighting function (eq. 7).
3.4 Frequency–magnitude distribution
We assume that a single universal frequency–magnitude distribution is valid for the entire study area, thus our model does not a priori account for spatial variability in the relative occurrence of small and large magnitude events as expressed by the b-value. However, this variability is introduced by the weighting function as explained above. Traditionally, methods consider the temporal variations in the completeness history within one volume or study area (Weichert 1980; Wiemer et al.2009b; Kijko & Smit 2012). Given the spatial variability of the completeness periods within the SHEEC catalogue, we extended the method by Wiemer et al. (2009b) to estimate the global a- and b-values considering the space–time complexity of the completeness using a maximum-likelihood approach. The temporal completeness history is available within 23 superzones (SHARE deliverable D3.2, Fig. 2). We combined the histories of all zones to estimate the best overall b- and a-values for the entire catalogue as: b = 0.9 and a = 5.87 (Fig. 6a). The magnitude range of our model is 4.5–8.6. The minimum magnitude of 4.5 was adopted as the lower magnitude of engineering interests while the 8.6 corresponds to the upper magnitude distribution across the SHARE target region. We quantify the uncertainty in the regional estimate of the a- and b-values with a bootstrap approach. To generate a bootstrap sample of the catalogue, we first select events in a given superzone, divide the catalogue in the completeness periods, bootstrap the catalogue within the completeness periods and then recreate the catalogue for the entire time span of the catalogue. This procedure ensures that the infrequent large events are well represented when bootstrapping. For the 1000 bootstrapped samples, we obtained quantiles for both, the b-value distribution (q0.25 = 0.89, q0.75 = 0.91) and a-value distribution (q0.25 = 5.85, q0.75 = 5.93; Fig. 6a).
Frequency–magnitude distribution (FMD) for the crustal part (a) and subduction part (b) of the SHEEC catalogue taking into account regional magnitude of completeness histories. The legend lists the a- and b-value together with the corresponding interquartile ranges derived from the 1000 bootstraps of the catalogue.
4 RESULTS
4.1 Earthquake rate forecast
We applied the method to crustal seismicity (depth ≤40 km) and subduction zone seismicity separately. The depth value for differentiating between deep and crustal seismicity is a compromise as the focal depth uncertainty in the catalogue is large in parts or very often not available; the latter is especially true for many larger historical events based on macroseismic intensities, but also for those located through the recordings of seismic networks. For both earthquake types, we estimated spatial probability densities and the overall a- and b-values of the truncated Gutenberg–Richter magnitude distribution. We first discuss the results of the crustal seismicity and then outline the adaption for the subduction seismicity.
4.1.1 Modeling crustal seismicity
We compared two different kernel-density estimation techniques for the seismicity component of our model: using a constant kernel size and a variable kernel size. We emphasize that the kernel size is not a spatial resolution parameter of past seismicity, but an expression of the degree of stationarity of seismicity. The kernel bandwidth optimization is based on retrospective forecast experiments; it consequently aims to estimate a spatial density that well describes the spatial distribution of future events. Such an optimization accounts for the prevalent observation that future events do not occur at the exact same location as past events. Our likelihood analysis supports the use of a constant kernel size considering target time spans of 5 yr as the likelihood gains per target event are slightly higher (insets Fig. 3). Nevertheless, we argue that the sensible use of a variable kernel size is more appropriate. The resulting density adverts to its advantages (Fig. 3b): On the one hand, a smaller bandwidth is more appropriate in high-density regions since the large number of earthquakes allows for a more accurate estimation of the density itself. On the other hand, a larger bandwidth is more appropriate in low-density regions where only a few events are observed (e.g. Northern Europe, Fig. 1a). The constant kernel spatial density results in a punctuated forecast, while the variable kernel one distributes event contributions to larger regions (e.g. Northern Europe, Fig. 1b). The latter resembles the arguments by Kafka (2007) and Liu et al. (2011), who argue that seismicity in intraplate regions may jump between fault system and the proximity to past events is in contrast to high-seismicity regions not the best guide for forecasts. The variable kernel technique is thus suitable for smoothing long tail distributions where under-smoothing in the tails is likely to cause difficulties (Silverman 1986). Our decision is supported by findings of Stock & Smith (2002) for New Zealand and the outcome of the RELM project (Schorlemmer et al.2010; Zechar et al.2013), where the forecast submitted by Helmstetter et al. (2007) showed the highest skill with respect to all other candidate models.
We converted the resulting spatial densities into spatial distributions of annual earthquake rates by scaling each spatial density with its corresponding activity, respecting the total productivity of the global Gutenberg–Richter magnitude distribution (Fig. 6a). Accordingly, the final annual rate of our forecast is purely based on an estimate from the earthquake catalogue data, whereas its spatial component incorporates contributions from both earthquake and fault moment-rate densities. The resulting earthquake rate model forecasts a total cumulative number (m ≥ 4.5) of 65.6 events yr−1 (cf. Fig. 6a). Fig. 7(a) shows the cumulative annual earthquake rate for all cells in our model space. The cumulative rates are dominated by the contributions from smaller magnitudes; thus, their spatial distribution is primarily controlled by the seismicity-based location probability density. Mapping incremental annual rates with increasing magnitude thresholds, the fault contribution increases and the spatial distribution of the rates is less affected by past earthquake occurrences (Fig. 7b). For the selected magnitude range 4.5 ≤ m ≤ 4.6, the spatial rate distribution resembles the one in Fig. 7(a), however, for the magnitude range 7.5 ≤ m ≤ 7.6, the seismically active faults in Greece and Turkey light up prominently compared to the other areas in Europe.
(a) Spatial distribution of cumulative annual earthquake rates (4.5 ≤ m ≤ 8.6). The total annual rate is 65.6 events yr−1 (cf. Fig. 6a). (b) Spatial rate distribution for selected incremental magnitudes. Note the increasing imprint of the fault contribution with increasing magnitude (cf. Fig. 5). (c) Resulting local departures from overall Gutenberg–Richter magnitude distribution (b = 0.9, cf. Fig. 6a).
Fig. 7(b) underlines that the final spatial component of our model is magnitude dependent (eq. 2). In regions of mapped faults, the contribution from past seismicity dominates the weighted location probability density for small events (4.5 ≤ m ≤ 5.5, 76 per cent seismicity +24 per cent faults), and the contribution from fault moment-rate points is predominant for large events (m ≥ 7.5, 20 per cent seismicity +80 per cent faults). The imposed functional form of the weighting function (Eq. 7) assures that local deviations from a plain Gutenberg–Richter distribution are minimized (Fig. 7c, see also Hiemer et al.2013, their fig. 6). We find that in particular the fast slipping faults deviate to smaller b-values while slow slipping faults show slightly increased values. The deviations are limited to a relatively small range of 0.2 (Fig. 7c).
4.1.2 Modelling seismicity in subduction zones
We applied the same approach to the subduction zone seismicity for events deeper than 40 km (Fig. 8). Due to the uncertainties in the focal depth, we did not differentiate between interface and intraslab seismicity, thus we assume that there is no difference in the spatial variation of these two types of earthquakes only due to the limitations in the given data set. We estimated spatial density probabilities based on the earthquake data and on the converted fault moment rate points, respectively (Fig. 9a). Smoothing was performed on the surface projection of the subduction zones considered (Calabrian Arc, Hellenic Arc and Cyprus Arc).
Input data for the subduction part of our model (top panel shows seismicity: event depth >40 km, bottom panel shows fault database: Calabrian Arc, Hellenic Arc and Cyprus Arc).
(a) Normalized spatial densities for the subduction zone model (top panel for seismicity, using a constant size kernel of 10 km; and bottom panel for the fault database, using a moment-weighted kernel). (b) Spatial distribution of cumulative annual earthquake rates (4.5 ≤ m ≤ 8.6). We used the same weighting function as for the crustal part of our model to weight the two spatial probability densities. The total annual rate is 8.1 events yr−1 (cf. Fig. 6b).
The seismicity-based density is less homogeneous than the subduction zone moment rate density, representing the spatial variability of seismicity and its apparent clustering. Seismicity density in the Calabrian Arc extends further to the northeast compared to the moment rate contribution. A similar feature is observed in the eastern Hellenic Arc, where the seismicity seems clustered from the eastern end of Crete to the southwestern border of Turkey. Considering the comparatively limited amount of data, we used the same weighting function as for the crustal part of our model for internal model consistency. We estimated the a- and b-values of a truncated Gutenberg–Richter magnitude distribution (a = 5.21, b = 0.95, and corresponding interquartile ranges from bootstrapping, Fig. 6b) for the complete part of the data set and scaled the weighted densities accordingly. Our resulting subduction model forecasts a total cumulative number (m ≥ 4.5) of 8.1 events yr−1. As highlighted by the densities (Fig. 9a) and the cumulative annual occurrence rates (Fig. 9b), the model's spatial component accentuates the location of observed seismicity that is laterally heterogeneous. The variability reduces for the larger magnitude events. Note that relatively large uncertainties due to the focal depth information remain for where events actually occur in the subduction zone. Resolving separate densities for inslab and interface events will only be possible for high-resolution catalogues including the characterization of each event with a moment tensor or a focal mechanism.
4.2 Consistency tests of earthquake rate forecasts
Procedures for testing and evaluating earthquake likelihood models have been established within CSEP (e.g. Jordan 2006; Zechar et al.2010b; Eberhard et al.2012), which is a community-supported infrastructure for conducting forecasting experiments in several regions around the world. Our model is expressed as a gridded earthquake rate forecast and suitable for such likelihood testing procedures: It specifies the annual expected number of earthquakes within non-overlapping spatial cells (256 210 cells, size of 0.1° × 0.1°), each of which is composed of 41 magnitude bins (4.5 ≤ m ≤ 8.6, 0.1 magnitude binning), and the model incorporates the assumption that the number of earthquakes in each bin is Poisson-distributed and independent of those in other bins.
We tested the consistency (1) between the total number of observed and forecasted events (CSEP: N-test), (2) between the observed and forecasted magnitude distribution (CSEP: M-test), (3) between the observed and forecasted spatial distribution (CSEP: S-test) and (4) between the observed and forecasted joint space–magnitude distribution (CSEP: conditional L-test, Werner et al.2010). These consistency tests are all governed by the same principle. One simulates many catalogues (in this study 10 000) that are consistent with the forecast and estimates the corresponding distribution of joint log-likelihoods (eq. 6). This distribution is then compared with the joint log-likelihood of the forecast given the observed data by computing a quantile score, that is the fraction of simulated likelihoods smaller or equal the observed one. The forecast is considered inconsistent with the observed data if the quantile score is lower than a critical significance value, which is commonly set to 0.05 (corresponds to 0.95 confidence in the results, for example Eberhard et al.2012). In case of a Poisson forecast, the N-test has an analytical solution and thus does not require simulations (Zechar et al.2010b). The objective of the S-test is to consider only the spatial distribution of the forecast and the observation. To isolate the spatial information, one computes the sum of all magnitude bins, and the resulting forecast sum is normalized so that it matches the observed number of target events. Similarly, the M-test only considers the magnitude distribution by summing over longitude–latitude spatial cells.
We tested four different earthquake rate forecasts: the SEIFA model presented, its single seismicity-based component (spatial density equals the seismicity density, SEI), its single fault-based component (spatial density equals the fault moment rate density, FA), and the classical ASM that is an alternative seismic source model in SHARE's logic tree for PSHA (Giardini et al.2013, http://www.efehr.org/). The first three only differ in their spatial earthquake rate distribution (implying identical results for N- and M-test), thus we used the S-test to detail their individual consistency.
We define target events as moment magnitudes m ≥ 5 crustal events (depth ≤ 40 km). We created two sub-catalogues from the SHEEC catalogue for retrospective likelihood testing, which span the time periods 1987–2007 (20 yr) and 1997–2007 (10 yr), respectively. The fact that the SHEEC catalogue ends in 2007 allows for pseudo-prospective consistency testing considering the time period of 2007–2013; these independent 6 yr of data were not incorporated in the initial model construction. For that purpose we used target events from both the CMT catalogue (Global Centroid-Moment-Tensor, http://www.globalcmt.org/CMTsearch.html) and the NEIC catalogue (National Earthquake Information Center, http://earthquake.usgs.gov/earthquakes/eqarchives/epic/) as these provide also moment magnitudes. The latter catalogues were declustered with the same window approach as the SHEEC catalogue prior to testing.
The total expected annual number of m ≥ 5 events within the study region equals 29.9 for the ASM and 23.3 for the SEIFA model. Considering the N-test we find that the ASM is inconsistent with the target events in all four investigated target catalogues (i.e. it overpredicts the observed total annual number, Table 1). However, the ASM shows a close match with the observed rates of the entire SHARE catalogue (Fig. 10). Notice that the model fits an a- and b-value in each area source, and each of those values were manually checked and adjusted to match the rates of observed events between m = 5.5 and 7.5. Thus the ASM is mainly constrained by the historical part of the catalogue with special consideration for a certain magnitude range, which explains the elevated rates in its overall summed frequency–magnitude distribution (Fig. 10). For the SEIFA model, the corresponding annual expected number of m ≥ 5 events is consistent with both, the pseudo-prospectively observed ones and the total observed number during the last 10 yr of the SHEEC catalogue (Table 1, Fig. 10). Considering the last 20 yr, however, the SEIFA model also overpredicts the total annual number (Table 1). With the M-test, we check for consistency of the observed and normalized modelled frequency–magnitude distribution so that the total number does not influence the result. We find that the magnitude dimension of both tested models are consistent with all four investigated target catalogues (see Table 2 for respective quantile scores). We obtained similar results for the conditional L-test, where we tested the joint magnitude–space dimensions normalized to the respective observed number of events (Table 3).
Frequency–magnitude distribution for SHARE's area source model (ASM) and smoothed seismicity-fault model (SEIFA). Both are consistent with observed magnitude distributions. The likelihood consistency tests were performed pseudo-prospectively using CMT/NEIC target events and retrospectively using SHEEC target events. We considered all m ≥ 5 earthquakes (depth ≤ 40 km) as target events.
Forecasted vs. observed number of m ≥ 5 target events (depth ≤ 40 km). Asterisks indicate failure of the N-test at a critical significance value of 0.05, that is inconsistency between the total forecast expectation and the number of target events. See Fig. 11 for the corresponding spatial distributions of target events.
| Observed/ | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| forecasted | (6 yr) | (6 yr) | (10 yr) | (20 yr) |
| SEIFA | 121/139.7 | 133/139.7 | 205/232.8 | 382/465.7* |
| ASM | 113/179.5* | 122/179.5* | 205/299.2* | 382/589.3* |
| Observed/ | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| forecasted | (6 yr) | (6 yr) | (10 yr) | (20 yr) |
| SEIFA | 121/139.7 | 133/139.7 | 205/232.8 | 382/465.7* |
| ASM | 113/179.5* | 122/179.5* | 205/299.2* | 382/589.3* |
Forecasted vs. observed number of m ≥ 5 target events (depth ≤ 40 km). Asterisks indicate failure of the N-test at a critical significance value of 0.05, that is inconsistency between the total forecast expectation and the number of target events. See Fig. 11 for the corresponding spatial distributions of target events.
| Observed/ | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| forecasted | (6 yr) | (6 yr) | (10 yr) | (20 yr) |
| SEIFA | 121/139.7 | 133/139.7 | 205/232.8 | 382/465.7* |
| ASM | 113/179.5* | 122/179.5* | 205/299.2* | 382/589.3* |
| Observed/ | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| forecasted | (6 yr) | (6 yr) | (10 yr) | (20 yr) |
| SEIFA | 121/139.7 | 133/139.7 | 205/232.8 | 382/465.7* |
| ASM | 113/179.5* | 122/179.5* | 205/299.2* | 382/589.3* |
M-test quantile scores (the fraction of simulated joint log-likelihoods smaller or equal the observed one). Both model's magnitude distributions are consistent with the observed data, because their quantile scores are greater than the significance value of 0.05. See Fig. 10 for the corresponding frequency–magnitude distributions of the models and the target events.
| M-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.88 | 0.64 | 0.63 | 0.18 |
| ASM | 0.87 | 0.47 | 0.78 | 0.30 |
| M-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.88 | 0.64 | 0.63 | 0.18 |
| ASM | 0.87 | 0.47 | 0.78 | 0.30 |
M-test quantile scores (the fraction of simulated joint log-likelihoods smaller or equal the observed one). Both model's magnitude distributions are consistent with the observed data, because their quantile scores are greater than the significance value of 0.05. See Fig. 10 for the corresponding frequency–magnitude distributions of the models and the target events.
| M-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.88 | 0.64 | 0.63 | 0.18 |
| ASM | 0.87 | 0.47 | 0.78 | 0.30 |
| M-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.88 | 0.64 | 0.63 | 0.18 |
| ASM | 0.87 | 0.47 | 0.78 | 0.30 |
Conditional L-test quantile scores (the fraction of simulated joint log-likelihoods smaller or equal the observed one). Both model's magnitude–space dimensions are consistent with the observed data (at the 95 per cent confidence level).
| L-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.64 | 0.63 | 0.44 | 0.65 |
| ASM | 0.34 | 0.68 | 0.47 | 0.82 |
| L-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.64 | 0.63 | 0.44 | 0.65 |
| ASM | 0.34 | 0.68 | 0.47 | 0.82 |
Conditional L-test quantile scores (the fraction of simulated joint log-likelihoods smaller or equal the observed one). Both model's magnitude–space dimensions are consistent with the observed data (at the 95 per cent confidence level).
| L-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.64 | 0.63 | 0.44 | 0.65 |
| ASM | 0.34 | 0.68 | 0.47 | 0.82 |
| L-test | CMT | NEIC | SHEEC | SHEEC |
|---|---|---|---|---|
| (6 yr) | (6 yr) | (10 yr) | (20 yr) | |
| SEIFA | 0.64 | 0.63 | 0.44 | 0.65 |
| ASM | 0.34 | 0.68 | 0.47 | 0.82 |
We devote more attention to the results of the consistency test for the spatial earthquake rate distribution, because the corresponding S-test allows for investigating the reliability of individual SEIFA model components (SEI, the purely seismicity-based, and FA, the purely fault-based end-member). We find that the spatial components of the ASM, SEIFA and SEI model are consistent with the spatial distribution of observed earthquakes applying the S-test (middle row panels in Fig. 11). The ASM model shows the highest joint log-likelihood per event given the observations. Note that for the retrospective cases (two right columns, Fig. 11) the target events were included when constructing the ASM.
Results for consistency test of spatial rate distribution of models (ASM, area source model; SEIFA, smoothed seismicity-fault model; SEI, smoothed seismicity end-member; FA, smoothed fault end-member). Each column corresponds to a different target catalogue using m ≥ 5 events (depth ≤ 40 km). The middle rows show consistency test results for all target events (total number N1, black circles and grey diamonds on top row map) and bottom rows for target events in the vicinity of mapped faults (number N2, black circles only). The inverted triangle represents the spatial log-likelihood/event of the observation given the forecast. The horizontal grey bars delineate the top 95 per cent of log-likelihoods from simulated catalogues that are consistent with the corresponding forecast. Non-filled triangles indicate a failure of the spatial consistency test, that is the respective forecast is inconsistent with the observed spatial distribution of target events.
To further detail differences, we created a subset of target events. We therefore considered only the target data that are located within the cells which were subjected to seismicity-fault density weighting, that is events that happened on or nearby a mapped fault (black circles in top row panels of Fig. 11). In these regions the ASM and SEIFA model are the only forecasts being fully consistent with all considered target catalogues (bottom row panels, Fig. 11). The target catalogues include surprises, i.e. events that did not happen on previously mapped faults and events that occurred in regions of low density of past-observed seismicity. Thus the individual components of the SEIFA model do not prove consistency when applying the S-test. The SEIFA model, however, accounts for such surprises, because it incorporates a weighted combination of both, the seismicity-based and fault-based location density. We investigated likelihood ratio maps for all model components with respect to a uniform forecast to underline our findings (Fig. 12). Such maps identify which spatial locations of earthquakes are less likely given a forecast compared to a uniform forecast (Werner et al.2010). The spatial inconsistency of the FA model is due to its strong concentration of event location density (compare Fig. 4 and Fig. 3b). The model suffers from losses in likelihood due to numerous off-fault events that are not compensated by equal or greater likelihood gains from well-located on-fault events. Note that we did not account for location uncertainties when evaluating the likelihoods of the individual target events.
Likelihood ratios between the (a) SEIFA, (b) SEI and (c) FA spatial forecast and a uniform forecast, respectively. All forecasts are normalized to sum to the 220 observed earthquakes (SHEEC 20 yr, see Fig. 11, bottom right-hand panel). Locations of earthquakes in blue areas are less likely given the forecast compared to the uniform forecast.
We performed the same consistency test for the subduction part of our model (Fig. 9b) and the model passed all outlined consistency tests at the 95 per cent confidence level. In comparison with the tests of the crustal part of our model, however, these tests are less reliable due to the limited number of target events, for example, 11 events within SHEEC's target time span of the last 10 yr and 30 events within its last 20 yr.
The consistency tests do not allow for addressing the question whether one model is better than the other. We applied the T-test (Rhoades et al.2011) to evaluate the significance of the observed spatial likelihood differences between the ASM and SEIFA model. The T-test is based on the Student's paired t-test to estimate whether the corrected average information gain per earthquake is significantly different from zero for a given pair of models or not. We find for all our considered target catalogues, that the SEIFA model has a positive information gain over the ASM. The likelihood gains are significant for the retrospective cases (filled circles, Fig. 13). We obtain the same results when using the W-test, which applies the Wilcoxon signed-rank test to the individual rate-corrected information gains and thus does not assume that they are normally distributed (as the T-test; Rhoades et al.2011).
Results for comparison T- and W-test of spatial rate distribution of models (ASM, area source model; SEIFA, smoothed seismicity-fault model). Each plot shows the mean and 95 per cent confidence interval of the information gain per earthquake of the SEIFA model over the ASM model for a different target catalogue. N1 is the total number of target events, and N2 is the number of target events in the vicinity of mapped faults (compare Fig. 11). Filled circles denote T-test significance at 95 per cent confidence level. The corresponding W-test significance is indicated by a small ‘W’ next to the average.
We emphasize that the crustal parts of both the ASM and the SEIFA model are consistent with the observed seismicity in the space and magnitude dimension considering target catalogues of different time spans and data sources. The SEIFA procedure highlights the importance of the instrumental part of the SHEEC catalogue, which is the reason why it successfully forecasts the total number of CMT and NEIC events. The ASM approach accentuates the historical part of the SHEEC catalogue, which results in a close agreement with the observed frequency–magnitude distribution of the entire SHEEC catalogue. Whether the emerged bulge in the magnitude range of m6 – m7.5 is an artifact or a feature of the Pan-European region remains to be understood. The testability of the models will enable to revisit the models and prospectively test their performance using much longer time spans within the European CSEP Testing Center.
5 CONCLUSIONS AND DISCUSSION
We presented a kernel smoothing approach applied to seismicity and fault moment release to model earthquake activity throughout Europe. The method is a zone-less alternative to other earthquake source models and has been used to capture epistemic uncertainties in the SHARE source model logic-tree. Our model extends the classical seismicity kernel-smoothing method to include the moment release information from faults. This extension addresses the strongest criticism directed to kernel-smoothing approaches: the disregard of slip rates and fault geometries. Furthermore, smoothed seismicity models are the antipode to the seismic gap model, in which large earthquakes are not expected where they had occurred in the past. Our approach includes slip rates on active faults and historically quiet faults. Note that we use this data component only for the estimation of the spatial distribution of earthquake rates. The overall productivity of our model is based on a global scaling law relation fit to the entire catalogue proposed here also for the first time. Regions, where no structures are mapped, are solely modelled by smoothing seismicity data, accounting for the present incomplete knowledge about the fault network.
Our approach is data-driven and aims to reduce the number of subjective judgments. The largest subjective influence arises from the determination of the frequency–magnitude distribution—we used a global approach for the crustal and the subduction zone separately to minimize its influence. Most influential in this respect is the determination of a regional space–time completeness history, which might be more difficult to solve in low than in high seismicity regions (Stirling & Gerstenberger 2010). We used retrospective forecast experiments to optimize the model parameters for the spatial components of our model: the optimal kernel size and the optimal function for weighting the seismicity-based and fault-based probability densities. For the magnitude dimension, we assumed that a truncated Gutenberg–Richter distribution is appropriate to describe the entire earthquake catalogue. The estimated a- and b-values account for regional differences in the temporal evolution of the magnitude of completeness.
Major improvements compared to the model by Hiemer et al. (2013) are: (1) we estimated the spatial seismicity density using optimized variable-sized kernels; (2) we used retrospective forecast experiments to optimize a weighting function between the seismicity- and fault-based spatial densities; (3) we estimated the a- and b-value of a truncated Gutenberg–Richter magnitude distribution based on a maximum likelihood approach that considers the spatial and temporal completeness history of the entire seismic catalogue using the approach by Wiemer et al. (2009b) for each single zone and (4) in addition we applied the technique for the first time to a subduction zone model, introducing spatially varying activity rates on a complex fault model.
Our results are affected by both the completeness and uncertainties of the input data. The main limitation is the intrinsic incompleteness of the fault database (e.g. Basili et al.2013) especially outside the area of well known plate boundaries (England & Jackson 2011). We used the upper 97.5 per cent of the fault-based density to represent areas of mapped faults (24 381 of all 265 120 cells, inset Fig. 4). We repeated our analysis using values of 90 per cent (12 599 cells), 95 per cent (18 248 cells) and 99 per cent (33 545 cells) to estimate the sensitivity of our results with respect to this choice. We found that the value does not alter our main findings presented in Figs 11 and 13: The forecasts are consistent with the target observations and they show a positive information gain over the ASM. Limitations due to data incompleteness are similarly true for the seismicity component, in particular for low-seismicity regions. The kernel-smoothing approach captures uncertainties in earthquake location, fault geometry and slip-rate estimate. Magnitude uncertainties have no impact on the kernel smoothing, because the kernel treats all magnitudes identically. However, they might lead to missing or additional events within completeness windows. Thus, the magnitude dimension of our model is affected by these uncertainties. The accuracy of the estimate of the regional a- and b-values of the entire SHEEC catalogue depends on the reliability of the previously conducted completeness-time assessment within SHARE's predefined superzones. The overall annual rate of our model is entirely driven by the a- and b-values estimate.
We estimated the optimal kernel bandwidth for the seismicity component of our model, but did not optimize the kernel bandwidth for the fault component individually. We optimized the weights between both densities assuming a weighting function's functional shape that aims at minimizing local deviations from a plain Gutenberg–Richter distribution. It would be insightful to estimate optimal kernel bandwidths and their weights jointly. Such a joint optimization should include constant, variable and moment-dependent kernels in order to reveal which approach yields the highest likelihood gain for which part of the data. The number of free parameters for the weighting function should be subjected to penalized-likelihood investigations as part of this optimization. Another important question for future studies is to determine whether including the historical part of the catalogue leads to a better forecast or not. We also used the entire declustered SHEEC catalogue to comply with the large-scale hazard community effort of the SHARE project. Wang et al. (2011) examined the importance of several assumptions and choices when constructing a smoothed seismicity forecast for California. They found that the inclusion of historical earthquakes failed to improve the forecast consistently.
The model concept is modular: the spatial densities and the magnitude density are exchangeable without changing the general flavour of the model. We used seismicity and geologic fault information for its spatial component, however, input from strain-rate models or recent geodetic observations could be implemented within the first step of the density estimation procedure. For each probability density function, a modeler can select the contributing data sets, as long as the choices are self-consistent. The time-independent implementation presented in this paper fits many requirements of time-independent probabilistic seismic hazard assessment by taking into account community-driven information, such as the regional completeness history of the cataloguer and the tectonic characterization of the study region. However, such assessments could be conducted independently (e.g. Hiemer et al.2013, for tectonic information weighting).
Our model is the first European zone-less forecast model that includes seismicity and fault information. We showed that it is a reliable and skillful alternative to the classical ASM; therefore we propose this to serve as a reasonable candidate reference model for the European region—as an earthquake rate forecast model as well as a seismic hazard model.
6 RESOURCES
All presented annual rate forecasts and their documentations are available at http://www.efehr.org/ (European Facility for Earthquake Hazard and Risk).
The authors would like to thank Max Werner and an anonymous reviewer for detailed constructive reviews, that helped us to improve many parts of the paper. We thank Damiano Monelli for feedback concerning SHARE's hazard computation. All figures were prepared with GMT version 4.2.0. This work has received funding from the EC-Research FP7-projects, SHARE, under grant agreement No. 226967 and NERA, under grant agreement No. 262330.
