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

A sequence of two major earthquakes, Mw 7.8 Pazarcik, and Mw 7.5 Elbistan, struck South-eastern Turkey in February 2023. The large magnitudes of the earthquakes and the short time between the two events raised questions about whether this sequence was an extremely rare disaster. Here, based on prior knowledge, we perform seismic hazard assessment for the region to estimate exceedance probabilities of observed magnitudes and ground motions. We discuss that many regional studies indicated the seismic gap in the area but with lower magnitude estimations. Observed ground motions generally agree with empirical models for the Pazarcik event. However, some records with high amplitudes exceed the highest observed amplitudes in an extensive database of shallow crustal earthquakes. We observe a notable trend of residuals for the Elbistan earthquake, leading to underestimation at long periods. We discuss potential advances in science for better characterization of such major earthquakes in the future.

the probability of exceedance values for the magnitude of the earthquakes and amplitudes of the observed ground motion levels.To obtain quantitative estimates of return periods of the fault ruptures and the observed ground motion levels, PSHA is performed for the locations of the seismic stations that recorded at least one of the two mainshocks.The recorded ground motions are also compared to predictive ground motion models (GMMs) to assess how extreme the observations are relative to the observed earthquake fault ruptures (i.e., the probability of observing ground motion levels at least as extreme as the recorded motions, conditional on the observed earthquake magnitudes and fault ruptures).

