Seismotectonic model and probabilistic seismic hazard assessment for Papua New Guinea

Abstract

Papua New Guinea (PNG) lies in a belt of intense tectonic activity that experiences high levels of seismicity. Although this seismicity poses significant risks to society, the Building Code of PNG and its underpinning seismic loading requirements have not been revised since 1982. This study aims to partially address this gap by updating the seismic zoning map on which the earthquake loading component of the building code is based. We performed a new probabilistic seismic hazard assessment for PNG using the OpenQuake software developed by the Global Earthquake Model Foundation (Pagani et al. in Seism Res Lett 85(3):692–702, 2014). Among other enhancements, for the first time together with background sources, individual fault sources are implemented to represent active major and microplate boundaries in the region to better constrain the earthquake-rate and seismic-source models. The seismic-source model also models intraslab, Wadati–Benioff zone seismicity in a more realistic way using a continuous slab volume to constrain the finite ruptures of such events. The results suggest a high level of hazard in the coastal areas of the Huon Peninsula and the New Britain–Bougainville region, and a relatively low level of hazard in the southwestern part of mainland PNG. In comparison with the seismic zonation map in the current design standard, it can be noted that the spatial distribution of seismic hazard used for building design does not match the bedrock hazard distribution of this study. In particular, the high seismic hazard of the Huon Peninsula in the revised assessment is not captured in the current building code of PNG.

1 Introduction

The territory of Papua New Guinea (PNG) is riven by the seismically active boundaries between three major tectonic plates and at least eight microplates (Baldwin et al. 2012; Koulali et al. 2015; Wallace et al. 2014). As a result, the region experiences significant seismic activity that includes large earthquakes, such as the 2018 Hela Province event (M7.5), which cause loss of life and widespread damage to buildings and infrastructure (McCue et al. 2018). In order to reduce societal and economic loss in future earthquakes, the likely ground motion from future seismic activity needs to be accommodated in PNG’s building code. The earthquake loading component of PNG’s current building regulations (Papua New Guinea National Standards Council 1983) is based on the work of Jury et al. (1982), which analysed the distribution of seismicity in the PNG region as known in the early 1980s. The 1982 seismic hazard map (Fig. 1) shows four zones within each of which the seismic hazard has a uniform value.

Fig. 1

The seismic zoning map of the national building code of PNG (Jury et al. 1982). It assumes that there is little variation in seismic hazard in the zones and the presented design PGA level for each zone is notionally the worst in the zone

A new seismic hazard model for PNG was produced in 2016 (PSHA16, Ghasemi et al. 2016), using probabilistic seismic hazard analysis (PSHA) and analysing a significantly larger and more complete earthquake database than was available in 1982. The PSHA16 results show a high level of seismic hazard on the Huon Peninsula and in the New Britain—southern New Ireland—Bougainville region, similar in part to the Jury et al. (1982) assessment. However, the PSHA16 study showed that the level of seismic hazard had in general been significantly underestimated in the 1982 study. In the present study we revise the PSHA16 model by implementing two key recommendations by Ghasemi et al. (2016): i.e., developing a more comprehensive earthquake catalogue for the region and including a fault–source model that makes use of recent Global Navigation Satellite System (GNSS) strain–rate studies (e.g. Koulali et al. 2015; Wallace et al. 2014) that estimate slip rates on major fault sources. To build the hazard model we also used two key global datasets that have become available recently: the latest versions of the ISC–GEM earthquake catalogue (version 6.0 of March 2019, an update of the catalogue described in Storchak et al. 2015) and a new subduction zone geometry model (Slab 2.0, see Hayes et al. 2018).

This manuscript describes the development of the 2019 PNG seismic hazard assessment (PSHA19) through: the development of a homogenous earthquake catalogue for the region; the development of seismic-source models for active faults and intraslab seismicity; and the selection and weighting of the ground-motion models (GMMs) used in the assessment. Finally, the updated hazard model is discussed in the context of the seismotectonic setting of PNG, and also in comparison with the seismic zonation map in the current design standards.

2 The seismotectonic setting of PNG

The tectonic setting of PNG is dominated by the northward movement of the Australian Plate since it’s separation from Antarctica in the late Cretaceous, and its convergence with the WNW movement of the Caroline and Pacific Plates to its north. The convergence of these three major plates is mediated by a combination of rotational and obliquely convergent motions of several microplates (Fig. 2). The complex interactions of these plates span the full complexity of tectonic processes on Earth (Baldwin et al. 2012) and there is pronounced seismic activity associated with subduction (both megathrust and intraslab), seafloor spreading and rifting, and collisional orogenesis.

Fig. 2

Tectonic setting of PNG. The boundaries of major tectonic plates (AU—Australian, PA—Pacific, CR—Caroline) and microplates or blocks (NBB—North Bismark, SBB—South Bismark, SSB—Solomon Sea, WLB—Woodlark, TBB—Trobriand, HLB—Highlands, PPB—Papuan Peninsula, ADB—Adelbert). The boundaries that are colored blue or black are considered as fault sources for the seismic hazard analysis, while boundaries colored red and gray are not. Toothed boundaries indicate convergence, with black-toothed boundaries indicating active subduction, gray uncertain activity, and blue active continental convergence. Microplate boundaries considered in the hazard analysis are: BTF—Bewani–Torricelli, SHTF—Southern Highlands Thrust, OSF—Owen Stanley, WBF—West Bismark, RMF—Ramu–Markham, WF—Weitin. Arrows show the direction of movement of the corresponding major plate (AU,CR and PA) in the ITRF2008 reference frame (Altamimi et al. 2011)

The complex tectonics of PNG have been the focus of many previous studies: its tectonic evolution (Baldwin et al. 2012; Davies 2012; Hall and Spakman 2003; Hamilton 1979; Holm et al. 2016; Pigram and Symonds 1991); its seismotectonics (Abers and McCaffrey 1988, 1994; Abers 1991; Abers et al. 1997; Kulig et al. 1993; Pegler et al. 1995); its active deformation (Abbott et al. 1997; Koulali et al. 2015; Wallace et al. 2004, 2014), and its crustal structure (Ferris et al. 2006; Stevens et al. 1998). Of the many studies of active tectonics in PNG, few if any agree on all the details of microplate configurations and the nature of the boundaries between them, and many disagree even on first-order features like the number of microplates. For example, although two of the most recent studies by Baldwin et al. (2012) and Koulali et al. (2015) agree in many details Baldwin et al. (2012); with microplates mostly based on Hill and Hall (2003) include a Sepik Block that is not present in Koulali et al. (2015).

3 Earthquake catalogue compilation and processing

A PSHA requires a reliable catalogue of past earthquakes from which robust estimates of earthquake recurrence parameters, and subsequently the expected level of ground-motions using ground-motion models (GMMs), can be made. Robust recurrence parameter estimates require a dataset that covers a wide range of magnitudes over a long period of time. To achieve this we unify the results from several available earthquake catalogues. Furthermore, because all modern GMMs are calibrated to the moment magnitude scale (M_W), we homogenize catalogue based on M_W to ensure consistency between predicted ground motions and earthquake magnitude recurrence rates.

We compiled a composite catalogue by merging three primary global datasets, i.e.:

• 1. The International Seismological Centre (ISC) catalogue (1900–2017), which is the most comprehensive database of earthquake locations and magnitudes for the region and time period of interest. It is a composite catalogue integrating solutions by global and regional partner agencies (Willemann and Storchak 2001).

• 2. The Global Centroid Moment Tensor (GCMT) catalogue (1976–2017), which is the most complete database of earthquakes with measured M_W and focal mechanism information for the region of interest (Ekström et al. 2012).

• 3. The joint ISC and Global Earthquake Model (ISC–GEM) catalogue (1904–2014), which reflects a significant enhancement of the ISC earthquake solutions for selected historical events (Storchak et al. 2015). The use of a consistent earthquake location procedure and all available phase arrival times to relocate each hypocenter, as well as its use of M_W magnitudes, are the key advantages of this catalogue.

We combined the data from these three catalogues into a single, unified catalogue, for which any catalogue entries having origin times separated by < 60 s and hypocenters separated by < 100 km in distance were assumed to refer to a single event. There was no need to refine the search criteria as both GCMT and ISC–GEM catalogues include only moderate to large events that are already listed in the ISC catalogue by other contributing agencies. The preferred location and magnitude for each event in the composite catalogue followed the selection hierarchies as detailed in Ghasemi et al. (2020) and listed in Table 1.

Table 1 Hierarchical criteria for selection of origin for homogenizing the catalogue

