1 Introduction

The Philippine archipelago is situated between the two converging major tectonic plates, the Sundaland-Eurasian plate moving southeastward and the Philippine Sea plate moving northwestward. These converging plates resulted in the complex tectonic features of the country characterized by two major seismotectonic regimes: the surrounding offshore subduction zones and the in-between crustal faults. These tectonic features make the Philippines vulnerable to high-magnitude seismic events (Bautista et al. 2001; Kreemer et al. 2003; Pailoplee and Boonchaluay 2016). For example, there are about 95 major earthquake events with a high moment magnitude (near or >7.0) recorded from 1960 to 2014, some of which have claimed a heavy toll on human life, and the economic and social resources of the country (Pailoplee and Boonchaluay 2016). Various areas of the Philippines, based on its earthquake distribution data (Bautista and Oike 2000; ISC 2022; PHIVOLCS 2022a, b; USGS 2022) experience different seismicity ranging from being aseismic zones to being high seismic hazard zones—areas which experienced up to nearly 8.0 surface-wave magnitude (Ms) earthquakes. As such, comprehensive seismic hazard assessment remains crucial to efficiently design earthquake-resilient structures and mitigate future seismic risks in the various respective locations of the country.

Several seismic hazard maps of the Philippines in terms of peak ground acceleration (PGA) have already been published (Acharya 1980; Su 1988; Molas and Yamazaki 1994; Thenhaus et al. 1994; Torregosa et al. 2002). The Philippine Institute of Volcanology and Seismology (PHIVOLCS) as the mandated government agency of the Philippines for monitoring, assessing, and mitigating risks and hazards related to earthquakes has published the most recent among them, the Philippine Earthquake Model (PEM) (PHIVOLCS 2017; Peñarubia et al. 2020). The common denominator among the previous seismic hazard maps of the country is the use of the Probabilistic Seismic Hazard Assessment (PSHA) method—the widely utilized method introduced initially and formalized by Cornell (1968). However, despite the wide application of PSHA, several works have been published expressing their critiques on the limitations of the PSHA method (Krinitzsky 1998; Castaños and Lomnitz 2002; Klügel 2007; Wang 2011; Kossobokov and Nekrasova 2012; Wyss et al. 2012; Panza et al. 2014; Mulargia et al. 2017).

These critiques on PSHA have demonstrated applicable limits in the physical and mathematical model of PSHA, as well as in its basic assumptions and application (Panza et al. 2012). An example critique of the PSHA is pointed out by Castaños and Lomnitz (2002), stating that the available historical seismicity data for the PSHA application are insufficient to develop sound and reliable statistical calculations. Similar commentary was arrived when analyzing one region of China with seismicity data from 780 BCE to 2015 by Zhang et al. (2021), wherein they concluded that the PSHA procedure with the available data is inadequate to derive reliable seismic hazard analysis. Klügel (2007) and Wang (2011) also demonstrated the incorrect treatment of uncertainty in the PSHA mathematical model. Moreover, Krinitzsky (1998) has shown the flaws of using the Gutenberg-Ritcher b-line in the PSHA probabilistic projection procedure. In practice too, PSHA results are not always realistic nor reliable (Mulargia et al. 2017). Many strong destructive earthquakes resulting in significant loss of human lives have occurred in regions categorized as low risk by PSHA hazard maps, especially the worldwide Global Seismic Hazard Assessment Program (GSHAP) maps (Kossobokov and Nekrasova 2012; Wyss et al. 2012).

With the critiques surrounding the current PSHA method, a relatively new multi-scenario and physics-based deterministic approach to seismic hazard assessment called the Neo-Deterministic Seismic Hazard Assessment (NDSHA) was introduced by Panza et al. (1996). In NDSHA, some of the limitations of PSHA are addressed. For instance, catalogue completeness at moderate-to-low magnitudes is not necessary for NDSHA, contrary to PSHA, since only potentially damaging events are considered (Panza et al. 2012). NDSHA uses a well-established physics-based model and considers upper-bound maximization in addressing uncertainties. Moreover, in practical application, NDSHA results in various countries remain largely unchallenged by real-world occurrences, with few exceptional events, e.g., the Morroco September 8, 2023 event which surpassed the expectation of NDSHA maps developed by Mourabit et al. (2013) and even of PSHA maps of Poggi et al. (2020). Indeed, NDSHA results have been corroborated in some other strong earthquakes (e.g., the Zemmouri-Boumerdes, Algeria, 2003 and Gujarat, India, 2011 events) where PSHA maps of GSHAP failed to predict (Zuccolo et al. 2011).

