Aftershock f or ecasts based on incomplete earthquake catalogues: ETASI model application to the 2023 SE T ¨urkiye earthquake sequence

The epidemic-type aftershock sequence (ETAS) model is the state-of-the-art approach for modelling shor t-ter m ear thquake clustering and is preferable for shor t-ter m aftershock forecasting. Ho wever , due to the large variability of different earthquake sequences, the model parameters must be adjusted to the local seismicity for accurate forecasting. Such an adjustment based on the ﬁrst aftershocks is hampered by the incompleteness of earthquake catalogues after a main-shock, which can be explained by a blind period of the seismic networks after each earthquake, during which smaller events with lower magnitudes cannot be detected. Assuming a constant blind time, direct relationships based only on this additional parameter can be established between the actual seismicity rate and magnitude distributions and those that can be detected. The ET AS-incomplete (ET ASI) model uses these relationships to estimate the true ETAS parameters and the catalogue incompleteness jointly. In this study, we apply the ETASI model to the SE T ¨urkiye earthquake sequence, consisting of a doublet of M 7.7 and M 7.6 earthquakes that occurred within less than half a day of each other on 6 Febr uar y 2023. We show that the ETASI model can explain the catalogue incompleteness and ﬁts the observed earthquake numbers and magnitudes well. A pseudo-prospective forecasting experiment shows that the daily number of detectable m ≥ 2 can be well predicted based on minimal and incomplete information from early aftershocks. Ho wever , the maximum magnitude ( M max ) of the next day’s aftershocks would have been overestimated due to the highl y v ariable b v alue within the sequence. Instead, using the regional b value estimated for 2000–2022 would have well predicted the observed M max values.


