Traditional seismic hazard analyses underestimate hazard levels when compared to observations from the 2023 Kahramanmaras earthquakes

On 6 February 2023, two large earthquakes (Mw 7.8 near Pazarcik and Mw 7.5 near Elbistan) ruptured multiple segments within the East Anatolian Fault Zone (EAFZ) and neighboring structures in Southeastern Turkey, a tectonically complex region at the junction of the Anatolian, Arabian, and African plates. Although instrumental-era major events were absent locally, historical evidence and active faulting indicated high seismic potential and a seismic gap near Kahramanmaras. Prior studies generally anticipated lower maximum magnitudes and treated key segments independently, leaving questions about the rarity of the observed multi-segment ruptures and the extremity of recorded ground motions. This study aims to provide a quantitative, non-speculative assessment of how rare the observed magnitudes and ground motions were by estimating exceedance probabilities via probabilistic seismic hazard assessment (PSHA) at recording sites and by comparing observations with ground motion models (GMMs).

Multiple works highlighted high seismic potential and a seismic gap at the southern EAFZ: historical large events (e.g., 1513 Turkoglu, 1822 Aleppo, 1872 Amik Lake, 1893 Malatya) and mapping/segmentation of EAFZ and DSFZ (Duman & Emre; Emre et al.). PSHA and source characterizations considered segment-specific and some multi-segment scenarios (e.g., Gülerce et al. considered multi-segment ruptures of Erkenek and Pazarcik and of both Surgu segments). GPS-derived slip rates constrained seismic potential (Aktuğ et al.), while recent seismotectonic analyses (Güvercin et al.) estimated maximum magnitudes and return periods for EAFZ segments (e.g., Pazarcik ~Mw 7.3, RP ~700 yr; Amanos ~Mw 7.4, RP ~900 yr). Most prior studies did not anticipate combined Amanos+Pazarcik rupture. Standard PSHA methodology and GMM development/assumptions are summarized in Baker et al. and related references, with ergodic vs non-ergodic modeling discussed in recent literature.

Recorded motions: Raw AFAD strong-motion accelerograms for the two mainshocks were curated from Gülerce et al., removing incomplete/noisy/unclear traces and processed per Akkar et al. Records within 100 km of the finite-fault plane and with known VS30 were retained, yielding 51 records for Mw7.8 and 25 for Mw7.5. Intensity measures computed were PGA and 5% damped RotD50 spectral accelerations at T=0.5 s and T=1.0 s; RotD50 was obtained by rotating orthogonal components and taking the median across orientations. PSHA framework: Using OpenQuake, site-specific hazard curves were computed for PGA, SA(0.5 s), and SA(1.0 s) at each recording station, based on the European Seismic Hazard Model 2020 (ESHM20) inputs. The hazard integral sums contributions across sources and ruptures; conditional exceedance probabilities use lognormal GMM residual assumptions. ESHM20 source model includes area and fault sources with a logic tree to capture epistemic uncertainties. For area sources, two MFD branches: double-truncated Gutenberg–Richter (TGR) with a/b values at 5th/50th/95th percentiles (weights 0.2/0.6/0.2) and a tapered Pareto (weight 0.4). Maximum magnitudes for area sources are set from catalog maxima plus 1σ (lower), +0.3 (median), +0.6 (upper). Fault-source Mmax is from Leonard’s scaling; slip rates (from EFSM20) feed activity rates via moment conservation; b-value is regional. Sources within 200 km of each site were included. The EAFZ is represented as a single long fault source (TRCF002); Surgu/Cardak are represented by TRCF03H, with area sources TRAS481, LBAS341, and TRAS434 overlapping the observed ruptures. Ground motion characterization: ESHM20 uses a scaled backbone approach with the partially non-ergodic Kotha et al. GMM to generate hazard curves. For residual analyses (period-dependent normalized residual ε), the NGA-West2 GMM of Boore et al. (BSSA14) was used because it includes events up to Mw 7.9, better matching the mainshocks. Residual and joint-exceedance analyses: Residuals ε over 0.01–10 s were computed relative to BSSA14, with ε interpreted via the standard normal exceedance probability. For broadband extremity assessment, joint exceedance probabilities at selected pairs of periods were approximated using the bivariate normal distribution with inter-period correlations from Baker & Jayaram (e.g., ρ for T=0.3–3.0 s and for T=0.5–1.0 s). Mean return periods were inferred as reciprocals of annual exceedance rates on site-specific hazard curves, with logic-tree median and 16th/84th ranges.