Recorded ground motions
The dataset of recorded ground motions for the two earthquakes is obtained from Gülerce et al. 10 , where they used the raw ground acceleration data published by the Disaster and Emergency Management Presidency of Turkey (AFAD).They eliminated potentially problematic records (i.e., records with incomplete trace, considerable noise content, or an unclear trace of the event) and processed the remaining records according to the procedure described in Akkar et al. 11 .Since this study focuses primarily on ground motions with relatively high amplitudes, only the records within 100 km of the fault plane are considered.In addition, the records at stations where the time-averaged shear-wave velocity within the upper 30 meters (V S30 ) is unknown are also eliminated, as this site proxy is a key input to GMMs.This resulted in 51 records for the Mw7.8 Pazarcik event and 25 records for the Mw7.5 Elbistan event.For all records, three ground motion intensity measures are calculated: peak ground acceleration (PGA) and 5% damped elastic spectral acceleration at periods of 0.5 and 1.0 seconds (SA T¼0:5s and SA T¼1:0s ).These three parameters are selected to roughly account for the variations in the frequency content of the ground motions.The RotD50 value for all intensity measures is calculated from the two orthogonal horizontal components.For two horizontal and orthogonal components, the RotD50 value is calculated by rotating the ground motion at small increments to get the acceleration at various orientations; calculating the intensity measure for each orientation; and calculating the median intensity across all orientations 12 .Figures 2 and 3 show the spatial distribution of the calculated intensity measures (i.e., PGA, SA T¼0:5s , and SA T¼1:0s ) for the Mw7.8 Pazarcik and Mw7.5 Elbistan earthquakes, respectively.The locations of the stations with respect to the fault planes are also presented as a zipped Keyhole Markup Language (.kmz) file, along with key information on the records and stations in Supplementary Data 1.
Probabilistic seismic hazard assessment for the region PSHA is a mathematical framework based on the total probability theorem for estimating the exceedance rates of intensity measure levels.This section provides a brief and simplified overview of the methodology without intending to provide a complete picture.A detailed explanation of the method, underlying assumptions, as well as recent developments can be found in Baker et al. 13 .
The mathematical equation for estimating the exceedance rate of a given intensity measure level (also called the seismic hazard integral) is expressed as follows [13][14][15] : where IM is the considered intensity measure (e.g., PGA), im is the considered level of IM (e.g., 1 g), λ IM > im ð Þis the occurrence rate of observing IM values greater than im, rup i;j is an individual fault rupture i at seismic source j, λðrup i;j Þ is the estimated rate of rupture i at the source j, and P IM > imjrup i;j ; site is the probability of observing an IM value greater than im at the considered site given that the considered rupture rup i;j takes place.By summing over all considered seismic sources in the region, and all potential ruptures in each seismic source, Eq. ( 1) estimates the total rate of exceeding im by considering the contributions from each seismic source in the region.Repeating Eq. ( 1) for many im values and plotting the exceedance rates against im values result in the most common output of PSHA, called the hazard curve.An example hazard curve for PGA, constructed for station 3125, located in Antakya, Hatay, is shown in Fig. 4.
Equation ( 1) can be separated into two components.The first one is the conditional probability of exceeding im, i.e., PðIM > imjrup i;j ; siteÞ, at the site of interest, given an earthquake rupture rup i;j .This part of the equation is most commonly calculated through empirical GMMs by assuming normally distributed model residuals 16 .Most empirical GMMs take the following simplified functional form to estimate the natural logarithm of intensity measures, lnIM 13 : where μ lnIM rup; site À Á is the median prediction, σ lnIM rup; site À Á is the total standard deviation of the model, and ϵ is a standard normal random variable with a zero mean and unit standard deviation, expressing the difference between observation and median prediction as the number of standard deviations.P IM > imjrup i;j ; site can then be calculated as: where ϕ is the cumulative distribution function of the standard normal distribution.Equation (3) can also be expressed as a  function of the normalized residual ϵ as: The second component is the occurrence rates of the considered earthquake ruptures, λðrup i;j Þ.An individual earthquake rupture can be generally defined with an estimate of the total released energy (seismic moment, M 0 , or moment magnitude, Mw) and the fault geometry.The rates of earthquakes exceeding a given magnitude are generally estimated with magnitudefrequency distributions (MFD).The most common MFDs used in practice are the double truncated Gutenberg-Richter 17 and characteristic earthquake models 18 , both of which are modifications to the original Gutenberg-Richter model 19 , which has the following form: where N M > m ð Þis the number of earthquakes with a magnitude M greater than m, and a and b are the parameters of the Gutenberg-Richter MFD.The a-value describes the overall seismic activity of the region, whereas the b-value describes the proportion of small-magnitude earthquakes to large-magnitude earthquakes.Both double truncated Gutenberg-Richter and characteristic earthquake models impose a maximum magnitude (M max ) to the classical Gutenberg-Richter relationship to eliminate very large magnitudes that are considered physically impossible.The characteristic earthquake model also assumes an increased probability of a so-called characteristic magnitude and magnitudes in its vicinity.The value of the characteristic magnitude is usually determined through historical evidence for the large-magnitude earthquakes at a given fault.On the other hand, a double truncated Gutenberg-Richter model introduces an upper bound magnitude (M max ) and a minimum considered magnitude (M min ) to the classical Gutenberg-Richter model while preserving the exponential decay.
This study uses publicly available inputs of the European Seismic Hazard Model 2020, ESHM20 20 to perform PSHA at each seismic station shown in Figs. 2 and 3. PSHA calculations are performed with the open-source OpenQuake platform 21 .Instead of using the already available seismic hazard curves computed in ESHM20, we compute site-specific curves in order not to be limited by the spatial resolution of the original model, which is designed for a continental scale.
The source model of ESHM20 includes area and fault source models for the Euro-Mediterranean region.To capture the epistemic uncertainty in the inputs of MFDs, the logic tree shown in  Fig. 5 is used.Here, in addition to the common double truncated Gutenberg-Richter MFD, a tapered Pareto distribution 22,23 with median a and b values is also used in a single branch with a weight of 0.4 for area sources.For the other branch of the area source model, a double truncated Gutenberg-Richter MFD is used with three further branches.These three branches combine the 5th, 50th (i.e., median), and 95th percentile values for a and b values with weights of 0.2, 0.6, and 0.2, respectively.The three a and b values are not exhaustively permutated as these combinations would correspond to extreme cases with non-trivial weights 20 .For the fault models, a single b value for the tectonic region is used, along with three different values for slip rate, as it was observed that the b-value uncertainties were less important in resulting MFDs when compared to the slip rate and M max uncertainties.The slip rates are compiled from the European Fault-Source Model 2020 24 .The slip rates are converted to activity rates assuming moment conservation.Maximum magnitudes for area sources in ESHM20 are estimated based on prior seismic activity.The lower bound of maximum magnitude (M max 1 ) is set equal to the highest observed magnitude for the zone, increased by one standard deviation.The  median (M max 2 ) and upper bound (M max 3 ) magnitudes are determined by adding a magnitude increment of 0.3 and 0.6, respectively.This value of 0.3 is the standard deviation of earthquake magnitude estimated from the entire unified catalog of ESHM20.The maximum magnitudes for fault sources are estimated with the fault scaling law of Leonard 25 .For each sitespecific PSHA calculation, seismic sources within a 200 km distance are considered.The locations of the area and fault sources in the source model of ESHM20 that coincide with the observed ruptures of the Kahramanmaras earthquakes are shown in Fig. 6, and the key parameters of these source models are shown in Table 1.In summary, the approach of ESHM20 for estimating the maximum magnitudes of seismic sources is based on a combination of prior major seismic activity, as well as the scaling of magnitudes with the fault rupture area 25 .For the area sources, maximum magnitudes are estimated considering the highest observed magnitudes in the region and the catalog uncertainty.
For the fault models, the maximum magnitudes are estimated based on the magnitude scaling with the fault rupture area.These two approaches are combined in the final seismic hazard estimates through the logic tree shown in Fig. 5 with a total of 21 branches.
The ground motion characterization of ESHM20 is based on the scaled backbone logic tree approach, where a single GMM is calibrated with various adjustment factors to capture epistemic uncertainties in the seismological properties of the region.For the backbone model, the partially non-ergodic GMM of Kotha et al. 26,27 is used, which has the following functional form for the median predictions: where ln IM e;l;s;r is the median prediction, and subscripts of e, l, s, and r stand for the event, tectonic locality, site, and region, respectively.In Eq. ( 6), lnμ describes the fixed-effect prediction without any adjustments for the event, tectonic locality, site, or region.δB 0 e;l , δL2L l , δS2S s , and δc 3;r are the random-effect adjustments for the event, tectonic locality, site, and region, respectively.The fixed-effect prediction, lnμ, is described as follows: where e 1 is the generic offset term, f R;g M w ; R JB À Á is the geometric spreading function, f R;a R JB À Á is the anelastic attenuation function, and f M M w À Á is the source function, representing the magnitude scaling of motion.R JB is the most common sourceto-site distance metric used in GMMs, defined as the distance to the surface projection of the fault.
The derivation of the logic tree used in the ground motion characterization of ESHM20 is described in Weatherill et al. 28 and Weatherill and Cotton 29 .Using described inputs of ESHM20, hazard curves for PGA, SA T¼0:5s , and SA T¼1:0s are generated for the seismic stations in our dataset with the Open-Quake Engine platform.