The resulting unified catalogue includes approximately 40,000 events in a geographical region bordered by 0°–12°S and 138°–160°E, covering the period 1900 to 2017 (Fig. 3). About 70% of events in the unified catalogue are those reported by the ISC, and a similar fraction of the magnitudes are expressed in terms of the m_b scale reported by either NEIC or ISC. By contrast, the preferred M_W measure is available for a relatively small number of events in the region, mostly for larger magnitudes (M_W > 5.0).

Fig. 3

Epicenters from the unified and homogenized earthquake catalogue of the PNG region covering the period 1900 to 2017

To homogenize the catalogue to the M_W magnitude scale, empirical models to estimate M_W converted from m_b or M_S were derived. Table 2 lists the derived empirical models using global as well as region-specific data. The models are evaluated using the general orthogonal regression (GOR) technique, which accounts for measurement errors in the m_b (or M_S) and M_W values, and provides a unique solution that can be used interchangeably between magnitude types. Below we outline the process of deriving a regression model for events reported in terms of m_b from NEIC and M_W from GCMT. We followed a similar process to derive the other empirical conversion models listed in Table 2.

Table 2 Hierarchical criteria for selection of magnitude for homogenizing the catalogue

Figure 4 shows the distribution of GCMT M_W versus NEIC m_b values at the global scale. It can be seen that there is a positive correlation between these magnitude measures, albeit with significant scatter and significant variation in the gradient of the data points, with possible “hinge points” at m_b values of about 4.5 and 6.5. This may suggest that a bilinear or exponential model may better explain the relationship between m_b and M_W. However, as demonstrated by Gasperini et al. (2013), the apparent flattening of the data distribution for m_b values lower than about 4.5 is an artefact which arises from the incompleteness of the M_W database for small earthquakes. In contrast, the other significant gradient change at m_b value of about 6.5 is due to the well-known m_b saturation at larger magnitudes (Kanamori and Anderson 1975). However, we note that there are few earthquakes in the M_W database with corresponding m_b values larger than about 6.5. This would limit the reliability of any fitted model in this range, regardless of their functional form. Furthermore, using nonlinear models would introduce large but poorly constrained correction factors in converting m_b values greater than about 6.5 to M_W (Fig. 4). This may bias the estimated annual rates of such large events that contribute significantly to the level of seismic hazard in PNG.

Fig. 4

Density (number of events per 0.1 magnitude intervals) of the events with reported GCMT M_W and NEIC mb values at the global scale. Empirical regression models between the magnitude scales are also displayed. Grey broken line represents the 1:1 relationship. Black dashed lines denote the boundaries given by one standard deviation of the linear model

In this study, given the above considerations, we chose a linear functional form to derive m_b–M_W empirical relationships (Fig. 4 and Table 2). Such models can be derived using either global or region-specific M_W datasets. We chose a model based on the global dataset as it is constrained by a more extensive dataset and the resulting model is not significantly different from a region-specific one. Figure 4 shows that the model fits the data reasonably well in the range of m_b 4.5–6.5. For data points with m_b lower than 4.5 or larger than 6.5, generally, the model under-predicts the observations. In practice, all of the events with m_b < 4.5 correspond to converted M_W values that are below the completeness level of the catalogue, and all events with m_b > 6.5 have either an observed M_S or M_W available, so there was no need to use M_W values converted from m_b in the range in which they are under-predicted. (Note that according to Table 2, M_S is preferred over m_b). Ghasemi et al. (2020) also explored the depth and time dependence of the magnitude conversion models used here, and concluded that the depth and time dependence was only marginally significant and therefore was not considered further.

In comparison with other similar magnitude conversions, the linear model of this study is in close agreement with linear model developed by Weatherill et al. (2016), produces steeper gradient in comparison with Scordilis (2006) and Das et al. (2011), and a gentler one in comparison with Di Giacomo et al. (2015) model. We did not explore in details the cause of the differences; however, the choice of regression technique (e.g. ordinary least-squares regression versus GOR), the criteria for selecting the dataset (e.g. magnitude truncation levels, time period of the dataset, etc.), and the tools implemented to perform the analysis are considered as the key contributors to the differences in the models.

PSHA methodology generally assumes a temporally homogeneous Poisson process to describe the occurrence of earthquakes, requiring that the inter-event times between successive earthquakes are statistically independent. In practice, to ensure the data fulfil this assumption we need to identify and remove dependent events (e.g., foreshocks and aftershocks) from the catalogue, which is known as “declustering”. There are numerous techniques of varying complexity used to decluster earthquake catalogues (see Stiphout et al. 2012, for an overview).

In this study, we considered two popular declustering algorithms referred to as mainshock-window and moving-window approaches, originally developed by Gardner and Knopoff (1974) and Musson (1999), respectively. Both techniques search for dependent events using magnitude-dependent time and distance windows, and the details of the implementations we used here are described by Weatherill et al. (2014). Several studies have modified the extent of spatial and temporal windows by considering the general characteristics of the seismicity in the regions of interest. In this study, we considered four time–distance window parameters as suggested by Gardner and Knopoff (1974) and Uhrhammer (1986) for California, Grunthal (as reported in Stiphout et al. 2012) for Europe, and Ghasemi et al. (2016) for PNG. Overall, we declustered the catalogue six times using either mainshock-window or moving-window algorithms in combination with different time and distance windows. Figure 5 shows the cumulative number of earthquakes with M_W⩾ 5.0 in the original (not declustered) catalogue in comparison with those in the declustered catalogues. There is a significant increase in the number of events since 1964 which is consistent among the original and declustered catalogues, which marks the advent of the ISC’s effort in 1964 to collect and analyse seismic data and publish a global earthquake catalogue. There are relatively few rapid changes in the cumulative number of events in the original catalogue. These temporary increases in the seismicity can be linked to the occurrence of dependent events that are typically removed from the original catalogue through the declustering process. However, we should note that for the majority of the “mainshocks”, the offspring events are likely to be masked by the high background seismicity of the region and cannot be visually detected in Fig. 5.

Fig. 5

Cumulative number of earthquakes with MW ⩾ 5.0 in the original (not declustered) catalogue and after declustering the catalogue following either mainshock-window or moving-window algorithms in combination with different time and distance windows

In Fig. 5, significant variations can be observed among declustered catalogues with a maximum of 79% and a minimum of 33% of events excluded from the original catalogue by applying the mainshock-window approach using Grunthal and Ghasemi windows, respectively. We should note that the results of the declustering process are sensitive to not only the selection of the time and distance windows but also to the algorithm used for declustering. As an example, it can be seen that the results of the moving-window approach using Uhrhammer and Grunthal windows are approximately similar with total removal of 70% and 67%, respectively, of events in the original catalogue. The opposite trend holds for the results of applying the same windows but using the mainshock-window algorithm. In this case, the Uhrhammer and Grunthal windows remove 56% and 79% of the events in the original catalogue, respectively. These different declustering approaches can result in different estimates for the a and b values in Gutenberg–Richter earthquake magnitude-recurrence law (Gutenberg and Richter 1944):

$$\text{log}_{10} \left( N \right) = a - bm$$

(1)

where N is the cumulative number of events having a magnitude > m. Ghasemi et al. (2020) show that, depending on the window and declustering algorithm used, the values of a and b estimated from the declustered catalogs can vary by as much as 35%, which will affect the estimated hazard levels in PNG.

We should emphasise that there is no “correct” way to decluster a catalogue. Catalogue declustering is not an “exact science”, as the same physical process rules the occurrence of dependent as well as mainshock events (Christophersen et al. 2011). In practice, this makes the comparison of the performance of different clustering methods a highly subjective process that even for a specific region varies from one cluster to another. In this study, as detailed in Ghasemi et al. (2020) we chose the mainshock-window approach using the Ghasemi et al. (2016) time and distance windows for declustering the compiled catalogue. This method in comparison with other selected declustering algorithms (Fig. 5), removes the minimum number of events while still producing a mainshock catalogue that is approximately Poissonian in time and space (Ghasemi et al. 2020).

To obtain robust earthquake recurrence statistics based on the declustered catalogue, we need to exclude earthquakes that are below the level of magnitude of completeness (M_C). The magnitude of completeness is defined as the minimum magnitude threshold above which all earthquakes in a space–time volume are reported in the catalogue (e.g. Wiemer and Wyss 2000). For a region of interest, this parameter may exhibit spatiotemporal dependency reflecting the evolution of seismic networks and advances in earthquake waveform processing techniques (Mignan and Woessner 2012). Ghasemi et al (2020) explored the temporal and spatial dependency of M_C parameter for the declustered PNG catalog by (a) examining the non-cumulative as well as the cumulative number of events as a function of time to determine when changes in network sensitivity occurred, and (b) comparing frequency–time plots for all of the area source zones in the region as defined in Ghasemi et al. (2016). These comparisons resulted in the three completeness models that are applicable to the entire region of interest and are illustrated in Fig. 6.