NDSHA was recently applied in various countries worldwide in place and in addition to PSHA including Egypt (Hassan et al. 2017), Greece (Moratto et al. 2007), India (Parvez et al. 2017), Sumatra, Indonesia (Irwandi 2015), Alborz Region, Iran (Rastgoo et al. 2018), and North Africa (Mourabit et al. 2013). Building upon the success of the application of the NDSHA method in different regions worldwide and considering the threat of future earthquake-related risks to the country, this paper introduced a regional-scale adoption of the NDSHA method in the Philippine archipelago with the specific objective of developing new seismic hazard maps of the country to supplement the existing PEM in comprehensively understanding its seismic hazard, corroborating these new seismic hazard maps with actual seismic event data, and subsequently identifying the high seismic hazard areas in the archipelago.

The rest of the paper is organized as follows. Section 2 details the methodology used in the application of NDSHA to the Philippines setting, together with the corresponding assumptions and input data needed in the computations. Section 3 discusses the preliminary results of the input data preparations and the results of the final computation. Specifically, the final computational results are new seismic hazard maps in terms of Peak Ground Displacement (PGD), Peak Ground Velocity (PGV), and Design Ground Acceleration (DGA) which have been corroborated by actual seismic events. Section 4 concludes the work.

2 Methodology and sources of data

2.1 NDSHA general procedure

The Neo-Deterministic Seismic Hazard Assessment (NDSHA) method considers the physical process of earthquake generation and wave propagation in a realistic medium and simulates realistic physics-based synthetic seismograms from which peak ground motion parameters can be extracted (Panza et al. 2001). This method can be applied on a regional scale where it only considers the seismic source characteristics and the path conditions of the seismic wave propagation up to the bedrock level or it may later be applied to a local scale where it also considers site-specific geotechnical conditions above the bedrock. In the regional-scale application such as used in this paper, the general procedure is outlined in Fig. 1. In this procedure, earthquake scenarios (represented by the seismic sources and their corresponding characteristics) and seismic wave path conditions (represented by the observation sites and their structural properties) need to be defined initially. Using modal summation and discrete wave number techniques, synthetic seismograms are simulated for each significant seismic-source-to-observation-site wave path (Panza 1985; Florsch et al. 1991; Pavlov 2009). Subsequently, peak ground motion parameters are extracted from these seismograms and are laid on maps to produce the NDSHA seismic hazard maps. For the in-depth discussions of the theories, developments, and applications of NDSHA, the readers are further referred to the works of Panza et al. (1996, 2001, 2002, 2012, 2013), Peresan et al. (2011), Magrin et al. (2016), Panza (2017), Rugarli et al. (2019), Panza and Bela (2020), Bela and Panza (2021), Mourabit et al. (2013), Hassan et al. (2017), Parvez et al. (2017), Rastgoo et al. (2018), Zhang et al. (2021), and Brandmayr et al. (2022).

Fig. 1
figure 1

General outline of regional-scale NDSHA application

Full size image

Generally, two groups of input datasets (see Fig. 2) are required to implement NDSHA: [i] the seismic sources (Fig. 2a) comprising of earthquake catalogues (seismicity), seismogenic zones, and earthquake focal mechanisms; and [ii] the regional structural models (Fig. 2b) comprising of the regional polygons and their corresponding one dimensional (1D) structural properties, to compute for the synthetic seismograms where ground motion parameters are extracted (Panza et al. 2001, 2012).

Fig. 2
figure 2

Outline of regional-scale NDSHA input data preparation

Full size image

2.2 Source data on seismicity

Earthquake catalogues represent important information that describes the historical seismicity of the investigated territory (Zhang et al. 2021). The catalogue typically contains information on the seismic events’ date (with exact time), epicenter coordinates (latitudes and longitudes), magnitude (with magnitude type in most cases), and focal depth. In this analysis, the sources of these data were: (1) the instrumental earthquake catalogue obtained from the PHIVOLCS (2022a, b) containing data recorded by seismic monitoring instruments from 1900 to 2022, and (2) the historical earthquake catalogue of Bautista and Oike (2000), the most comprehensive study on historical earthquakes in the Philippines, which contained inferred data on earthquakes happened before the advent of seismic monitoring instruments from 1600 to 1900. Specifically, the considered earthquakes were only those located between 116°E and 129°E longitude and between 3°N and 23°N latitude, defining the extreme coastal boundaries of the Philippines with an allowance of around 150 km due to the possible effect of earthquakes within such range.