Results
Magnitude-frequency distributions.The fault source model of ESHM20 is not finely segmented unlike most regional seismic source characterization studies (e.g., Emre et al. 3 ).Hence, the entire EAFZ is modeled as a single fault source (TRCF002) with a total length of 576 km which almost completely encompass the observed fault rupture of 350 km for the Mw7.8 Pazarcik earthquake.The Mw7.5 Elbistan earthquake occurred on a relatively less studied fault structure which is potentially mapped with lower accuracy compared to the EAFZ, resulting in a moderate disagreement between the fault model of ESHM20 and the observed rupture (Fig. 6).In particular, the northeast extension of the fault zone (Surgu-1) is not mapped in the fault source model of ESHM20 within the TRCF03H fault.However, the total length of the TRCF03H and the observed rupture are quite similar as the TRCF03H fault further extends towards the southwest.When the agreement between the area source model of ESHM20 and the observed fault ruptures are investigated, it is observed that the rupture of the Mw7.8 earthquake is not confined to a single area source, but it has substantial portions in two: LBAS341 and TRAS481.Similarly, the observed rupture of the Mw7.5 Elbistan earthquake corresponds to two area sources, TRAS434 and TRAS481, although most of the rupture is within TRAS434.
Initially, MFDs for the sources of ESHM20 are investigated to estimate the probability of observing earthquakes with magnitudes larger than or equal to 7.8 and 7.5 at sources coinciding with the ruptures of Pazarcik and Elbistan earthquakes, respectively.MFDs corresponding to different logic tree branches (Fig. 5) are combined by calculating a weighted average with the same weights from the logic tree.The estimated annual rates are converted to mean return periods for easier interpretation, as shown in Fig. 7.
For the Mw7.8 Pazarcik earthquake, the mean return period is estimated as about 1000 years with the fault source model, TRCF002, which resembles the observed fault rupture very similarly towards the southern end (Fig. 6).The two area sources in the same region, TRAS481, and LBAS341, yield mean return period estimates of about 4000 and 3500 years, respectively.For the Mw7.5 Elbistan earthquake, the mean return period is estimated as about 7800 years with the corresponding fault source, TRCF03H, and 2600 years, with the corresponding area source, TRAS434.