I N T RO D U C T I O N
As observed for the M 7.7 mainshock that occurred on 6 February 2023 in SE T ürkiye, thousands of aftershocks usually follow a large earthquake (mainshock) closely in space and time.The largest of these aftershocks can be as destructive or deadly as the mainshock or even worse.Therefore, in the immediate aftermath of a mainshock, accurate predictions of the expected frequency and magnitude of aftershocks are essential for planning emergency decisions and risk mitigation in the affected area.In general, the aftershock occurrence is well described by the Omori-Utsu (OU) decay for the earthquake rate, R = K 0 /( c + t ) −p , with parameters K 0 , c and p , and the Gutenberg-Richter (GR) distribution for the earthquake magnitudes (Utsu et al. 1995 ).Generic parameters can be used to get a first estimate of the expected activity (Page et al. 2016 ), but the decay parameters are not universal and are known to v ary widel y between dif ferent mainshocks.Therefore, the estimation of sequence-specific parameters is preferable, but this task is complicated by (i) large aftershocks triggering their own aftershock sequences and by (ii) the incompleteness of earthquake catalogues after mainshocks when only the largest aftershocks are detected, which is statistically observed (Kagan 2004 ;Omi et al. 2013 ) and firmly proven by detailed analysis of waveform data (Peng et al. 2006(Peng et al. , 2007 ; ;Enescu et al. 2007Enescu et al. , 2009 ) ).
The first problem can be addressed by the well-known epidemictype aftershock sequence (ETAS) model which considers secondary aftershock triggering and describes the actual earthquake rate as the sum of a background rate ( r ) and the rates of ongoing OUtype aftershock sequences from past earthquakes (Ogata 1988 ).The ETAS model assumes that the catalogue is complete for earthquakes of magnitude m ≥ M c and gives the earthquake rate of m ≥ M c events by where the index i refers to past events at times t i with magnitudes m i ≥ M c (Ogata 1988 ).Here, the proportionality factor of the OU decay function scales exponentially with the trigger magnitude (parameters K and α), in agreement with observations (Utsu et al. 1995 ;Hainzl & Marsan 2008 ).For a given earthquake sequence, the five ETAS parameters ( r , K , α, c , p ) are generally estimated by the maximum likelihood method.Ho wever , the ETAS estimates can be strongly biased by shorttime aftershock incompleteness (STAI).For example, the empirical relation of Helmstetter et al. ( 2006 ) predicts that m = 4 aftershocks are not fully recorded within the first 2 hours after a M = 7.7 earthquake, such as the 2023 T ürkiye event, while m = 3 events are only partially recorded, if at all, within the first 2 d.In contrast, the regional completeness magnitude before the mainshock was less than 2. This strong increase in completeness magnitude can be explained by the increased seismic noise due to the coda waves of the mainshock and its high number of aftershocks (Kagan 2004 ;Hainzl 2016a , b ;de Arcangelis et al. 2018 ).To estimate sequence-specific parameters from early aftershocks with eq. ( 1), M c must be set to a high value, resulting in a small sample size and therefore uncertain parameters, or by trying to replenish the missing data (Zhuang et al. 2017 ;Zheng et al. 2021 ).An alternative approach is to explicitly model the time-dependent incompleteness of the catalogue and use the incomplete magnitude range as well.
The frequency-magnitude distribution of recorded earthquakes can usually be well described by the model of Ogata & Katsura ( 1993 ), hereafter referred to as OK93.It is given by the GR distribution, N ( m ) = 10 a − bm , multiplied by the detection probability, which is expressed by a cumulative normal distribution with mean μ and standard deviation σ .The time evolution of both model parameters can be first estimated from the observational data using smoothness constraints (see details in Sec. 3 ).Then the detection probability can be used to estimate unbiased OU and ETAS parameters.Omi et al. ( 2013Omi et al. ( , 2014Omi et al. ( , 2015 ) ) have shown in a series of papers that this approach works mostly w ell.How ever, it does not provide point information because it involves smoothing μ and σ in time.Fur ther more, it cannot estimate the inherent incompleteness of future catalogue data, which may be important for prospective forecasts that are tested with recorded data.
A direct dependence of the detection probability on the earthquake rate has recently been proposed by Hainzl ( 2016a ), based on the simple assumption that the local seismic network cannot detect smaller-magnitude events within a certain period after an earthquake, the so-called blind time T b .Taking this additional parameter into account, Hainzl ( 2016b ) showed that STAI explains the apparent dependence of the OU-parameter c value on the mainshock magnitude.Mizrahi et al. ( 2021 ) first implemented a slightly different rate-dependent detection probability function in ETAS, while Hainzl ( 2021 ) extended ETAS with the pre viousl y deri ved detection probability function, hereafter called ETASI model, where the parameter T b and thus the detection probability is estimated simultaneously with the other ETAS parameters, without smoothing.The ETASI approach is thus much simpler and straightforward, but it relies on the validity of the blind-time approach (see description in Section 2 ).
In this paper, we explore the potential of the method for shortterm aftershock forecasting using the exemplary data of the 2023 T ürkiye sequence, described in Section 4 .We focus only on the forecasts of the number of earthquakes and their magnitudes, ignoring their spatial distribution in order to simplify the presentation and to av oid unresolved prob lems of the space-time ETAS model related to the use of isotropic kernels (Hainzl et al. 2008 ).The results are presented and discussed in Sections 5 and 6 .

E TA S I M O D E L
In the standard ETAS approach, the earthquake catalogue is assumed to be complete for m ≥ M c at any time, and R 0 (eq. 1 ) is fitted directly to the observed data.In contrast, ETASI considers the STAI of earthquake catalogues in the model fit.Assuming (i) a blind time T b which follows each earthquake during which events smaller than that of the earthquake cannot be detected and (ii) a GR distribution for the earthquake magnitudes, the detection probability p d for an earthquake with magnitude m at time t is given by where R 0 ( t ) is the true earthquake rate of m ≥ M c events (Hainzl 2016b ).With this detection probability, both the detectable (apparent) rate R and the probability density function of detectable magnitudes f ( m , t ) can be determined anal yticall y; in particular, the detectable rate is given by (3) (Hainzl 2016b(Hainzl , 2021 ) ).
The parameters of the underlying ETAS model can be estimated by the log-likelihood function LL , which is given by 2021 ).Unlike the log-likelihood function of the standard ETAS model, the LL optimization of the b value and the ETAS parameters cannot be separated in the ETASI model, because f depends not only on the magnitude but also on R 0 ( t ) and consequently on the ETAS parameters.In addition to the standard parameters ( b , μ, K , α, c , p ), the LL -value depends on the blind time T b .Therefore, the maximum likelihood search involves seven instead of six free parameters.Note that ETASI can be easily extended to a space-time ETAS model as described in Grimm et al. ( 2022 ).
Based on the estimated ETASI parameters, the catalogue completeness can be e v aluated b y inverting eq. ( 2 ) for the magnitude which is detectable with a given probability p d at time t , namely In practice, one may be interested in the completeness level defined by the magnitudes detectable at 50 per cent or 90 per cent; thus one uses eq. ( 6) with p d = 0.5 or 0.9, respecti vel y.