Since the NDSHA application only considers potentially damaging earthquakes, this study only considered earthquakes with a magnitude of at least 5.0, which is conventionally taken as the magnitude of the least potentially damaging earthquake (D’Amico et al. 1999) and with a depth of at most 40 km since strong ground motion is mainly controlled by shallow sources (Vaccari et al. 1990). The 40 km maximum focal depth was considered since the maximum crustal thickness in the Philippines is around 33 km according to the Philcrust3.0 model (Parcutela et al. 2020) with the allowance for possible error in the focal depth measurement. Moreover, only seismic events in the works of Bautista and Oike (2000) with good quality reports (Quality A and B as defined in Baustista and Oike 2000) were considered. Aside from being poor-quality reports, the neglected events from Bautista and Oike (2000) mostly have a magnitude below 5.0.

2.3 Source data on seismogenic zones

In NDSHA computation, it is imperative to identify seismogenic zones, which are seismic source areas characterized by assumed uniformity in tectonic behavior. With the flexibility of NDSHA, the method also allows the inclusion of seismogenic nodes as a supplement to seismogenic zones. Seismogenic nodes are areas anticipated to become potent sources of high-magnitude earthquakes in the future but with usually no historical seismicity (Zuccolo et al. 2011). Nevertheless, in our analysis, we focused solely on seismogenic zones. This decision was made due to the necessity for extensive separate geological, geotectonic, and morpho-structural investigations to identify seismogenic nodes in the country (Gorshkov et al. 2003), which were beyond the scope of our analysis. Seismogenic zones were specifically delineated in our analysis based on (1) the spatial distribution of the seismicity of the country using the considered earthquake catalogue (PHIVOLCS 2022a, b; Bautista and Oike 2000) and the other available earthquake catalogues from seismic monitoring agencies, e.g., ISC (2022) and USGS (2022), and (2) the tectonic features of the archipelago considering specifically the major active faults and subduction zones within its vicinity. Each zone was assigned with a single representative focal mechanism (i.e., the centroid moment tensor (CMT) information specifically the earthquake’s strike, dip, and rake angle values that describe earthquake source deformation). This representative value was the focal mechanism of either the strongest event or the best obtained event within the zone extracted from the Global Centroid Moment Tensor (GCMT) database—the most comprehensive database for global earthquakes’ CMTs (Ekström et al. 2012).

2.4 Discretization, smoothing, and finalizing of seismic sources

In the computation, seismic point sources were defined through the discretization and smoothing of the historical seismicity (earthquake catalogue) and its subsequent delimiting within the seismogenic zones. In the discretization process, the area was divided into 0.2° × 0.2° grid cells from which the earthquake with maximum magnitude within each cell was chosen as representative magnitude and was placed in the cell center (see Fig. 3a). In the subsequent smoothing process, the magnitude of discretized seismicity within a central cell was compared with its adjacent cells, and the higher magnitude was chosen. This smoothing procedure serves to accommodate the spatial extension of seismic sources, compensate for possible catalogue incompleteness, and mitigate errors in earthquake localizations (Panza et al. 1990). The smoothing process can be done using different radii, e.g., 1, 2, and 3 as illustrated in Fig. 3b, c, and d, respectively. For a detailed description of the smoothing process, the readers are referred to Panza et al. (2001). In our analysis, a radius of 3 cells was adopted, aligning with methodologies applied in other NDSHA studies conducted globally (see e.g., Moratto et al. 2007; Irwandi 2015; Rastgoo et al. 2018; Zhang et al. 2021). Then, the smoothed point sources were delimited within the seismogenic zones. Meaning, only those cells with centers contained within the seismogenic zones were considered seismic point sources. In the cases where the cell within seismogenic zones had no specified magnitude or registered below 5.0, a magnitude of 5.0 was automatically assigned. This decision was based on the hypothesis that within identified seismogenic zones, the potential for significant earthquakes capable of causing damage remains consistent (Zuccolo et al. 2011). Moreover, a magnitude of 5.0 is conventionally taken as the least potentially damaging earthquake (D’Amico et al. 1999).