Mean return periods of observed ground motion amplitudes.
For each station in the compiled dataset, mean return periods for the observed levels of three parameters, PGA, SA T¼0:5s , and SA T¼1:0s are calculated from the hazard curves calculated in the previous section.Calculation of the mean return period for the observed PGA level of 0.94 g at station 3125 located in Antakya, Hatay, is illustrated in Fig. 4. Spatial distributions of estimated mean return periods are shown in Figs. 8 and 9 for the Mw7.8 Pazarcik and Mw7.5 Elbistan earthquakes, respectively.Estimated return periods and corresponding observed ground motion levels are shown in Fig. 10 for both earthquakes.The maximum estimated return periods for exceeding the observed ground motion levels during the Mw7.8 Pazarcik earthquake are approximately 6900, 4100, and 6900 years for PGA, SA T¼0:5s , and SA T¼1:0s , respectively.Maximum estimates for the Mw7.5 Elbistan earthquake are 1200, 1700, and 6900 years for PGA, SA T¼0:5s , and SA T¼1:0s , respectively.Figure 10 also displays the 16 th and 84 th percentile return period estimates among all individual logic tree branches.Here, it is observed that the estimates are highly variable and sensitive to logic tree branches.The variability is most prominent for very high return periods, especially for the upper limit of the confidence interval.
Comparison of recorded ground motions to GMMs.The recorded ground motions can also be compared to GMMs with inputs constrained to represent the observed fault rupture and subsurface conditions at the site of interest.Although the GMM of Kotha et al. 26,27 was used for the PSHA in the previous section, the GMM of Boore et al. was preferred for the residual analyses performed in this section, as the dataset of Boore et al. 30 contains earthquakes with magnitudes as large as 7.9.However, the maximum magnitude earthquake used in the derivation of the Kotha et al. 26,27 model was 7.4, which is lower than the magnitudes of both earthquakes analyzed in this study.The residuals at vibration periods between 0.01 and 10 seconds for the ground motions of both earthquakes are shown in Fig. 11 according to Boore et al. 30 model.Since it is a standard normal variable, the ϵ values can also be interpreted as the probability of exceeding the observed levels of ground motion, conditional on the observed earthquake rupture and the known subsurface conditions (i.e., V S30 ).The corresponding percentile values for ϵ values are also given in parentheses in Fig. 11.The residuals for the recorded motion at station 3135, located in Arsuz, Hatay, are observed to be generally high across all periods, with a maximum value of about 3 (corresponding to a probability of exceedance of only 0.1%) at a period of 0.3 seconds.A more in-depth investigation of this record is provided in Fig. 12.Here, the elastic response spectra of the observed motion are compared to the 475-and 2475-year return period design spectra described in the Turkish Building Earthquake Code 31 and Boore et al. 30 .It is observed that the recorded motion generally exceeds both code provisions and model predictions in almost the entire period range.Figure 12 also presents ground acceleration and velocity series, as well as the Fourier spectrum of the recorded motion.The station is in the forward directivity region and located very close to the fault rupture (Fig. 12b), therefore, the motion is expected to be considerably affected by the finite-fault effects.The high amplitudes and relatively short duration of the motion also support the potential forward directivity effects.However, when the velocity trace of the record is investigated (Fig. 12e), no clear pulse-like features are observed.
One peculiar thing about this ground motion record is its broadband nature, with very high amplitudes across all vibration periods.The probability of observing such broadband motions can be approximated through the multivariate normal distribution for two or more vibration periods, as spectral accelerations at multiple periods are well approximated by a joint normal distribution 32 .For simplicity, two periods are selected for this analysis as 0.3 and 3.0 seconds to represent shorter and longer periods, respectively.The period of 0.3 seconds is selected as the period of maximum residual, and 3.0 seconds is selected as a period of high positive residual, which is also further away from 0.3 seconds so that the correlation between the two vectors of residuals is not very high.The observed residuals at these two periods are 2.96 and 1.21 (Fig. 12).The correlations between spectral acceleration values at different vibration periods are estimated with the correlation model of Baker and Jayaram 33 .This model is preferred even though it was derived from an older version of the NGA database, as the correlation estimates are shown to be insensitive to the choice of GMM 33 .This model predicts a correlation of 0.2535 for GMM residuals at periods of 0.3 and 3.0 seconds.Using this information, the probability of observing a ground motion that is at least as extreme as the observed motion at periods of 0.3 and 3.0 seconds is calculated as follows: Pr ϵ T¼0:3s ≥ 2:96 ϵ T¼3:0s ≥ 1:21 ρ ϵ T¼0:3s ;ϵ T¼3:0s ¼ 0:2535 For completeness, the same probability is calculated for periods of 0.5 and 1.0 seconds, as these periods were used prior to summarize the observed ground motions.The correlation between the residuals for these two periods is 0.749, calculated with the model of Baker and Jayaram 33 .The observed residuals at 0.5 and 1.0 seconds are 2.19 and 1.70, respectively.Consequently, the probability of observing a ground motion that is at least as extreme as the observed motion at periods of 0.5 and 1.0 seconds is calculated as follows: Pr ϵ T¼0:5s ≥ 2:19 ϵ T¼1:0s ≥ 1:70 ρ ϵ T¼0:5s ;ϵ T¼1:0s ¼ 0:749 These low probabilities indicate that the observed high amplitudes are not localized peaks, but the recorded motion at this station is considerably higher than expectations within the entire period range.The low probabilities and high residuals calculated in this section might indicate that the observed ground motions are much higher than the expected levels for such largemagnitude earthquakes.However, the employed GMMs are derived from datasets that lack records from such largemagnitude earthquakes recorded very close to the fault ruptures.In fact, the records of this earthquake sequence will likely lead to updates of existing empirical models, especially for constraining the scaling of motions with larger magnitudes at short source-tosite distances.Fig. 7 Mean return periods for the observed earthquake magnitudes of 7.8 and 7.5.a Mean return period of a Mw7.8 earthquake at the sources coinciding with the rupture of the Pazarcik earthquake.b Mean return period of a Mw7.5 earthquake at the sources coinciding with the rupture of the Elbistan earthquake.Mean return periods are calculated for both area and fault sources of ESHM20.In both figures, the individual return periods estimated for a single branch of the logic tree are also plotted against magnitude.
Comparison of recorded ground motions to previous recordings.This section compares the recorded ground motions from the two earthquakes to some of the most extreme ground motions in the NGA-West2 database 34 , containing more than 20,000 records from 599 shallow crustal earthquakes in active tectonic regimes.Records from the NGA-West2 database with the highest spectral acceleration values at discrete periods between 0.01 and 10 seconds are identified for comparison.This selection criterion resulted in 19 unique ground motion records from nine earthquakes, as shown in Table 2.A comparison of the response spectra of these 19 records against five of the records from 6 February 2023 is shown in Fig. 13.Here, four records are selected from the Mw7.8 Pazarcik earthquake at stations 3126, 3138, 3139, and 3135.In addition, one record from the Mw7.5 Elbistan earthquake at station 4612 is also included in the comparison.It is observed that the selected five records are generally comparable to the maximum values in the entire NGA-West2 database.Particularly, SA T¼1:5 value at station 3138 (1.20 g) is slightly higher than the maximum value in the NGA-West2 database (1.17 g).Similarly, SA T¼1:8 value at station 3139 (1.12 g) is shown to be higher than the maximum value in the NGA-West2 database (1.08 g).