Fig. 6

Density of earthquake magnitudes with time for the catalogue declustered using the Ghasemi window. Solid lines show the candidate M_C models

The candidate models use a single-corner (i.e. M_W⩾ 5.0 since 1980), two-corners (i.e. M_W⩾ 7.5 since 1900, M_W⩾ 5.0 since 1980), or four-corners (i.e. M_W⩾ 7.5 since 1900, M_W⩾ 6.4 since 1925, M_W⩾ 5.6 since 1964, and M_W⩾ 5.0 since 1980). Ghasemi et al. (2020) found that, for the same declustering algorithm, the differences in Gutenberg–Richter parameters by applying different Mc models is less than 5%. We chose the four-corner M_C model as the most appropriate model which agrees well with the observations (see Fig. 6) while keeping the maximum number of events to obtain robust earthquake recurrence statistics.

Finally, in order to develop earthquake recurrence models for the megathrust and intraslab components of the seismotectonic model described in Sect. 4, it was important to classify earthquakes in the catalogue as being either shallow crustal, megathrust or intraslab events. We followed Ghasemi et al. (2020) in classifying as a megathrust earthquake any event occurring within a depth threshold of 50 km for the New Guinea Trench and 40 km for the New Britain Trench (both eastern and western segments), and whose hypocentre was within a Euclidean distance of 20 km from the respective subduction megathrust. All other earthquakes were either considered to be shallow crustal events if their depths were 40 km, or intraslab events if their depths were greater than 40 km. Shallow crustal events were either ascribed to the modelled crustal faults or as part of “background”, spatially distributed seismicity, according to the magnitude-based criteria described in Sect. 4.1.

4 Seismotectonic source model

Unlike previous PSHA study in PNG (PSHA16) that is exclusively point source based, we developed a seismotectonic model that combines the active faults and spatially distributed sources. Furthermore, as outlined below, this source model incorporates earthquake rate information from both the instrumental earthquake catalogue and geological and geodetic studies of neotectonic faults to characterize earthquake sources in the region.

4.1 Fault sources

Mainland PNG formed the leading edge of the Australian Plate after it separated form Antarctica in the late Cretaceous. As the Australian Plate moved northward during the Cenozoic, it created and collided with island arcs, and these have left a record in both a crustal fabric riven by several major active faults, and in active subduction along and offshore its northern margin. These accommodate the relative movement of the various microplates depicted in Fig. 2, and are therefore potential sources of large earthquakes.

Below we consider each of these faults, grouped in terms of subduction and crustal faults. Ghasemi et al. (2020) have reviewed in detail the geological evidence for the locations of the active fault traces and faulting regimes, as well as the earthquake history, and we only summarise these details here. The results of the GNSS studies of Wallace et al. (2014) and Koulali et al. (2015) are used to choose appropriate faulting parameters, which are listed in Table 3. These fault sources are intended to account for all earthquakes of magnitude 6 and above considered in the PSHA presented below, except for the Aure Fold-and-Thrust Belt (AFTB) and Bismark Sea Seismic Lineation (BSSL; th e hatched areas in Figs. 7 and 8). The AFTB and BSSL are zones of diffuse seismicity, indicating that internal deformation of the northern part of the Papuan Peninsula Block and relative motion between the Adelbert and North Bismark Blocks, respectively, are accommodated by complex fault systems instead of by one main fault. For this reason, these are the only zones where the shallow distributed seismicity is taken to include earthquakes with magnitude > 6.

Table 3 Fault parameters for active plate and block boundaries. AL and YC refer to Anderson and Luco (1983) and Youngs and Coppersmith (1985) magnitude–frequency distribution models

Fig. 7

Fault sources and historical earthquakes on mainland Papua New Guinea. As in \* MERGEFORMAT Fig. 2, AU and CR denote the Australian and Caroline Plates, respectively, while microplates are: HLB—Highlands Block, ADB—Adelbert Block, SBB—South Bismark Block, NBB—North Bismark Block, PPB—Papuan Peninsula Block, NNGB—North New Guinea Block. In addition to the New Guinea Trench megathrust, the major crustal faults indicated are: BTF—Bewani–Torricelli Fault, WBF—West Bismark Fault, RMF—Ramu Markham Fault and SHTF—Southern Highlands Thrust Fault. AFTB denotes the western margin of the Aure Fold-and-Thrust Belt. The megathrust faults depicted as black-toothed curves and the crustal faults plotted in blue are all considered as fault earthquake sources, while the crustal faults plotted in red are not. Also displayed are post-1990 Global CMT solutions for Mw ≥ 7, hypocenters from the ISC–GEM catalog with Mw ≥ 6.5 with a “≈” to indicate those which were tsunamigenic, and events from the GEM Historical Earthquake Archive (Albini et al. 2014), all with depth ≤ 40 km. Green, magenta, and orange toothed curves denote the western, central and eastern sections, respectively, of the Ramu–Markham Fault. The hatched areas denote the Aure Fold-and-Thrust Belt (AFTB) and Bismark Sea Seismic Lineation (BSSL), where distributed seismicity accounts for earthquakes with magnitude greater than 6

Fig. 8

Fault sources and historical earthquakes on mainland Papua New Guinea. As in \* MERGEFORMAT Fig. 2, AU denotes the Australian Plates, and the microplates are: SBB—South Bismark Block, PPB—Papuan Peninsula Block, WLB—Woodlark Block, TBB—Trobriand Block. In addition to the New Britain Trench and Trobriand Trough megathrusts, the major crustal faults indicated are: OSF—Owen-Stanley Fault, NF—Nubara Fault, RMF—Ramu–Markham Fault and WF—Weitin Fault. The megathrust faults depicted as black-toothed curves and the crustal faults plotted in blue are all considered as fault earthquake sources, while the crustal faults plotted in red are not. Also displayed are post-1990 Global CMT solutions for Mw 7.0, hypocenters from the ISC–GEM catalog with Mw 6.0 (with a “≈” to indicate those which were tsunamigenic), and events from the GEM Historical Earthquake catalog, all with depth ≤ 40 km. For the OSF, red indicates pure thrust (normal) on the western(eastern) segments, magenta oblique thrust (normal) on the central western (eastern) segments, and blue indicates pure sinistral slip. Also indicated is the rotation pole and rate of the AU/TBB Euler pole. The hatched area denotes the Aure Fold-and-Thrust Belt (AFTB), where distributed seismicity accounts for earthquakes with magnitude greater than 6

4.1.1 Subduction megathrusts

The New Guinea Trench is the only subduction megathrust affecting mainland PNG, and accommodates convergence between the Caroline Plate and North Bismark Block on the one hand, and the North New Guinea Block on the other (CR, NBB and NNGB, respectively, Fig. 7). It extends for almost 1000 km, from 134.5°E offshore of the northern coast of West Papua, Indonesia, to where at 143.5°E it appears to merge with onshore thrust faulting that forms the SW edge of the Adelbert Block (ADB in Fig. 7). The ISC–GEM catalogue includes 10 earthquakes with 7.0 ≤ Mw ≤ 7.8 and depth ≤ 40 km that occurred near the eastern (PNG) side of the New Guinea Trench, but farther west it has experienced large historical earthquakes—an M_W 7.9 event in 1914 (Okal 1999) and an Mw 8.2 event in 1996 (Henry and Das 2002). We have assigned the New Guinea Trench a maximum magnitude of 9.0, since its 1000-km length could host an event of this size (e.g., Strasser et al. 2010). We have also set 90 mm/years oblique convergence rate for this trench with a coupling factor of only 10% based on Koulali et al. (2015).

The New Britain Trench (NBT) accommodates convergence between the Solomon Sea and North Bismark Blocks (SSB and NBB, respectively in Fig. 8), with northward subduction of the SSB beneath the island of New Britain. The NBT is the most active of all the earthquake sources in PNG, accounting for four of the six PNG earthquakes of magnitude 8 and above in the ISC–GEM catalogue. All geodetic studies of PNG (Koulali et al. 2015; Tregoning et al. 1998; Wallace et al. 2004, 2014) find convergence rates increase from 80 to 90 mm/year off the western tip of New Britain at 148°E to over 150 mm/year at 153°E, where the NBT axis changes strike from WSW–WNE to NW–SE. We break the NBT into west and east segments at this change in strike, and adopt convergence rate of 132 mm/year from Koulali et al. (2015). Because none of the geodetic studies constrain the locking of the NBT megathrust, here we adopt a conservative approach and assume it is fully locked. If we assume the largest NBT earthquake would rupture the full 600-km length of the NBT, the scaling relation of Strasser et al. (2010) implies a mean M_W value of 8.7. However, we are guided by the conclusion of Kagan and Jackson (2013) study of global megathrust earthquake that “magnitude 9 earthquakes can be expected in any major subduction zone”; we therefore adopt M_W 9 as the maximum earthquake magnitude for the NBT.