Fig. 3
figure 3

Sample illustration of the smoothing process using different radii

Full size image

Furthermore, the seismic point sources’ focal depth in the analysis was considered to be a function of magnitude (i.e., 10 km for M < 7.0 and 15 km for M ≥ 7.0) to simplify the computations, to account for the magnitude–depth relationship demonstrated in the statistical properties of shallow earthquake occurrences (Caputo et al. 1973; Molchan et al. 1997), and to account for the large errors generally affecting the focal depth for historical events reported in the earthquake catalogues (Panza et al. 2001).

2.5 Regional structural models

The second group of input datasets for NDSHA computation consists of regional structural models (Fig. 2b). These refer to the average physical properties of several flat anelastic layers of the earth from the bedrock of the crust down to the upper mantle defined in terms of thickness, density, P- and S-wave velocities, and their corresponding attenuation Q values. These structural models stand as the source-site path of the seismic wave propagation. For this part of our analysis, the Philippine archipelago was geographically divided into 1° × 1° areas, with adjustments to fit the complex physical topography of the Philippines (see also Parvez et al. 2017), each area having its defined one-dimensional (1D) representative structural model. The upper structural properties down to a depth of about 100 km were specifically extracted from the Litho1.0 from Pasyanos et al. (2014), the most recent crustal and lithospheric earth model. The attenuation Q values for P-waves, however, were calculated using the following equation: Qp = (9/4) Qs, (Anderson 1989; Stein and Wysession 2003). For deeper extension of structural models down to about 1100 km, which is necessary for NDSHA computation of synthetic seismograms, the most recent standard global reference structures, specifically the AK135 model (Kennett et al. 1995) for densities and waves velocities and the QRFSI12 model (Dalton et al. 2008) for attenuation Q values were considered. In the analysis, the considered bedrock has at least 1500 m/s shear wave velocity, that is hard rock (Soil Class A) in the National Structural Code of the Philippines (NSCP) 7th edition (ASEP 2016). Thus, it is worth noting that a layer of sediment with thickness varying from 0.15 to 5.40 km in Litho1.0 (Pasyanos et al. 2014) was neglected in this NDSHA study because of having lower shear wave velocity ranging from 300 to 900 m/s.

2.6 Simulation of synthetic seismograms and development of the new hazard maps of the Philippines

Synthetic seismograms were computed on every observation site within the studied territory considering the seismic wave generation from the defined various point sources and its propagation through the source-site medium with defined structural properties. Similar computation processes to the well-established procedures in literature (see Panza et al. 2001, 2012; Rugarli et al. 2019) were used in the analysis. The chosen observation sites for the computation were the nodes of the regular grids used in the discretization where the cell centers that fell within the seismogenic zones were also assigned as the seismic point sources as described earlier (see Sect. 2.4).

Synthetic seismogram simulation includes the complex computations of the modal summation technique (Panza 1985; Florsch et al. 1991) and recently of the discrete wave number technique (Pavlov 2009). After applying these computations in our analysis, each observation site had several seismograms (tens, hundreds, or even thousands) from which the seismogram with the maximum parameters was selected and associated with the site. In particular, seismograms were generated using a 1 Hz cut-off frequency, aligning with the resolution of commonly available data sources (Panza et al. 1996). This frequency cut-off was deemed adequate for producing seismograms with dependable peak ground parameters, such as peak ground displacement (PGD) and peak ground velocity (PGV). However, it should be noted that this frequency cut-off may not be sufficient for accurately estimating ground acceleration (Panza et al. 1999). The resulting PGD and PGV values were then superimposed onto the Philippine map, contributing to the development of new seismic hazard maps for the country.

In this study, the NDSHA acceleration parameter, known as the Design Ground Acceleration (DGA), was derived by extending the computation results to higher frequencies, specifically up to a period near 0 s, using an appropriate design response spectrum (for further details, refer to Panza et al. 2001). Our analysis utilized the normalized design response spectrum from the Uniform Building Code (UBC) 1988 edition (International Conference of Building Officials 1988), specifically for Soil Type 1. This choice was made because the UBC 1988 served as the reference code for the old NSCP, and Soil Type 1 encompasses bedrock with an S-wave velocity >0.8 km/s. Subsequently, the new seismic hazard maps describing the maximum ground shaking at the bedrock were developed by overlaying the obtained DGA parameters, together with the PGD and PGV parameters, at each node across the Philippine map.