E S T I M AT I O N O F T H E C ATA L O G U E C O M P L E T E N E S S B Y O K 9 3
To assess the ability of ETASI to correctly model catalogue completeness, we use the independent OK93 model (Ogata & Katsura 1993 ).Similar to ETASI, a GR distribution of earthquakes is assumed in OK93, but the detection function is not related to the earthquake rate and is characterized by independent parameters.In particular, the detection probability is modelled by a cumulative normal (Gaussian) distribution with mean value μ and standard deviation σ (Ringdal 1975 ), and the probability density distribution for the magnitudes becomes where C is a normalization constant.Here, μ is the magnitude value at which 50 per cent of earthquakes are detected, and σ refers to the range of magnitudes in which earthquakes are only partially detected.
To characterize the non-stationarity of the magnitude distribution, we assume time-dependent parameters b ( t ), μ( t ) and σ ( t ).In particular, we use β = ln (10) b and assume that the parameters ln ( β), ln ( μ) and ln ( σ ) are expressed by flexible functions of time t using cubic B-spline functions.For a given data set ( t i , m i with i = 1, . . .N ), we consider the log-likelihood function given by using the shortcuts β i ≡ β( t i ), μ i ≡ μ( t i ), σ i ≡ σ ( t i ).To account for smoothness in time, we consider a penalized function for the smoothness constraints with corresponds to the logarithm of a posterior function, where the smoothness constraints correspond to the logarithm of a Gaussian prior function since the functions β( t ), μ( t ) and σ ( t ) have linear coefficients of cubic B-spline functions (Ogata & Katsura 1993 ).In order to optimize the strengths of the constraints, we maximize the integration of the posterior function with respect to the six weights w β 1 , w β 2 , w μ 1 , w μ 2 , w σ 1 and w σ 2 , which is equi v alent to minimizing the ABIC (Akaike 1980 ). Finall y, b y fixing the optimal weights, we calculate the maximum a posterior solution of the coefficient parameters of the penalized log-likelihood function in eq. ( 10).See Ogata & Katsura ( 1993 ) for a more detailed calculation procedure.
Here, the penalty applied to the so-called natural time (event index), which advances by one between two successive events, may work better than the penalty applied to real time in cases where the occurrence times are highly clustered.

D ATA
The earthquake catalogue is provided by the Disaster and Emergency Management Organization of T ürkiye (AFAD).The catalogue is based on automatic processing, partially manually revised, of waveform data from a permanent network.An in-depth analysis of the AFAD catalogue characteristics is provided by C ¸ıvgın & Scordilis ( 2019 ).For this study, earthquakes were selected in the box defined by longitude 33.5-41.0• and latitude 34.0-40.0• , centred on the 2023 sequence.The region was chosen to be broad enough to include most of the aftershock activity, but to exclude larger clusters of independent events at the edges of the region.For the ETASI and OK93 fits, we used data recorded between 1 January and 8 March 2023.For visual comparison and analysis of the background b value, we also used the earthquake catalogue from the same region recorded between the beginning of 2000 and the end of 2022.Fig. 1 shows the temporal and spatial distribution of the data.In Fig. 1 (a), the magnitudes are plotted against the occurrence times of the recorded earthquakes.The chosen cut-off magnitude M c = 1.95 for the ETASI model application is indicated by the horizontal dotted line, which corresponds to m ≥ 2 events due to the magnitude binning of 0.1.Although the catalogue is complete for m ≥ 2 events according to the OK93 model before the mainshock, it is incomplete for the first few days after the mainshock.Ignoring this incompleteness may lead to a significant underestimation of the aftershock productivity of the largest earthquakes due to an unrealistically small α value (Kagan 2004 ;Hainzl et al. 2013 ).In Fig. 1 (b), the spatial distribution of the fitted m ≥ 2 earthquakes is shown by blue symbols, indicating that the activity is well separated in space from other activity.

R E S U LT S
The study addresses two questions, namely (i) how well ETASI can fit the complete sequence and (ii) how good 1-d forecasts are based on parameter estimates from incomplete data of preceding early aftershocks.To answer both questions, we fitted the m ≥ 2 earthquakes from 1/1/2023, that is T 1 = −35.8d relative to the mainshock, where we ignored potential aftershock-triggering by preceding events in 2022, as there is no apparent clustering at the beginning of the sequence (see Fig. 1 a).

