1 Introduction

Probabilistic seismic hazard analysis (PSHA) is usually carried out based on the framework described in Cornell (1968) and Esteva (1968). The analysis requires consideration of the seismic source model, magnitude-recurrence relation, and attenuation relation (that is, the ground motion prediction equation or ground motion model). The framework can be illustrated in a logic tree if combinations of the source zone, magnitude-recurrence relation, and ground motion model are considered, where the estimation of the seismic hazard is performed by folding the tree by applying the total probability theory (McGuire 2004). This can be carried out by using the numerical integration technique (McGuire 2004) or simulation technique (Musson 2000; Hong et al. 2006). The seismic hazard assessment results are important for seismic structural reliability analysis and the selection of optimal seismic design (Rosenblueth 1976; Goda and Hong 2006; Hong and Hong 2007).

Many advances have been reported on the assignment of the source model, modeling and characterizing magnitude-recurrence relations, and the development of ground motion models. The advances include a better understanding of the earthquake occurrence mechanism, the ground motion characteristics, and the statistical techniques used to develop the source model, the frequency distribution of the magnitude and the occurrence rate, and the ground motion models (Gerstenberger et al. 2020). The history of the application of the PSHA framework to estimate Canadian seismic hazard maps (CanSHMs), which were implemented in Canadian structural design codes, was presented in Adams (2019) (see also Hong et al. 2006). According to Adams (2019), Natural Resources Canada (NRCan) (and its predecessors) has been responsible for six generations of seismic hazard maps for Canada since 1953. The mapped Canadian seismic hazard formed the basis for the seismic design provisions recommended in the National Building Code of Canada (NBCC). These studies pointed out that the last three generations of CanSHMs are developed based on a source model that consists of a set of delineated source zones. The source model defined by a set of source zones (Kolaj et al. 2019; Kolaj et al. 2020) was used to develop the 6th generation of CanSHM.

The boundary of each source zone could be determined using existing geological or geophysical knowledge. Each source zone is assumed to have similar seismicity characteristics. However, defining the boundary may introduce uncertainty and could potentially involve a circular argument (Hong et al. 2006). The use of spatially smoothed seismic source models (that is, seismic source model obtained based on the historical earthquake catalogue and spatial smoothing techniques) was considered by several researchers (Kagan and Jackson 1994; Frankel 1995; Woo 1996; Stock and Smith 2002; Hong et al. 2006; Xu and Gao 2012; Feng and Hong 2021). The smoothed seismic source model is obtained by spatially smoothing the historical earthquake occurrence using a kernel function with a fixed or adaptive bandwidth. The smoothing of the spatial occurrence could be carried out by neglecting or considering the earthquake magnitude dependency (Frankel 1995; Woo 1996).

By considering the ground motion models used to develop the 4th generation of CanSHM, Hong et al. (2006) showed that the degree of smoothing affects the estimated seismic hazard. The smoothing in their study was carried out using the Voronoi polygon or the spatial kernel smoothing with a constant bandwidth, and the seismic hazard assessment was carried out only for two regions in Canada. Systematic quantification of differences in the estimated and mapped seismic hazard at Canadian sites that are obtained based on the delineated seismic source model and spatially smoothed seismic source models, however, is unavailable, although the quantification is valuable to identify epistemic uncertainty in the estimated seismic hazard. One may argue that the use of the spatial smoothing approach could not incorporate geological and tectonic features while the use of the delineated source modeling approach is often related to the subjective definition of the seismic source zones. However, the latter was refuted by arguing that a circular argument exists between the selection of seismicity patterns and the shape of tectonic plates (Hong et al. 2006). Therefore, both modeling approaches are valid for seismic hazard modeling and a comparison of the results obtained from these approaches is necessary to quantify their differences.

Note that the ground motion models used for the mapping of the 5th and 6th generations of CanSHMs (denoted as CanSHM5 and CanSHM6) differ. The ground motion models used for CanSHM6 were summarized in Kolaj et al. (2020). They also indicated that different software was used for the seismic hazard assessment for CanSHM5 and CanSHM6.

The main objectives of the present study are (1) to map the Canadian seismic hazard by applying the spatially smoothed source models and (2) to quantify the differences between CanSHM6 and that obtained by using the spatially smoothed source models. For consistency, in all cases, the ground motion models used to develop CanSHM6 are adopted in the present study for the seismic hazard assessment. We used two approaches to develop the spatially smoothed source models. The first one uses the kernel smoothing technique with a constant bandwidth and the second one uses the adaptive kernel smoothing technique. In both cases, the smoothing was carried out by considering the earthquake magnitude dependency. The considered historical earthquake catalogue, kernel function, and estimation of the bandwidth are described in the following section. This is followed by the description of the estimation and mapping of the seismic hazard and the comparison of the estimated seismic hazard.

2 Application of Spatial Smoothing Techniques to Develop Smoothed Seismic Source Models

In this section, we summarize the historical earthquake catalogue used in the present study. Based on such a catalogue, we develop the spatially smoothed seismic source model by using fixed as well as adaptive kernel smoothing techniques. Moreover, we compare the magnitude-recurrence relation obtained based on the spatially smoothed source models to that obtained based on the delineated seismic source model that was used for the development of CanSHM6.

2.1 Historical Earthquake Catalogue and Its Completeness Assessment

The historical earthquake catalogue for developing the spatially smoothed source model in the present study is from the 2011 Canadian Composite Seismicity Catalogue (CCSC).Footnote 1 The catalogue consists of the 2011 CCSC West Catalogue and the 2011 CCSC East Catalogue. The 2011 CCSC West Catalogue contains a list of 34,934 historical seismic events that occurred from January 1700 to December 2010, where the epicentral location of each event was within the latitude and longitude of 43.5° N–75° N and 160° W–110° W. The 2011 CCSC East Catalogue contains a list of 11,893 events that occurred from January 1534 to December 2010, where the epicentral location of each event was within the latitude and longitude range of 35° N–80° N and 110° W–45° W. The consideration of this catalogue rather than a more up-to-date catalogue is justified since Adams et al. (2019) indicated that the (delineated) seismic source model used in CanSHM6 remains mostly the same as that used in CanSHM5, which was developed by considering the seismic events that occurred before 2011 (Halchuk et al. 2015), and that the estimated seismic hazard is relatively insensitive to the additional seismic events since the development of the source model used in CanSHM5. In other words, the use of CCSC facilitates the comparison of the assessed seismic hazard obtained by using the delineated seismic source model and the spatially smoothed source model, which are developed based on the same historical seismic catalogue.