The Trobriand Trough marks the subduction of SSB beneath the Papuan Peninsula (Fig. 8). The existence of a distinct Trobriand Block (TBB in Figs. 2 and 8) and activity on the Trobriand Trough are contentious issues in the tectonics of PNG. Most studies argue that the Trobriand Trough is no longer active, since there is no well-defined pattern of seismicity that might indicate subduction of the SSB beneath the Papuan Peninsula (Hall and Spakman 2003). Furthermore, the GNSS studies of Wallace et al. (2004, 2014) and Koulali et al. (2015), all make the case that the GNSS data require little or no convergence along the Torbriand Trough. However, since the entire SSB is submarine and therefore not accessible to direct observation, there appears no way for GNSS to constrain Trobriand Trough convergence. Our model for a separate Trobriand Block bounded in the north by the Trobriand Trough, follows that suggested by Kington and Goodliffe (2008) and Cameron (2014). Our motivation for including a nominal 10 mm/year convergence, with seismic coupling 0.5, is driven mainly by the Letz et al. (2016) study of a major historical earthquake that appears to have occurred on or near the Trobriand Trough in 1895. This event generated a tsunami with initial receding or “draw-down” polarity in the Trobriand Islands, indicating it was a thrust earthquake.

For the New Guinea and New Britain megathrust fault sources, Ghasemi et al. (2020) showed that the Anderson and Luco (1983) exponential magnitude–frequency distribution (MFD) fits the catalogue seismicity rates better than Youngs and Coppersmith (1985) hybrid exponential/characteristic earthquake MFD. However, we believe it would be wrong to conclude that the purely exponential model is a better representation of the actual earthquake activity than the hybrid model, because of two important sources of bias: (1) the simple approach described at the end of Sect. 3 for classifying events in the catalogue as megathrust earthquakes fails to screen out all low-magnitude intraslab or shallow crustal events, and (2) the limited duration of the catalogue fails to account for large-magnitude megathrust events with very long recurrence times. Both of these factors inevitably bias the megathrust-classified seismicity away from the hybrid model, which includes more large-magnitude at the expense of low-magnitude events. Therefore, despite the better fit of the exponential model to the megathrust-classified catalogue seismicity, we have chosen as a megathrust earthquake recurrence model a logic tree comprising two branches: 1) the hybrid exponential/characteristic MFD of Youngs and Coppersmith (1985), with weight 0.6, and 2) the exponential MFD of Anderson and Luco (1983), with weight 0.4. For the Trobriand Trough, the lack of any recorded seismicity suggests that the hybrid MFD of Youngs and Coppersmith (1985) may be more appropriate.

4.1.2 Crustal faults

The crustal faults considered as earthquake sources in this study are the Bewani-Toricelli, Ramu–Markham, Southern Highlands Thrust, Owen-Stanley and Weitin Faults (BTF, RMF, SHTF, OSF, and WF, respectively, see Figs. 7 and 8).

From west to east, the BTF, RMF and OSF correspond to the boundary between the North New Guinea, Adelbert, South Bismark and Trobriand Blocks, respectively, to the north and the Highlands and Papuan Peninsula Blocks to the south. These generally mark the boundaries of various episodes in which oceanic lithosphere was obducted and accreted to the northern margin of the Australian plate. From west to east, this system of faults undergoes a pronounced transition in slip rate and faulting regime.

The GNSS observations of Koulali et al. (2015) establish that the present-day oblique convergence along the northern coast of mainland PNG is partially partitioned between the New Guinea Trench and the BTF, with both accommodating oblique sinistral motion of about 12 mm/year. The GNSS observations are not dense enough to determine exactly which fault strand is active, and so the BTF indicated in Fig. 7 follows Koulali et al. (2015) in indicating a single, approximate fault trace that follows the axis of the Bewani–Torricelli Mountain Range. Although few recorded earthquakes can be confidently ascribed to the BTF, Johnson and Molnar (1972) identify two strike-slip events in 1962 and 1968 consistent with sinistral motion on the BTF. A similar sense of motion was also inferred from dammed lakes and landslides reported in the aftermath of an earthquake in 1935 (Bird 1981; Dow 1977) to which Everingham (1974) assigns a magnitude of 7.0. As indicated in Table 3, we have characterized the earthquake activity of the BTF with the hybrid characteristic MFD of Youngs and Coppersmith (1985), having M_max of 7.5, a slip rate of 12 mm/year, and a seismic coupling of 0.5.

The RMF (Fig. 7) marks the collision of the Adelbert and Finisterre Arcs with northern New Guinea that began in the Late Miocene (Abbott et al. 1994; Hill and Hall 2003), accommodating convergent motion as these two arcs were accreted to the PNG mainland. We have followed Holm et al. (2016) and Koulali et al. (2015) in depicting the RMF as connecting with the New Britain Trench off the eastern tip of the Huon Peninsula, and continuing westward to the intersection of the BTF with the New Guinea Trench. It should be noted however, that the RMF’s expression on land south of the Finisterre Range is believed to have no association with the Solomon Sea Plate subducting at the New Britain Trench. The westernmost Solomon Sea Plate is folded along an WNW–ESE axis aligned with and sinking beneath mainland PNG’s Ramu–Markham Valley, with a doubly-vergent Wadati–Benioff zone at depths greater than 45 km that is distinct from the shallower seismicity of the RMF (Pegler et al. 1995).

Koulali et al. (2015) showed that convergence along the RMF is highly variable, fluctuating between 5 and 17 mm/year west of 146.5° longitude, where some share of the convergence is accommodated by the Aure Fold-and-Thrust Belt and the Southern Highlands Thrust Fault (AFTB and SHTF, respectively, in Fig. 7), to as much as 49 mm/year along the southern margin of the Huon Peninsula at 146.5°–148.0° longitude. In the fault source model we use for seismic hazard assessment, we have simplified Koulali et al. (2015) slip model for the Ramu–Markham fault to consist of three segments—western, central and eastern segments (see Fig. 2) having slip rates of 10, 17 and 47 mm/year, respectively. The RMF is expressed on land near the city of Lae in out-of-sequence thrusts, known as the Wongat and Gain Thrusts that dip northward from surface traces just north of both Lae and the Ramu–Markham Valley. Abbott et al. (1994) have determined that the Wongat Thrust accommodated a minimum of 15 km of horizontal shortening before active convergence shifted around 117 ka to the Gain Thrust, which has since accommodated 7.4 km of shortening. This agrees with Koulali et al. (2015) GNSS study of present-day motion, which found that the GPS data were best fit with convergence of about > 40 mm/year confined to the Gain Thrust.

For the central and eastern sections of the RMF, we follow the study of Stevens et al. (1998), in assigning seismicity to a steep ramp with a dip of 40° that descends from the surface trace of the Gain Thrust to a depth of 15 km, where it intersects a low-angle detachment of 12° dip extending 50 km north of the valley. Stevens et al. (1998) established this geometry from an analysis of seismic waveforms and co-seismic displacement measurements of four large (Mw 6.3 to 6.9) earthquakes that occurred in 1993, showing the largest ruptured both a steeply dipping shallow ramp and a low-angle, mid-crustal detachment, while another two ruptured the low-angle detachment only and the other had rupture confined to the steep ramp. For this detachment, we assigned a seismic coupling of 1.0 as suggested by Koulali et al. (2015) with an M_max of 8.0 since its shallow dip accommodates a potentially large seismogenic zone.

The OSF (Figs. 2 and 8) originated as megathrust associated with northward Paleogene subduction of the Australian Plate, over which the ophiolites were obducted during arc-continent collision (Daczko et al. 2011; Davies and Jaques 1984). In this study we have followed Koulali et al. (2015) and Wallace et al. (2014) in taking the OSF to stretch over 500 km along the Papuan Peninsula, from its intersection with the RMF near Lae, all the way down the Papuan Peninsula, and extending off the northern shore of the peninsula’s eastern tip. Although the Koulali et al. (2015) block model joins the three blocks WLB, SSB and TBB together as Woodlark Block, the only relative motion our model implies among these three blocks is a nominal 1 cm/year convergence at the TBB/SSB boundary. To a first approximation we can therefore apply Koulali et al. (2015) rotation pole for the Woodlark Block with respect to the PPB to all of WBL, SSB and TBB. Because this pole lies only 20 km from the OSF, with respect to PPB, TBB rotates anticlockwise in a tight arc. This would.

