Probabilistic seismic hazard analysis in low-seismicity regions: an investigation of sensitivity with a focus on Finland

We investigate probabilistic seismic hazard analysis (PSHA) in low-seismicity regions in which epistemic uncertainties are largely due to the sparsity of data, with a focus on Finland, northern Europe. We investigate the sensitivity of site-specific PSHA outcomes to different choices of basic input parameters, starting from preexisting PSHA models of the nuclear licensees in the country, without producing a final hazard curve. The outcome shows that the parameters and models needed to estimate future seismicity rates from actual observations, in particular the b value, seismicity rates, and the largest possible magnitude, Mmax, as well as the median ground-motion prediction equation, play significant roles. The sensitivity also depends on the spectral frequency; for example, the effect of Mmax is significant especially for a low-frequency hazard at annual frequency of exceedance 10−5 but more moderate for peak ground acceleration. The delineation of seismic source zones (SSZs) remains ambiguous in regions of low seismicity. This, combined with the dominance of the host SSZ and its seismicity parameters, may have a substantial impact on the outcome. Our results are quantitatively applicable to Finland, but may also be of relevance to other low-seismicity regions in Europe and elsewhere. For future work we recommend the exploration of PSHA sensitivity with focus on the host SSZ with its immediate vicinity and the b value around the site of interest.