For each event, the catalogue provides the occurrence time (year/month/day/hour/minute), epicentral location (latitude and longitude), focal depth (km), magnitudes (for example, body wave magnitude Mb, Nuttli magnitude MN, local magnitude ML, surface wave magnitude MS, coda magnitude MC, duration magnitude MD, and moment magnitude M), and the original catalogues that reported the event. In the catalogue, M is used as the most preferred magnitude, and the reported Modified Mercalli Intensity (MMI) is adopted to estimate M for early earthquake events due to the unavailability of seismic instruments. M is used for the present study. The historical seismic events that occurred within the spatial coverage, including those that occurred outside of Canada but within about 500 km distance of the Canadian border, are employed to develop the spatially smoothed source model and shown in Fig. 1.

Fig. 1
figure 1

Spatial distribution of historical earthquake events with M ≥ 3.0, where the red and blue color symbols denote the events reported in the East Catalogue and West Catalogue, respectively

Full size image

A detailed explanation of CCSC was given by Fereidoni et al. (2012). The number of induced seismic events in this catalogue is small and with a small magnitude. Since according to Halchuk et al. (2015), the removal of such events is difficult if not impossible, no attempt was made to remove such events. A preliminary process for the catalogue is carried out in the present study. The earthquake events for the 2011 CCSC East Catalogue and 2011 CCSC West Catalogue with occurrence time earlier than 1660 and 1860 are removed based on the first year of complete coverage for time-magnitude categories (Hong et al. 2006). Note that the historical earthquake catalogue is associated with unequal observation periods for different earthquake magnitude intervals. This affects the evaluation of the magnitude-recurrence relation (Weichert 1980). However, the unequal observation periods are unknown and must be assessed through the completeness analysis (Albarello et al. 2001). For the assessment, consider that the catalogue is reported up to the most recent time TR, and the earliest time when an event with M greater than or equal to a threshold magnitude MT is TMT,S. Let ∆Tmax = TR − TMT,S denote the maximum observation period for the considered MT. The completeness analysis of the catalogue for events with magnitude greater than Mmin,j since or after the time TC,j, for j = 1, …, NC, is to be carried out, where NC is the total number of lowest magnitude cases considered for completeness analysis, and TC,j ∈ [TR, TMT,S). Given Mmin,j, ΔTC,j = TR − TC,j is divided into 2N elementary non-overlapping subintervals of equal duration δt = ΔTC,j/(2N), and that a comparison of the earthquake occurrence rates in the i-th and (N + i)-th intervals, denoted as λi and λN+i, is carried out. According to Albarello et al. (2001), the probability of completeness of the catalogue for a given value of TC,j (that is, within the specific time span ΔTC,j), \(P_{j} (C\left| {T_{C,j} } \right.)\), is given by:

$$P_{j} \left( {C\left| {T_{C,j} } \right.} \right) = \frac{{\Delta T_{C,j} }}{{\phi \Delta T_{\max } }} \times \left( {0.5} \right)^{{N^{\prime}}} \times \sum\limits_{k = m}^{{N^{\prime}}} {\frac{{N^{\prime}!}}{{k!(N^{\prime} - k)!}}} ,$$
(1)

where ϕ is a constant equal to or greater than 1; \(N^{\prime}\) denotes the total number of cases where \(\lambda_{i} \ne \lambda_{N + i}\) for i = 1, …, N; m represents the observed number of cases where \(\lambda_{i} < \lambda_{N + i}\). An estimate of the p-quantile of TC,j, \(t_{C,j,p}\), is obtained by solving:

$$\sum\limits_{k = 1}^{{K_{p} }} {P_{j}^{*} \left( {C\left| {t_{C,j,p} } \right.} \right)} = p,\;K_{p} \in \left[ {1,N_{j,p} } \right],$$
(2)

where

$$P_{j}^{*} \left( {C\left| {t_{C,j,p} } \right.} \right) = P_{j} {{\left( {C\left| {t_{C,j,p} } \right.} \right)} \mathord{\left/ {\vphantom {{\left( {C\left| {t_{C,j,p} } \right.} \right)} {\sum\limits_{k = 1}^{{N_{j,p} }} {P_{j} \left( {C\left| {t_{C,j,p} } \right.} \right)} }}} \right. \kern-0pt} {\sum\limits_{k = 1}^{{N_{j,p} }} {P_{j} \left( {C\left| {t_{C,j,p} } \right.} \right)} }},$$
(3)

in which \(N_{j,p}\) is the total number of likely values \(t_{C,j,p}\) to be considered. \(t_{C,j,p}\) is considered as the point estimate of TC,j. The inter-quantile range, \(|t_{C,j,0.25} - t_{C,j,0.75} |\), is used as a measure of the uncertainty in TC,j. It is further assumed that \(T_{C,j} > T_{C,k}\) for Mmin,j < Mmin,k. For the numerical analysis to be carried out, N is set equal to 120, ϕ is taken equal to 1.5, the values of p equal to 0.25, 0.5, and 0.75 are considered (Feng et al. 2020), and M ranging from 3 to 7.25 and 3 to 8.0 is considered for completeness analysis of 2011 CCSC East Catalogue and 2011 CCSC West Catalogue, respectively. Based on these considerations, the obtained \(t_{C,j,p}\) for the employed historical seismic events for the CCSC East Catalogue and the CCSC West Catalogue are shown in Fig. 2a and b.