In the analysis, to reduce computational effort, the source-receiver distance (epicentral distance) was kept below an upper threshold, which was taken to be a function of the magnitude associated with the source (Panza et al. 2001). Specifically, source-receiver distance of 150, 200, 400, and 800 km for M < 6, 6 ≤ M < 7, 7 ≤ M < 8, and M ≥ 8 events, respectively, was adopted similar to Parvez et al. (2017) and Hassan et al. (2017). Moreover, lateral heterogeneities of the wave paths were considered in a rough way. Meaning, that if the source-site path crossed one or more boundaries between adjacent structural models, the signal was computed assuming the model of the site as representative of the whole path (Panza et al. 2001).

3 Results and discussion

3.1 Seismicity and seismogenic zones

The combined seismicity of 1622 events from historical and instrumental earthquake catalogues is illustrated in Fig. 4. Specifically, the combined seismicity was composed of 129 events from the historical earthquake catalogues of Bautista and Oike (2000) and 1493 events from the instrumental earthquake catalogue of PHIVOLCS (2022a, b). In terms of magnitude distribution, the combined seismicity had 63 events with a magnitude of 7 and above, 275 events with a magnitude of 6–6.9, and 1284 events with a magnitude of 5–5.9. It is also noticeable that a significant number of these events occurred in the subduction zones surrounding the Philippines, particularly (1) the major Philippine Trench (subduction zone of the Philippine Sea Plate) along the eastern coasts of the archipelago from the Visayas islands down to the south of Mindanao Island (see the location of seismogenic zone 18 on Fig. 5) and (2) the Manila Trench (subduction zone of the West Philippine Sea Basin) along the Western coast of the Luzon Island (see the location of seismogenic zones 3 and 5 on Fig. 5). Furthermore, the narrow island of Palawan Province in the westernmost part of the Philippines connecting to the Northern Borneo of Malaysia is noticeably an aseismic zone with no significant nearby seismic events based on our considered catalogues.

Fig. 4
figure 4

Seismicity of the Philippines considered in this work

Full size image
Fig. 5
figure 5

The delineated seismogenic zones of the Philippines with the corresponding focal mechanisms

Full size image

The proposed delineated 30 seismogenic zones with their corresponding chosen focal mechanisms used in the analysis are shown in Fig. 5. The focal mechanism which corresponds to the possible fault mechanism of the corresponding seismogenic zones is represented in the illustration (Fig. 5) by the black-and-white ‘beachballs.’ Most of the seismogenic zones within the archipelago had strike-slip focal mechanisms (i.e., beachballs showing the crossing of the black and white colors), especially those within the Philippine Fault system from the central Luzon area down to the Davao region in the south (that is the seismogenic zones 9, 13, 20, and 23). While the coastal seismogenic zones where the subduction zones were located (seismogenic zone 18 for the Philippine Trench, seismogenic zones 5 and 3 for Manila Trench, and seismogenic zone 25 for Cotabato Trench) showed either normal or reverse focal mechanisms (i.e., beachballs not showing the crossing of the black and white colors). These chosen representative fault mechanisms were the focal mechanisms of the maximum earthquake that had occurred historically within the region or at least the best-obtained earthquake event.

3.2 Result of discretization, smoothing, and delimitation

Figures 6, 7, and 8 show the results of the discretization process, the smoothing process, and the delimitation of the seismic sources’ spatial coverage with the boundaries of the seismogenic zones, respectively. In Fig. 6, each cell with recorded historical earthquakes is represented by one seismic source placed on the cell center, while those cells with no recorded seismicity do not have seismic source representation (empty cells). Most seismogenic zones had several empty cells. With the subsequent smoothing process, these empty cells were mostly filled with seismic source representation which even expanded past the seismogenic zones (Fig. 7). However, only seismic sources within the polygons of the seismogenic zones were considered. Hence, seismic sources’ spatial coverage was delimited as shown in Fig. 8.

Fig. 6
figure 6

The result of the discretization process of the seismicity of the Philippines

Full size image
Fig. 7
figure 7

The result of the smoothing process of the seismicity of the Philippines

Full size image
Fig. 8
figure 8

The seismic point sources after discretization, smoothing, and delimitation with seismogenic zones