result in a pronounced along the strike change in both slip rate and faulting regime along the OSF. Fault motion is dominated by extension off the eastern the tip of the Papua Peninsula, is mostly sinistral in the central part of the fault, and transitions to thrust in the west as it nears the intersection with the Ramu–Markham fault (Koulali et al. 2015; Wallace et al. 2014). To simulate this we have assigned different combinations of thrust, sinistral, and oblique slip as indicated in Table 3 on the different fault segments shown in Fig. 8.

Along the northern margin of the Australian Plate, the Highlands Block (HBL) exhibits geologically recent (5 Ma, see Hill and Gleadow 1989) deformation along a fold-and-thrust belt variously referred to as the Papuan Fold Belt (Dow 1977; Hill and Gleadow 1989; Ripper and McCue 1983), Highlands Fold-and-Thrust Belt (Koulali et al. 2015; Wallace et al. 2014) and the (eastern) New Guinea Fold-and-Thrust Belt (Abers and McCaffrey 1988). The most recent deformation has occurred on the fold-and-thrust belt’s southern margin, which has also been identified as a zone of high seismicity (Abers and McCaffrey 1988; Ghasemi et al. 2016; Ripper and McCue 1983), with the largest events being a Mw 6.9 in 1922, a Mw 7.2 in 1976, and a Mw 7.5 event in 2018 (Fig. 2). We again follow Koulali et al. (2015) in ascribing these large earthquakes to a Southern Highlands Thrust Fault (SHTF in Fig. 2) which is the southernmost expression of the Highlands Fold-and-Thrust Belt. From their analysis of GNSS crustal velocities, Koulali et al. (2015) estimate about 10 mm/year of convergent motion along the SHTF, and we assume the earthquake activity to be described by a hybrid characteristic MFD (Youngs and Coppersmith 1985) with a maximum magnitude of Mw 7.5 and a seismic coupling of 0.5 (Table 3).

We also included in our source model the Weitin Fault (WF) at the eastern boundary between the NBB and the SBB where it traverses New Ireland (see Fig. 8). This is a major sinistral strike-slip fault (Lindley 2006) that has experienced very large earthquakes, in particular an event in 2000 that may have been as large as Mw 8.2 (Park and Mori 2007), although the ISC-GEM catalogue indicates Mw 8.0 for this event. Because the SSB undergoes an extremely rapid (7.8°/Myr) clockwise rotation (Wallace et al. 2004), the left-lateral motion along this fault is very high, i.e. 144 mm/year (Koulali et al. 2015). It is also associated with a high level of seismicity; hence, we have assigned a seismic coupling of 1.0 to this fault (Table 3, see also Fig. 9).

Fig. 9

Observed cumulative number of earthquakes compared with cumulative seismicity calculated for the branches of the recurrence model logic tree corresponding to 25% variations in slip rate and changes in maximum magnitude (M_max) of 0.125 and 0.250 for the (a) Ramu-Markham, and (b) Weitin Fault, respectively. The inset maps show the location of the fault surface traces (red solid curves), and for the Ramu-Markham Fault contours at 18 and 34 km depth (red dashed curves). Also plotted in the inset maps are epicentres of M ≥ 6.0 earthquakes assigned to each fault (black circles), and those not occurring on the plotted faults (grey circles)

The fault parameters used to calculate MFDs for all crustal fault and megathrust sources considered in this study are indicated in Table 3. In addition, we have allowed for epistemic uncertainty in M_max and slip rate, by considering a logic tree comprising branches with perturbations of ± 25% in slip rate and ± 0.25 units in M_max for all sources except the RMF, which has only ± 0.125 unit variation in M_max (since it has a historical event of magnitude 8.0, we took this as the lowest possible M_max, and allowed for M_max as high as 8.25). The logic tree weights assigned were 0.25 for the mean model whose parameters are given in Table 3, 0.125 for the four branches with a perturbation in either slip rate or Mmax, and 0.0625 for the four branches with perturbation in both slip rate and M_max. In Fig. 9, we compare the observed cumulative number of earthquakes for the catalogue with those calculates for different logic tree branches for the RMF and WF. Note, however, that in the PSHA calculation the logic tree was collapsed in order to achieve better computational efficiency. As can be seen in Fig. 9, the observed cumulative seismicity for the RMF and WF falls within the range spanned by the logic tree branches.

4.2 Distributed sources

We modelled shallow crustal earthquakes with depth lower than 40.0 km and magnitude lower than 6.0 by a spatially distributed earthquake source model. Note that these events are not considered in the “fault sources” model. We followed Frankel (1995) approach to develop and characterise the distributed sources model. Here, seismic sources are modelled by a set of point sources located on a regular grid with the spacing of 0.1° × 0.1°. Each point source generates earthquakes with a focal depth distribution following the probability mass function with values of 0.25, 0.5, and 0.25 for depths at 10.0 km, 20.0 km, and 30.0 km, respectively. The focal mechanisms of the point sources vary spatially, reflecting the representative mechanism of each of the area source zones as assigned by Ghasemi et al. (2016). Each point source follows a Gutenberg–Richter recurrence law to produce earthquakes with magnitudes in the range of Mw 5.0–6.0. For those sources within the Bismarck Seismic Sea Lineation and Aure Fold-and-Thrust Belt zones (Figs. 2, 7), we increased the magnitude upper bound to 7.2 as these zones of distributed faulting were not included in the “fault sources” model. The upper bound limit is set considering the historical seismicity within the aforementioned zones.

To estimate the Gutenberg–Richter recurrence parameters (a- and b-values), we used a subset of the declustered catalogue that includes events classified as shallow crustal, i.e. events with depth ≤ 40 km, and shortest distance to the New Britain and New Guinea subduction interfaces more than 20 km. Figure 10 shows the magnitude–frequency distribution (MFD) of shallow crustal earthquakes in the catalogue. Based on the observed MFD, recurrence parameters were estimated using the maximum likelihood method of Weichert (1980) with a 0.2 magnitude-unit interval. It can be seen that the data closely follow the Gutenberg–Richter distribution with a b-value of 1.023 ± 0.024. To explore the spatial variation of b-value, we estimated recurrence parameters for individual area-source zones defined by Ghasemi et al. (2016). The results suggest that a fixed b-value of 1.023, the value calculated for the entire region, can be used for the distributed point sources, as the estimated b-values for 16 (out of 18) area zones did not differ significantly from the one calculated for the entire region.

Fig. 10

Magnitude–frequency distribution of shallow crustal earthquakes in the PNG region. Solid and dashed lines represent the fitted Gutenberg–Richter model

To calculate the annual rate of earthquakes greater than M_W = 0.0 (from Eq. 1, $N_{0} = 10^{a}$), first, we defined a grid with spacing of 0.1° in latitude and longitude and point sources located at the centre of the grid cells. The observed annual rate of earthquakes are then calculated based on number of events within each grid cell and using a Gutenberg–Richter distribution with a b-value of 1.023. We also considered temporal variation in completeness by normalizing the observed rate in each magnitude bin according to the completeness duration of the specific bin, as described in Woo (1996). Finally, for each grid cell the observed annual rate of earthquakes is smoothed using an isotropic Gaussian kernel with fixed correlation distance of 30 km. This subjective value is justified as the resulting smoothed seismicity maps based on this value better represent the general seismicity patterns across the region. Figure 11 shows a map of the smoothed Gutenberg–Richter a-values across the PNG region.

Fig. 11

Spatially smoothed annual rates of shallow crustal events in the region assuming a fixed Gutenberg–Richter b-value of 1.023

4.3 Intraslab sources

The complex subduction environment that surrounds Papua New Guinea, both to the west and the east, means that the contribution of intraslab earthquakes to the seismic hazard in the archipelago is significant. Given the complexity of the slab geometries in the region, adopting a smoothed seismicity model similar to that used for shallow distributed sources may mix earthquakes from different slabs, so a more refined approach to the characterisation of intraslab seismicity is needed. In this study, the intraslab seismic sources are represented by a set of finite ruptures within the continuous slab volume following the methodology described in Weatherill et al. (2017).