Discussion
The Kahramanmaras earthquake sequence on 6 February 2023 initiated many discussions on whether these earthquakes were anticipated or not.This study offers an investigation of the observed fault ruptures and recorded ground motions to estimate how extreme the two earthquakes are, based on the prior scientific knowledge of the region.The seismic potential of the region was well known, but the maximum magnitude estimates were generally lower than the observed magnitudes of 7.8 and 7.5.In particular, Amanos and Pazarcik segments of EAFZ (Fig. 1b) were generally considered independent fault segments, and a combined rupture of the two segments was not considered in previous regional seismic hazard studies.However, a recent seismic hazard model, ESHM20, modeled most of EAFZ as a single fault source with a length of more than 500 kilometers.
The use of fault sources for seismic hazard assessment is generally preferred over area sources when the seismotectonic structure is well known, and records of prior seismic activity allow precise modeling of fault locations and MFDs.This was indeed the case for the EAFZ, as many studies focused on its geometry, segmentation, slip rate, and seismic potential [3][4][5][6][7] .Conversely, area sources are generally preferred where the fault structures are not very well known due to a lack of previous studies and historical large earthquakes in the vicinity.The fault source model yielded a mean return period of about 1000 years for the observed magnitude of 7.8 at EAFZ, whereas the two corresponding area sources yielded estimates between 3500-4000 years.For the Surgu and Cardak faults, the fault source model estimates a return period of 7800 years, whereas the area source model estimates 2600 years.In summary, the fault source model yields a lower estimate for the EAFZ, while the area source model yields a lower estimate for the Surgu and Cardak faults.This might indicate that fault sources work better in well-studied regions such as EAFZ.On the other hand, lack of information and higher epistemic uncertainty on fault geometry and seismicity at some sources might not allow accurate modeling and constraint of the source as a fault source model.This might be the case for the Surgu and Cardak faults, where the fault source model gives a very high mean return period of 7800 years.The maximum amplitudes are generally recorded in Hatay during the Mw7.8 Pazarcik earthquake.Accordingly, the highest return periods are also estimated for Hatay, as the ground motions recorded at many stations in Hatay exceeded mean return periods of 1000 years.The high amplitudes in Hatay could be attributed to rupture directivity effects, as most of the province is in the forward-directivity region.However, the concentration of high amplitude records in Hatay could also be related to the denser seismic network coverage, as the highest number of stations recording the Mw7.8 Pazarcik earthquake are also in Hatay in our dataset with 22 stations.After Hatay, the highest number of stations are in Kahramanmaras, with nine records, and there are no other provinces with more than five records in the dataset.
The Mw7.5 Elbistan earthquake was not as densely recorded as the Mw7.8 Pazarcik earthquake, with 25 records within 100 km rupture distance retained in our dataset.The highest amplitudes for this event are recorded in station 4612, in the Goksun district of Kahramanmaras, with PGA, SA T¼0:5s and SA T¼1:0s values of 0.58, 1.02, and 0.61 g, respectively.The mean return period estimates for these amplitudes are 750, 830, and 910 years for PGA, SA T¼0:5s and SA T¼1:0s , respectively, for this station.Interestingly, these are not the highest return period estimates for this earthquake.For PGA, the highest return period is estimated for station 4406, located in Malatya, with 1250 years for the observed value of 0.44 g.For SA T¼0:5s and SA T¼1:0s , highest return periods are estimated for station 3802 located in Kayseri for the observed values of 0.39 and 0.38 g, with 1600 and 6900 years, respectively.This discrepancy is due to Kayseri being located in a relatively  low-hazard region.Still, station 3802 recorded considerably high amplitudes, especially at longer periods.The standard residential buildings in Turkey are constructed for seismic loads corresponding to a return period of 475 years, determined according to TBEC 31 based on location and seismic site class.This return period is determined by assuming an average service time of 50 years for a building and aiming for a 10% probability of exceedance in those 50 years, assuming a Poisson distribution for exceedance rates.Similarly, return periods of 1000, 2500, 5000, and 10,000 years roughly correspond to the probability of exceedance values of 5%, 2%, 1%, and 0.05%, respectively, in a 50-year service time.The maximum return period estimates for the observed ground motions in these two earthquakes are about 7000 years, corresponding to a probability of exceedance value of 0.7% in 50 years.This implies that this earthquake sequence simply resulted in much higher demands on structures compared to the levels they are designed for.However, these very high return periods (or low probabilities) should be interpreted very carefully and should not be thought to justify the observed widespread damage and losses.This study focused mostly on the most extreme motions recorded during the earthquake sequence.These few records can actually be considered local peaks of the spatial distribution of the motion.Therefore, the majority of the earthquake-affected region is actually thought to have experienced much lower ground shaking levels, which can also be inferred from the spatial distribution of peak ground motion parameters shown in Figs. 2 and 3.
When the behavior of residual spectral accelerations is investigated, it is observed that there is no apparent trend for the residuals of the Mw7.8 Pazarcik earthquake, meaning that  the observed ground motions are generally in agreement with empirical models.However, some records have very high residuals across various periods (e.g., on stations 3126, 3129, 3135).For the Mw7.5 earthquake, a positive trend against the vibration period is observed for residuals, indicating substantial underestimation at longer periods.
Two stations in the compiled dataset for this study recorded spectral accelerations higher than the maximum values in the NGA-West2 database at periods of 1.5 and 1.8 seconds, roughly corresponding to the first mode periods of mid to high-rise reinforced concrete buildings.Since near-fault recordings of large-magnitude earthquakes are rare in recorded ground motion datasets, physics-based ground motion modeling might be preferred over empirical modeling for large earthquakes.Physicsbased ground motion modeling might allow more accurate modeling of important physical phenomena, such as the effects of fault geometry, rupture directivity, rupture velocity, and longperiod amplifications due to deep basins.
Recent advances in non-ergodic ground motion modeling are also promising tools for better ground motion estimation.In an ergodic process, the spatial variation of the variable is equivalent to its temporal variation, meaning that a model derived from observations at spatially distributed locations can approximate the distribution over time at a single site.This assumption has to be made in generic ground motion modeling and PSHA due to the limited data in similar conditions repeating over time.However, with the increasing density of ground motion networks and number of observations, it is becoming more and more viable to develop non-ergodic GMMs to account for the site-and region-specific characteristics of ground motions 35,36 .These non-ergodic GMMs have the potential to estimate ground motions more accurately, especially with the increasing site-and region-specific ground motion data, by accounting for repeatable source, path, and site effects.However, most ground motion databases still lack near-fault records of large-magnitude events, which are generally the most critical scenarios for the seismic hazard assessment of sites located near active faults.Therefore, care must be taken to extrapolate these models to extreme conditions, such as imposing physical constraints for near-source magnitude scaling 35 .The large number of near-source observations from this earthquake sequence will likely lead to better characterization of ground motions from large magnitude earthquakes in ergodic ground motion modeling, as well as the estimation of ergodic terms for partially non-ergodic models.These potential improvements in ground motion modeling could yield much more robust hazard estimates in the future.
This study investigates the Mw7.8 Pazarcik and Mw7.5 Elbistan earthquakes and their ground motions separately.However, the Pazarcik earthquake most likely triggered the Elbistan earthquake, which occurred only nine hours later.Therefore, these two earthquakes are not statistically independent due to this triggering.Common practice PSHA does not consider triggering of secondary mainshocks or short-term clustering of large magnitude earthquakes.An active area of research, physics-based simulations of earthquake sequences, is a promising tool for modeling the complex behavior of fault ruptures, which can impose physical constraints on short-term clustering of large magnitude events and multi-segment fault ruptures in complex settings [37][38][39][40][41][42][43][44][45] .Incorporation of such models into PSHA might lead to much more accurate earthquake recurrence estimates.
One final note is necessary about the very high return periods estimated in this study, both for the earthquake recurrences and observed ground motion amplitudes.),.9][50][51][52] ), these validations are bound to be for short time windows.In particular, Beauval et al. 48alculated the minimum required observation window to accurately calculate the mean return period of a Poisson process.They showed that a minimum observation window of 12,000 years is required for a return period of 475 years to ensure a coefficient of variation of 20% or less.Return periods of 3000 and 10,000 years require minimum time windows of 75,000 and 250,000 years, respectively, for the same threshold of 20% coefficient of variation.Consequently, the very high return periods estimated in this study contain substantial uncertainty and are likely sensitive to changes with improved modeling and longer instrumentation windows in the future.The return period estimates are also shown to be highly sensitive to input combinations in hazard calculations (see Fig. 7 and Fig. 10), highlighting the importance of accounting for the epistemic uncertainty in seismic hazard studies through the careful construction of logic trees.