Fit to the whole sequence
The maximum likelihood estimates (MLE) for the recorded m ≥ 2 earthquakes in 2023 ( T 1 = −35.8d and T 2 = 30 d) using the ETAS and ETASI models are given in Table 1 .The ETASI model provides a significantly better fit in terms of the corrected information gain per earthquake (IGPEc), which considers the additional free parameter of the ETASI model.The ETAS model naturally provides biased estimates with a poor fit because, as discussed in Section 2 , it uses incomplete aftershock data immediately following the mainshock.The estimated parameters are also significantl y dif ferent between ET AS and ET ASI; in particular, the α value increases strongly when the ETASI model is applied.While the ETAS model obtains a low value of 0.30, the ETASI model obtains a value of 1.24.A similar trend is observed for the b value, which increases from 0.56 to 0.87.Both observations are consistent with those observed for mainshock-aftershock sequences in California (Hainzl 2021 ).Similarly, the estimated blind time of T b = 162 s is in the range of the estimated values in California.The model fits are shown in Fig. 2 as a function of logarithmic time after the mainshock.In particular, the model rates are compared with the event rates of the recorded m ≥ 2 earthquakes in Fig. 2 (a).The plot shows a surprising result from the ETAS model.Ignoring STAI, it expects an increasing aftershock rate immediately after the mainshock despite the Omori-Utsu decay of direct aftershocks triggered by the mainshock.This is due to the low value of α, which gives significant weight to the triggering of secondary aftershocks (aftershocks triggered by aftershocks); α = 0.3 means that eight magnitude 4.7 aftershocks together trigger, on average, the same number of aftershocks as the magnitude 7.7 mainshock, but they are typically 1000 times more frequent (assuming b = 1).Thus, the cumulative triggering of the first aftershocks exceeds those of the mainshock leading to an increasing rate.In contrast, Table 1.Parameters that maximize the log-likelihood value ( LL ) for the 9438 m ≥ 2.0 earthquakes in 2023, that is between T 1 = −35.8and T 2 = 30 d relative to the mainshock.The column IGPEc refers to the corrected information gain per earthquake relative to the ETAS model, which is determined by the difference of the corrected Akaike information criterion (AICc) between the ETASI and the ETAS model, normalized by twice the event number N , where AICc = −2 LL + 2 N p + ( N p + 1) / ( N − N p − 1) with N p being the number of fit parameters (Rhoades et al. 2014 )  the ETASI model assumes a much higher α value.The true rate of m ≥ 2 earthquakes is estimated as the dashed line, which is about a factor of 500 higher than the observed rate immediately after the mainshock.Ho wever , considering the detection probability function (eq. 2 ) related to the detection limits, the detectable (apparent) ETASI rate R ( t ) agrees very well with the observed event rate.
The magnitude versus time plot of the recorded events (Fig. 2 b) shows that many small magnitude aftershocks immediately after the largest earthquakes were missed.This deficiency is compared with ETASI and OK93 estimates.The red and magenta curves show the magnitude value that can be detected at a given time with a probability of 50 and 90 per cent, respecti vel y, according to eq. ( 6).The 50 per cent cur ve ag rees well with the independent Downloaded from https://academic.oup.com/gji/article/236/3/1609/7512212 by Bibliothek des Wissenschaftsparks Albert Einstein user on 10 April 2024 OK93 estimate of the μ value (green curve), which also refers to a detection probability of 50 per cent.Both results agree within the uncertainties of μ with two exceptions: (i) For a short period after the largest aftershock, OK93 gives a higher completeness magnitude than ETASI and (ii) at the time of the fourth largest event in the sequence, the m = 6.4 earthquake that occurred about 14.7 d after the mainshock, ETASI expects a significant jump in incompleteness, while OK93 gives only a moderate increase.
The estimated catalogue completeness is also compared with the empirical completeness relation, M c ( m , t ) = m − 4.5 − 0.75log ( t ), using the mainshock magnitude m = 7.7 (dashed line).Although this relation has only been calibrated to mainshocks in California (Helmstetter et al. 2006 ), the 50 per cent curve approximately follows it until the occurrence of the largest aftershock.Ho wever , with the occurrence of the largest aftershock, the completeness magnitude jumps up and then decays faster than expected from the empirical relationship.Fur ther more, the 90 per cent probability cur ve is al wa ys about one magnitude larger than the 50 per cent curve, indicating that many earthquakes are still missed between the two curves.Thus, it is crucial to use a time-dependent detection probability function instead of a time-dependent threshold value such as the empirical function of Helmstetter et al. ( 2006 ).
The ETASI model also predicts a time-dependence of the frequency-magnitude distribution of the catalo gue e vents according to eq. ( 5).In contrast, the standard ETAS approach implicitly assumes a complete catalogue with a constant b value, where the mean magnitude is constant over time.Fig. 2 c shows the comparison of the predicted mean magnitude of the ETASI model (curve) with the observed value for non-ov erlapping consecutiv e bins of 25 m ≥ 2 earthquakes (crosses).The mean magnitude of the observed aftershocks is strongly time-dependent, with the highest values immediately after the mainshock and the largest aftershock, as visually expected from panel (b).A further jump in the mean magnitude is observed at the time of the M 6.4, which occurred about 14.7 d after the mainshock.The ETASI model well fits the observed trend.The OK93 model also follows the observed data very closel y, onl y missing the peak after the M 6.4 e vent, probabl y due to the smoothness constraints.
We compared the log-likelihood values for the observed magnitudes using (i) a constant GR distribution with b = 0.56 (Table 1 ), (ii) eq. ( 5) for ETASI with constant parameters (Table 1 ) but timevarying rate R 0 ( t ), and (iii) eq. ( 8) for OK93 with time-varying b , μ and σ .The results in Table 2 show that both OK93 and ETASI clearly outperform GR with an Information Gain Per Event (IGPE) of 0.17.Note that the applied IGPE value is not corrected for the number of model parameters.The OK93 model, which has many free parameters and flexibility in fitting the time e volution, gi ves the best result, but only with a small IGPE of 0.002 with respect to ET ASI.The ET ASI has much less freedom to model the evolving incompleteness with time because only one additional parameter ( T b ) is introduced to model the time evolution of the recorded earthquake rate and magnitude distribution simultaneously.Thus, the simultaneous good fit of the earthquake rates and magnitudes demonstrates the consistency of the model with the data.
A more detailed analysis of the time-dependent frequencymagnitude distribution is shown in Fig. 3 , where the histograms of the observed magnitudes are compared with the ETASI and OK93 results in non-overlapping successive time windows, each containing 500 m ≥ 2 aftershocks, except for the first interval containing the less numerous foreshocks.Here, the ETASI curves are calculated by summing eq. ( 5) evaluated at the occurrence times of the obser ved ear thquakes, and the shaded area refers to the corresponding 90 per cent confidence interval.The OK93 curves are computed similarly using eq.( 8).Overall, the model results are in good agreement with the observ ations.Howe ver, for earl y aftershocks, the ETASI slope at high magnitudes tends to underestimate the slope of the large-magnitude tail, which is related to an underestimation of the b value, as shown next.
Fig. 4 (a) shows the b value for the preceding activity that occurred in the same region between 2000 and 2022 as a function of the cut-off magnitude using MLE (Aki 1965 ).The b estimate becomes stable within its uncertainties for M c > 3.5 and scatters around 1.07.Ho wever , the OK93 model indicates that the b value is strongly time-dependent.In 2023, OK93 gives b = 0.76 ± 0.07 for the earthquakes that occurred before the mainshock and b values of about 1.2 for the early aftershocks, decaying to about 0.85 for the later aftershocks (see Fig. 4 b).The b value estimated by OK93 for aftershocks occurring 10 d or later after the mainshock agrees with the constant value estimated by ETASI.
Using OK93, we observe a 50 per cent increase in the b value coinciding with the mainshock event.This coseismic rise in b value generally aligns with findings by Gulia et al. ( 2018 ), who, in studying 58 aftershock sequences, noted a 20-30 per cent increase in b value following mainshocks.Changes in b values are commonly interpreted as indicators of stress variations, as rock experiments have demonstrated a decrease in b values of acoustic emissions during rock fracturing with heightened differential stress (Scholz 1968(Scholz , 2015 ) ).Assuming a consistent mechanism for both mainshocks and aftershocks, the average shear stress in the aftershock zone decreases due to the mainshock (Gulia et al. 2018 ;Sharma et al. 2023 ).When combined with the laboratory-deriv ed inv erse relationship between stress change and b value, this might be mistakenly used to explain the increased b value in aftershocks.Ho wever , despite the decrease in average shear stress, aftershocks occur in areas where stress increases locally, necessitating a positive correlation between b values and stress changes to explain the observation.Such a positive relation was observed by Sharma et al. ( 2023 ), who conducted a comprehensi ve anal ysis of the b -v alue dependence on stress changes induced by mainshock events across 127 mainshock-aftershock sequences, suggesting a crucial role of structural heterogeneity and strength v ariations.Ne vertheless, the rapid decay of the b value in our case remains puzzling.Although previous studies have reported decays on a timescale of years (Tormann et al. 215 ;Gulia et al. 2018 ), the b value, in our case, returned almost to its pre-mainshock value ( b ≈ 0.8) within 10 d.
The time dependence of the b value is not considered by the ETASI model, which assumes a constant b for simplicity.The significantly larger b value for early aftershocks (according to OK93) explains the ETASI misfit of the magnitude tail of these events shown in Fig. 3 .However, the simultaneous estimation of the timedependent b value and the catalogue completeness is challenging and may only be reliable with sufficient data.
For forecasts based on limited input data, assuming a constant value of b may be the only way to obtain a robust estimate.For this reason, we tested estimates assuming constant b from the beginning of 2023 until the time t indicated on the x -axis of Fig. 4 (b).The ET ASI estimate (ET ASI( t )) is very close to the value estimated for the whole period (dashed red line).We also applied the b -positive method, which w as pre viousl y shown to be a robust b -value estimator for catalogues with time-varying completeness (van der Elst 2021 ).In this case, only the difference m between the magnitudes of subsequent events is considered, and the b value is estimated 0.56 −6897.90 − 0 .17OK93 eq. ( 8) b ( t ) −5295.0 0.17 0 .002ETASI eq. ( 5) 0.87 −5311.90.17 0   systematically increases with increasing m c .For m > 0.7, the estimates are slightly below 1.07, while for m > 1.0, they are slightly above 1.07.The strong dependence of the b -positive results on m c indicates that a stable estimation of b is difficult for this data set.
The results show that the true b value was time-dependent and w as, for earl y aftershocks, similar to the b v alue of the background activity in the period 2000-2022, and thus significantly larger than the value estimated by ETASI assuming a constant b for the preshocks and aftershocks.