Fig. 2
figure 2

Curves corresponding to p = 0.25, 0.50, and 0.75 for selected events by considering a the CCSC 2011 East Catalogue and b the CCSC 2011 West Catalogue

Full size image

Many source zones are included in the source model for developing CanSHM6 (Kolaj et al. 2019). In each source zone, the earthquake is assigned to one of the five earthquake types: (1) stable crustal earthquake, (2) active shallow crustal earthquake, (3) subduction inslab earthquake with focal depth greater than 30 km (which is referred to as subduction inslab30), (4) subduction inslab earthquake with a focal depth greater than about 55 km (which is referred to as Subduction Inslab55), and (5) subduction interface earthquake; these earthquake types will be simply referred to as type 1 to 5, respectively. By considering these, each of the historical seismic events in the catalogues is assigned an earthquake type according to its geographical location (that is, latitude, longitude, and focal depth) and magnitude. The value of TC,j, \(t_{C,j,p}\), for an earthquake event with an assigned earthquake type of a magnitude greater than Mmin,j is shown in Fig. 2.

Since no interface earthquake could be identified from the historical catalogue based on the adopted criteria for the geographical location and magnitude, the spatial smoothing is carried out only based on the remaining four types of earthquakes. Note that the use of p = 0.5 was by Albarello et al. (2001). By adopting such a value, the earthquake events identified based on the completeness analysis (that is, earthquakes within the completed catalogue) are shown in Fig. 3. A comparison of the results shown in Figs. 1 and 3 indicates that the number of seismic events shown in Fig. 3 is much less than that shown in Fig. 1 and that the seismic activity varies geographically. The earthquakes within the completed catalogue for p = 0.5 is used in this study for developing the spatially smoothed source models.

Fig. 3
figure 3

Historical seismic events that are identified based on the earthquakes within the completed catalogue for p = 0.5

Full size image

2.2 Spatially Smoothed Seismic Source Models

In the present study, both the fixed kernel smoothing technique and adaptive kernel smoothing technique are considered to develop the spatially smoothed seismicity. The area within Canada is covered by a grid system with 0.5° × 0.5° (latitude × longitude) cells. For the fixed kernel smoothing technique, the smoothing procedure given by Woo (1996) is followed. Woo (1996) used the following kernel function \(K({\mathbf{M}}_{i,j} ,{\varvec{x}} - {\varvec{x}}_{i,j} )\) (Vere-Jones 1992):

$$K({\mathbf{M}}_{i,j} ,{\varvec{x}} - {\varvec{x}}_{i,j} ) = \frac{\alpha - 1}{{\pi [H({\mathbf{M}}_{i,j} )]^{2} }} \times \left( {1 + \frac{{\left| {{\varvec{x}} - {\varvec{x}}_{i,j} } \right|^{2} }}{{[H({\mathbf{M}}_{i,j} )]^{2} }}} \right)^{ - \alpha } ,$$
(4)

where \({\varvec{x}}\) is a point of interest (for example, the center of a cell), \({\mathbf{M}}_{i,j}\) is the earthquake magnitude of the i-th earthquake type and the j-th seismic event whose location is denoted as \({\varvec{x}}_{i,j}\), \(\left| {{\varvec{x}} - {\varvec{x}}_{i,j} } \right|^{2}\) represents the distance from \({\varvec{x}}\) to \({\varvec{x}}_{i,j}\), and α is taken as equal to 1.8 (Feng et al. 2020). In Eq. 4, H(Mi,j) is the bandwidth function for the i-th earthquake type that is given by:

$$H({\mathbf{M}}_{i,j} ) = a_{1,i} \exp (a_{2,i} {\mathbf{M}}_{i,j} ),$$
(5)

where a1,i and a2,i are the model parameters. These model parameters are estimated based on regression analysis by using one of the following three options and the earthquakes within the completed catalogue with M ≥ 3.0, which are grouped in bins according to the value of earthquake magnitude where the bin width equal to 0.25 is used in all cases. Option 1 follows the procedures suggested by Molina et al. (2001), where the earthquake events are grouped according to the number of bins of earthquake magnitude, and the pair of samples for each bin is formed by using the magnitude at the center of the bin and the shortest spatial distance between any pair of events within the bin. The samples are presented in Fig. 4 and are used for regression analysis based on Eq. 5. The fitted curve is also shown in Fig. 4. Rather than carrying out regression analysis based on Option 1, in Option 2, the regression is carried out by using the average of the shortest distance from a given event to any other events within the bin (Zuccolo et al. 2013). The samples as well as the fitted curve based on this option are also shown in Fig. 4. In Option 3, as suggested by Goda et al. (2013), the magnitude of each seismic event within a bin and its shortest distance to any event within the bin is used as a paired sample. The collection of the samples, as well as the fitted curve based on the parametric model shown in Eq. 5 for the considered seismic events within the completed catalogue, are shown in Fig. 4. In general, the fitted curves based on Option 3 are between the curves determined based on Options 1 and 2. Therefore, the fitted model based on Option 3 is considered in the following for the numerical analysis.

Fig. 4
figure 4

Samples and fitted curves are defined by Eq. 4 by considering three different options

Full size image

Before carrying out the spatial smoothing for the development of the seismic source models, we plot the historical seismic events within the completed catalogue with the magnitude M greater than Mmin = 4.8 as shown in Fig. 5. This Mmin = 4.8 is the same as that adopted for the development of CanSHM6 (Kolaj et al. 2020). Based on the adopted kernel function shown in Eq. 4 and the bandwidth shown in Eq. 5, the smoothed magnitude-recurrence relation based on the fixed kernel smoothing technique for the k-th cell with magnitude greater than M by considering the i-th earthquake type, denoted as λF,i,k(M) per year and per km2, is evaluated according to:

$$\lambda_{F,i,k} ({\mathbf{M}}) = \sum\limits_{j = 1}^{{N_{T} }} {I({\mathbf{M}}_{i,j} \ge {\mathbf{M}})I_{E} (E_{j} = i)K({\mathbf{M}}_{i,j} ,{\varvec{x}}_{k} - {\varvec{x}}_{i,j} )/\Delta t_{C,j,0.5} } ,$$
(6)

where \(I({\mathbf{M}}_{i,j} \ge {\mathbf{M}})\) is an index function that equals 1 for Mi,j ≥ M and zero otherwise, Ej denotes the earthquake type for the j-th event, which can be 1, 2, 3, or 4, \(I_{E} (E_{j} = i)\) is an index function that equals 1 for \(E_{j} = i\) and zero otherwise, NT is the total number of earthquakes in the completed catalogue considered, \({\varvec{x}}_{k}\) is the center of the k-th cell, and \(\Delta t_{C,j,0.5} = T_{R} - t_{C,j,0.5}\). The sum of λF,i,k(Mmin) for i = 1 to 4, denoted as λF,k(Mmin), represents the occurrence rate by considering all four earthquake types.

Fig. 5
figure 5

Historical earthquake events within the completed catalogue with M ≥ 4.8

Full size image

The adaptive kernel smoothing technique is also considered. For the adaptive kernel smoothing technique, the local bandwidth for the k-th cell with a magnitude greater than M is:

$$c_{i,k} ({\mathbf{M}}) = \left\{ {{{\left( {\prod\limits_{m = 1}^{{n_{c,i} }} {\lambda_{F,i,m} ({\mathbf{M}})} } \right)^{{1/n_{c,i} }} } \mathord{\left/ {\vphantom {{\left( {\prod\limits_{m = 1}^{{n_{c,i} }} {\lambda_{F,i,m} ({\mathbf{M}})} } \right)^{{1/n_{c,i} }} } {\lambda_{F,i,k} ({\mathbf{M}})}}} \right. \kern-0pt} {\lambda_{F,i,k} ({\mathbf{M}})}}} \right\}^{\gamma } ,$$
(7)

where nc,i is the total number of the cells that cover the area within Canada for i-th earthquake type and γ is a model parameter that is taken as equal to 0.5 (Silverman 1986). By using the estimated \(c_{i,k} ({\mathbf{M}})\), the magnitude-recurrence relation for the i-th earthquake type and the k-th cell, λA,i,k(M), is given:

$$\lambda_{A,i,k} ({\mathbf{M}}) = \sum\limits_{j = 1}^{{N_{T} }} {\frac{{I({\mathbf{M}}_{i,j} \ge {\mathbf{M}})I_{E} (E_{j} = i)}}{{\Delta t_{C,j,0.5} }}\frac{\alpha - 1}{{\pi [c_{i,k} ({\mathbf{M}}) \times H({\mathbf{M}}_{i,j} )]^{2} }} \times \left( {1 + \frac{{\left| {{\varvec{x}}_{k} - {\varvec{x}}_{{{(}i,j{)}}} } \right|^{2} }}{{[c_{i,k} ({\mathbf{M}}) \times H({\mathbf{M}}_{i,j} )]^{2} }}} \right)^{ - \alpha } .}$$
(8)

The sum of λA,i,k(Mmin) for i = 1 to 4, denoted as λA,k(Mmin), represents the occurrence rate by considering all four earthquake types that are obtained based on the adaptive kernel smoothing. Based on the above, the calculated λF,k(Mmin) and λA,k(Mmin) are shown in Fig. 6a and b, respectively (log represents the logarithm with base 10). Figure 6 shows that λA,k(Mmin) is more concentrated than λF,k(Mmin) for regions with concentrated historical events that occurred in eastern Canada, and a smoother seismicity is observed for western Canada. Both λF,k(Mmin) and λA,k(Mmin) follow the spatial distribution of the historical earthquake events, and λF,k(Mmin) for cells in the Cascadia subduction area is greater than λA,k(Mmin). The parameters used for spatial smoothing are given in Table 1. These parameters were determined based on the procedure given by Goda et al. (2013). In the same table, we also included the estimated total occurrence rate by considering each earthquake type for M ≥ 4.8. The obtained rate based on the fixed kernel smoothing, \(\lambda_{F,k}^{T} (4.8)\), is the same as that obtained by using the adaptive kernel smoothing, \(\lambda_{A,k}^{T} (4.8)\).

Fig. 6
figure 6

The plotted values are the rates per cell (0.5° × 0.5°) rather than per unit area: a λF,K (Mmin) and b λA,k (Mmin) based on Mmin = 4.8. The squares shown in the plots indicate the location of the identified sites to be considered

Full size image
Table 1 Bandwidth parameters used for smoothing and total occurrence rate by considering each earthquake type
Full size table

To appreciate the magnitude-recurrence relation, we plot the estimated λF,k(M) and λA,k(M) in Fig. 7 for 15 selected representative sites that are identified in Fig. 6 (squares). The results presented in Fig. 7 indicate that the trend observed from λF,k(M) agrees well with that from λA,k(M). Overall, λA,k(M) for the considered sites is greater than λF,k(M). For four sites, Vancouver, Kelowna, Whitehorse, and Destruction Bay, there are appreciable differences between λF,k(M) and λA,k(M). A further inspection of the results indicates that such a discrepancy is due to the use of a narrower bandwidth in the adaptive kernel smoothing as compared to that in the fixed kernel smoothing.

Fig. 7
figure 7

Estimated λF,k(M) and λA,k(M) per 0.5° × 0.5° cell based on Mmin = 4.8 for 15 representative sites. The square, cross, and triangle are used in the plots to indicate Mobs for earthquake types 1, 2, and 4. FKSM and AKSM denote fixed kernel smoothing technique and adaptive kernel smoothing technique, respectively

Full size image