The slabs extending from both the New Guinea and the New Britain trenches are initially defined using the Slab 2.0 dataset (Hayes et al. 2018; Fig. 12). However, the poorly defined deeper seismogenic structure beneath Papua New Guinea that appears to connect the two subduction systems is instead modelled preferably as a deep volumetric source. Here no interface is determined, however, an intraslab geometry is defined separately from the New Britain trench and constructed in the same manner as that of the main subduction intraslab sources described in the following.

Fig. 12

Geographical extents of the slab surfaces used in the construction of the intraslab sources represented by contours for the New Guinea subduction (west) and New Britain subduction (east). The original Slab 2.0 (Hayes et al. 2018) slab extents are shaded

4.3.1 Intraslab source geometry

As many modern subduction intraslab ground-motion models require representations of the earthquake rupture as a finite fault, it is necessary to devise a method that can allow the construction of such faults whose area and rupture aspect ratio scale with magnitude, albeit respecting the seismogenic thickness of the slab. Such ruptures should not propagate into the mantle above or below the slab limits. It is also recognised that the orientation of the fault rupture will change within the slab depth as the driving stress transition from those associated with crustal flexure and compression at shallower depths to extension from gravitational slab pull further in the Wadati–Benioff zone (Astiz et al. 1988).

The intraslab geometries is constructed as follows. Initially the Slab 2.0 geometry is used to represent the upper surface of the slab volume, and this is rendered into a regularly spaced mesh. A lower surface is then defined by projecting the upper surface downward along a normal vector by the slab thickness, which for the current study is assumed to be approximately 40 km. Within these surfaces multiple layers of points are rendered parallel to the two surfaces, which form the centroids of the ruptures, or the hypocentres if the ground-motion model requires hypocentral distances. The orientations of ruptures are specified for each user-defined depth interval in terms of a relative azimuth and inclination with respect to the corresponding point on the upper surface. One is therefore able to control the dip, orientation and rake of the ruptures according to depth, thus better reproducing the variations in focal mechanism typically found within Wadati–Benioff zones.

For each earthquake the corresponding upper and lower seismogenic depths are defined based on the projected intersection of the ruptures with the upper and lower seismogenic surfaces of the slab, which ensures the ruptures are constrained to the slab volume in the vertical dimension. A single magnitude–frequency distribution is defined for the entire slab volume and the activity rates of the slab distributed evenly among the individual point sources. Finally a magnitude to rupture-area scaling relation and prior aspect ratio are defined, which are then used inside the PSHA calculation to generate finite ruptures for each source according to the magnitude. These finite ruptures are constrained within the slab volume, and allowed to propagate horizontally exceeding their prior assumed aspect ratio once the down-dip width of the ruptures reach the seismogenic thickness. For application to the New Guinea and New Britain subductions, the magnitude to rupture-area scaling relation of Strasser et al. (2010) is adopted, and prior aspect ratios set to 1.0. At depths shallower than 100 km the ruptures are inclined at 20˚ with respect to the local inclination of the upper surface but with strike aligned parallel to the subduction zone. The style of faulting in the upper 40 km is set to reverse, corresponding to flexure or compressive stresses in the shallowest part of the plate (e.g. Astiz et al. 1988), while for depths between 40 and 80 km a mixture of normal and reverse faulting is considered, and for depths greater than 80 km all ruptures are considered to be normal. None of the GMMs selected in the subsequent section explicitly consider style-of-faulting as a parameter for intraslab earthquakes, so this latter choice will not impact on the hazard and is assigned purely for consistency with the regional tectonics.

4.3.2 Intraslab recurrence

For each of the subducting slab volumes described previously, with upper surfaces illustrated in Fig. 12, a corresponding earthquake sub-catalogue is extracted from the compiled regional catalogue. As only a relatively minor proportion of the earthquakes in the regional catalogue contain moment tensors, identification of intraslab earthquakes is done primarily on the basis of proximity of the hypocentre to the slab volume. Earthquakes whose epicentres fall within the surface projection (plus a 20 km buffer) of the slab geometry, and depths within the upper and lower surfaces, are retained for the calculation of the intraslab magnitude–frequency distribution.

Using the regional estimates of completeness, a double-truncated Gutenberg–Richter model is fit to the observed seismicity in each intraslab source using the Weichert (1980) approach. Both zones are found to have similar b-values of 1.238 ± 0.11 for the New Guinea slab, and 1.251 ± 0.04 for the New Britain slab, with rates above M_W 5.0 of 2.26 ± 0.156 per year and 18.5 ± 0.155 per year respectively. Joint optimisation of the fit of the truncated Gutenberg–Richter magnitude–frequency distribution to the observed completeness-adjusted activity rates indicates that the best fitting upper bounds are M_W 7.7 for the New Guinea trench and M_W 8.0 for the New Britain trench, which appears to be consistent with other intraslab source models as well as the largest observed intraslab earthquakes around the globe (Astiz et al. 1988; Strasser et al. 2010). As it must be acknowledged that the observed number of earthquakes above M_W 7.0 is generally small in both in-slab zones, there is little constraint on the shape of the magnitude frequency distribution at larger magnitudes. The Bayesian M_MAX estimator of Kijko (2004) is therefore also applied in order to provide an independent constraint using the fitted b-values and their respective uncertainties. For the New Guinea slab, with a maximum observed magnitude of M_W 6.9, this produces an estimate of M_MAX 7.22 $\pm$ 0.34, whilst for the New Britain slab, whose maximum observed magnitude is M_W 7.3, the resulting M_MAX is 7.43  $\pm$  0.12. In spite of the large number of events used to constrain the b-value, we find the Bayesian M_MAX estimators underestimate the rates of events for magnitudes close to the maximum observed due to the tapering of the truncated Gutenberg-Richter model. In addition, such values would appear decidedly non-conservative given the physical size of the slab and the precedents for larger in-slab earthquakes elsewhere in the world. We have therefore retained the previous estimates of M_W 7.7 and M_W 8.0 for application in the PSHA to the New Guinea and New Britain slabs respectively.

5 Ground-motion models

The selection of GMMs for use in PSHA is one of the most challenging aspects of model development in regions where little or no strong motion data is available. Given difficulties in characterising appropriate GMMs for a range of tectonic region types in regions with few instrumental records, it is common-practice to select multiple GMMs within a weighted logic-tree to capture the epistemic uncertainty in ground motions. In this study we use pre-selection criteria to identify candidate models (Bommer et al. 2010; Cotton et al. 2006); then, we can explore the resulting distribution of median ground-motion values to identify GMMs that we believe cover the centre and range of epistemic uncertainty. This last point is assessed in two ways: by trellis plots that capture the distribution of GMM scaling with magnitude, distance, and spectral period, and by using a more information-theoretic perspective using Sammon’s mapping to interpret the model space occupied by a suite of candidate GMMs. For more detailed information on the theory and interpretation of Sammon’s maps the reader is referred to Scherbaum et al. (2010) and its references.

In the following we present in more detail the selection and ranking of GMMs for active shallow crustal earthquakes. For other tectonic region types considered in this study, i.e. subduction interface and subduction intraslab earthquakes, we followed a similar approach and the results are presented in Ghasemi et al. (2020).

Shallow crustal earthquakes in tectonically active regions are well-represented by the available GMMs in the scientific literature. With the development of models from extensive databases of strong motions (e.g. Akkar et al. 2014, Ancheta et al. 2014), a greater number of candidate models can be identified that pass the pre-selection criteria. The tectonics of shallow seismicity in the region described above highlight the need to consider GMMs that can apply to different faulting regimes, and that can scale appropriately in the near-fault region. Searching the database of ground motion models available in the scientific literature we find a total of 29 candidate GMMs; however, this number includes cases of multiple models from the same literature source to account for regionalisation as well as a number of calibrations of global (i.e., ergodic) models to the local source scaling or attenuation properties of some regions (i.e., partially non-ergodic). Since we know very little about the attenuation or stress-drop scaling of seismicity in Papua New Guinea, we cannot necessarily exclude the possibility that such characteristics may be analogous to other regions of the globe for which global GMMs have been calibrated; hence we consider both the original model and its local calibration.

For the Sammon’s maps we define the high dimensional space from scenario ground motions for magnitudes 4.5 ≤ M_W ≤ 8.0 in increments of 0.25 magnitude units, and twelve distances in the range 15 ≤ R_JB (km) ≤ 250. Calculated ground-motions used to generate the Sammon’s maps assume a vertically dipping strike-slip scenario to minimise the influence of near-fault hanging wall effects in the models and differences in finite rupture distances. The same general trends were observed for a dipping reverse-faulting scenario but are omitted here for brevity. As the hazard model is constructed for the rock site class, the reference site condition in these comparisons is fixed at V_S30 = 760 m/s.