Introduction
The safety of critical infrastructure, such as nuclear power plants (NPPs), is sustained by several defense systems, including protection against major natural phenomena that pose an external hazard to the NPP. Resilience to strong vibratory ground motions triggered by natural earthquakes is evaluated by seismic risk analysis, which is composed of three main parts: assessment of seismic hazard, exposure, and vulnerability (e.g., Scawthorn 2008).
We deal with probabilistic seismic hazard analysis (PSHA), a deductive methodology used to estimate the probability that a threshold value of a selected ground-motion parameter will be exceeded at a target site or region during a given period, as outlined by Cornell (1968). It is the de facto standard procedure for seismic hazard analysis worldwide (e.g., Budnitz et al. 1997aBudnitz et al. , 1997bSolomos et al. 2008;Hanks et al. 2009;USNRC 2012). Sitespecific PSHA, relevant to NPP safety, aims to quantify uncertainties about the location, magnitude, and resulting ground motion of future earthquakes and combine them into a ground-motion intensity distribution at the site of interest (e.g., Reiter 1990). The primary outcomes of site-specific PSHAs are seismic hazard curves and uniform hazard spectra. The vulnerabilities of NPPs are analyzed considering the estimated seismic hazard levels in order to assess risks and to mitigate damage from future earthquakes.
Since its development in the early 1960s, following largely from the collaboration between Allin Cornell and Luis Esteva and others (McGuire 2008), PSHA has become the tool to guide engineering design and decision-making. Testing the sensitivity of PSHA results to the different models and input parameters remains an inherent part of elucidating the methodology (Grünthal and Wahlström 2001;Beauval and Scotti 2004;Petersen et al. 2004;Sabetta et al. 2005;Avital et al. 2018, among others). PSHA has also been criticized for failing to map high-hazard regions in which subsequent devastating earthquakes have occurred (e.g., Stein et al. 2012) and for ignoring earthquake physics (Mulargia et al. 2017). Moreover, conducting PSHA in essence means taking seismicity models at tectonically active plate boundaries and stable continental regions (SCRs) be similar. However, the hazard associated with identified SCR faults may be overestimated and underrated elsewhere (England and Jackson 2011). It has been proposed that large SCR earthquakes may be non-recurrent phenomena, so the concepts of recurrence or fault slip rate would not apply (Calais et al. 2016).
The quantitative assessment of seismic hazard is challenging in low-seismicity regions such as Finland, northern Europe, where the level of natural seismicity is low, with no earthquake known to have exceeded moment magnitude (M w )4.5 in more than three centuries (Ahjos and Uski 1992 and subsequent;Fig. 1). Small earthquakes dominate the available earthquake catalog, and the b value is used to extend the observed magnitude ranges to non-observed ranges. Identification and removal of non-seismic events and assessment of completeness levels becomes critical. Much of the observed seismicity comes from diffuse sources, meaning that it is difficult to discern an obvious correlation between the pattern of seismicity and geological structures (IAEA 2016(IAEA , 2022. The complete absence of strongmotion recordings from natural earthquakes complicates construction of local GMPEs; typically GMPEs are borrowed from tectonically similar regions. The logic tree must be populated by seismic source zone (SSZ) and GMPE models that capture epistemic uncertainty, but construction of such models on the basis of sparse data may remain insufficient (e.g., Bommer 2012). Differences in PSHA outcomes can also follow from procedural differences in how the experts interact with the analysts and with each other, in how they share or do not share data, or as a result of a lack of feedback and proper communication between experts and analysts (Budnitz et al. 1997a(Budnitz et al. , 1997b. Finland has an active nuclear program, which covers all stages of NPP development including planning, commissioning, operation, and decommissioning. Site-specific seismic hazard analysis has been mandatory for NPPs since 1988. For older NPPs with narrow seismic margins, changes in the calculated seismic hazard may have large effects on the respective risk. PSHA has been carried out retrospectively for these NPPs (e.g., Varpasuo Fig. 1 Seismicity of the target region with three seismic source zone (SSZ) models. Blue dots indicate earthquakes that occurred between 2000 and 2020 and for which reliable focal depths are available. Gray dots denote earthquakes that occurred between 1970 and 1999 (years inclusive). The black squares denote the sites of interest for probabilistic seismic hazard analysis (H = Hanhikivi; L = Loviisa; O = Olkiluoto). The retained SSZs were A and B for the L site and B and C for the O site (green). The SSZs D, E, F (blue) and I, J, K (red) were retained from two alternative SSZ zonings for the H site. The earthquake data were retrieved from Ahjos and Uski 1992 and subsequent. The geologic map is from Nironen et al. (2016) © Geological Survey of Finland. Abbreviations are as follows: FI(Finland), SE(Sweden), RU(Russia), Ku(Kuusamo), GoB(Gulf of Bothnia), LBB(Ladoga-Bothnian Bay zone). Pink color inside zone 'A' is the onshore portion of the Vyborg rapakivi granite batholith inside Finland. LBB divides the Karelian province in the northeast from the Svecofennian province in the southwest and Puttonen 1985), followed by risk evaluation (Varpasuo et al. 1996). The nuclear licensees or license applicants are responsible for seismic hazard analysis and determination of the design basis earthquake for their facility. The analyses are submitted to the Radiation and Nuclear Safety Authority (STUK) for review and acceptance every 10 years (YVL 2019). Since seismic design is not required for residential and commercial buildings in the country, the seismic hazard and risk communities focus on the nuclear industry. National traits such as a taciturn communication culture combined with a desire to avoid conflict have been identified as potentially problematic in the decision-making processes related to nuclear safety (OECD 2019a). The comparative investigation by the Organisation for Economic Co-operation and Development (OECD) Nuclear Energy Agency (OECD 2019b) focuses on PSHAs of many European countries and SCR regions elsewhere, such as Sweden, the Czech Republic, Korea, and eastern Canada.
Motivated by complexities of decision-making based on a variety of hazard results, STUK launched a new project to obtain additional insight into the PSHA process for the review of future PSHAs in Finland. The new calculations were aimed to examine the sensitivity of PSHA results to different choices of input, using seismicity data available in Finland. The project did not aim to produce complete PSHA outcomes to be imposed on regulations or future work as such. However, the sensitivity calculations help to explore which annual frequency of exceedance (AFE) should be adopted and from which hazard curve the ground-motion value should be read in the future. The target PSHA result in Finland is traditionally the AFE 10 −5 median hazard curve, with peak ground acceleration (PGA) of 0.1 g as the minimum. In the new project, calculations were also carried out for AFE 10 −7 , which is the basis for the design-extension condition (DEC) loads, i.e., extremely rare scenarios, in which the design goals for the nuclear plant are basic physical integrity.
In the following, we outline the seismotectonics of the target area, describe the execution of the SENSEI (SENsitivity study of SEIsmic hazard prediction in Finland) project, present selected PSHA results with relevance to Finland, adjacent areas, and possibly also to other low-seismicity regions, and discuss the results and the scope of this approach for seismic hazard analysis.

Seismotectonics and seismic source zones
The target region, Finland, is situated within the SCR of the Fennoscandian Shield, where the surface area covers some of Earth's most ancient continental crust from Precambrian time. The Fennoscandian Shield formed by the progressive accretion of large crustal units of different ages and exhibits a geochronological zonation in which the oldest ages are found in the northeast (Gaál and Gorbatschev 1987). In the Finnish territory, the crust consolidated during the Archean and the Proterozoic (Fig. 1). The Karelian province to the northeast formed largely during the Lopian orogenies (2.9-2.6 Ga) and the Svecofennian province to the southwest formed during the Svecofennian orogeny (1.9-1.75 Ga). The northwest-southeast-trending Ladoga-Bothnia Bay (LBB) zone is a major tectonic structure that divides the two provinces. Vaasjoki and Sakko (1988) explained the evolution of this zone by a progressive accretion of crust from the northwest to the southeast.
The Fennoscandian Shield features an overall low level of natural seismicity with areas of enhanced microearthquake activity and seismic quiescence. Enhanced seismicity occurs in particular along the eastern coast of Sweden, in the Gulf of Bothnia bordered by Finland to the east and Sweden to the west, in Lapland and Kuusamo (Fig. 1). The marked northeastern-southwestern seismicity trends have been associated with postglacial faults and the western flank of the Gulf of Bothnia (e.g., Olesen et al. 2013;Lindblom et al. 2015). In the Finnish territory, the most active northeast-southwest trending zone extends from Kuusamo over the border to Russia. Reliable focal depth estimates are available since 2000 and show a wide range of depths from very shallow down to 40 km and even deeper (e.g., Uski et al. 2012). On average, seismicity is shallower in the south than in the north of the country.
NPP-related PSHAs are based on seismicity generated by the present stress regime, so the effects of glaciations are not considered. The first entry of the earthquake catalog available dates back to 1375 CE, but data are relatively sparse until the systematic collection of observations using macroseismology since the late 1800s and the instrumental monitoring of local seismic events since the mid-1950s (Ahjos and Uski 1992 and subsequent). Figure 1 shows the SSZs that host the NPPs and their nearest SSZ. The SSZs have primarily been delineated using spatially distinct seismicity patterns. The spatial variation of seismicity is expected to reflect local changes, i.e., seismic stress released in the current stress field, which favors preexisting zones of weakness, ductile shear zones, and other weak zones (Korja and Kosonen 2015). The major LBB structure does not stand out as distinct to seismicity. SSZ 'A' is the geologically uniform Vyborg rapakivi granite batholith that is discussed later.

The SENSEI project: objective and workflow
The SENSEI project was conducted in 2019-2020 under the auspices by STUK, with reporting extended to 2021. The unprecedented feature was that SENSEI did not aim to provide new, complete PSHA outcomes to be imposed by the regulator. With the help of new sensitivity calculations, STUK aimed to obtain further insight into the PSHA for future PSHA reviews. The focus was on exploring the sensitivity of PSHA in the country. Also, SENSEI gave an opportunity to compare PSHA for the sites of interest to the nuclear industry in Finland: Loviisa (L) hosting the oldest NPP from the 1970s, Olkiluoto (O) with two older units and a new one, as well as a repository for spent nuclear fuel, and an NPP considered at Hanhikivi (H). The background seismicity differs between the northern site and the two southern sites (Fig. 1).
The project combined an expert group (EG) and a calculation group (CG) with expertise on regional seismology, geology and geophysics, earthquake catalogs, existing PSHAs, earthquake engineering, and seismic qualification of NPPs. The work consisted of a structured sequence of discussions between the EG, the CG, and STUK, based on the site-specific PSHAs conducted in the previous years. An initial list of issues to be explored was compiled from the questions raised by the EG at the kickoff meeting on 5-6 March 2019, based on reviews of the preexisting PSHA documentation. The growing list was modified, prioritized, and converted to work tasks at each meeting. A detailed description of the project execution is provided by Mäntyniemi et al. (2021).
The starting points of the sensitivity investigation were the PSHA models developed previously for the three nuclear sites. Their results were given as median PGA hazard curves and the associated response spectra at AFE 10 −5 in accordance with the national guidelines (YVL 2019). The starting point for the new calculations was obtained from the preexisting results by a so-called pruning process, in which the logic-tree branches with the largest influence on the outcome were retained. The pruning process preserved the core of the models without introducing significant changes to the results. Table 1 outlines how the preexisting logic trees with 144 branches were reduced to four branches. Regarding SSZs, for example, the main hazard contributor was found to be the SSZ hosting the facility, with the nearest SSZs acting as secondary contributors (cf. Figure 1). The dominance of the host SSZ and its vicinity can also be noticed in the PSHA conducted in the SCR of South Africa (Bommer et al. 2013 Figs. 2, 3, 5;Bommer et al. 2015 Fig. 18). When the more distant SSZs were removed, their joint contribution was lost, and when a logic-tree branch was removed, its associated weight was redistributed among the remaining branches.
Only limited review and reevaluation of the Gutenberg-Richter (GR) parameters a GR and b GR was possible, so values from previous PSHAs were used (Table 1).
The sites of interest have very hard (base) rock (VHR) conditions, with shear-wave velocity V s > 3000 m/s, and the focus is on the high-frequency content of the ground motions. The spectral frequencies used throughout the work were selected as 1, 5, and 25 Hz in addition to PGA (100 Hz). The selection was partly based on the old spectra of national guidelines.
The large number of ground-motion prediction equations (GMPEs) available worldwide (Douglas 2021) declines rapidly when selecting equations suited for SCRs. Two GMPEs were selected based on mutual consent between the EG, CG, and STUK: Next Generation Attenuation for Central and Eastern North America (NGA-East; Goulet et al. 2018, The NGA-East attenuation model is a suite of median predictions coupled with a set of uncertainty (randomness, σ) estimates. It has been developed for hard rock (V s > 3000 m/s) in central and eastern North America, which is tectonically and geologically an acceptable analog to Finland and adjacent areas. NGA-East models were used with the weighted average mean prediction in the pruned logic-tree models. The weighting of the median predictions in the logic tree followed the original proposal by Goulet et al. (2018). The weights of the 17 branches were spectral-frequency dependent, based on the weighting factors in their Tables 9.2a and b. The σ model used with NGA-East was the composite ergodic sigma (σ) model for Central-North America (CENA). The model coefficients (σ1, σ2, σ3 and σ4) for the central, high and low prediction of σ were taken from Goulet et al. (2018;Table 11-24). The coefficients were used to calculate the magnitude and spectral-frequency dependent σ of the model, using Eq. 11-14 from Goulet et al. (2018). Fenno-G16 is an adaptation of the G16 GMPE by Graizer (2016). It is partly based on data recorded at seismic stations in Finland, valid for magnitudes M w ≥ 2, and thus suited to be a complementary GMPE. A significant number of exploratory sensitivity calculations were run using the four spectral frequencies, starting from the pruned logic trees. The effects of the minimum magnitude (M min ), GR parameters a GR and b GR , depth distribution, GMPE, and maximum magnitude were sequentially quantified (Table S1 in the Supplementary information shows the scheme for site L). In the final models of the sensitivity investigation, Fig. 3 The final logic tree for the three sites of interest with branching for seismic source zoning, maximum magnitude, the Gutenberg-Richter parameters a GR and b GR (b + with standard deviation added and b − subtracted), the ground-motion prediction equation (GMPE) median prediction, and the GMPE randomness (σ) estimate. The GMPE nodes refer to the 17 mean predictions of the NGA-East GMPE, and σ refers to its three ergodic σ models (Goulet et al. 2018Youngs et al. 2021). The seismic source zoning models H_3 and L_2 are modifications made during the SENSEI project. The weights of the logic-tree branches are shown except for the GMPEs which depend on spectral frequency and were taken from Table 9-2 of Goulet et al. (2018) a new logic-tree structure, based on discussions during the project, was used for each site. The NGA-East GMPE was used with all its 17 median predictions and three aleatory variability estimates described by the ergodic sigma model. With ergodicity it is assumed that the ground-motion variability at the sites is equal to the ground-motion variability present in a global dataset (Goulet et al. 2018;Al Atik et al. 2010). Table S2 in the Supplementary information summarizes the assumptions and input values of the calculations.

Results
The results from the pruned logic tree and from the final logic tree are investigated separately.

Calculations using the pruned logic tree
In the single-branch calculations, the effect of change of a single input parameter was tested. They show the change on the outcome separately. Comparison was with the baseline models of the pre-SENSEI projects at AFE 10 −5 . In the case of site L, the shape of the baseline spectrum deviates from that of the new calculations, especially at 5 Hz and 25 Hz (Fig. 2a). The lowest acceleration values are related to the NGA-East weighted average GMPE, while the Fenno-G16 GMPE gives higher values, which was also noted during its development (Fülöp et al. 2020). An increase in the largest possible magnitude (M max ) from 6.5 to 7.5 had less effect on the output than an increase in M min from 2 to 4 when using the Fenno-G16 GMPE; however, they have a contrasting effect, because increasing M min decreases the hazard, also shown, for example, by Bender and Campbell (1989). A relatively low M min is justifiable due to the vulnerabilities of the older plants in the country, which were not designed to withstand earthquake ground motions.
The former delineation of the host SSZ of site L covered the geologically uniform Vyborg rapakivi granite batholith, which extends to offshore areas and over the border to Russia (SSZ 'A' in Fig. 1) and implied that future seismicity can occur anywhere with equal probability. The uneven spatial distribution of observed seismicity prompted to new delineations and subsequent sensitivity tests. Modifying the host SSZ had an anticipated large impact on the output, because the observed seismicity was distributed over approximately one-third of the former SSZ.
For site O, the shape of the baseline spectrum similarly deviates from the new calculations at 5 and 25 Hz, and the lowest acceleration values are associated with the NGA-East GMPE (Fig. 2b). A large effect followed from incorporating new data into the distribution of focal depths, according to which earthquakes occur in the uppermost 13 km of the crust in the southern part of the country, rather than in the entire crust. In contrast, the focal depths determined for the vicinity of site H imply an overall seismogenic crust, and the new depth distribution 0-26 km had a lesser impact (Fig. 2c). The data do not suffice for tailoring site-specific depth distributions. In all the diagrams in Fig. 2, the lowest output values shown are for the NGA-East weighted average GMPE with M max 5.5. When M max was increased to 6.5, the results increased substantially.

Final sensitivity analysis
The final logic trees used had 4590, 2295, and 6885 branches for sites L, O, and H, respectively (Fig. 3). Each logic tree had five branching levels, and the different numbers of branches follow solely from the number of alternative SSZ zonings available for each site. The M max distribution was updated from the median of M w 5.5 used previously to a median of M w 6.64, based on the Bayesian method (Johnston et al. 1994). The highest M max value used was M w 7.77, which is close to the proposal by Vanneste et al. (2016).
The hazard distribution among the logic-tree branches was assessed qualitatively for PGA and 1 Hz at AFE 10 −5 , i.e., at the highest and lowest frequency investigated (Fig. 4). The high-frequency plots (Fig. 4a-c) are more similar to each other than to the low-frequency plots (Fig. 4d-f). The branches generating a higher hazard concentrate at certain nodes that vary with the spectral frequency. A higher value of hazard is generated for low values of b GR (e.g., Fig. 4a, d). The hazard shows a significant variation with alternative SSZ zonings for sites L and H and for given branches of the NGA-East GMPE (e.g., 16 and 17). M max plays a less significant role in controlling the hazard for PGA but features more prominently at low frequency. Figure 5 shows the shape of the spectrum for site L, peaking at 20-25 Hz. The spectral-peak to-PGA ratio obtained was in the range of 2.30-2.33 and 2.05-2.17 for the mean uniform hazard spectra at 10 −5 and 10 −7 AFE, respectively. The sites have VHR conditions, and focus is on the high-frequency content of the ground motions. The magnitude-distance disaggregation of the results shows the expected trend of significant contribution to the hazard originating from the near-distance range at AFEs 10 −5 at the three sites (Fig. 6). However, also longer distances contribute to the hazard at 1 Hz and AFE 10 −5 . In particular, the contribution of intermediate-magnitude earthquakes at short distance is significant. Site H is close to the seismically active eastern coast of Sweden (Fig. 1), which contribute to the hazard. Figure 7 shows the PGA hazard curves for the entire AFE range with the mean of all logic-tree branches and the mean of the subsets branched at each branching level for M max (Fig. 7a-c), GR parameters ( Fig. 7d-f), and GMPE median predictions (Fig. 7g-i). Results corresponding to M max show the smallest dispersion, and the different NGA-East GMPEs a larger dispersion around the global mean, especially at low AFEs. On the other hand, the GR parameters have a general effect on the entire range of AFEs and primarily depend on the target location. A larger dispersion of the GR parameters is seen at sites L and O due to the very low number of earthquakes in their vicinity. In the case of L, low-magnitude earthquake sequences further complicate the determination of seismicity rates.
In the nuclear field, with focus on low AFEs, the key PSHA models and parameters may be different from those related to residential and commercial buildings that is AFE 10 −2 to 10 −3 (EN 1998(EN -1:2004. A key outcome is the significant influence of the GMPE for low AFEs. This reinforces the idea of using GMPEs with a quantified estimate of median prediction uncertainty and the disadvantage when relying on a small number of GMPEs (e.g., Fig. 7g-i;Sabetta et al. 2005). The effect of σ is also significant at low AFEs when reasonable σ values (e.g., > 0.6) are used but not as much as that of GMPE median predictions. Generally, uncertainties tend to increase the spread of the hazard for low AFEs.
The PGA and 1 Hz hazard values, which control the design basis earthquake in the Finnish guidelines (YVL 2019), are plotted in Fig. 8. The figure has been derived from 1 3 Fig. 4 Influence of logic-tree nodes on a-c peak ground acceleration (PGA) and d-f 1 Hz frequency output for 10 −5 annual frequency of exceedance at sites a, d Loviisa (L), b, e Olkiluoto (O), and c, f Hanhikivi (H). Line tone represents the hazard on the branch normalized to the maximum value of all branches (i.e., scale is from 0 white to SA max /PGA max black. PGA min , PGA max , and the mean PGA values are given)

Fig. 5
Uniform hazard spectrum in the SCR very hard rock setting of Finland: Example from site L (Loviisa) for frequencies considered by the Next Generation Attenuation East ground-motion prediction equation (Goulet et al. 2018Youngs et al. 2021). The blue lines denote selected frequencies 1 Hz, 5 Hz, 25 Hz, and 100 Hz (PGA). The continuous lines denote the mean hazard and the dashed lines the median hazard. Spectral acceleration is given in fractions of g. The annual frequency of exceedance is 10 −5 . The corresponding figures for Olkiluoto and Hanhikivi are given in Figure S1 in the Supplementary information Fig. 6 Disaggregation results by magnitude and hypocenter distance for the final logic tree. Example for peak ground acceleration (PGA) at Loviisa (a), Olkiluoto (b), and Hanhikivi (c) for annual frequency of exceedance (AFE) 10 −5 . Distances are from 0 km in 20-km increments, and magnitudes are from M w 4 in unit increments; blue bars are for the range of M w 4 to M w 5, orange for M w 5 to M w 6, gray for M w 6 to M w 7 and yellow for M w > 7. The z-axis shows the contribution to the total hazard of unit-width magnitude bins. Disaggregation plots for PGA for AFE 10 −7 and 1 Hz spectral acceleration for AFEs 10 −5 and 10 −7 are given in Figure S2 in the Supplementary information the hazard curves of Fig. 7 for 10 −5 AFE. The values at each branching level represent the spread of hazard observed at that level. The most important parameters influencing the hazard are the GMPE median prediction and the GR parameters. M max plays a prominent role for the low frequency. SSZ models also have a visible effect, especially for site L.

Discussion
SCRs with low tectonic loading rates and low levels of seismicity limit the knowledge of the occurrence of large earthquakes. They provide less data to PSHA than plate boundaries due to an absence of information about active fault geometries, fault slip rates, lock widths, directivity, associated maximum magnitudes, etc. It is standard practice to delineate SSZs instead of active faults when there is an apparent lack of correlation between geological structure and seismic activity, as already remarked on by Cornell (1968), but the delineation process remains largely subjective. The level of detail included may lead to SSZs with too small amounts of data for proper statistical analysis (e.g., De Crook et al. 1989), but formally enlarging SSZs to expand the data set may lead to unreliable results (e.g., De Crook and Egozcue 1992). Musson (1997) advocated designing alternative zonings based on geology and seismicity, and opposed focusing on major features only, since moderate earthquakes can break relatively short fault segments. SSZ models have been designed and existing models modified in the framework of Pan-European hazard mapping projects based on the common criteria for geological structures of different scales, tectonics, and seismicity (e.g., Arvidsson and Grünthal 2009). Musson (2000) pointed out the benefits of testing the SSZ model against the assumption of equal probability of earthquake occurrence: homogeneity of the spatial distribution of seismicity within each polygon can be demonstrated by statistical testing. Combined with the dominance of the host SSZ, the absence of guidelines on SSZ delineation means that large uncertainty may be incorporated into PSHA in low-seismicity regions. It is thus reasonable that SSZ models continue to be included in sensitivity analyses, although there are other important aspects of PSHA, such as the probabilistic computations (Beauval and Scotti 2004).
In the case of site L, the host SSZ stands out as the geologically uniform Vyborg rapakivi granite batholith, but the observed seismicity is strongly clustered in space, also confirmed by testing. The modification of the host SSZ attempted to obtain a more evenly distributed seismicity, which came with a significant effect on the outcome (Fig. 2a). This SSZ also illustrates well the trade-off between obtaining a statistically sufficient amount of data and including local features, such as shallow swarm-type seismic activity. SSZs ignoring the entire batholith have also been delineated in the past. Moreover, the relation between the observed low-magnitude swarms and potentially larger earthquakes is unknown but very relevant to seismic hazard at site L. The effect of M max is significant especially for low-frequency hazards at AFE 10 −5 , while for PGA the effect is more moderate (Fig. 8).
In a pre-SENSEI PSHA project, two SSZ models for site H were delineated independently by two groups based on careful scrutiny of geological, geodetic, geophysical, and seismological data (Korja and Kosonen 2015). Despite such efforts, the models were supplemented by a third one in the SENSEI project, mainly to compare seismicity rates of the previous models with those of an overlapping, preexisting SSZ model designed for sites L and O (L-O model). One SSZ of the L-O model was incorporated into other compatible SSZs in the original models for site H, and an additional logic-tree branch accounting for the uncertainty of the original host SSZ was considered. This is illustrative of the ambiguity in interpreting the distribution of past earthquakes.

3
The delineation of SSZs tends to have a direct bearing on the GR parameter values. In the case of sparse data, even the different interpretations of historical earthquakes may affect them (Grünthal and Bosse 1997). The trade-off has been tackled by determining prior b GR values for superzones, and updating them if supported by local data, when mapping Pan-European seismic hazard (e.g., Woessner et al. 2015). In the SENSEI project, it was noticed how the GR parameters have a significant impact on hazard at AFE 10 −5 (Fig. 7). One would expect the effect to diminish for higher AFEs, which fall into the range of actual observations; however, the uncertainty of the GR parameters may also preclude stabilization of the hazard for high AFEs (Figs. 7g,h). The scarce seismicity data in the southern part of the country lead to higher uncertainty of the GR parameters, creating a large spread of hazards for high AFEs. The more abundant data in the northern part result in a narrower spread of the hazard for site H (Fig. 7i).
From the tornado plots at AFE 10 −5 (Fig. 8), the median prediction of NGA-East emerged as a strong parameter influencing the level of hazard. This highlights the need to develop epistemic uncertainty estimates for GMPEs or include a sufficient number of GMPEs in the logic tree to reduce bias (Sabetta et al. 2005;Bommer and Scherbaum 2008). The selected GMPEs have more significant impact than weighting of the GMPE branches in the logic tree, as long as the weights do not favor any given GMPE (Sabetta et al. 2005).The GMPE sensitivity calculations reconfirmed that any GMPE must have a well-quantified median-prediction variability, or σ, which can only be achieved by calibration using a large set of ground motions. The importance of the GMPE prediction variability with lower AFE is shown in Fig. 7(g, h, i). The effect of the median GMPE prediction and σ is quantified for AFE 10 −5 PGA in Fig. 8.
Preliminary comparisons showed that Fenno-G16 is aligned with the high predictions of the NGA-East suite of GMPEs for the key range of magnitudes M w 4-6; it also has σ in the range of 0.85-0.9 similar to the range of high ergodic σ for NGA-East (0.686-0.901). In contrast, the low and central σ of NGA-East are in the range of 0.415-0.687 and 0.529-0.792, respectively. Fenno-G16 is compatible with some higher branches of the NGA-East suite, and gives a larger median prediction than the weighed-average of NGA-East. Intriguingly, the difference is more significant for the range of M w 4-5, in which Fenno-G16 is better constrained with data, and decreases for M w > = 6 (see Fig. 13 of Fülöp et al. 2020). This could indicate that attenuation in Fennoscandia is even lower that in CENA, where the data for NGA-East originates.
The attenuation-related epistemic uncertainty should be quantified and included in the GMPE used. For example, Toro (2006) pointed to the epistemic uncertainty in the median ground motion (σ μ ) as a major sensitivity factor in PSHA: σ μ plays an important role in regions with few recordings from damaging earthquakes and for evaluation of critical facilities requiring hazard at AFE 10 −4 or lower. Here, the effect of σ μ is explicitly included with the use of the 17 branches of the NGA-East GMPE and was found to be dominating Fig. 7 Influence of input parameters on the peak ground acceleration hazard curves for sites Loviisa (L), Olkiluoto (O), and Hanhikivi (H) using the final logic tree. The effect of the a, b, c maximum magnitude (M max ), d, e, f Gutenberg-Richter a GR and b GR parameters, and g, h, i ground-motion prediction equation (GMPE) median prediction. The gray lines are hazard curves from all branches of the logic tree. Orange, green, and blue lines are the means of all branches crossing the M max , a GR and b GR , and GMPE nodes in Fig. 4. The number of colored lines corresponds to the number of nodes in the branch. The black line is the overall mean hazard together with a GR and b GR (Fig. 8). The finding is in line with those in regions of moderate seismicity in Europe (SWISSNUCLEAR 2004). The target PSHA result in Finland is traditionally the AFE 10 −5 median hazard curve. Current national guidelines stipulate a minimum threshold of 0.1 g PGA median hazard to be accounted for in the safety studies in the country (YVL 2019). The AFE 10 −5 mean hazards are shown in Fig. 8, resulting in higher hazard values than the previous outputs, in which only the H site exceeded the 0.1 g threshold. Nuclear regulators interchangeably use median and mean hazard as a basis for decision-making, gravitating toward the mean hazard (OECD 2019b). Mean hazard curves appear to be less stable and sensitive to the most Branching levels for the zoning maps (red), maximum magnitude (M max ) (orange), Gutenberg-Richter a GR and b GR parameters (green), ground-motion prediction equation (GMPE) median predictions (blue), and GMPE randomness (σ) (black). The order of the branching levels was kept the same as in all previous panels; thus, the typical tornado shape is not apparent. The square sizes are proportional to the weights of the logic-tree branches included; the weights at each branching level add to unity. The black vertical lines show the mean hazard severe alternative models included in the analysis, especially for very low annual frequencies of exceedance. On the other hand, mean hazard is more consistent with risk estimation (e.g., McGuire et al. 2005). The remarkable stability of the median hazard was noticed during the pruning process (Table 1), when the preexisting logic tree could be reduced to a few branches preserving the median hazard. This stability also means that the epistemic uncertainties incorporated into the logic trees played only a small role in the original models. The additional arguments concerning inconsistency with fragilities in risk calculations (McGuire et al. 2005) presented a compelling case to focus on the mean hazard outputs in this investigation. We argue that the SENSEI project format can be seen as a possible general framework for providing guidance on how to conduct a PSHA in a low-seismicity country with a restricted expert pool. A well-known example of a body providing both technical and procedural guidance on PSHAs is the Senior Seismic Hazard Analysis Committee (SSHAC; Budnitz et al. 1997aBudnitz et al. , 1997b. A system of rating different sophistication levels for SSHAC PSHAs has been proposed (e.g., Hanks et al. 2009). High-level SSHAC projects have also been executed in SCRs, such as the Central and Eastern United States (USNRC 2012) and South Africa (Bommer et al. 2015). The SSHAC procedure has also influenced other methodologies to PSHA (e.g., Aldama-Bustos et al. 2019;Tromans et al. 2019).
The SENSEI project was markedly distinct from any other PSHA in Finland, or any SSHAC-level procedure, because it did not aim to arrive at final results. However, some parallels nevertheless suggest themselves: the SENSEI project involved participants focusing on different aspects of the PSHA and seismic risk analyses. The EG members were unaffiliated with the national stakeholders and did not participate in any previous PSHA in the country. It represented a partnership approach, in which all participants contributed their expertise and shared knowledge throughout the project. The members of the expert group were unaffiliated with the national stakeholders. The SENSEI project was a costeffective way to explore, transfer, and increase knowledge. The higher-level SSHAC procedures are highly structured, with the project participants having clearly defined roles, whereas the SENSEI project was more like a forum for discussion steered by the feedback from the expert group.
In a low-seismicity region, the general sparsity of empirical data makes capturing the related epistemic uncertainties very challenging and relevant. The models are typically logic trees, in which the knowledge and uncertainties of key parameter values and conceptual models are incorporated. Substantial aid can be obtained in the process from the SENSEI type of approach.

Conclusion
The SENSEI project consisted of a sensitivity investigation of PSHA for the three sites of interest for the nuclear field in Finland. The project scrutinized the existing consensus and, at the same time, the expert interactions allowed exploration of different alternatives. No final PSHA outcome had to be provided. The project format has clear advantages for countries with low seismicity and restricted expert pools focusing on seismic hazard.
The sensitivity investigation focused on the basic input parameters of the PSHA. The main sources of uncertainty have been identified as the GR parameters and the epistemic uncertainty of the GMPE, both amplified by data scarcity. The GR parameter estimation is linked to the SSZ delineation. With the GMPE used, the epistemic uncertainty related to attenuation should be quantifiable and included in the PSHA model. On the basis of the outcomes of the present investigation, we recommend the exploration of PSHA sensitivity with focus on the host SSZ with its immediate vicinity and the b value around the site of interest.