The maximum magnitude, denoted as Mobs, in the estimated λF,k(M) and λA,k(M), as shown in Fig. 7 for each site, is directly controlled by the observed maximum magnitude for the earthquakes in the completed earthquake catalogue. They are equal to 7.40, 8.45, 6.10, and 7.00, for the earthquake type 1 to 4 (that is, eastern stable crustal earthquake, western active shallow crustal earthquake, subduction inslab earthquake with focal depth greater than 30 km, and subduction inslab earthquake with focal depth greater than about 55 km), respectively. These values can be smaller than the central values of the maximum magnitude adopted in developing CanSHM6, which are 7.8, 8.27, 7.00, and 7.60, for the earthquake type 1 to 4, respectively. Therefore, as a conservative measure, we consider that the maximum magnitude, Mmax, values for the four types of earthquakes are equal to the maximum of the mentioned two sets of values. In addition, if the maximum comes from Mobs, an additional 0.1 is included as a conservative measure. That is, Mmax equal to 7.80, 8.55, 7.00, and 7.60 for the four types of earthquakes are employed in the PSHA carried out in the following.

For cases where Mobs is less than Mmax for the i-th earthquake type, the magnitude-recurrence relation for M greater than Mobs is modeled using:

$$\lambda_{F,i,k} ({\mathbf{M}}) = \lambda_{F,i,k} ({\mathbf{M}}_{\text{obs}} )\frac{{\exp ( - \beta_{k} {\mathbf{M}}) - \exp ( - \beta_{k} {\mathbf{M}}_{{{\text{max}}}} )}}{{\exp ( - \beta_{k} {\mathbf{M}}_{\text{obs}} ) - \exp ( - \beta_{k} {\mathbf{M}}_{\max } )}},$$
(9)

in which the regression analysis is employed to estimate βk using the magnitude-recurrence relation for M up to Mobs that is determined for the k-th cell based on the fixed kernel smoothing technique. Similarly, λA,i,k(M) for M ≥ Mobs, is defined by using Eq. 9, except that the subscript F is replaced by A and βk is determined based on the magnitude-recurrence relation determined for the k-th cell based on the adaptive kernel smoothing technique.

Since the assessed seismic hazard for a site of interest x is mainly dominated by the earthquakes that occurred at a distance to this site smaller than 250 km (Goda and Hong 2009), it is instructive to estimate and compare the seismic occurrence rate in a circle centered at the site of interest with a radius equal to 250 km by considering different seismic source models. We evaluate such a rate, denoted as λ250km(M), by integrating the magnitude-recurrence rate within the circle for the site of interest by considering the source model adopted in developing CanSHM6 (which is referred to as CanSSM6), the source model developed by the fixed kernel smoothing technique (FKSM), and the source model developed by the adaptive kernel smoothing technique (AKSM). The obtained values of λ250km(M) for a few selected sites are shown in Fig. 8. In general, the rate shown in Fig. 8 obtained by considering FKSM is similar to that obtained by considering AKSM. However, a relatively large discrepancy is observed for Whitehorse, Kelowna, Calgary, and Edmonton. This is attributed again to the bandwidth used in these two smoothing techniques.

Fig. 8
figure 8

λ250km(M) for 15 representative sites based on fixed kernel smoothing technique (FKSM), adaptive kernel smoothing technique (AKSM), and the regional source model for CanSSM6. Western Canada model (WCM), Eastern Arctic model (EAM), and Southeastern Canada model (SCM) are used for CanSSM6. H2, HY, and R2 are the sub-models for the regional source model

Full size image

Note that CanSSM6 consists of three regional seismic source models that cover the regions in Canada (Kolaj et al. 2020). These three regional models are the Western Canada model, the Eastern Arctic model with two sub-models (H2 and R2 with weights equal to 0.6 and 0.4, respectively), and the Southeastern Canada model with three sub-models (H2, HY, and R2 with weights equal to 0.4, 0.4, and 0.2, respectively). For each sub-model, the nine magnitude-recurrence relations (with three sets of magnitude-occurrence rate models in combination with three values of Mmax), each with its associated probability, are considered. By folding the logical tree representing the source models, the overall magnitude-recurrence relation was developed and used as the input to develop CanSHM6 by using OpenQuake (Pagani et al. 2014). We use such an overall magnitude-recurrence relation to evaluate λ250km(M) for the considered sites of interest as shown in Fig. 8. The comparison shown in Fig. 8 indicates that the estimated λ250km(M) for sites located in western Canada by using FKSM and AKSM is relatively greater than that obtained based on the seismic source model adopted for the development of CanSHM6. This is especially the case for M greater than about 6.5. For sites located in eastern Canada, the estimated λ250km(M) by different methods are more comparable. Therefore, it is expected that the estimated hazard by using different seismic source models is likely to result in different seismic hazard estimates. The quantification of the differences in the estimated seismic hazard is presented in the following sections.

3 Seismic Hazard Assessment

First, we describe the ground motion models (GMMs) that were used for CanSHM6 and adopted in the present study. We then present the adopted simulation procedure used for seismic hazard assessment based on the delineated or spatially smoothed seismic source models. We estimate the seismic hazard as well as the uniform hazard spectrum for a few selected sites by using different seismic source models.

3.1 Ground Motion Models Used for the Development of CanSHM6

Before carrying out PSHA by using the developed FKSM and AKSM, we briefly summarize the GMMs adopted for developing CanSHM6 (Kolaj et al. 2019). For eastern stable shallow crustal earthquakes, two sets of GMMs are adopted, each with a probability of 0.5. The first set of GMMs was given by Atkinson and Adams (2013), hereafter AA13. This set was developed for eastern stable crustal earthquakes; it included three GMMs (or three GMM curves) to account for the epistemic uncertainty. These three were termed “upper,” “central,” and “lower” GMMs with associated weights equal to 0.25, 0.50, and 0.25, respectively. The central GMM in AA13 is constructed by taking the arithmetic average of the five candidate GMMs: one from Pezeshk et al. (2011), one from Atkinson and Boore (2011), one from Atkinson (2008) and Atkinson and Boore (2011), and two from Silva et al. (2002) for eastern North America. The upper and lower GMM curves are obtained by ± the standard deviation (in the logarithmic scale), each with a weight of 0.25, and the central GMM represents the median curve with a weight of 0.50. The second set, which is based on the Central and Eastern North America project, was given by Goulet et al. (2017), hereafter Gea17. The Gea17 is obtained by weighing a set of 13 candidate GMMs by using weights that depend on the natural vibration period Tn; the standard deviation for Gea17 is the same as that for AA13. Note that Kolaj et al. (2020) considered that the standard deviation for Gea17 is the same as that for AA13. Although it would be interesting to investigate the impact of the assumed standard deviation on the estimated seismic hazard, this is beyond the scope of the present study.