Full size image

3.3 Structural properties

For the second input data of this regional-scale NDSHA study, 84 regionalized polygons were produced (Fig. 9), starting from polygon 1 of the western tip of Palawan up to polygon 84 of the southeastern tip of Davao peninsula. Each of these polygons had its defined own structural properties. These polygons were all originally 1° × 1° square areas—the same resolution of the Litho1.0 of Pasyanos et al. (2014), the source of the shallow structural properties in this work—covering all the Philippine territory including the off-land sea territories. However, those that did not cover any island (or inland territory) were excluded and some polygons were adjusted to cover a significant portion of inland only, since seismic hazards in this work were assessed specifically on areas of possible human habitation.

Fig. 9
figure 9

The 84 structural polygons considered in the study

Full size image

3.4 Result of the NDSHA computation

The NDSHA computation resulted in six (6) different maps that illustrate the spatial distribution of the maximum ground motion parameters (i.e., PGD, PGV, and DGA) of the Philippines, as shown in Figs. 10, 11, 12, 13, 14, and 15. PGD maps (Figs. 10 and 11) are expressed in terms of cm, PGV maps (Figs. 12 and 13) are expressed in terms of cm/s, and DGA maps (Figs. 14 and 15) are expressed in terms of “g” ratio, where “g” is the acceleration due to gravity with an approximate value of 980.665 cm/s2. Figures 10, 12, and 14 are related to horizontal ground motion, while Figs. 11, 13, and 15 are related to vertical ground motion. The resulting vertical ground motion parameters represent a unique attribute of NDSHA. Other standard methods (e.g. PSHA) rely on the assumption that the amplitude of the vertical component of strong motion, which could be defined as a fraction of the horizontal one, generally could be less than the horizontal components (Hassan et al. 2017). However, this is not necessarily true for high-frequency ground motion in the near-source condition (Reiter 1991; Shrestha 2009), since directivity, propagation effect, and local site condition may combine and produce a dominant vertical component (Hassan et al. 2017).

Fig. 10
figure 10

Horizontal PGD map

Full size image
Fig. 11
figure 11

Vertical PGD map

Full size image
Fig. 12
figure 12

Horizontal PGV map

Full size image
Fig. 13
figure 13

Vertical PGV map

Full size image
Fig. 14
figure 14

Horizontal DGA map

Full size image
Fig. 15
figure 15

Vertical DGA map

Full size image

It is important to note that the values in the produced regional-scale hazard maps are presented as discrete ranges of values in geometrical progression close to 2, consistent with the discrete graduation of macroseismic intensity (Panza et al. 2001). The factor of 2 is comparable with the real resolving power of available global data sources (see, e.g., Cancani 1904; Lliboutry 2000). Moreover, this approach helps to account for the inherent uncertainties and limitations of NDSHA, which are characteristic of any physical earthquake model (Panza 2017). More specific hazard estimates can be obtained at a local scale through specialized studies, such as those conducted by Rugarli et al., (2019). This is particularly important as local site-specific conditions above the bedrock level might cause possible amplification or de-amplification of ground motions.

The result of this NDSHA provides more ground motion parameters than in the PEM of PHIVOLCS (2017) which has only horizontal PGA values. The newly developed maps (Figs. 10, 11, 12, 13, 14, and 15) have both horizontal and vertical ground motion parameters for displacement, velocity, and acceleration. The PGV parameter is considered to be a preferable representative measure of earthquake intensity, which is directly connected with energy demand (Kouteva-Guentcheva 2010). The horizontal DGA parameter from these maps is comparable to the PGA values of the PSHA maps of the PEM. The horizontal DGA values and the horizontal PGV values can be related to the seismic acceleration coefficient, Ca, and seismic velocity coefficient, Cv, in the NSCP 7th edition (ASEP 2016) used in the seismic design of regular structures in the Philippines. The outcomes of NDSHA offer a wealth of earthquake parameters crucial for design considerations. Key ground motion parameters pertinent to seismic engineering or design can be easily extracted from the database of synthetic seismograms generated during the analysis. These synthetic seismograms retain all necessary data, including peak values, ensuring their accessibility for further utilization and analysis (Romanelli et al. 2010). Parameters such as the period of the peaks, the spectral acceleration at different spectral periods, and even the complete time series are readily available. Structural designers and engineers can leverage these parameters as seismic inputs, particularly in the design of seismic-resilient infrastructures (Panza et al. 2012). However, since the regional-scale NDSHA application only considers the ground motion on the bedrock level, a significant layer of sediment with varying thickness up to 5.40 km was neglected in the analysis, so the discrete values of these maps (Figs. 10, 11, 12, 13, 14, and 15) still need adjustment to suggest exact values with the further application of local-scale NDSHA.