Magnitude recurrence: For the Mw7.8 Pazarcik rupture, the fault source TRCF002 implies a mean return period (MRP) ~1000 years; overlapping area sources TRAS481 and LBAS341 imply MRPs ~4000 and ~3500 years, respectively. For the Mw7.5 Elbistan rupture, the fault source TRCF03H gives MRP ~7800 years, whereas area source TRAS434 gives MRP ~2600 years. Ground motion exceedance: For Mw7.8, maximum MRPs across stations for exceeding observed motions are approximately 6900 years (PGA), 4100 years (SA 0.5 s), and 6900 years (SA 1.0 s). For Mw7.5, maxima are ~1200 years (PGA), ~1700 years (SA 0.5 s), and ~6900 years (SA 1.0 s). Station-specific highlights: station 3125, Antakya, observed PGA 0.94 g (illustrated hazard curve); station 4612 (Goksun) recorded PGA 0.58 g, SA0.5 1.02 g, SA1.0 0.61 g with MRPs ~750, 830, and 910 years; station 4406 (Malatya) had the highest PGA MRP ~1250 years for PGA 0.44 g; station 3802 (Kayseri) had the highest SA MRPs for SA0.5=0.39 g (1600 years) and SA1.0=0.38 g (6900 years). Residuals and extremity: For Mw7.8, overall residuals show no strong trend, but some sites (e.g., 3126, 3129, 3135) have very large positive ε across periods; station 3135 reached ε≈3 at T=0.3 s (exceedance probability ~0.1%). For Mw7.5, residuals trend upward with period, indicating underestimation at longer periods by BSSA14. Joint exceedance at station 3135 indicates extreme broadband motion: Pr(ε0.3≥2.96, ε3.0≥1.21 | ρ=0.2535)≈0.0005 (~1/2000); Pr(ε0.5≥2.19, ε1.0≥1.70 | ρ=0.749)≈0.0086 (~1/116). Comparisons with NGA-West2 show some 6 Feb 2023 records meet or exceed historical maxima: SA1.5 at station 3138 is 1.20 g (exceeding NGA-West2 max 1.17 g); SA1.0 at station 3139 is 1.12 g (exceeding NGA-West2 max 1.08 g). Highest amplitudes and MRPs cluster in Hatay for Mw7.8, consistent with forward-directivity and dense station coverage.

Quantitative recurrence estimates indicate that a Mw7.8 event on the EAFZ is not extraordinarily rare under a fault-source model (~1000-year MRP), whereas area-source models imply longer MRPs; for the Surgu/Cardak structures the opposite holds, likely reflecting greater epistemic uncertainty in geometry and rates. Ground motions, especially in Hatay during Mw7.8, reached levels corresponding to MRPs up to ~7000 years (≈0.7% probability in 50 years), well beyond standard building design targets (475-year MRP), underscoring that observed demands can greatly exceed code levels locally. While Pazarcik motions generally align with GMMs, several sites showed broadband, extreme residuals; Elbistan residuals increase with period, suggesting underestimation at long periods by current ergodic models. These findings support that traditional seismic hazard analyses, constrained by segment-level magnitude limits and ergodic GMMs developed without abundant near-fault large-event data, can underestimate hazard levels, particularly for multi-segment ruptures and long-period motions. The study highlights paths forward: physics-based ground motion modeling for large, near-fault events to capture rupture geometry, directivity, rupture velocity, and basin effects; development and use of non-ergodic/partially non-ergodic GMMs leveraging growing data to account for repeatable source, path, and site effects; and integration of physics-based earthquake sequence simulators to represent triggering and clustering, not captured in standard PSHA. The sensitivity of MRP estimates to logic-tree branches and the practical difficulty of validating very long return periods further emphasize the need to quantify epistemic uncertainties carefully.

By performing site-specific PSHA with ESHM20 inputs and comparing near-source recordings from the 6 February 2023 Mw7.8 and Mw7.5 earthquakes to empirical GMMs, the study quantifies exceedance probabilities for the observed magnitudes and ground motions. Results show that while the seismic potential of the region was recognized, traditional analyses tended to underestimate hazard levels for combined, multi-segment ruptures and at longer periods, with several recorded spectra among the most extreme ever observed. The work underscores that some observed motions correspond to multi-millennial return periods, far exceeding standard design levels. Future research should: (1) incorporate physics-based rupture and wave propagation models into hazard and risk assessments for large, near-fault events; (2) expand and utilize non-ergodic GMMs with constraints for near-source magnitude scaling; (3) integrate physics-based earthquake sequence simulators to treat triggering and clustering; and (4) refine source models and logic trees with improved fault characterization and richer datasets, reducing epistemic uncertainty and improving robustness of hazard estimates.

High MRPs inferred for magnitudes and motions carry substantial epistemic and aleatory uncertainty and are sensitive to logic-tree choices; practical validation is limited because extremely long observation windows (tens to hundreds of millennia) would be required to test such return periods. GMMs used for hazard and residual analyses are extrapolated to magnitudes and near-fault distances not well represented in existing databases, especially at large magnitudes and long periods. The ESHM20 fault representation (e.g., EAFZ as a single long source) may not fully capture detailed segmentation and multi-fault interactions; Surgu/Cardak geometry is less certain. Standard PSHA does not model short-term triggering and clustering between mainshocks. Station distribution and site-characterization uncertainties (e.g., VS30 proxies) may influence site-specific estimates.