The four sets of GMMs given by the NGA-West2 project (Abrahamson et al. 2014; Boore et al. 2014; Campbell and Bozorgina 2014; Chiou and Youngs 2014), hereafter ASK14, BSSA14, CB14, and CY14, each with a weight of 0.25, are employed for western active shallow crustal earthquakes.

The GMMs given by Atkinson and Boore (2003), Abrahamson et al. (2016), García et al. (2005), and Zhao et al. (2006), hereafter AB03, Aea16, Gea05, and Zea06, are adopted for subduction inslab earthquakes; and the GMMs given by Abrahamson et al. (2016), Atkinson and Macias (2009), Ghofrani and Atkinson (2014), and Zhao et al. (2006), hereafter Aea15, AM09, GA14, and Zea06, are considered for subduction interface earthquakes. The weights assigned for each of the four GMMs for subduction inslab are 0.25. The same is done for the GMMs for subduction interface earthquakes.

The predicted median of the peak ground acceleration (PGA), spectral acceleration (SA) at the natural vibration period, Tn = 0.2, 0.5, and 2 s for VS30 = 450 m/s based on the adopted GMMs for CanSHM6 (Kolaj et al. 2020) is illustrated in Fig. 9. In the figure, δ is the fault dip, ϕ is the strike angle, λ is the rake angle, LR is the rupture length, WR is the rupture width, Ra = LR/WR, HSU is the upper seismogenic depth, HSL is the lower seismogenic depth, and Fh is the focal depth. The geometric representation of these parameters is depicted in Fig. 10.

Fig. 9
figure 9figure 9

Ground motion models (GMMs) used to estimate the seismic hazard by considering peak ground acceleration (PGA), and spectral acceleration (SA) at Tn = 0.2, 0.5, and 2 s for the four considered earthquake types

Full size image
Fig. 10
figure 10

Representation of the earthquake rupture plane

Full size image

For the predicted PGA and SA values presented in the figure for a given scenario event defined by M, δ, ϕ, λ, Ra, HSU, and HSL, the rupture area, Arup, is calculated using the following rupture scaling relationship:

$$\log A_{{{\text{rup}}}} = b_{1} + b_{2} \times {\mathbf{M}}.$$
(10)

Values of b1 = − 4.366 and b2 = 1 are considered for eastern stable crustal earthquake ruptures according to the US Nuclear Regulatory Commission (2012); the relation in such a case is referred to as CEUS2011. Values of b1 = − 3.42, − 3.99, and − 2.87; b2 = 0.9, 0.98, and 0.82 for the strike-slip, reverse, and normal faults, respectively, are considered according to Wells and Coppersmith (1994); the relationship in such a case is referred to as WC1994. Both CEUS2011 and WC1994 are used in Kolaj et al. (2020). Based on the calculated Arup, \(W_{R} = \sqrt {A_{{{\text{rup}}}} /R_{a} }\) and \(L_{R} = A_{{{\text{rup}}}} /W_{R}\). Note that the parameters selected for the plots shown in Fig. 9 are within the ranges considered in Kolaj et al. (2020).

Figure 9 shows that the relative variability for the predicted median ground motion is increased with the increased Tn. Also, an increased variability for the predicted PGA and SA by using different GMMs for a scenario event is observed as the distance to the rupture surface decreases and the magnitude increases.

3.2 Procedure Used for Probabilistic Seismic Hazard Analysis (PSHA)

The GMMs, which were described in the previous section and adopted for developing CanSHM6, are used in the following for the PSHA by considering FKSM and AKSM. The procedure for carrying out PSHA is based on the simulation steps given by Hong et al. (2006). For comparison purposes, we also considered CanSSM6. In this case, rather than using OpenQuake, we also use the same simulation procedure for consistency. The major steps for the simulation procedure include: (1) sample earthquake occurrence according to the adopted stochastic model for each earthquake type (for example, Poisson process); (2) generate the magnitude based on the considered earthquake type and magnitude-recurrence relation for each sampled seismic event; (3) sample rupture geometry parameters for each seismic event based on the corresponding probability distribution models as describe in CanSHM6 (see Tables 2 and 3); (4) evaluate the ground motion measures (that is, PGA and SA) at a site of interest by using the adopted GMMs; and (5) extract the annual maximum of PGA or SA to form their empirical probability distributions that can be used to estimate the T-year return period value of PGA, AT, and of SA, ST(Tn). In all cases, the analysis is carried out using the adopted GMMs for VS30 = 450 m/s, which was considered for CanSHM6.

Table 2 Distribution models for ϕ, δ, and λ (in degrees) for seismic events with different tectonic settings for spatially smoothed seismicity model [values used for CanSHM6 are from Kolaj et al. (2020)]
Full size table
Table 3 Distribution of Hhypo, [HSU, HSL] (value unit: km), and rupture scaling relationship adopted for the considered tectonic regions [values used for CanSHM6 are from Kolaj et al. (2020)]
Full size table

3.3 Comparison of the Estimated Site-Specific Seismic Hazard

