3.1. Methodology and Definition of Input Parameters
The first numerical formulation of PSHA was developed in 1968 by [
6], and the methodology has been followed by many studies such as [
56,
57,
58,
59]. In PSHA, all potential earthquakes are interpreted, taking into account their resulting ground motion levels and the likelihood of their occurrences. Instead of searching for a worst-case scenario, PSHA estimates the probability of different levels of ground motion intensity being exceeded within a certain time frame and at specific locations. In order to provide an overview of the PSHA process and its outcomes, a framework is depicted in
Figure 4, and the main steps of PSHA used in this study are outlined as follows:
- (a)
Definition of seismicity: The first step of PSHA requires a collection of events occurring within a specific timeframe and location, along with their associated information, such as magnitude, time, location, depth, etc. These events may include both instrumental and historical records.
In this stage, seismic activity data are gathered from several earthquake catalogs, including Instituto Português do Mar e da Atmosfera (IPMA) [
60], Instituto Geográfico Nacional (IGN) [
61], United States Geological Survey (USGS) [
34], International Seismological Centre (ISC) [
62], National Earthquake Information Center (NEIC) [
63], and Global Centroid Moment Tensor (GCMT) [
64]. IPMA [
60], considering the beginning of the instrumental era for Portugal (since 1961). Due to varied magnitude scales across earthquake catalogs, one of the challenges is in the process of homogenization of earthquakes. The IGN [
61] in Spain employs different magnitude types, further complicating catalog construction. A study by [
65] divides the earthquake catalog into three sub-periods, each linked to specific magnitude types, including L
g magnitude and body wave magnitude
mbLg(MMS),
mbLg(L), and
mb(VC). For this study, the M
w scale is chosen as the primary magnitude scale, with all data converted accordingly using formulas adapted from [
65] (see
Table 2 and
Table 3).
The earthquake catalogs from USGS [
34], ISC [
62], NEIC [
63], and GCMT [
64] is harmonized using the provided equations as given in
Table 3. For the magnitude scales, including surface wave magnitude and local magnitude (
and
, the following formulas are employed [
66,
67]:
The conversion of magnitude methodology aligns with a prior PSHA study [
45] for mainland Portugal. In summary, this approach ensures a standardized seismic catalog for PSHA in mainland Portugal, aligning with established practices. Herein, it is important to note that a minimum threshold magnitude of 3.5 has been assigned for the PSHA framework.
Figure 5 presents 3047 instrumental records corresponding to all earthquakes that occurred between 1961 and 2023, along with historical records before 1961. It is worth noting that the literature on historical earthquakes in Portugal predominantly provides information in terms of intensity level.
Once different magnitude scales are converted into
, a composite earthquake catalog is created by eliminating duplicates from different earthquake catalogs. Herein, it is important to note that different earthquake catalogs may provide the same earthquake with varying locations and magnitudes due to several factors related to the data collected by seismographs and the methodologies used to process that data. The variations in reported earthquake locations and magnitudes across different catalogs are a result of differences in seismograph network density and distribution, data processing algorithms, velocity models, magnitude estimation methods, and the availability of data. These differences highlight the importance of using multiple sources to get a comprehensive understanding of seismic events. The algorithm developed for this study prioritizes the selection of the
scale for earthquakes when it is available. In cases where
is not provided, the algorithm converts other magnitude scales to
using the previously mentioned empirical equation. To determine the final magnitude, the algorithm selects the highest converted value. For the earthquake’s location, the algorithm evaluates the latitude and longitude coordinates provided by different catalogs. When the coordinates from different catalogs are in close proximity, the algorithm assigns the earthquake’s location based on the most frequently repeated coordinates among the various catalogs. The next step involves performing a declustering analysis, which consists of two steps: pre-processing and magnitude-dependent filtering. The pre-processing step focuses on eliminating fore/aftershocks by applying a declustering process based on time and space windows determined by the methodology proposed by [
68]. This methodology incorporates time and space window approaches derived from previous studies [
69,
70]. The second step involves considering all events with magnitudes larger than 6.0 as main shocks. It should be noted that a certain level of uncertainty is involved in distinguishing between fore/aftershocks using the approach developed by [
68], which identifies 2049 mainshocks. Finally,
Figure 6 presents all the mainshocks in the study area.
Historic seismic data often suffer from incomplete records due to limited observational technology and sparse data collection practices. This lack of comprehensive seismicity information can lead to gaps in understanding the true earthquake activity of a region. Stepp’s completeness test [
71], developed in 1972, plays a crucial role in addressing these limitations by assessing the reliability and thoroughness of earthquake catalogs. By analyzing the magnitude–frequency distribution of recorded earthquakes, Stepp’s method [
71] identifies the point at which the catalog becomes complete for earthquakes of a certain size. This allows researchers to validate whether the historical data accurately reflects the region’s seismic activity or if smaller events may be underrepresented. Consequently, Stepp’s test helps adjust for historical data gaps and improves the reliability of seismic hazard assessments. For this reason, to determine the completeness level of the earthquake catalog, the algorithm by [
71] is adopted. The main goal of the completeness test is to determine whether the catalog includes all significant earthquakes within a given time period and region. To do this, the test analyzes the magnitude–frequency distribution of earthquakes, looking at how many earthquakes of different sizes are recorded. By plotting these data, the test identifies a completeness threshold (the smallest earthquake size that is consistently recorded). This threshold helps to ensure that the catalog is accurate and that all relevant earthquakes are included. As a result of this approach, the findings reveal that the number of earthquakes ranging from 3.5 to 4.0 magnitude exhibited a noticeable peak after 2001, while the 4.0–4.5 magnitude range experienced similar behavior in 1998. Furthermore, the cumulative number of earthquakes in the 4.5–4.9 and 5.0–5.4 magnitude ranges has shown consistent increases since 1985 and 1973, respectively. Conversely, the remaining magnitude classes do not display any discernible peak response, suggesting they are relatively complete datasets. This is also depicted in the graph provided by
Figure 7.
The mean rate of occurrence is defined by the Poisson process as follows:
where
,
,…
are the number of earthquakes per unit time interval (
).
represents the mean rate of occurrence, which is computed through dividing the number of earthquakes by the interval years.
The variance of the mean rate of occurrence,
is evaluated by the following equation based on the study of [
71]:
The unit time interval of one year provides the standard deviation of the estimate of the mean
, as follows:
where
represents sample length in years.
Finally,
Figure 7 presents the stability of the mean rate of occurrence for each magnitude class. The y axis refers to
(standard deviation of the estimated mean of the cumulative number of events), which is calculated using Equation (6). The graph includes fitting lines that are plotted based on the defined years for each magnitude class’s level of completeness. The mean rate of the magnitude class within that time interval is fixed to establish the trendline.
Figure 8 shows that the magnitude classes of 3.5–3.9 and 4.0–4.5 begin to align with the trendline within the last 30 years. Similarly, the magnitude class of 4.5–4.9 starts conforming to the fit line within the last 40 years. The magnitude class of 5.0–5.4 shows a fit with trend line 4 over the past 50 years. On the other hand, the remaining magnitude classes, including 5.5–6.0, 6.0–6.5, 6.5–6.9, and 7, are considered complete and maintain consistency over the past 59 years. It should be noted that the magnitude class ≥ 7 overlaps with the 6.5–6.9, which implies the same behavior.
- (b)
Definition of seismic sources: Identifying all potential seismic sources within a chosen site, which can be in the form of point, area, or line sources (faults). The goal is to characterize their mechanism types, slip rates, dip angles, and other relevant attributes. Additionally, the process involves eliminating inactive seismic sources, known as “dead” seismic sources, which are no longer capable of generating tectonic activity within the study area.
For the scope of PSHA, a comprehensive literature review is conducted for the mainland of Portugal, and several seismic sources, including area source models and line source models, are found in the literature [
15,
41,
72,
73]. These source models partially or entirely cover the study area. Therefore, seismogenic zones that fall entirely within the study area are directly taken for this research. In contrast, those that partially cover or exceed the study area boundaries are adjusted based on seismic activity to determine each activity parameter.
A total of 16 zones are involved in the SHARE source model [
19,
20,
31], which is plotted in
Figure 9a. The required adjustments, including rescaling of Source Nos 262, 255, 254, 250, 254, 253, 248, 247, 244, 243, and 0, were made in this model to make the Share model [
19,
20,
31] compatible with the study area. The final form of the SHARE Model [
19,
20,
31] consists of 15 areal sources, as shown in
Figure 9b, which depicts that most of the seismicity occupies zones 253, 247, and 254.
The study of [
41] identified and correlated some significant historical and instrumental earthquakes with the geographical and tectonic features of the region. The line source model is shown in
Figure 10a. These faults and fault traces are also described in the QAFI Database [
73] with related information about their associated parameters: geometric, kinematic, quaternary activity, slip rate, paleo earthquakes, and seismic parameters. For this study, the line source models by [
41,
73] were adapted. Also,
Figure 10b shows the line sources model with assumed potential area sources (PASes) within the region along with the past earthquakes. A PAS is employed in PSHA for regions with seismic activity without specific information about faults or tectonic regimes. A polygon is modeled based on the area’s observed seismic activity and considered a potential area source in such cases. This study defines six PASes within the region where no active faults have been identified. While modeling these polygons, the mechanism types of the seismogenic sources are considered to be aligned with the tectonic regime.
To sum up, two types of seismic source models are employed in PSHA analysis. The first model is the areal model, adapted from the area sources model. This model assesses seismic hazards from distributed sources over an area. The second model is the line sources combined with the defined six PASes. This model accounts for seismic hazards associated with specific fault lines with their characteristics and for unspecified regions without information regarding faults. Combining these two models achieves a more comprehensive and accurate assessment of seismic hazard, considering both distributed area sources and specific fault lines within the region. In order to provide a more detailed account of the geological and physical characteristics of the source models,
Table 4 and
Table 5 are provided. The source segmentation proposed by [
46] is adopted for identifying the area source model regarding source mechanism. Historical and instrumental records are used to correlate with the area sources. Herein an additional note is that as part of the evaluation process in PSHA, seismic sources located within a 300 km radius of the chosen site are considered, which is proposed by [
1]. Therefore, the area covers a radius of 300 km, which is shown in
Figure 11. These specific sources are highlighted in bold in
Table 4 and
Table 5. Among these sources, Source No. 253 stands out as it experienced the Lisbon earthquake in 1755 and recorded a maximum magnitude of 7.8 in 1969. This high activity rate and its involvement in devastating events designate Source No 253 as an active tectonic feature. In the region, HO and SWIM-1 are the closest faults to the potential epicenter of the significant Lisbon earthquake. Notably, the most prominent instrumental record with a magnitude of 7.8 also occurred near these faults (see
Figure 12). Another observation is that HO has the reverse type of mechanism, whereas SWIM-1 has a strike-slip mechanism despite their coinciding location, as shown in
Figure 11. This variation in fault mechanism can have significant implications for seismic activity in the region. According to the methodology of Wells and Coppersmith, 1994 [
74] utilized in the study, both Fault SWIM-1 and Fault HO have the potential to generate earthquakes with a maximum magnitude of up to 8.5. The approach by [
74] considers the fault dataset, including the fault mechanism and their rupture length, to estimate the potential earthquake magnitude. The fact that both faults can generate high-magnitude earthquakes further emphasizes their significance in SHA for the area. The coexistence of reverse and strike-slip mechanisms at the exact location might imply complex tectonic processes in the region, warranting careful consideration in assessing potential seismic risks. By incorporating this information into the seismic hazard analysis, the study ensures a more comprehensive understanding of the potential earthquake scenarios in the region.
Finally, to choose an appropriate maximum magnitude for use in PSHA for the associated faults or area sources, the study compares the instrumental, historical, and W&C magnitudes (only for the line source model) [
74]. The highest value among these three estimates is selected as the input for PSHA, considering it the potential maximum magnitude for SHA.
- (c)
Description of the magnitude recurrence relationship: This task involves the definition of a mathematical or statistical model used to evaluate the likelihood of earthquakes of different magnitude levels occurring over a specific time.
The magnitude recurrence relationship, also called magnitude–recurrence distribution or Gutenberg–Richter law introduced by [
75], is a key concept that describes the relationship between the magnitude of earthquakes and their frequency of occurrence within a period of time and region for each seismic source in PSHA. The Gutenberg–Richter law [
75] assumes that the distribution of earthquake magnitudes in a given region and period follows an exponential pattern. This relationship is mathematically represented as follows:
where
is the logarithm of the number of earthquakes exceeding a certain magnitude m within a specific time period and a specific region,
and
are constants that depend on the seismicity characteristics of the region, and m is the magnitude of the earthquake. The “a” parameter describes the expected number of earthquakes of magnitude equal to or larger than m. The parameter “b” is a slope that represents the rate at which the number of earthquakes decreases with increasing magnitude in the magnitude–frequency distribution.
In PSHA, the minimum value of magnitude, denoted as
, represents the expected intensity level or magnitude that could potentially cause damage. For this particular study, it is assumed that
is equal to 3.5. On the other hand, there is an upper limit, denoted as
, which represents the maximum magnitude earthquake that is likely to occur in the region of interest. This upper limit can be determined either by analyzing historical earthquake catalogs or by estimating it through empirical relationships between rupture dimensions and magnitude. Because of the constraints imposed by both the lower and upper magnitude limits, the probability density function of earthquake magnitudes can be expressed as shown below [
76]:
where
is a seismotectonic parameter that is utilized to describe the relative frequency of large earthquakes in comparison to smaller ones. It helps to characterize the distribution of earthquake magnitudes. The parameter
serves as a normalizing constant, which adjusts the cumulative distribution function in such a way that it reaches unity at the upper magnitude limit
.
To reduce uncertainty, instead of relying on an estimated value for the seismicity parameter λ, the average annual number of earthquakes with a magnitude equal to or higher than
is considered. The
for each seismic source zone is determined using two alternative algorithms: the least squares regression (LSR) and the maximum likelihood (MLH) approaches. The LSR aims to minimize the sum of squares of the differences between observed and estimated values, while the MLH seeks to maximize the likelihood of the estimated function fitting the observed data. The
value of the incomplete catalog is assessed based on the approach of Aki-Utsu [
77] as follows:
where,
is the average magnitude, and
corresponds to the minimum magnitude of the incomplete catalog.
The
values corresponding to the artificially completed earthquake catalogs are calculated by taking into consideration the time interval, which can be explained by the difference in the observed mean activity rate (λ). The algorithm for the MLH approach is proposed by [
78], which estimates
values by dividing an incomplete earthquake catalog into sub-catalogs to define the level of completeness of each catalog separately (i.e., sub-catalog 1, sub-catalog 2, …, sub-catalog n), as shown below:
where
is the Aki-Utsu [
77] estimator and
denotes the rate of events that have a magnitude equal to or greater than the level of completeness which is calculated as follows:
where
corresponds to the total number of events at a specific magnitude level for sub-catalog
i, and in the denominator, the total number of events across all magnitude levels is defined, symbolized by
.
Finally,
Table 6 and
Table 7 present the results in terms of activity rates of the area sources model and line sources model, respectively, for both complete and incomplete catalogs. Note that the activity rates labeled as “corrected” in
Table 6 and
Table 7 refer to values calculated using only the complete earthquake catalog obtained through the methodology described by [
71], as outlined previously. Despite the LSR method mostly tending to estimate higher activity parameters than MLH, the parameters from both approaches generally show a close agreement for the incomplete earthquake catalog. LSR may lead to larger activity parameters compared with MLH, particularly if there are large earthquakes, which can disproportionately affect the estimated parameters. However, for the case of the completed catalog, a noticeable variation is observed between the methods, along with larger values compared to the estimators for the incomplete earthquake catalog. According to the results obtained from both MLH and LSR methods, it is evident that the activity parameters of the line source model tend to overestimate the seismicity in the region, while the parameters from the area model lead to underestimation of seismicity in the region. For instance, sources of 247, 253, and 254 display higher λ values, implying that the past seismicity has a major composition of small earthquakes. On the other hand, the presence of lower
values such as in sources 242, 0, 243, and 262, indicated a larger proportion of large earthquakes relative to small earthquakes in the region. In the case of the line source model, most of the faults exhibit lower values of λ and
that can be attributed to the occurrence of larger magnitude earthquakes in close proximity to the faults, while only a few smaller magnitude earthquakes are observed.
- (d)
Selection of Ground Motion Models (GMMs): GMMs predict the intensity of shaking during an earthquake based on various factors, such as earthquake size and distance from the site. Choosing models that are well suited for the specific region improves the accuracy of the predictions, as models are often developed based on regional seismic data and conditions. It is important to review existing literature, compare different models, and seek expert recommendations to select the most suitable ones. This step involves choosing GMMs by considering various factors such as the region’s tectonic setting, geological conditions, and the available data on ground motion recordings to predict the potential ground shaking level at a particular site.
Numerous researchers have developed GMMs for the purpose of estimating ground motion intensity on a global scale [
27,
79,
80,
81,
82,
83,
84]. The process of selecting the most suitable GMMs for the mainland of Portugal involves several important steps. These steps include identifying seismic sources relevant to the region, considering the soil conditions and insights from local experts, and adhering to past studies that were conducted for the same region. The selection of GMMs for this study was meticulously considered to ensure their relevance to the seismic characteristics of Lisbon. According to recommendations by [
18], who conducted PSHA for Portugal, the model of [
27] is suggested for seismic sources located in inland areas with rock sites. This GMM is widely used across Europe and is based on data from highly active regions like Italy, Greece, and Turkey, as demonstrated by [
17]. Additionally, local experts for mainland Portugal endorse using the GMM of [
29] for SHA in the region. The mentioned GMMs are employed in the calculation process of PSHA to estimate PGA values at selected sites for this study. The reason for their utilization is their suitability and compatibility with the seismic characteristics of the region under investigation.
- (e)
Identification of Soil Class: This step involves the classification and characterization of different types of soil materials in the region. The seismic response of the ground under seismic load can differ according to the type of soil where seismic waves pass. Local soil characteristics play a critical role in how seismic waves are amplified. Soft soils can significantly amplify shaking compared to firmer soils or bedrock, affecting the intensity of ground motion at the surface. Incorporating site-specific soil data, such as soil density, stiffness, and layering, allows for adjustments to the GMMs to account for these amplification effects. It is noted that local soil characteristics are incorporated as the Vs30 value, which is used as an input for the GMMs.
For the soil model, the study conducted by [
85] is adopted. The study focuses on the LTV region, including MAL, which is densely populated with approximately 3.5 million inhabitants and has a history of devastating historical earthquakes resulting in significant economic and human losses. In response to this concern, the researchers investigated creating detailed average shear wave velocity in the top 30 m (Vs30) and soil classification maps for the LTV region. They achieved this by utilizing in situ shear-wave velocity measurements, along with P- and S-wave seismic velocities from seismic refraction data and some cross-hole datasets. Additionally, the mapping process incorporated lithostratigraphic studies and analyses of boreholes drilled for water supply and geotechnical investigations. They identified 82 soil profiles within Lisbon in terms of their geological, lithological features, age, and soil classes according to Eurocode 8 [
32].
Figure 12 shows the Vs30 map of the surrounding region of LTV, developed using the natural neighborhood interpolation method in ArcGIS PRO [
86]. The Vs30 map of the region depicts that the resolution of values is higher in the western part of LTV, which can be attributed to exhibiting a more precise spatial distribution of sites with known Vs30 in the region. Overall, the region is predominantly occupied by B soil type according to Eurocode 8 [
32]. In addition, the location of the National Museum of Archeology (MNA) in Monastery of Jerónimos, a monument with great importance in Portuguese history, is shown in
Figure 13 to assess its soil class and then calculate its hazard level later in PSHA.
- (f)
Utilized Logic Tree Algorithm:
Figure 14 shows the utilized logic tree algorithm, consisting of five branches: seismic zonation, earthquake catalog, seismicity models, Mmax, and GMMs. The first branch of the framework displays two seismogenic models, which are employed to assess seismicity levels within the region. Furthermore, in the context of PSHA, the analysis exclusively focuses on mainshocks, as indicated in the catalog branch. This is because the rate of ground motion intensity measure exceedance at the specific location is determined by the largest magnitude of an earthquake within each cluster of events, as explained by [
6]. Subsequently, two statistical techniques are employed to estimate seismicity parameters for input into the PSHA. The initial method is the MLH, which seeks to maximize the likelihood of the estimated function aligning with the observed data. The second approach is the LSR, which aims to minimize the sum of squares of the discrepancies between the observed and estimated values. The next step involves determining the maximum magnitude of an earthquake for each seismic source, which can be achieved through instrumental or historical records. In the case of the line source model, the estimation by [
75] is employed, taking into account the rupture length of the faults. The highest obtained value is then considered as M
max within the logic tree framework. Finally, the appropriate GMMs for Portugal, namely those proposed by refs. [
27,
29], are employed in PSHA.