Fig. 1
Fig. 1 Comparison of the observed fault ruptures to the active fault map of Turkey.a Surface projection of fault ruptures for 6 February earthquakes, along with the simplified traces of EAFZ and NAFZ b Surface projection of fault ruptures for 6 February earthquakes in comparison to segmentation of active fault map of Turkey 4 .The surface projections of the fault ruptures are obtained from the finite-fault models of the United States Geological Survey 8 .A blue rectangle in a indicates the selected frame in b.The estimated maximum magnitude for each segment is shown in parentheses in the legend of b.

Fig. 2
Fig. 2 Maps of recorded ground motion intensities for the Mw7.8 Pazarcik earthquake.Maps show intensities in terms of a peak ground acceleration, PGA, b spectral acceleration at a period of 0.5 s, SA T¼0:5s , c spectral acceleration at a period of 1.0 s, SA T¼1:0s .

Fig. 3 Fig. 4
Fig. 3 Maps of recorded ground motion intensities for the Mw7.5 Elbistan earthquake.Maps show intensities in terms of a peak ground acceleration, PGA, b spectral acceleration at a period of 0.5 s, SA T¼0:5s , c spectral acceleration at a period of 1.0 s, SA T¼1:0s .

Fig. 5 Fig. 6
Fig. 5 The source model logic tree of ESHM20 for shallow crustal earthquakes.The subscripts 1, 2, and 3 represent the lower bound, median, and upper bound values for the inputs, respectively.TGR stands for double truncated Gutenberg Richter MFD, and SR stands for slip rate.b tecno is the median b value for the tectonic region.