Seismic hazard assessment is carried out by using the procedure, the source models (CanSSM6, FKSM, and AKSM), and the GMMs described and summarized in the previous sections. The consideration of CanSSM6 is for comparison purposes. For the analysis, one million years of seismic activities are simulated. Based on the analysis results, the empirical probability distribution of the annual maximum PGA in the upper tail region is shown in Fig. 11. The figure shows that the probability distribution of the annual maximum PGA in the upper tail region can be adequately modeled as a lognormally distributed random variable. The fact that the distributions in each plot do not have the same slope indicates that the magnitude of the uncertainty of the annual maximum PGA in the upper tail depends on the considered seismic source model. In several cases, the distributions are separate, suggesting that the estimated return period value of the annual maximum PGA for a site differs by considering different source models. Depending on the considered site, the large difference in the distribution could be between those obtained by using CanSSM6 and FKSM, by using CanSSM6 and AKSM, or by using FKSM and AKSM. In addition, a simple analysis indicates that the value of the coefficient of variation based on the fitted distribution to the upper tail region varies spatially. For example, it is about 2 for Vancouver and Victoria; ranges from about 5 to 9 for Ottawa, Montreal, and Quebec, and is greater than about 10 for Edmonton, Winnipeg, Toronto, Halifax, and St. John’s. The observed large spatial variation is consistent with that indicated in Hong et al. (2006). Although the spatially varying coefficient of variation could have a significant impact on the development of reliability-consistent or risk-target seismic design maps for Canadian sites (Hong and Hong 2006), the investigation of this topic is beyond the scope of the present study.

Fig. 11
figure 11

Empirical probability distributions of annual maximum peak ground acceleration (PGA) for the upper tail region for 23 locations

Full size image

To better appreciate the difference in the T-year return period of the annual maximum PGA, AT, we estimate A2475 (that is, quantile with an exceedance probability of 2% in 50 years) from the empirical distributions shown in Fig. 11. We show the estimated A2475 in Table 4. We include the values reported by Kolaj et al. (2020) (that is, from CanSHM6) in the table as well for comparison purposes. In addition, we also plot the empirical distribution for annual maximum SA at the vibration period Tn, denoted as S(Tn), for Tn = 0.2 and 2.0 s. Since the observed trends are similar to those shown in Fig. 11, they are not included (to save space). However, the estimated T-year return period values of S(Tn), ST(Tn), from the empirical distribution are shown in Table 4 for T = 2475 years, and Tn = 0.2 and 2.0 s. Similarly, for comparison purposes, we include the reported values of S2475(Tn) by Kolaj et al. (2020) in Table 4.

Table 4 Estimated peak ground acceleration (PGA) and spectral acceleration (SA) with Tn = 0.2 and 2 s, PE = 2% for selected sites (unit: g)
Full size table

To quantify the differences in the predicted AT and ST(Tn) shown in Table 4, we define the relative difference \(r_{A/B} = (y_{A} - y_{B} )/y_{B}\), where A and B represent CanSSM6, FKSM, AKSM, or CanSHM6, respectively, and yA represents AT or ST(Tn) by using the source model A, and yB is defined similarly. By considering the PGA, \(r_{{\text{CanSHM6/CanSSM6}}}\) for the sites with AT greater than 0.1 g ranges from − 7.9 to 7.8% with a mean of practically zero. These values become − 4.6%, 7.6%, and 1.8% by considering ST(0.2), and − 1.4%, 10.7%, and 2.7% by considering ST(2.0) if ST(Tn) greater than 0.1 g is considered. This comparison indicates that the results obtained by our approach and OpenQuake agree well.

By considering again AT or ST(Tn) (for FKSM) greater than 0.1 g, the obtained minimum, maximum and mean of \(r_{CanSHM6/FKSM}\) equal − 114.3%, 80.4%, and 22.6% for PGA, − 100.1%, 79.7%, and 19.6% for ST(0.2), and − 90.3%, 80.4%, and − 4.5% for ST(2.0). Similarly, the obtained minimum, maximum, and mean of \(r_{{\text{CanSHM6/AKSM}}}\) equal − 56.0%, 78.4%, and 18.5% for PGA, − 222.5%, 71.9%, and 7.0% for ST(0.2), and − 34.7%, 65.9%, and − 0.2% for ST(2.0) if AT or ST(Tn) (for AKSM) greater than 0.1 g is considered. These calculated relative differences reflect the observed difference in the source models.

By considering again AT or ST(Tn) (for AKSM) greater than 0.1 g, the calculated \(r_{{\text{FKSM/AKSM}}}\) equal − 107.4%, 30.4%, and − 8.8% for PGA, − 280.6%, 27.9%, and − 20.8% for ST(0.2), and − 8.1%, 29.2%, and 3.2% for ST(2.0). These results indicate that the obtained seismic hazard by using FKSM and AKSM, on average, are consistent although in some cases there are significant differences. Again, these observed differences reflect the difference in the obtained spatially smoothed source models.

Rather than showing the estimated AT or ST(Tn), we plot the estimated uniform hazard spectra (UHS) in Fig. 12 for the 23 selected sites by considering T = 2475 years. In all cases, the trends of the UHS are consistent for considered source models and CanSHM6. Since the logarithmic scale is used for the vertical axis of the plot, the differences in the estimated UHS by using different source models appear to be relatively small. The UHS from CanSHM6 and the estimated UHS by considering CanSSM6 in the present study agree very well. Also, the estimated UHS by considering FKSM and AKSM are consistent.

Fig. 12
figure 12

Uniform hazard spectra (UHS) based on CanSHM6, and those estimated in the present study by considering CanSSM6, the fixed kernel smoothing technique (FKSM), and the adaptive kernel smoothing technique (AKSM)

Full size image

4 Seismic Hazard Maps Obtained by Considering Different Seismic Source Models

To map the seismic hazard, we estimate AT or ST(Tn) by following the procedure used in the previous section but for each cell of 0.5° × 0.5° in a grid system covering the regions in Canada. In all cases, we consider VS30 = 450 m/s, Tn = 0.2, 0.5, 1, 2, 5, and 10 s, and each of the three source models (CanSSM6, FKSM, AKSM). The obtained A2475 (that is, S2475(0)) and S2475(0.2) are shown in Fig. 13. Note that although we have plotted and inspected the maps for other Tn values, they are not included here (to save space).