From both the trellis plots and the Sammon’s maps (Figs. 13, 14, respectively) a greater degree of convergence can be seen in the intermediate and longer-period spectral accelerations (i.e., 0.5 ≤ T (s) ≤ 5) than at shorter periods. The Sammon’s maps here are revealing, as they indicate a slight separation between two groups of models: one that includes Campbell and Bozorgnia (2014), McVerry et al. (2006), and Zhao et al. (2016b) (highlighted for the PGA and Sa (0.2 s) case via the ellipses shown in Fig. 14), and one that encompasses the remaining models. Given the position of Zhao et al. (2016b) in the trellis plots, these may be representative of models showing higher short-period accelerations, which may be a consequence of the implicit stress-drop scaling within the model and/or the local site characteristics.

Fig. 13

Trellis plots showing the scaling of the candidate active shallow crustal GMMs (grey lines) with period for three different magnitude scenarios: MW 5.0, 6.5, 7.5 (top, middle, and bottom rows respectively) and three different distances: RJB 20 km, 50 km

Fig. 14

Sammon’s maps for active shallow crustal GMMs for PGA (top left); Sa (0.2 s) (top right); Sa (1.0 s) (bottom left); and Sa (2.0 s) (bottom right). The separated group of GMMs for short periods (PGA and Sa (0.2 s)) are marked by the red ellipses

For the active shallow crustal environment, three ground motion models are selected for inclusion in the logic tree: Cauzzi et al. (2015), Chiou and Youngs (2014), and Zhao et al. (2016b). The latter represents the case in which the seismological characteristics of shallow crustal earthquakes in PNG are more directly analogous to those of Japan, from which the data set of Zhao et al. (2016b) is taken. This corresponds to the higher short-period ground-motion end of the epistemic range, which converges with the general body of the models at higher magnitudes. The Zhao et al. (2016b) model is assigned a weight of 0.4. The remaining body of the distribution is well represented by Chiou and Youngs (2014) and Cauzzi et al. (2015). Though, both models are fit to global databases of strong motions, Cauzzi et al. (2015) is significantly skewed toward Japan. In spite of this it is found to diverge from Zhao et al. (2016b) at short periods, which may reflect some differences the influence of the different site characterisation and modelling assumptions between the two models. Despite the overlap between Cauzzi et al. (2015) and Cauzzi et al. (2015) in terms of the source region, we believe that the contrasting modelling assumptions and parameterisations represent an appropriate degree of epistemic uncertainty in the ground motion. Therefore, both models are retained in the logic tree.

All Three models also differ in the characterisation of near-fault hanging wall effects, with Cauzzi et al. (2015) adopting rupture distance (R_RUP) as the preferred metric but with no explicit hanging wall amplification, while Chiou and Youngs (2014) include an explicit hanging wall term. Both models are broadband, meaning that (with the inclusion of Zhao et al. 2016b) seismic hazard over a spectral range of 0.01 to 5 s can be constrained. The weightings for Chiou and Youngs (2014) and Cauzzi et al. (2015) are split evenly with each assigned a weight of 0.3.

The influence of volcanism in the ground motion attenuation is pertinent for Papua New Guinea given the presence of many active volcanoes in New Britain and northern Papua New Guinea. Of the selected GMMs, Zhao et al. (2016b) explicitly accounts for the influence of the travel path through regions of low Q associated to volcanic zones via the use of the volcanic distance term, R_VOLC. In application, however, this term is difficult to utilize, requiring not only a prior definition of the low Q zones but also the computational tracing of the travel path through the zones; something rarely implemented in most PSHA softwares including OpenQuake Engine. In contrast to Japan, however, in which the belt of active volcanism is located well onshore, bisecting the country down its central mountain belt, in Papua New Guinea and New Britain however, Quaternary volcanism is mostly constrained to offshore islands and the extreme northern coastlines. In the case of earthquakes from shallow onshore sources, we would expect a negligible influence from volcanism on resulting ground shaking. For the northward dipping New Britain trench, the travel path for in-slab events affecting the northern coast of New Britain is predominantly through subducting oceanic crust rather than low attenuation mantle. For the New Guinea slab, although the subducting slab dips below the western Papua New Guinea directly there is a notable lack over Quaternary volcanism in the region. Given the tectonic and geological complexity of the region, combined with an absence of data to verify the influence of volcanism in seismic attenuation, it is decided to effectively turn off the volcanic distance term of Zhao et al. (2016b) (i.e. set R_VOLC = 0 km), rather than unduly attenuate the strong shaking too strongly. This is naturally a conservative assumption, but we believe it is appropriate for the application to seismic design.

Following a similar use of trellis plots and Sammon’s mapping for exploration of the model space of candidate GMMs, Ghasemi et al. (2020) show how to select and weight subsets of potential GMMs for the other tectonic region types considered in this study, corresponding to subduction interface and subduction intraslab earthquakes. This allows the complete ground-motion characterisation model for the PSHA19 to be defined as depicted in Fig. 15. While there is limited (to no) data from which to undertake quantitative assessments of the relative performance of individual GMMs, this logic tree aims to capture the range of plausible ground motions that could be observed from future large earthquakes in the PNG region.

Fig. 15

Ground-motion logic tree for application in Papua New Guinea considering the three tectonic region types: active shallow crust, subduction interface and subduction intraslab

6 Hazard results

For a specific site, PSHA combines the modelled uncertainties in the earthquake location, magnitude, and ground-motion level to derive the seismic hazard curve, i.e. the probability of exceeding given ground-motion level over a specified time period. We computed hazard curves for a grid of sites with spacing of 0.1° in latitude and longitude. The calculations were for bedrock sites with average shear-wave velocity in the upper 30 m of 760 m/s. For each site, a series of hazard curves were computed for spectral accelerations of periods in the range of 0.0 to 4.0 s. From the computed hazard curves, we generated a series of seismic hazard maps for the selected spectral acceleration ordinates for a 10% and 2% probability of exceedance in 50 years. These exceedance probabilities refer to 1/475 and 1/2475 annual exceedance probability (AEP), respectively. Figure 16, as an example, shows the mean PGA hazard map for PGA and for 10% probability of exceedance in 50 years. Overall, it can be seen that the level of seismic hazard is mainly controlled by the modelled fault sources representing the key tectonic boundaries in the input source model (Fig. 2). It reaches the highest values in the coastal areas of the Huon Peninsula and New Britain/Bougainville regions. The Huon Peninsula accommodates high level of seismicity associated with the Ramu–Markham Fault zone and the deeper Solomon Sea plate (Pegler et al. 1995). Similarly, the New Britain/Bougainville region is the host of a significant number of moderate-to-large earthquakes linked with the New Britain Trench. As in Fig. 1, the lowest level of hazard is calculated in the southern part of the New Guinea Highlands Block. This block is assumed to be part of the stable Australian plate experiencing a low level of local seismicity.

Fig. 16

Spatial distribution of mean PGA values for 10% probability of exceedance in 50 years. The location of major cities of PNG are shown by black circles

There is significant uncertainty associated with the mean values presented in Fig. 16. Figure 17 shows the PGA hazard maps derived from 15 and 85th percentile hazard curves. Notably, the PGA level varies in the range of 0.5–1.6 g in the high hazard part of the New Britain/Bougainville coastal region. Note that Fig. 17 only depicts the epistemic uncertainty of PGA estimates due to variability in the GMMs as captured by selecting multiple GMMs in a logic-tree framework. Such uncertainty is further illustrated in Fig. 18, which shows the ratio of the PGA hazard of the Britain/Bougainville coastal region calculated using the Zhao et al. (2016a) GMM relative to the Atkinson and Boore (2003) GMM. It can be seen that the Zhao et al. (2016a) GMM results in significantly higher hazard across the region with factors up to 75%. This comparison shows the variability among the selected GMMs and the vital need to develop a database of recorded ground-motions in PNG. Such a database could be used to develop PNG-specific GMMs and also guide the process of selection and ranking of GMMs for future PSHAs.

Fig. 17

1/475 AEP PGA hazard maps corresponding to the 15th percentile (left) and 85th percentile (right).

Fig. 18

The ratio of 10% probability of exceedance in 50-year PGA using the subduction interface model of Zhao et al (2016b) and Atkinson and Boore (2003)