Pseudo-pr ospectiv e forecast experiment
To analyse the possibility of forecasting the number of aftershocks and their maximum magnitude based on limited and incomplete catalogue infor mation, we perfor med a retrospective forecasting test.Specifically, we computed forecasts for the next day starting at different start times T s , using earthquakes up to T s as input for parameter estimation.The first forecast started immediately after the mainshock, using parameters estimated for m ≥ 2 earthquakes that occurred between one month before and after the M L 6.8 Elazi g earthquake occurred on 24 January 2020, in the same spatial region (see the yellow point in Fig. 1 ).All subsequent forecasts use only activity within 2023 as input for the parameter calibration to allow for more localized and sequence-specific parameter estimates.
The second forecast starts immediately after the largest aftershock (0.38 d after the mainshock) and the third after 0.5 d.Subsequent forecasts start at one to 30 d after the mainshock, in 1-d increments.
Table 3 and Fig. 5 show the parameters resulting from the MLE of the ETASI model for the different starting times.In Fig. 5 , the result for the initial estimation is shown in grey, while the results for the data fit after the 2023 mainshock are colour-coded according to T s .The latter parameter estimates are found to be quite stable from the start, despite the high level of incompleteness of the input data.The coefficient of variation (the ratio between the standard deviation and the mean) is less than 6 per cent for r , α, p , T b and b , and larger only for K (21 per cent) and c (39 per cent).Ho wever , due to the positive correlation between c and p and the anticorrelation between K and α, the temporal decay (normalized to the first 100 d, Fig. 5 b) and especially the aftershock productivity as a function of the earthquake magnitude (Fig. 5 c) are very stable.The estimate for the 2020 Elazi g sequence differs mainly in its lower scaling of the aftershock productivity with earthquake magnitude, corresponding to a lower α value of 1.05 compared to the values ranging between 1.2 and 1.3 for the 2023 data.
Our results show that the blind time T b is rather stable and estimated in the range 142-163 s.Thereb y, T b slightl y increases with T s with the smallest T b for the smallest data sets, which include the largest events ( T s = 0.38 and 0.5 d).Therefore, the result does not indicate an increase of T b with the event magnitude, which might be explained by the fact that detection algorithms usually involve a fixed time window; such as in the case of the classical shor t-ter m-average to long-term-average ratio (ST A/LT A) picker, which typically uses LTA windows around 30 s (Earle & Shearer 1994 ).During this period, the arri v al of seismic waves related to a preceding e vent pre vents the picking of a smaller magnitude event within the STA window.The distance-dependent delayed arri v al of different phases of the preceding events and path effects (scattering of the wave field) can further extend the ef fecti ve blind time.These aspects do not directly depend on the magnitude and will set the lower limit of the blind time for small events.A magnitude dependence of the blind time might be related to the source duration.Considering the scaling of the rupture length L with magnitude and typical rupture velocities, L = 10 −2.44 + 0.59 M (Wells & Coppersmith 1994 ), subsurface rupture length, all mechanisms) and assuming a rupture propagation of 2-3 km s −1 , this leads to a source duration of 40-60 s for M = 7.7, about 5 s for M = 6, and about 1 s for an M = 5 event.Thus, the effectiv e e xtension of the blind time due to the magnitude effect may only be rele v ant for M > 7 events, and a constant blind time is a good overall approximation.
Each 1-d forecast is made b y e v aluating 1000 Monte Carlo simulations using the estimated ETASI parameters and the preceding activity in 2023 as inputs.Specifically, we use the inverse transform method for our Monte Carlo simulations of the ETAS model, where interevent times and magnitudes are randomly chosen from inhomogeneous Poisson and GR distributions (see, e.g.Felzer et al. 2002 ).
Here, the magnitude of each sampled earthquake in the forecast period is randomly selected from the doubly truncated GR distribution in the range [1.95, 8.0] with given b .Because the simulations are compared with incomplete observ ations, onl y the detectable events are filtered in the output catalogues.Events following any larger event within the estimated value of T b are removed according to the ETASI assumption.In particular, for each simulated event ( t i , m i ), all subsequent events with m j < m i and t i < t j < t i + T b are marked as unobservable and finally removed.The filtered output catalogues are then e v aluated regarding the number of m ≥ 2 events and the maximum magnitude in the forecast period.
Due to the aforementioned problems in estimating b , alternative forecasts are performed.The first one uses the b value estimated by ETASI for each input data set.In contrast, the preferred second approach assumes a constant b = 1.07 for all forecasts, that is the b value estimated for the background activity between 2000 and 2022.We also analysed the forecasts using the b -positive estimates of early aftershocks within [0, T s ] for T s > 0 or the background b value for T s = 0.For each forecast, the mean value and the 90 per cent confidence interval of the simulations are then compared with the obser vations concer ning the number of m ≥ 2 earthquakes and their maximum magnitude.To account for parameter uncertainties of the ETAS model, we used Bayesian inference employing Markov Chain Monte Carlo (MCMC) for sampling the posterior parameter distribution (see e.g.Shcherbakov et al. 2019 ).Here, we applied the Metropolis-Hastings algorithm with flat priors and Gaussian proposals with a standard deviation of 0.05 of the maximum loglikelihood value.For each forward simulation, a set of parameters is randomly selected from the full posterior distribution.Ho wever , the results indicate that the uncertainties are mainly driven by the aleatoric uncertainties (intrinsic randomness of the magnitudes and trigger times), and the epistemic uncertainties are rather small; that is, Monte Carlo simulations with MLE lead to confidence intervals, which are only slightly smaller.
Fig. 6 shows the results of the forecasting experiment.The predicted numbers of detectable m ≥ 2 aftershocks are found to be al wa ys close to the real ones, even for the first forecasts where ETASI parameters are only estimated based on background events or pre-shocks with no or very limited and incomplete aftershock information.The forecasts with the different b value approaches yield similar results, and the forecasts with the background b value or b -positive estimates seem to be only slightly better.Ho wever , this changes significantly for the maximum magnitude ( M max ) forecast.
Here, the forecasts based on the underestimated b values generally yield significantly overestimated M max values.In contrast, the forecasts with fixed b = 1.07 and the similar forecasts with b -positive estimates agree with the observation in 29 out of 32 cases within the 90 per cent confidence interval.Note that three outliers are, on av erage, statistically e xpected in the case of the chosen confidence interv al. Howe ver, the mismatch of M max forecast for the first day may also be related to a sampling bias since the occurrence of the largest aftershock with a magnitude difference of only 0.1 is known to be rare.For the first day (and two other days), the forecast with underestimated b values fits better.