One notable observation on these maps (Figs. 10, 11, 12, 13, 14, and 15) is that the ground motion values of the Palawan and Sulu archipelago in the eastern part of the country are relatively lower than the ground motion values of the rest of the country. This is a well-established knowledge regarding the country’s seismicity as conformed by the referenced seismic map of the Philippines in the National Structural Code of the Philippines (NSCP) 7th edition (ASEP 2016).

From these maps (Figs. 10, 11, 12, 13, 14, and 15), high-level hazard areas in the country with maximal values for DGA in the range of 0.6–1.2 g, PGV in the range of 60–120 cm/s, and PGD in the range of 30–60 cm were identified. They are located around five areas as re-illustrated in Fig. 16. The first area includes the Central Luzon region and some parts of the Ilocos Region, Cordillera Administrative Region, Cagayan Valley Region and Calabarzon Region. This area is known for two of the most devastating earthquakes in Philippine history, the 1990 Central Luzon earthquake and the pre-instrumental 1645 Central Luzon Earthquake reported during the Spanish period (Punongbayan et al. 1992). The second area includes the southeastern part of the Davao region in Mindanao which is known for its proximity to the seismically active southern part of the offshore Philippine Subduction Zone (PSZ) and the southern segment of the inland Philippine fault system (PFS). The third area includes a large portion of the Soccsksargen region in western Mindanao specifically the provinces of Sultan Kudarat, South Cotabato, Maguindanao del Sur, and the western part of Sarangani. This area is known for the destructive offshore 2002 Palimbang earthquake and the epicentral location of the 1976 Moro Gulf Earthquake. The fourth area includes the Western Visayas Region which is the location of the 1948 M8.3 Antique Earthquake and other notable earthquakes. The fifth area is the northern shore of Cagayan Valley including the surrounding areas and nearby islands which is the epicentral location of a pre-instrumental 1619 offshore earthquake with a magnitude around 8.00. Another notable area with high DGA value is the surrounding areas of the epicentral location of the recent July 2022 Luzon earthquake with 7.0 magnitude. This region, however, has lower PGV and PGD values.

Fig. 16
figure 16

Locations of high-level hazard areas

Full size image

3.5 Validation

To verify the validity of the resulting NDSHA seismic hazard maps, these were compared with actual seismic events. Six recorded actual earthquakes with observed vertical and horizontal PGA data provided by PHIVOLCS were used. Table 1 shows the comparison. From the table, the observed PGA’s are all lower than the upper limit of the calculated DGA class ranges except for one event, the 22 August 2017 Leyte earthquake. This suggests that historical events corroborate the new NDSHA earthquake hazard maps. The notably high PGA values of the Leyte earthquake on August 22, 2017, which surpassed the calculated horizontal DGA and substantially exceeded the calculated vertical DGA, may be attributed to two primary reasons. Firstly, the earthquake’s epicenter lies near the liquefiable areas within the province, as indicated by the liquefaction susceptibility map of the Philippines (PHILVOLCS 2018). Such proximity to liquefiable zones could significantly amplify ground motion, thus accounting for the observed disparities. However, it is important to note that the regional-scale NDSHA maps provide estimates of ground motions solely at the bedrock level, without factoring in specific site conditions above the bedrock. Secondly, there is a possibility that the raw data provided by PHILVOLCS (2022c) may require adjustments. This conjecture arises from the fact that earthquakes with Intensity VI should have lower PGA values than earthquakes with Intensity VII. Moreover, the observed PGA’s of the two Intensity VI events in Table 1 (i.e. the 5 March 2017 Surigao del Norte earthquake and the 22 August 2017 Leyte earthquake) exceeded significantly the equivalent PGA range of events with Intensity VI in Modified Mercalli Intensity (MMI) Scale and PHIVOLCS Earthquake Intensity Scale (PEIS) which is equal to 0.092–0.18 g (Wald et al. 1999; see also Inoue et al. 2015 for PEIS and MMI relationship).