To explore the type, size, and spatial location of earthquakes with significant contributions to the seismic hazard, hazard deaggregation analysis was conducted for 24 major cities in PNG. While beyond the scope of this manuscript, such analyses assist with identifying design earthquakes for major population centres in PNG. The full set of results are not shown here, however they generally identify the nearby “fault sources” as the key contributing sources to the calculated level of seismic hazard. For example, the deaggregation analysis for Mendi (Fig. 19) suggests that large earthquakes with magnitude modal value of Mw 7.3 occurring on the Southern Highlands Thrust Fault (SHTF) are contributing the most to the hazard. As expected, the SHTF in our source model is the closest active tectonic feature to Mendi. It is also considered as the causative fault of the recent 2018 Mw 7.4 earthquake that caused significant loss of lives and damage in the Highlands region, including Mendi.

Fig. 19

Seismic hazard deaggregation for Mendi for 10% probability of exceedance in 50-year PGA. The solid line shows the Southern Highlands Thrust Fault trace

Note that the results do not rule out a hazard contribution from other source types, such as distributed and intraslab sources. For example, the deaggregation of seismic hazard for Bulolo City (Fig. 16) indicates a considerable contribution (20%) from intraslab sources to the overall level of seismic hazard at this location. The geographic deaggregation of intraslab sources shows that mainly large earthquakes with magnitude and distance modal values of 7.0 and 105, respectively control the contribution from intraslab sources to the seismic hazard at Bulolo City. These events are located on the north–northwest side of Bulolo City occurring in the west segment of the New Britain slab. The size and location of the controlling earthquakes coincides reasonably well with those of the recent, 2018 MW 7.2 event that post-dated the catalogue used for this study. Such examples demonstrate the value of deaggregation analysis for identifying controlling earthquakes and developing realistic earthquake scenarios for preparedness and impact mitigation activities in PNG.

In comparison with the Ghasemi et al. (2016), generally the results of this study are higher at the location and vicinity of the fault sources modelled in this study (Fig. 2, and Section 0), and lower elsewhere. Such discrepancies reflect the differences of these studies in the modelling of the earthquake sources in the region as well as selection of ground-motion models. As an example, Fig. 20 compares the computed uniform hazard spectra for Lae with that from Ghasemi et al. (2016). In addition, the uniform hazard spectrum calculated using the source model developed in this study in combination with the ground-motion models from Ghasemi et al. (2016) is also displayed. It can be seen that at periods longer than ~ 0.5 s, the differences between this study and Ghasemi et al. (2016) are mainly controlled by the selection of different ground-motion models. However, at shorter periods (T <  ~ 0.5 s) differences in the earthquake source models become the key factor in explaining the discrepancies among the mentioned studies. To further clarify this point, we carried out the disaggregation analysis for Lae at period of PGA using earthquake source models developed in this study and Ghasemi et al. (2016). To remove the impact of using different ground-motion models on the results, we used both source models in combination with ground-motion models selected in this study to perform the disaggregation analysis. The results (not presented here) indicated that, unlike this study, in Ghasemi et al. (2016) there is a significant contribution from deep earthquake sources linked to the New Britain subduction. It should be noted that in 2016 study, we modelled inslab events as different layers of point sources, however in this study, as discussed before, we implemented a more refined approach to the characterisation of inslab seismicity in the region.

Fig. 20

Comparison of the computed uniform hazard spectrum for Lae (labelled as PNG19) with that from Ghasemi et al. (2016) (labelled as PNG16). The uniform hazard spectrum calculated using the source model developed in this study in combination with the ground-motion models selected for PNG16 model is also displayed

It is difficult to compare our results with those of Jury et al. (1982), on which the current PNG building code is based (Fig. 1), because the methods and presentation of the results are very different. Such comparison would be meaningful only by assuming that the presented 1/475 AEP PGA values, similar to the Jury et al. (1982) results, represent seismic design levels for building construction in PNG. If this assumption holds, Fig. 21 shows the comparison between the seismic design levels calculated in this study versus those of Jury et al. (1982). The results are displayed as the ratio between mean 1/475 AEP PGA values and those of Jury et al. (1982). It can be seen that the design levels of the current building code vary significantly, depending on the location, relative to the present study. The increase of the hazard level by a factor of 143% at Lae City in the Huon Peninsula region is of particular interest to this study since Lae is PNG’s industrial capital and hosts considerable human and economic activity.

Fig. 21

The ratio of the results of this study for 10% probability of exceedance in 50-year PGA relative to the assumed design PGA levels of the current building code of PNG

In Jury et al. (1982), Lae is located within Zone II with PGA design level of 0.3 g, considerably lower than that of the New Britain/Bougainville region (Zone I, Fig. 1). By contrast, the results of this study indicate that the PGA level of 0.73 g at Lae that is comparable with the high level of hazard in the New Britain/Bougainville region. Note that the seismic zoning map of Jury et al. (1982) relies merely on the distribution of large earthquakes (M > 6.0) in the period 1900 to 1978, justifying the geometries of the zones represented in Fig. 1. Naturally in this study, in comparison with Jury et al. (1982), we compiled a more comprehensive catalogue with a lower-magnitude threshold to model earthquake occurrence rates in the region. We also included fault sources in the hazard model to better constrain the hazard level at the vicinity of active tectonic features. The deaggregation results for Lae City (Fig. 22) suggests that large earthquakes (with modal values of magnitude Mw 8.0, and distance 15 km) mainly contribute to the estimated high level of 1/475 AEP PGA. The controlling earthquakes are located on the north side of the city and occur on the Ramu–Markham Fault (RMF) system. The occurrence of such large earthquakes on this fault is quite plausible given the results of the recent GNSS studies and records of historical earthquakes. The recent GNSS studies indicate that a significant portion of the RMF appears to be locked and accommodating up to 61 mm/year of convergence between the New Guinea Highlands and the Australian plates (Koulali et al. 2015; Wallace et al. 2014). Furthermore, as previously mentioned, there are at least five historical earthquakes in the late 19th and early twentieth century on the RMF with magnitudes approximately in the range 7.0 to 8.0.

Fig. 22

Seismic hazard deaggregation for Lae for 10% probability of exceedance in 50-year PGA. The solid line shows the Ramu–Markham Fault trace

7 Conclusion

A PSHA was performed to underpin a revision of the earthquake loading component of the Building Codes of PNG, which is presently based on the seismic zoning map of Jury et al. (1982). The results of this study indicate significant differences between the updated national seismic hazard assessment with that used in the Building Codes of PNG. The updated assessment offers many important advances over its predecessor. It is based on a modern probabilistic hazard framework and considers an earthquake catalogue augmented with an additional four decades-worth of data. It also considers advances in ground-motion modelling through the use of multiple GMMs. This has enabled a more comprehensive seismic hazard assessment and its uncertainties. In particular, for the region of interest, we compiled a homogeneous earthquake catalogue based on M_W covering the period 1900–2017. We also applied qualitative and quantitative tests to verify the methods used for processing the earthquake catalogue. In the hazard input model, the individual fault sources are implemented representing active major and microplate boundaries to better constrain the expected levels and patterns of ground motion for any given site. In this model, the intraslab sources are also represented in a realistic way by using the continuous slab volume to constrain the finite ruptures of such events.

We suggest that future improvements in the national seismic hazard assessment of PNG can derive from several sources. For the short term, these are: (1) an improved earthquake catalogue, augmented with both more recent earthquakes and more information for historical events, and a better classification scheme (e.g. distinguishing megathrust from intraslab events); (2) accounting for additional epistemic uncertainty associated with alternative seismotectonic models and scaling relations; (3) more robust selection of GMMs based on data from PNG earthquakes.

In the longer term, further improvements in seismic hazard assessment in PNG will rely on a fully functional seismic network, which is expected to yield a significant number of locally recorded earthquakes that will improve the delineation of seismic zones and the accuracy of the Gutenberg–Richter MFDs. While in the short term data from this network will be crucial for compiling a dataset of ground motions to enable the selection and ranking of GMMs, in the long term these data will be critical for developing PNG-specific GMMs to reduce the epistemic uncertainties associated with future national hazard assessments. However, the maintenance of seismic networks in PNG faces several challenges including: adverse environmental factors such as heat and humidity; unreliable/intermittent communications networks and; the need to provide a high level of security for network infrastructure.

The tectonic model of the region is quite complex and not definitive, with even the number of micro-plates uncertain, not to mention their boundaries and Euler vectors. Several of the components of our seismotectonic model, in particular the SHTF and the WF, may involve deformation over a network of faults rather than the simplified single-fault models we have presented here. In the long term it will be essential to continue and improve on previous GPS campaign measurements, and these efforts would benefit from the establishment of a network of permanent, continuous GPS (cGPS) stations. This is especially important in places like the Ramu-Markham Fault and the Southern Highlands Thrust, where high earthquake activity is proximate to important infrastructure and populations centres.