Fig. 8
Fig. 8 Maps showing the calculated mean return periods for the observed ground motion levels during the Mw7.8 Pazarcik Earthquake.Maps show return periods of a PGA, b SA T¼0:5s , c SA T¼1:0s for the Mw7.8 Pazarcik earthquake.

Fig. 9
Fig. 9 Maps showing the calculated mean return periods for the observed ground motion levels during the Mw7.5 Elbistan earthquake.Maps show return periods of a PGA, b SA T¼0:5s , c SA T¼1:0s for the Mw7.5 Elbistan earthquake.

Fig. 10
Fig. 10 Estimated return periods plotted against the observed ground motion levels.For a PGA, b SA T¼0:5s , c SA T¼1:0s .Stations at which some of the highest amplitudes are recorded or highest return periods are estimated are annotated.Red scatters show the data for the Mw7.8 Pazarcik earthquake, and blue scatters show the data for the Mw7.5 Elbistan earthquake.The lines passing through the scatters indicate the 16th and 84th percentile estimates.

Fig. 11
Fig. 11 The normalized residual spectral acceleration (ϵ) plotted against vibration periods between 0.01 and 10 seconds.For a Mw7.8 Pazarcik earthquake, b Mw7.5 Elbistan earthquake.An individual record with very high ϵ values, recorded at station 3135 for the Mw7.8 Pazarcik earthquake is plotted in blue.The corresponding probability of exceedance values for ϵ values are given in parentheses.