Fig. 13
figure 13

Seismic hazard map in terms of A2475 and S2475(0.2): CanSSM6, fixed kernel smoothing technique (FKSM), and adaptive kernel smoothing technique (AKSM) are used as the source model for the maps shown in the first to third rows, respectively

Full size image

An inspection of the results shown in Fig. 13 indicates that, in general, the spatial trends of the seismic hazard maps obtained by considering different source models are similar. This is expected as the comparison shown in the previous sections indicates that, on average, the estimated AT or ST(Tn) by considering different source models are in good agreement. To quantify the differences between different maps, we calculate the average of the ratio of ST(Tn) obtained by considering FKSM to that obtained by considering CanSSM6. The obtained values for Tn = 0, 0.2, 0.5, 1, 2, 5, and 10 s are 0.98, 0.95, 0.97, 0.98, 1.03, 1.15, and 1.21, respectively. This indicates that, on average, the use of CanSSM6 could overestimate or underestimate the seismic hazard as compared to that obtained by using FKSM. We also calculate the average of the ratio of ST(Tn) obtained by considering AKSM to that obtained by considering CanSSM6. The obtained values for Tn = 0, 0.2, 0.5, 1, 2, 5, and 10 s are 0.95, 0.95, 0.93, 0.93, 0.95, 1.02, and 1.08, respectively. In this case, the use of CanSSM6 could overestimate the seismic hazard as compared to that obtained by using AKSM, depending on the considered vibration period. Given the absence of compelling justifications to dismiss the spatially smoothed source models, one may consider multiple seismic source models as part of the logic tree for the seismic hazard assessment. In such a case, one may assign a weight of 0.5 for CanSSM6, 0.25 for FKSM, and 0.25 for AKSM to develop the combined source model. By adopting such a combined source model and repeating the seismic hazard mapping analysis, the obtained maps in terms of A2475 and S2475(0.2) are also shown in Fig. 14. The calculated average of the ratio of ST(Tn) obtained by considering the combined source model to that obtained by considering CanSSM6 is 0.98, 0.97, 0.97, 0.98, 0.99, 1.04, and 1.07, for Tn = 0, 0.2, 0.5, 1, 2, 5, and 10 s, respectively. The average of the ratio indicates that the estimated seismic hazard by considering the combined source model becomes relatively consistent. But this consistency, which is based on average behavior, should not be extended to an individual site.

Fig. 14
figure 14

Seismic hazard map in terms of A2475 and S2475(0.2) by considering the combined source model, that is, the combination of CanSSM6, fixed kernel smoothing technique (FKSM), and adaptive kernel smoothing technique (AKSM) with weights equal to 0.5, 0.25, and 0.25, respectively

Full size image

It must be emphasized that in recent years, the assessment of the seismic hazard due to the induced seismic events caused by hydraulic fracturing was also carried out (Atkinson et al. 2015; Atkinson et al. 2016). The results indicate that the magnitude-recurrence relation for the induced seismicity and their ground motion attenuation characteristics are different from the natural seismic events. The studies further showed that the seismic hazard due to the induced seismic events in Alberta, Canada, could be greater than that due to natural hazards. This can be explained by noting that the seismic hazard due to natural seismic events in Alberta is very small. Such seismic hazard is not included in Figs. 11, 12, 13, and 14.

5 Conclusion

The current Canadian seismic hazard maps are developed by considering a delineated source model (that is, CanSSM6). These maps provide information about the T-year return period values of annual maximum spectral acceleration (SA) at the natural vibration period Tn, referred to as ST(Tn), as well as the annual maximum peak ground acceleration (PGA), denoted as AT, represented as ST(Tn).

In the present study, we developed two new source models, namely FKSM and AKSM, using the historical earthquake catalogue and employing the fixed kernel smoothing technique and adaptive kernel smoothing technique, respectively. For the development, numerical analyses were conducted to assess catalogue completeness and determine the appropriate bandwidth for the smoothing process. Using the newly developed source models, we performed site-specific seismic hazard assessments and created seismic hazard maps. These assessments and mappings took into consideration the different source models, namely CanSSM6, FKSM, and AKSM. The analysis results yielded the following findings:

  1. (1)

    The magnitude-recurrence rate for a specific site or area of interest varies depending on the chosen source model, although in general, they exhibit similarities. This discrepancy suggests that the magnitude-recurrence obtained from CanSSM6 may not align entirely with that derived from the historical catalogue for a particular site or region of interest. Consequently, this difference is likely to be reflected in the assessed seismic hazard in terms of PGA and SA.

  2. (2)

    The upper tail regions of probability distributions for PGA or SA can often be effectively modeled as lognormally distributed random variables. This observed behavior remains consistent irrespective of the chosen source model.

  3. (3)

    In most cases, if the assessed intensities of AT and ST(Tn) for a site using either source model are large, the obtained seismic hazards by using different source models tend to agree well. However, it is observed that the assessed AT and ST(Tn) for a site may be significantly large when using one model, while they can be comparatively small when another source model is utilized.

  4. (4)

    The application of FKSM generally yields smoother seismic hazard contours compared to those obtained through the use of AKSM.

  5. (5)

    On average, the use of CanSSM6 may underestimate the seismic hazard in comparison to the results obtained by using FKSM. However, the use of CanSSM6 can either underestimate or overestimate the seismic hazard relative to AKSM, depending on the specific vibration period considered.

  6. (6)

    Given the absence of compelling justifications to dismiss the spatially smoothed source models, it is suggested that multiple seismic source models can be considered as part of the logic tree for seismic hazard assessment. An example of such a combined source model, which assigns weights of 0.5 to CanSSM6, 0.25 to FKSM, and 0.25 to AKSM, is employed to map the seismic hazard. On average, the results obtained by using CanSSM6 and the combined source models become relatively consistent. It must be emphasized that this consistency, which is based on average behavior, should not be extended to an individual site.