Table 1 Comparison of actual events’ observed PGAs (PHIVOLCS 2022c) and their corresponding calculated DGAs
Full size table

Assuming no data adjustment was needed, the student’s t test was performed to statistically compare the mean values of the observed raw data of actual seismic events and their corresponding calculated data. The comparison was specifically between the observed PGA’s (horizontal and vertical) and the upper limits of the corresponding calculated DGA’s. The hypothesis posited that the observed PGA is greater than the upper limit of the calculated DGA. The Student’s t test result in Table 2 shows that the p values are closer to unity, indicating that the observed horizontal and vertical PGA values of the actual events are not significantly greater than the upper limit of the class range of the calculated DGA’s of NDSHA. In other words, the actual ground accelerations corresponding to the major seismic events in Table 1 have not been statistically greater than the upper limit of ground accelerations calculated by the NDSHA in this work. Thus, the NDSHA maps (Figs. 10, 11, 12, 13, 14, and 15) have been corroborated by the observed seismic events’ data.

Table 2 Result of the statistical comparison using student’s T Test between the observed PGA of actual events and the upper limit of the corresponding calculated DGAs on their locations
Full size table

A recent significant earthquake event, registering a magnitude (Mw) of 7.4, struck Surigao Del Sur, Mindanao Philippines, with its epicenter located at 8.44°N, 126.59°E, on December 2, 2023 (Llamas et al. 2024). This event presented an invaluable opportunity to further corroborate the newly developed maps. This event had a maximum reported intensity of PEIS VII (Destructive) equivalent to MMI VII (Llamas et al. 2024). According to the correlation established by Wald et al. (1999), Intensity VII corresponds to a horizontal PGA range of 0.18–0.34 g and a horizontal PGV range of 16–31 cm/s. Conversely, our developed maps indicate DGA and PGV values ranging from 0.30 to 0.60 g and 30 to 60 cm/s, respectively, in areas surrounding the epicenter. Notably, both the maximum calculated DGA and PGV values encompass the actual felt intensity of the event, providing strong evidence for the validity of the newly developed seismic hazard maps.

Comparison of the newly developed NDSHA maps with the existing PEM—developed using PSHA—is not discussed in this paper since the comparison will not make one validate the other. Our goal is to supplement the existing PEM with the new earthquake hazard maps of the Philippines using NDSHA, thus providing a more comprehensive depiction of the country’s seismic hazard landscape. Moreover, the differences in input data and assumptions (e.g. the period of the used seismicity data and the assumed characteristic of the soil where the hazards are estimated) for both PSHA and NDSHA should be considered for comparison, requiring better data availability.

4 Conclusion

New developed seismic hazard maps of the Philippines in terms of Peak Ground Displacement (PGD), Peak Ground Velocity (PGV), and Design Ground Acceleration (DGA) using the approach of regional-scale NDSHA by utilizing available data in the literature and repositories of government and non-government institutions (i.e., PHIVOLCS, USGS, and ISC) are presented in this paper. These maps show estimates of ground motions over the Philippine archipelago on the bedrock levels considering the earthquake wave generation and travel from the hypocenter of various simulated earthquake scenarios to the bedrock level of the observation sites. This NDHSA analysis was able to identify five high seismic hazard areas in the country with high PGD, PGV, and DGA values, specifically in the areas of Central Luzon Region, Davao Region, Soccsksargen Region, Western Visayas Region, and the northern shore of Cagayan Valley—areas known to have significant seismic events. Statistical validation and comparison were done to validate the resulting values of these NDSHA maps with acceptable results. These maps could deepen our understanding of the seismic hazard landscape of the country and could supplement the probabilistic perspective of the PEM. Furthermore, they opened potential avenues for future research works on further application of the NDSHA methodology within the Philippines. One promising avenue involves refining the methodology through the implementation of a more nuanced local-scale NSDHA. This entails incorporating site-specific soil conditions above the bedrock to accurately assess the potential amplification of de-amplification effects. Moreover, the utilization of NDSHA could augment the investigation, analysis, and design of critical infrastructures throughout the country. However, the successful realization of these pursuits hinges upon access to more detailed geological and geophysical data specific to the archipelago, including localized structural models having higher resolution. Such enriched datasets would enhance the efficacy and precision of future seismic hazard assessments and engineering endeavors in the country.