Fig. 12
Fig.12Investigation of the ground motion recorded at station 3135, located in Arsuz, Hatay.a 5% damped elastic response spectrum of the ground motion.The response spectrum is compared to the 475-and 2475-year return period design spectra defined in Turkish Building Earthquake Code31 , and the GMM of Boore et al.30 .b The location of station 3135 with respect to the fault rupture of the Mw7.8 Pazarcik earthquake.c Recorded ground accelerations at station 3135 for East-West and North-South components.Cumulative Arias intensity, i.e., Husid plot, is shown with the dashed line.d Fourier amplitude spectra of the recorded ground motions for East-West and North-South components.e Recorded ground velocity at station 3135 for East-West and North-South components.

Table 1
Key parameters of the area and fault sources in the source model of the European Seismic Hazard Model, ESHM20, that corresponds to the fault ruptures of the Mw7.8 Pazarcik and Mw7.5 Elbistan earthquakes.

Table 2 A
summary of ground motion records with the highest spectral acceleration values in the NGA-West2 database between discrete periods between 0.01 and 10 seconds.Fig.13Comparison of the 5% damped elastic response spectra of five records from 6 February 2023 Kahramanmaras earthquakes against the records with the highest amplitudes in the NGA-West2 database.Records at stations 3126, 3138, 3139, and 3135 are from the Mw7.8 Pazarcik earthquake.The record at station 4612 is from the Mw7.5 Elbistan earthquake.For the NGA-West2 database, records with maximum spectral acceleration values at discrete periods between 0.01 and 10 seconds are selected.