D I S C U S S I O N A N D C O N C L U S I O N
Early aftershock forecasting can be critical for risk reduction and response efforts.Due to the stochastic nature of the earthquake nucleation and propagation, statistical forecasts that provide probabilities and expectation values for earthquake occurrences are the best we can do so far.The standard ETAS model accounts for OUtype aftershock decay and secondary aftershock triggering, where the model parameters must be adjusted by data fits.Ho wever , earthquake catalogues are highly incomplete in the first hours to days after a large event.Therefore, it is important to address the issue of catalogue incompleteness when fitting model parameters.
Using a time-dependent threshold value, such as the empirical relationship derived for California (Helmstetter et al. 2006 ), is not recommended due to the influence of STAI even for larger events.Earthquakes are also missed in the upper magnitude range as indicated by the high level of the μ + 2 σ curve (OK93, Fig. 1 a) and the M c, 90% curve (ETASI, Fig. 2 b).Thus, parameter estimations must take into account the detection probability as a function of the earthquake magnitudes.The most sophisticated approach commonly used is first to estimate the time dependence of the detection probability using OK93 and then fit the OU or ETAS model (Omi et al. 2013(Omi et al. , 2016 ; ;Page et al. 2016 ).Ho wever , ETASI provides a simpler, straightforward and closed-form formulation that ef fecti vel y accounts for catalogue incompleteness.ETASI assumes a blind time, denoted by T b .While this assumption may not be appropriate for semi-manual earthquake catalogues that lack consistent processing over time, it may be appropriate for automatically processed catalogues, such as those available in near-real time after a mainshock.
Application of ETASI to the 2023 earthquake sequence in T ürkiye demonstrates its ability to accurately fit the data in terms of the number and mean magnitudes of the aftershocks (Fig. 2 ).The temporal evolution of the full magnitude distribution is also generally well reproduced, indicating that the rate-dependent detection probability is appropriate.Only small deviations in the large magnitude tail indicate an underestimation of the b value for the early aftershocks.Our analysis with the OK93 model confirms that the b v alue w as strongl y time-dependent, with v alues less than 0.8 before the pre-shocks, then jumping to values of about 1.0 and 1.2 after the mainshock and the largest aftershock, respecti vel y.Then, within the first 10 d of the aftershock sequence, the b value returned almost to the pre-mainshock level, with values scattering around 0.85 after that.The latter value is consistent with the ETASI estimate, assuming a constant b value.
Our retrospective forecast experiment shows that the daily number of m ≥ 2 earthquakes can be well predicted based on the limited and highly incomplete information from early aftershocks.However, the maximum magnitude of the earthquakes within the next 24 hr would have been overestimated using the ETASI-estimated b value.This overestimation is probably related to the time-dependent and significantly increased b values for early aftershocks.The M maxforecasts are significantly improved when using the b value estimated for the regional earthquake that occurred in the preceding 2000-2022 period or, similarly, using b -positive estimates with high m c for early aftershocks.Thus, it is crucial to use an appropriate estimation of the b value to accurately predict M max , and a promising future approach would be to combine ETASI with a timedependent estimation of b to improve the forecasting capabilities further.

A C K N O W L E D G M E N T S
We are grateful to Robert Shcherbakov and the two other anonymous re vie wers for their thoughtful comments and suggestions, which helped us to improve this paper.SH benefited from the RISE project funded by the European Unions 2020 research and innov ation pro g ramme under g rant ag reement number 8211185, and TK and YO were supported by MEXT Project for Seismology towards Research Innovation with Data of Earthquake (STAR-E) grant number JPJ010217.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The earthquake data can be downloaded from https://deprem.afad.gov.tr/event-catalogue .We downloaded the 2023 data on 8 March 2023 and the 2000-2022 catalogue on 21 April 2023 in csv format.

Figure 1 .
Figure 1.Observed seismicity: (a) Magnitudes versus time of earthquakes recorded in 2023 in the region defined by the boundaries of the map (b).The solid and dashed curves refer to the detection probability of the OK93 model (eq.7 ) using μ and μ + 2 σ , and the horizontal dotted line refers to the applied magnitude cut-off of M c = 1.95 for ETASI.In (b), m ≥ M c earthquakes are plotted, with events occurring between 2000 and 2022 shown in grey and 2023 events are shown in blue.In both plots, the symbol size is scaled by the earthquake magnitude.The yellow circle marks the epicentre of the M L 6.8 Elazi g earthquake in 2020.

Figure 2 .
Figure 2. Model fits as a function of logarithmic time after the mainshock: (a) Earthquake rate of m ≥ 2 earthquakes, where the crosses refer to the observed rates and the lines indicate the mean rates for the optimized models: standard ETAS (solid grey), detectable ETASI ( R , solid red) and actual ETASI ( R 0 , dashed red).(b) The detected event magnitudes along with the completeness levels estimated by OK93 and ETASI, where the solid curves (green: OK93 and red: ETASI) refer to the magnitudes detected with a probability of 50 per cent, while the dashed red curve refers to a detection probability of 90 per cent according to ETASI.Also shown is the empirical completeness-relation of Helmstetter et al. ( 2006 ) (black dashed line).(c) Mean earthquake magnitudes: Observations computed in non-overlapping bins of 25 events are shown by error bars, where the horizontal bar refers to the time period of the events used, and the vertical bar refers to ±1 standard deviation assuming an exponential function with the observed mean value.The OK93 estimate is shown by the green curve (using eq. 8 ), while the solid red curve shows the mean value related to eq. (5).

Figure 3 .
Figure 3.Comparison of the observed (bars) and expected (lines) frequency-magnitude distributions in the subsequent periods indicated in the upper right of each panel, each containing 500 m ≥ 2 aftershocks (except for the foreshock window).The green curves refer to the OK93 model, while the red curves show the ETASI results with the 90 per cent confidence interval indicated by the grey shaded area.

Figure 4 .
Figure 4. b Value: (a) MLE of b as a function of the cut-off for the regional earthquakes between 2000 and 2022.The grey horizontal line refers to b = 1.07, and the blue and grey error bars refer to the traditional and b -positive methods, respectively; both are slightly offset for better visibility and lead to similar results.(b) MLE of ETASI for b using the recorded m ≥ 2 earthquakes up to the time t indicated on the x -axis (ETASI( t ), red dashed line) or the entire sequence (ETASI, red dotted line).The results for b -positive for m ≥ 2 aftershocks up to time t are also shown for three different thresholds m c .Unlike the other results, the OK93 model estimates are localized in time and shown by the green curve with its 2 σ confidence interval.The OK93 result indicates that the value of b was close to the background level in the first few days and then decayed to the value estimated by the ETASI model.onlyfor those that are positi ve.Theoreticall y, this estimator is unbiased for m ≥ 0 even for large variations in the completeness magnitudes if (i) the completeness magnitude does not increase between earthquakes and (ii) no earthquakes are recorded below the completeness magnitude.Ho wever , since a large fraction of events

Figure 5 .
Figure 5. (a) ETASI parameters as a function of the forecast start time T s , where the estimates are based on MLE using the input catalogue up to T s .(b) Normalized Omori-Utsu function and (c) the productivity relationship associated with these parameters, where colours refer to time as defined in (a).

Figure 6 .
Figure 6.Results of the pseudo-prospective forecast experiment: (a) number and (b) maximum magnitude of m ≥ 2 earthquakes within the next day as a function of the forecast starting time T s .In both panels, the black points refer to the observations and the curves with their 90 per cent-confidence intervals refer to the forecasts.Grey results refer to simulations based on the b value estimated by ETASI, while simulation results with the b value estimated from background events are shown in red.Dashed lines refer to the mean forecasts using b -positive estimates of early aftershocks within [0, T s ] for T s > 0 and the background b value for T s = 0. .

Table 3 .
ETASI parameters used for the daily forecasts start at time T s relative to the mainshock, where the MLEs have been performed for the period [ T 1 , T 2 ] consisting of N recorded m ≥ 2 events.