Earthquake forecasting from paleoseismic records

Domain: earth_science
Published: 2024-03-02
Editorial outcome: minor_revision

Abstract
Forecasting large earthquakes along active faults is of critical importance for seismic hazard assessment. Statistical models of recurrence intervals based on compilations of paleoseismic data provide a forecasting tool. Here we compare five models and use Bayesian model-averaging to produce time-dependent, probabilistic forecasts of large earthquakes along 93 fault segments worldwide. This approach allows better use of the measurement errors associated with paleoseismic records and accounts for the uncertainty around model choice. Our results indicate that although the majority of fault segments (65/93) in the catalogue favour a single best model, 28 benefit from a model-averaging approach. We provide earthquake rupture probabilities for the next 50 years and forecast the occurrence times of the next rupture for all the fault segments. Our findings suggest that there is no universal model for large earthquake recurrence, and an ensemble forecasting approach is desirable when dealing with paleoseismic records with few data points and large measurement errors.

Introduction
Large earthquakes are one of the most devastating natural hazards, not only because they cause a significant number of casualties, but also because they lead to widespread infrastructure damage. Researchers from multi-disciplinary backgrounds have been endeavouring to understand their underlying mechanisms and forecast their occurrences. Elastic rebound theory1 forms the basis of the standard earthquake cycle of strain accumulation and release. This theoretical basis, combined with analysis of long-term earthquake records, support the proposition that the recurrence times of large earthquakes along the same fault can be modelled with renewal processes2–6. A renewal process is a statistical model that describes event occurrences in time, treating each new occurrence as a renewal in which the system is reset. It assumes that the times between events (the inter-event time, such as the time between two consecutive earthquake occurrences) are independent and identically distributed. Probabilistic Seismic Hazard Analysis7 frequently uses the Poisson process, which is a special type of renewal process that describes events effectively occurring randomly in time. Recent analysis of global paleoseismic records nevertheless suggests that large earthquake recurrence on individual fault segments is typically more periodic than expected from a Poisson process8–10. Such periodic patterns can be captured by other renewal processes, such as the Gamma renewal process. An alternative to renewal processes for modelling paleoseismic records uses a Long-Term Fault Memory Model11, which assumes that the timing of future events is dependent on not only the time elapsed since the most recent event (as in renewal models) but also the previous inter-event times. However, when considering global dataset and synthetic earthquake records there appears to be no significant correlation or anti-correlation between successive inter-event times for the vast majority of the earthquake records9,12. Therefore, in this study we only focus on renewal processes.

Previous earthquake forecasts using paleoseismic data were mainly based on one single best-model selected from a number of candidate models according to an information criterion, such as the Akaike Information Criterion (AIC13). This approach ignores the fact that there can be model uncertainty, i.e. uncertainty as to which model is best. If the AIC values for several models are not very different, it may be better to use a combination of the forecasts from those models, rather than rely solely on the forecast from the model with the lowest AIC. In real-world applications such model uncertainty may have a dominant effect over uncertainty in the estimates of the parameters in each model14,15. This becomes a serious limitation when the sample size is relatively small (only 14 out of 93 earthquake records considered in this study have more than 10 events). Failing to acknowledge model uncertainty by selecting a single best model may produce erroneous or unrealistically precise forecasts.

In this study, we use five candidate renewal processes to investigate the recurrence patterns of large earthquakes using paleoseismic records from 93 worldwide fault segments previously compiled in Griffin et al.9 and supplemented with additional records from other publications (See Data availability section). We use Bayesian model-averaging to carry out probabilistic forecasts of future large earthquakes to account for model uncertainty. We compare the forecasts from each single model with that from model averaging. We find that there is no single best model that universally describes the recurrence of large earthquakes for the 93 fault segments considered here, nor for fault segments with the same faulting styles, from the same tectonic region, or even within the same fault system. We provide the distribution of the probability that at least one large earthquake will occur in the next 50 years along each fault segment. We also carry out a leave-one-out test to compare the performance of the model-averaged forecast and the Poisson forecast, as the latter is frequently used in national hazard models7.

Results
No universal model for large earthquakes
Modelling of large-earthquake recurrence times from paleoseismic records typically either uses the BPT renewal process, because of its physical explanation of the earthquake process4,6,16, or selects the best model among a few candidate renewal processes, based on an information criterion3,17. As discussed above, inference based on a single best model ignores potential model uncertainty.

In order to minimise the impact of model uncertainty we use a model-averaging approach and consider geological and historical records of large-earthquake occurrence times from 93 worldwide fault segments (Table 1). Here, we consider a “fault segment” as a section of a fault that is recognised as being geometrically distinct, and that is likely to define the boundaries of at least some earthquake ruptures on the fault. Some models, such as the third Uniform California Earthquake Rupture Forecast (UCERF316) and the 2022 New Zealand National Seismic Hazard Model18, relax the strict fault segmentation assumption in their earthquake probability model component by considering interactions between different fault subsections. Nevertheless, for each individual fault subsection the fundamental elastic-rebound theory part of the model uses a Brownian passage-time (BPT) renewal model. Here we focus on the elastic-rebound theory part of the earthquake probability model component, and thus do not consider interactions between fault segments. Once this layer of modelling is improved, interactions between subsections can be considered in order to forecast large earthquakes that may rupture multiple fault subsections and/or faults. Improved forecasts from this layer will strengthen a holistic hazard model that includes fault models, deformation models, earthquake rate models, and earthquake probability models.

We select fault segments with records of at least five large earthquakes in order to be able to fit two-parameter models to the inter-event times. We only consider the occurrence times of earthquakes that left detectable geological evidence9. Earthquakes that leave a geological signature are assumed to be large and significant for seismic hazard, but we do not explicitly consider their magnitudes, or the magnitude distribution (i.e. characteristic vs Gutenberg-Richter) in our model. The number of large earthquakes in the paleoseismic record of each fault segment is small, with a maximum of 35 events for the Chile Megathrust19, only 3 fault segments having more than 20 events, and 14 having more than 10 events. The measurement errors associated with most dated earthquake ages are large, resulting in 1\sigma uncertainties of 25±12% for inter-event times in the data considered here. To maximise the reliability from the dataset of all event ages and to make good use of the measurement errors, we simulate 100 Monte Carlo (MC) samples from the empirical distribution of the occurrence times for each fault segment provided in the literature. If any earthquake has more than one dated age published in the literature, we take the average of the different age distributions and sample from this averaged distribution.

The 100 MC samples of occurrence times for each fault segment form each dataset. We fit five different renewal processes to each dataset, including the Poisson process, and the Gamma, Weibull, BPT and lognormal renewal processes (see Methods). In each model, we allow the parameters to vary for different MC samples from the same fault, which captures the similarities between the 100 MC samples. We use a Markov Chain Monte Carlo (MCMC) algorithm to generate samples from the joint posterior distribution of all the model parameters for each model. These posterior samples of parameters are then used to simulate future earthquake occurrence times from each fault, thereby providing us with forecasts from each model.

For each model fitted to each fault segment, we calculate the Watanabe-Akaike Information Criterion (WAIC20), a measure of the prediction performance of a model in Bayesian analysis. For model-averaging, we calculate the WAIC weight15 for each model for all 93 fault segments (see Methods), as shown in Fig. 1a. Model-averaged forecasts of future earthquakes for each fault are obtained by combining the forecasts from each individual model using these WAIC weights (see Methods). For a fault segment with any single model having weight ≥0.95, the model-averaged forecast is very close to the forecast from that model (i.e., the single-best model), as we would expect. The Weibull model best fits 41 (44.1%) fault segments, the Gamma model best fits 5 (5.4%) fault segments, the BPT model best fits 4 (4.3%) fault segments, and the lognormal model best fits 15 (16.1%) fault segments. The Poisson model has no weight greater than 0.7. The remaining 28 (30.1%) fault segments have weights <0.95 for each single model, suggesting that model-averaging will be most beneficial for those fault segments.

WAIC weights provide a numerical comparison of the amount by which a model is better at prediction than another, as they show how much weight should be given to the prediction from each of these models when calculating a model-averaged prediction. Based on the WAIC weights in Fig. 1a, we can see that the predictive performance of the Weibull model is not uniformly better than the others. For 33 fault segments, the WAIC weight for the Weibull model is close to 0, which suggests that for these fault segments the estimated predictive accuracy of the other models is much better than that of the Weibull model. Among the 11 San Andreas Fault segments we find three different best-fit models. Although they all suggest that every segment exhibits quasi-periodic behaviour (see discussion later), segments with larger sample sizes were best fit by a Weibull model. It is unclear if this is simply an outcome of sampling or due to real differences in recurrence behaviour between the different segments, including how neighbouring fault segments interact. Variability in observed earthquake recurrence behaviour at paleoseismic sites on the San Andreas Fault has been proposed to be (at least partially) due to overlap of ruptures occurring on neighbouring segments21 and this has been supported by studies using earthquake simulators22. In contrast, relatively strong quasi-periodic recurrence on the Alpine Fault has been attributed to its geometric simplicity and relative isolation from other faults23, although its fault geometry has also been invoked to explain the variability in earthquake inter-event times24. The persistence (or otherwise) of rupture barriers between segments25 may also be a significant factor controlling the distribution of inter-event times observed on a fault segment. Therefore, because of the limitations of the available data, even within a single well-studied fault system, we cannot use a universal single best model, and it is not clear that one exists.

Best-model frequencies in Fig. 1b, c indicate that the Weibull renewal process is more likely to best reproduce recurrence times for fault segments with longer records compared to the BPT renewal process. The fault segments that have the Weibull renewal process as the single-best model have on average about 1.59 (95% CI 1.10–2.47) times more large earthquakes than those for which the BPT renewal process is the single-best model. The fault segments that have the BPT renewal process as the single-best model all have fewer than 10 events in the paleoseismic records. When we removed the last event in each record, in order to carry out retrospective forecasts (see section Assessment of Prediction Error), the fault segments with more than 15 events all have the Weibull renewal process as the single-best model; only one fault segment with 12 events has the BPT renewal process as the single-best model. It is therefore difficult to determine whether the dependence of the favoured model on sample size is due to real differences in the best model or simply an outcome of the small amount of data.

Previous studies26–28 have shown that the standard deviations of the scaled (divided by the mean value) inter-event times along several fault segments appear to be constant. For each of the 93 fault segments, we calculated the scaled inter-event times for each MC sample, and then reported the median and the 2.5% and 97.5% quantiles of the standard deviations of the scaled inter-event times calculated for the 100 MC samples. Fault segments with higher rates of earthquake occurrence appear to have smaller standard deviations of the scaled inter-event times (Supplementary Fig. 1 (e)). The median standard deviations of the scaled inter-event times for the 41 fault segments that were best fit by a Weibull model are all smaller than 0.8 except for one case. In contrast, more than half of the median standard deviations of the scaled inter-event times for the 15 fault segments that were best fit by a lognormal model are larger than 0.8 (Supplementary Fig. 1a, d). The large median standard deviations are related to fault segments that have inter-event times with long-tailed distributions which result in a large spread of data values. These long-tailed distributions are best fit by a lognormal renewal process.

Probabilistic forecasting
The estimated probability that at least one large earthquake will occur in the next 50 years along each fault segment is shown in Fig. 2 and Table 1. The highest probability of a large earthquake in the next 50 years appear to be along the Carrizo and Wrightwood segments of the San Andreas Fault, and the Mentawai Segment of the Sumatra megathrust, with median probabilities over 50%. The upper limit of the 95% credible interval (CI) of this probability exceeds 75% for the Wrightwood segment of the San Andreas Fault. These results are consistent with previous suggestions that these faults are in or approaching the beginning of a new seismic cycle29–33. Other sites with high probabilities of rupture (>20%) in the next 50 years include the Alpine Fault, the Chile Margin, Hayward Fault, most other segments of the San Andreas Fault, the San Jacinto Fault, and certain segments of the Dead Sea Transform and North Anatolian Fault (Table 1). If one fault segment has a high probability of rupture, it may affect the probability of rupture at a neighbouring fault segment. This is not built in our methods for the global data analysis, but can be considered for a particular local fault system. The lowest probabilities of a large earthquake in the next 50 years are along the Kusaihu segment of the East Kunlun Fault, the Futugawa Fault, the Eastern Piedmont segment of the Helanshan Fault, and the Laohushan segment of the Qilianshan Fault, with the upper limit of the 95% credible interval of these probabilities all being lower than 0.001%. On average, faults along a plate boundary are about 32 (95% CI 9–125) times more likely to have a large earthquake in the next 50 years than intraplate faults.

Model-averaging and best-model approaches can give significantly different forecast probabilities for an earthquake occurring within the next 50 years. The maximum difference in the forecast probability between the two approaches is about 14% for the Thousand Palms segment of the San Andreas Fault, where the model-averaging approach gives a higher probability than the best-model approach (Table 1). Having said that, for about 90% of the fault segments, this difference is less than 0.1%. When considering the forecast occurrence times, for the 28 fault segments where the single-best model has a WAIC weight less than 0.95, the two approaches give quite different results. The maximum difference between the credible intervals obtained using the model-averaging and single-best model approaches is over 10,000 years, and the difference is over 50 years for 43% of the 28 fault segments that benefited from a model-averaging approach. Figure 3 shows the forecast occurrence time of the next large earthquake along each fault using the model-averaging approach, while Table 1 lists the probabilities of a large earthquake occurring within the next 50 years. A comparison of the forecasts from different models is in the supplement (Supplementary Figs. 2–4).

A quantitative comparison of our probabilistic forecasts with those published in other studies is challenging (Supplementary Data 1). Reports of the mean recurrence interval estimates are equivalent to forecasting the next occurrence time using a Poisson process (this comparison is in Supplementary Figs. 2–4). A meaningful comparison is difficult when the forecast of the probability of the next earthquake occurrence is specified in terms of a start date and a fixed time period, as these vary between studies. For example, UCERF3 provides 30 year probabilities from 2014 CE aggregated by parent fault section16, whereas we provide 50 year probabilities from 2022 CE on each fault segment. The UCERF3 forecasts for the San Andreas fault segments all have similar mean values and overlapping ranges, whereas our forecast probabilities show greater variability between segments, with several non-overlapping credible intervals (Supplementary Fig. 5). We found previous forecasts for another eight fault segments in our study. All except two of our forecasts overlap within uncertainties with previous values, although our uncertainty bounds are typically narrower (Supplementary Fig. 5). Having said that, UCERF3 provides the ranges of forecast probabilities rather than 95% confidence intervals.

Clustering or periodic behaviour
For a Gamma or Weibull renewal process, if the estimated shape parameter \alpha < 1, then the process tends to have clustering behaviour, while if \alpha > 1, the process tends to have quasi-periodic behaviour, with larger \alpha suggesting more periodic behaviour. For a BPT renewal process with probability density function as defined in Eq (3), smaller \beta (for \beta < 1) suggests a more symmetrical probability density function of inter-event times and hence more periodic behaviour. A lognormal renewal process describes quasi-periodic behaviour, with a smaller standard deviation \sigma corresponding to more periodic behaviour.

Based on the parameter estimates from the Gamma and Weibull renewal processes, five fault segments appear to show clustering behaviour, with the upper 95% credible limit of the shape parameter from each of these models being less than 1. These are Cadell, Dunstan, Lake Edgar, Solitario Canyon, and Waitangi, all of which have low earthquake occurrence rates. For Waitangi, the BPT renewal process has a WAIC weight over 0.95, and the estimates of \beta are over 2.5, confirming a clustering behaviour. For Lake Edgar, the Gamma renewal process has a WAIC weight over 0.95. Six fault segments appear to demonstrate near Poisson behaviour, Dead Sea Beteiha, Dead Sea Qatar, Dead Sea Taybeh, Langshan Piedmont Xibulong East, Reelfoot, and Wharekuri, with the 95% credible interval of the shape parameter for both the Gamma and Weibull model containing 1. The remaining 82 fault segments show quasi-periodic behaviour. These results are consistent with a previous study9 that used different tests to check quasi-periodic recurrence behaviour.

A regression model for the relationship between the shape parameter of a Weibull renewal process and the earthquake rate, tectonic setting, faulting type, and the number of earthquakes of each of the 93 fault segments (Fig. 4) suggests that when the earthquake rate increases 10 fold (e.g., from 0.0001 to 0.001), the shape parameter of the Weibull renewal process increases by about 24% (95% CI: 5%–45%). This suggests that fault segments with higher earthquake rates tend to have more periodic behaviour. The shape parameter of the Weibull renewal process for fault segments located in a stable continental intraplate setting is about 87% (95% CI: 51%–148%) of that for fault segments located at or near a plate boundary. Although the 95% CI is wide and covers 100%, the posterior density plot (Supplementary Fig. 6) suggests that there is a high probability (0.7) that the latter may be more periodic than the former. There are only five fault segments from a stable continental intraplate setting, so to reach a more robust conclusion, more data from intraplate fault segments are needed. These findings are consistent with and nuancing those from past studies9,10. The shape parameter of the Weibull renewal process for fault segments located in an active intraplate setting (predominantly faults in central China) is about 1.4 (95% CI: 1.1–1.9) times that for fault segments located at or near plate boundary, suggesting that the former appear to be more periodic than the latter. The shape parameter of the Weibull renewal process for fault segments located in a subduction region is about 1.9 (95% CI: 1.1–3.3) times that for fault segments located at or near plate boundary, suggesting that the former appear to be more periodic than the latter. On average, the Weibull shape parameter for reverse faults is about 69% (95% CI: 46%–104%) of that for normal faults, suggesting that large earthquakes recur more periodically on normal faults than on reverse faults. Note that obtaining long earthquake records for reverse faults is often more difficult than for normal faults, due to progressive burial of the evidence for previous earthquakes. With more data becoming available in the future, one could investigate whether variations in strain rate or kinematics have a significant impact on record length and event preservation.

The likelihood of systematically missing very short inter-event times in the paleoseismic records and hence biasing our analysis is low. If fault re-rupture commonly occurred shortly after previous earthquakes, then we would expect to see this frequently in the historical record of surface rupturing earthquakes, which we largely do not34. While preservation of paleoseismic evidence is an important consideration in the interpretation of any paleoseismic record, the characteristics of the geomorphic and geological setting, and the relative rate of tectonic processes to non-tectonic geomorphic processes, will control whether evidence of past earthquakes is preserved35. Ideally, paleoseismic studies consider these factors in site selection and interpretation; while it is acknowledged that missed events are a possibility, it is not clear that this should lead to a systematic bias in our statistics towards or away from periodicity.

Assessment of prediction error
For each fault, we removed the last event in the record in order to carry out retrospective forecasts using both the model-averaging and single-best model approaches. The single-best model remained the same as that for the full dataset for 49% of the fault segments (46 out of the 93). It appears that the single-best model is more likely to change for fault segments with fewer recorded earthquakes. The number of earthquakes along fault segments for which the single-best model changed is about 27% (95% CI: 11%–39%) fewer than for those fault segments for which the single-best model remained the same. This demonstrates the large model uncertainties for paleoearthquake data, again suggesting that a model-averaging approach is preferable.

Figure 5 shows the 95% credible intervals of the forecast of the last earthquake occurrence time, with 0 representing the mean of the recorded last earthquake occurrence time. Out of the 93 fault segments, the model-averaged forecast successfully covered 79 of the mean true occurrence times. The forecasts from the Poisson process successfully covered 89 of the mean true occurrence times, which at first glance may suggest that it outperforms the model-averaged forecasts. However, it does this by having much wider credible intervals (on average twice as wide as the credible interval from the model-averaged forecast); i.e. it has much more uncertainty in the forecast. When examined in more detail, we see that model-averaged forecasts routinely outperformed Poisson process forecasts, provided that there were sufficient events left in the record. Specifically, about 83% (77 of the 93 fault segments) of the model-averaged forecasts have much smaller mean squared errors (MSEs) than the forecasts from the Poisson process (Supplementary Fig. 7). MSE is the average squared difference between the forecast value and the true value, which is equal to the sum of the variance and the bias squared, and provides a measure of the trade-off between accuracy and precision. For about half of the fault segments (45 of the 93), the MSEs of the forecasts from the Poisson process are more than twice of that from the model-averaged forecasts (Supplementary Fig. 8). The 14 fault segments for which the model-averaged forecast 95% credible interval didn’t cover the true mean were characterised by few events being left in the record: 6 had only 4 events left in the record and thus too few to fit models with more than two parameters, while a further 7 fault segments had fewer than 7 events. Even though in some situations with a small number of events in the record, the less informative, more uncertain Poisson-based forecasts seem to cover the true value, the majority of fault segments with small numbers of events are still better represented by the model-averaging approach (e.g., 23 out of the 30 fault segments which had 4 events in the retrospective forecasts have model-averaged forecasts with much smaller MSE than the Poisson forecasts). It is anticipated that for most hazard modelling purposes the smaller errors associated with the model-averaged forecasts favour their use. Having said this, the Poisson process may still be a valuable model when limited data are available, which is the case for many fault segments that are not included in this study because they have less than five events in the record.

The MSEs of the retrospective forecasts from the single-best model approach are very close to those from the model-averaging approach for the majority of the fault segments (Supplementary Fig. 9), all within two times relative difference. Model-averaging with WAIC weights is not usually designed to achieve a better MSE than the single best-model approach. However, when there is some uncertainty as to the best model, model-averaging outperforms a single best-model primarily in terms of better representing all the uncertainties.

The Bayesian model-averaging approach presented here explicitly considers model uncertainty based on the data and associated measurement errors, rather than relying on selection of a best model. Retrospective testing shows that model-averaging provides more informative and accurate forecasts compared with a single-best model approach or assuming a Poisson process (i.e., random earthquake recurrence). The earthquake probabilities presented in this study also provide a testable hypothesis of future earthquake occurrence that can be evaluated at a global scale.

Methods
Models
For each fault segment, we obtain 100 sequences of Monte Carlo (MC) samples for the ages of the sequence of large earthquakes in the paleoseismic record (see Data and Resources section). Each sequence of MC samples is then considered a realisation of the occurrence times of large earthquakes along that fault9,36. We denote them by t k0 < t k1 < ... < t kN i ≤ T, where k=1, 2, …, 100 denotes the kth MC sample, N i denotes the number of earthquakes in the record for the ith fault segment with i=1, 2, …, 93, and T denotes the censoring time which we take as the year 2022. The inter-event times for the kth MC sample of the ith fault are then x k1 = t k1 - t k0 , ..., x kN i = t kN i - t k(N i-1) .

For each earthquake record, we fit the following five models. The first is a Poisson process with occurrence rate Z k \lambda for the kth sequence of MC samples, where Z k is a random variable (k=1, ⋯, 100) capturing the similarities between the 100 MC samples for each fault. The second model is a Gamma renewal process with the inter-event times for the kth sequence of MC samples following a Gamma distribution with probability density function
f(x; \alpha, \lambda, Z k, Y k) = \frac{\Gamma(Z k \alpha) (Y k \lambda) Z k \alpha x Z k \alpha-1 exp{-Y k \lambda x} }{g} (1)
where Z k \alpha and Y k \lambda are shape and rate parameters, with Z k and Y k being random variables (k=1, ⋯, 100) capturing the similarities between the 100 MC samples for each fault. The third model is a Weibull renewal process with inter-event times for the kth sequence of MC samples following a Weibull distribution with probability density function
f(x; \alpha, \lambda, Z k, Y k) = Z k \alpha (Y k \lambda) Z k \alpha x Z k \alpha-1 exp{- (Y k \lambda x) Z k \alpha} (2)
where Z k \alpha and Y k \lambda are shape and rate parameters. The fourth model is a Brownian Passage-Time (BPT, also called inverse Gaussian) renewal process with inter-event times for the kth sequence of MC samples following a BPT distribution with probability density function
f(x; \mu, \beta, Z k, Y k) = \sqrt{\frac{Y k \beta}{2 \pi (Z k \mu) 2 x 3}} exp{- \frac{(Z k x - \mu) 2 Y k \beta}{2 \mu 2 x}} (3)
where Z k \mu and Y k \beta are the mean and coefficient of variation of the distribution. The fifth model is a lognormal renewal process with inter-event times for the kth sequence of MC samples following a lognormal distribution with mean Z k+ \mu and standard deviation Y k \sigma, both on the log scale.

Estimation
Given the kth MC sample of the earthquake occurrence times t k0, t k1, ..., t kN i along the ith fault, with final censoring time T, the likelihood of the kth MC sample for each model is
L(\theta k; t k1, ..., t kN i, T) = (1 - F(T - t kN i; \theta k)) \prod j=1 N i f(t kj - t k(j-1); \theta k), (4)
where k=1, ⋯, 100, \theta k=(\lambda, Z k) for the Poisson process, \theta k=(\alpha, \lambda, Z k, Y k) for the Gamma and Weibull renewal processes, \theta k=(\mu, \beta, Z k, Y k) for the BPT renewal process, and \theta k=(\mu, \sigma, Z k, Y k) for the lognormal renewal process.

The different MC samples from the same fault should have similar recurrence patterns. To reflect this, we assume that both Y k and Z k follow a distribution with mean 1, i.e.,
Y k \sim Gamma(1=\sigma 2 Y, 1=\sigma 2 Y), Z k \sim Gamma(1=\sigma 2 Z, 1=\sigma 2 Z), (5)
where \sigma Y and \sigma Z are the standard deviations of Y k and Z k, respectively.

A Markov Chain Monte Carlo (MCMC) algorithm generates samples from the joint posterior distribution of \theta k, \sigma Y and \sigma Z given the occurrence times t k0, t k1, ..., t kN i from the kth MC sample and the censored time T, using software JAGS and the R2jags package in R37. We use three chains, half-normal priors for \alpha and \beta, and weakly informative half-t prior distributions for the variance parameters \sigma Y and \sigma Z 38, i.e.
\alpha \sim N(0, 100 2)T(0,), \lambda \sim N(0, 100 2)T(0,),
\mu \sim N(0, 100 2)T(0,), \beta \sim dt(0, 0:04, 3)T(0,),
\sigma \sim dt(0, 0:04, 3)T(0,), \sigma Y \sim dt(0, 0:04, 3)T(0,),
\sigma Z \sim dt(0, 0:04, 3)T(0,):

For the MCMC algorithm, we use three chains with 5,010,000 iterations, discarding the first 10,000 iterations as burn-in, and use a thinning rate of 1000. The scale reduction factors of the Gelman-Rubin convergence diagnostic are all less than 1.02, indicating convergence39,40.

Model-averaged forecasts
To calculate Bayesian model-averaged forecasts, we combine the posterior distributions of the forecast under each model using model weights. We use prediction-based Bayesian model-averaging (PBMA)15 with the model weights calculated using the Watanabe-Akaike Information Criterion (WAIC20). This is far less sensitive to the priors for the parameters than classical Bayesian model-averaging (CBMA), which uses posterior model probabilities15. PBMA is sometimes referred to as Bayesian model combination41. Unlike CBMA, it does not involve the assumption that one of the models is true. The WAIC for model k is calculated as
WAIC k = -2 \sum i=1 n log(p(y i | y, k)) + 2p k (6)
where p(y i | y, k) is the pointwise posterior predictive density from model k, which can be estimated using the mean of the posterior MCMC sample of p(y i | \theta k, y, k) for model k; and p k is a correction for overfitting. A common choice for p k is
p k = \sum i=1 n var log p(y i | \theta k, y, k), (7)
where each term in the summation can be estimated by taking the variance of the posterior MCMC sample of log p(y i | \theta k, y, k) for model k. The WAIC weight for model k is given by
p(k | y) / exp(-(WAIC k - m i in WAIC i) = 2) (8)
WAIC is a prediction-based criterion, analogous to AIC in the non-Bayesian setting, and use of Eq. (8) to define the model weights is motivated by the form of AIC weights42. Alternative approaches to prediction-based model averaging in seismology have been proposed43,44. WAIC is known, in large samples, of being equivalent to Bayesian leave-one-out cross-validation (B-LOO)45. When B-LOO is used in model averaging it is known as Bayesian stacking, and has the useful property that, for large samples, it leads to the best linear combination of the posterior distributions of the forecasts from each model, whereas CBMA will lead to use of the posterior distribution of the forecast from the single best model41–43 (which is why some authors refer to CBMA as a tool for model selection41). We would therefore expect WAIC weights to provide a close-to-optimal linear combination of the posterior distributions of the forecasts from each model, whilst being much less computationally-intensive than Bayesian stacking.

For each model, we can obtain a posterior MCMC sample of the forecast quantity: either the forecast occurrence time of the next large earthquake or the forecast probability of at least one large earthquake occurring in the next 50 years along the specified fault segment. These forecasts are conditioned on the fact that there was no large earthquake between the last large earthquake occurrence time and the year 2022. The five posterior MCMC samples of the forecast quantity were combined into one model-averaged posterior sample by randomly taking the value from one of the five posterior samples at each iteration. The probability weights for the random sampling are the WAIC weights for the corresponding five models.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have proposed a Bayesian model-averaging approach for forecasfing large earthquakes along acfive faults. They emphasise the importance of model uncertainty in the convenfional single-model approaches and ensemble all the candidate models with their weights derived from the model selecfion criterion WAIC to mifigate the effect of model uncertainty on forecasts.\nSelecfion of recurrent-fime distribufion has long been controversial for long-term earthquake forecasts. The authors argue that no single model consistently provides the best fit across the global catalogs, or even within the individual fault systems. Their approach presents a plausible solufion in this situafion. Addifionally, the authors discussed the periodicity of earthquake recurrence and revealed through regression analysis that the degree of periodicity depends not only on the occurrence rate (as discussed in the previous studies) but also on the faulfing type and the tectonic region.\nThe authors provide the probability forecasts in the next 50 years for the global 93 fault segments. Their retrospecfive analysis shows good accuracy of the proposed approach for most of the fault segments.\nThe paper is well organised and insighfful. However, I'd like to provide feedback on several specific points:\nL86: The authors describe that the AIC can’t provide guidance as to how much befter a parficular model is compared to other models. However, since the AIC is an esfimator of the predicfion error of a fifted model in terms of the log-likelihood of a new observafion, it can indeed assess the relafive predicfive competence of different candidate models.\nL148: Though the authors have revealed relafionships between the number of events and the single-best model in their catalogs, they do not provide discussion on characterisfic behavior of earthquake recurrence in each single model. For example, calculafing the normalised recurrence intervals, which are divided by the mean recurrence intervals of the respecfive fault segments, grouping them by the selected single-best models, and comparing the histograms between the single-best models might help us to observe some characterisfic recurrence behavior in each model. In this case, I suggest referring the following paper:\nNishenko, S. P., & Buland, R. (1987). A generic recurrence interval distribufion for earthquake forecasfing. Bullefin of the Seismological Society of America, 77(4), 1382-1399.\nL207-222: The ranges of the CIs in this paragraph appear too short compared to the sample size of 93 and the dispersion of the data points in Figure 4. I suspect that this might be because all the 100 MC samples and/or 5000 MCMC samples for each segment are used as the data points in evaluafing the CIs. If this is the case, it is not appropriate to discuss the stafisfical significance using the CIs from these leveraged samples.\nL244: From Figure S4, the authors have shown that the forecast values from the model-averaging approach have much smaller mean square errors (MSEs) than those from the Poisson process. However, it is difficult to see the relafive improvement in accuracy from this figure because the fault segments have a wide range of scales in their recurrence fimes. Therefore, I suggest adding a figure and comments on the relafive differences in MSEs between the model-averaging approach and the Poisson process.\nL252: favours -> favour ?\nL265: The accuracy of forecasts for the single-best model approach does not seem to be discussed in the retrospecfive tests in Secfion 5.\nL292: The sentence ends with ‘..’.", "reviewer_2": "The paper aims at introducing a new methodology to esfimate large earthquake occurrence on faults. Although the topic is of great interest and the paper contains a very interesfing approach, I have some concerns that should be carefully addressed; some of them might be crifical.\nBelow I list my main concerns\n- This is likely the most crifical point. In the \"Introducfion\" secfion, the authors did not quote recent bibliography on the earthquake recurrence models on faults, such as the recent Neely et al (SRL, 2023) and the work made in UCERF3 (Field et al.", "reviewer_2015": "Long-Term Time-Dependent\nProbabilities for the Third Uniform California Earthquake Rupture Forecst(UCERF3), BSSA).\nInterestingly, Field et al (2015) discuss the impact on the earthquake occurrence models (in\nthis case, the Brownian Passage Time also used here) when some assumptions are relaxed, such\nas relaxing a strict faults segmentation. In this paper, the authors are assuming a strict faults\nsegmentation which is, at best, debatable (so many recent examples, including Kaikoura eqk).\nIf this assumption is relaxed, the models used by the authors may be incoherent and/or\ninconsistent. What is the impact of relaxing faults segmentation, and what is the reliability of\nsuch an assumption? These are important issues that should be carefully and satisfactorily\naddressed in the paper.\nReply: We thank the reviewer for the suggestion of the two references and for giving us the\nopportunity to clarify the aim of our research study. We have added these two references in the\nintroduction and added discussions around the reviewer’s comments. Please see Lines 57-64,\n88-99, 171-182, 216-218 in the revised manuscript. We also provide some further discussions\nbelow.\nWe address each of the two studies mentioned by the reviewer individually.\n1. The Long-Term Fault Memory Model proposed by Neely et al. (2023) assumes that\npartial strain release during earthquakes can lead to accumulated strain on a fault, and\nthat therefore the timing of future events is dependent not only on the time elapsed since\nthe most recent event (as in renewal models) but also previous inter-event times, which\nmay or may not lead to an accumulation of strain that brings forward the timing of the\nnext rupture. This model presents an interesting hypothesis, and the authors applied the\nmodel to the Mojave section of the San Andreas fault in California. However, Griffin\net al. (2020) used the ‘memory coefficient’ of Goh and Barabasi (2008) to characterise\nthe autocorrelation of successive inter-event times for 80 long-term earthquake records.\nFor the vast majority of earthquake records there was no significant correlation or anti-\ncorrelation between the length of successive inter-event times. In addition, Kempf and\nMoernaut (2021) generated synthetic earthquake records using Monte-Carlo sampling\nfrom typically parameterised renewal models and showed that despite these being\nindependent events it is possible to generate records that show positive or negative\nmemory coefficients. These authors argued therefore that caution should be practised\nwith respect to interpreting memory from short paleoseismic records. We therefore did\nnot consider the Long-Term Fault Memory model as a generic model for global large\nearthquake recurrence, although we do not discount the applicability of the model in\nsome particular situations.\n2. The study by Field et al. (2015), of which some aspects relevant to this study are\ndetailed in Field (2015), discusses the use of renewal models in the context of\ndevelopment of unsegmented fault and rupture models in the UCERF-3 model, which\nis developed for California earthquake rupture forecast.\nUCERF-3 consists of four components, namely the fault model, the deformation model,\nthe earthquake rate model, and underlying these three, the earthquake probability\nmodel. In their earthquake probability model, they relaxed the strict fault segmentation\nassumption in UCERF-2 by subdividing the main fault sections into many fine\nsubsections and considering interactions between different subsections. For each\nindividual fault subsection, although interactions between fault subsections are\nincluded, the fundamental elastic-rebound theory part of the model is a BPT renewal\nmodel.\nWe aim to use a global record to find a better approach for the elastic-rebound theory\npart of the earthquake probability model component, that is, a better way to model the\ninter-arrival times which have been frequently modelled using either the BPT renewal\nprocess, because of its physical explanation of the earthquake process (Matthews et al.,\n2002; Nomura et al., 2011; Field et al., 2015), or the best model among a few candidate\nrenewal processes, based on an information criterion (Ogata, 1999; Rhoades et al.,\n2011). Once this layer of modelling is built, one can then add interactions between small\nsubsections to forecast large earthquakes along specific individual faults. Improving the\nforecast from this layer will result in a holistic hazard model that includes fault models,\ndeformation models, earthquake rate models, and earthquake probability models. We\nuse the global data set to demonstrate that there is no universal model and the BPT\nrenewal process is not necessarily the best model. Choosing a different renewal model\nor choosing to average over several renewal processes may produce a better forecast.\nThe fundamental question is around fault segmentation, and in particular for the San\nAndreas which a) has many paleoseismic sites along different fault segments; and b)\nhas been the subject of several other studies to estimate time-dependent earthquake\nprobabilities.\nIt is also noted that our study attempts to present a globally applicable method, and that\nin some more well-studied regions such as the San Andreas Fault System it may be\nappropriate to develop fault-specific approaches such as that presented by Field et al\n(2015) that build in additional information about fault segmentation and geometry. It is\nnot our intention to repeat such studies. That said, the question of segmentation and co-\nrupture of multiple fault segments in less frequent, larger magnitude events vs rupture\nof individual segments in more frequent, smaller magnitude events is significant from\na seismic hazard point of view. The challenge arises because, apart from a small number\nof historical events where we have information regarding whether co-rupture of\nneighbouring fault segments occurred, uncertainties in the paleoseismic record cannot\nresolve whether two fault segments ruptured in the one event, or in two separate events,\neven when dating uncertainties overlap. That is, two separate events a day apart on\nneighbouring fault segments will look very different to a single event rupturing both\nsegments from a hazard point of view, yet these two possibilities will be\nindistinguishable by comparison of paleoseismic records from the two segments. While\nmany studies assume overlapping dating uncertainties to indicate co-rupture (e.g. the\n‘stringing pearls’ approach; Biasi and Weldon 2009), these are subject to a number of\nsubjective decisions. For example, in Biasi and Weldon 2009 (p478) overlapping\nruptures are arbitrarily discarded based on a series of rules, such as if the median is\nmore than 40 years from the mean, or the width of the combined age distribution is\ngreater than 90 years, and several other criteria. While these all may or may not be\nconsidered reasonable criteria to impose, we cannot independently determine from the\ndata alone whether events at adjacent paleoseismic sites that have overlapping age\ndistributions are the same event or different events. At best, paleoseismic records can\nexclude co-rupture in the case that there is negligible overlap in age uncertainties\nbetween sites. Some inference on the likelihood of co-rupture can be gained from\nanalysis of historical surface rupturing events (e.g. Biasi and Wesnousky 2021; Walsh\net al 2023), however the problem of combining data from multiple paleosesimic sites\nis still largely intractable without recourse to some kind of expert judgement.\nIn the Field et al (2015) approach, time-dependent probabilities are assigned to different\nruptures (i.e. different combinations of fault segments) by taking averages of the\nrelevant parameters (mean inter-event time and time of most recent event) over the\nparticipating fault segments. To the best of our knowledge, the relative likelihood of\ntwo adjacent segments rupturing together in one event vs independently in two events\nis not constrained by the paleoseismic data. Note that the value for the aperiodicty is\nassumed (P552 of Field, 2015), while the paleoseismic data is used only to constrain\nthe mean recurrence intervals. The model also ‘assumes the availability of a long-term\nearthquake-rate model, which gives the frequency of all possible superseismogenic\nruptures on a fault or fault system’ (Field et al., 2015 p515). This seems to be quite\ndifferent to what we are doing here, where we are fitting renewal models to the paleo\ndata directly at each site. The method developed in our study could be used as an input\nto a UCERF3 style model, in terms of providing more accurate constraints on the mean\nrecurrence intervals and aperiodicty, but it is beyond the scope of this study to\nimplement a UCERF3 approach.\n- I am not a paleoseismologist, but my feeling is that the data that can be collected from a fault\nmay be filtered. In essence, if two large earthquakes occur nearby in time, they cannot be\ndistinguished (in many cases). This means that the observations have a low-frequency filter,\nand removing short inter-event times the data appear more periodic than what they are. Maybe\na short discussion on the time resolution of the catalogs may be worthwhile.\nReply: Thank you for raising this issue. Please see our detailed response below. We have\nincluded a short discussion in the revised manuscript on the time resolution of the catalogs as\nsuggested by the reviewer. See Lines 285-295.\n1. Perhaps the most commonly used event marker in paleoseismic studies is colluvium\nformed due to collapse of a fault scarp. Where there is a component of dip-slip, fault\nscarps are often unstable, and if collapse does not occur co-seismically (as is often the\ncase), it is typically seen within days to months, and occasionally a few years after an\nearthquake (e.g. many of the 2016 Kaikoura earthquake scarps degraded rapidly, e.g.\nKearse et al., 2018). Strike-slip earthquakes often lead to vertical cracking at the surface\nof the fault plane, and again these are unstable structures that typically collapse and/or\ninfill quickly. Future events that re-rupture the same fault segment will typically fault\nthe deposited colluvium, before deposition of a new (unfaulted) colluvial unit above\nthe previous unit. In this situation, exposure of the faulted and unfaulted colluvial units\nis an indicator of two separate events, and identification of multiple events is dependent\nentirely on the stratigraphic and structural relationships of the units, and not the length\nof time between their formation. i.e., the stratigraphy will look approximately the same\nregardless of the inter-event time, barring perhaps some inter-event soil formation or\nother deposition, depending on the environment, in the case of a long inter-event time.\nNote that the identification of the two events based on stratigraphic relationships is a\nseparate question as to whether current dating techniques can resolve their absolute\ntiming. While preservation of paleoseismic evidence is an important consideration in\nthe interpretation of any paleoseismic record, the characteristics of the geomorphic and\ngeological setting, and the relative rate of tectonic processes to non-tectonic\ngeomorphic processes, will control whether evidence of past earthquakes is preserved\n(McCalpin and Nelson 2009). Ideally, paleoseismic studies consider these factors in\nsite selection and interpretation; while it is acknowledged that missed events are a\npossibility, it is not clear that this should lead to a systematic bias in our statistics\ntowards or away from periodicity.\n2. As discussed by Williams (2022), if fault re-rupture commonly occurred shortly after\nprevious earthquakes, then we would expect to see this frequently in the historical\nrecord of surface rupturing earthquakes, which we largely do not. To quote from that\npaper:\n“It seems important to recognize, however, that such short‐term re‐rupturing of\nindividual faults or fault sections is almost entirely absent from the historical and\nmodern records of large earthquakes. In California, for example, even the most famous\nexamples of short‐term re‐rupture of the same fault section had IETs of 40 yr\n(Imperial fault, 1940–1979; southern San Andreas Mojave section, 1812–1857;\nnorthern San Andreas Santa Cruz Mountain section, 1838–1890–1906). It is im∼portant\nto note that in the case of the San Andreas Mojave and Santa Cruz Mountain sections,\nthese re‐rupture events were in fact captured in the paleoseismic records of the affected\narea (Weldon et al., 2004; Streig et al., 2014). Were short term re‐ruptures as common\nglobally as predicted by truly random recurrence, it seems reasonable that at least\nseveral examples would have been observed, particularly given the preponderance of\nactive seismogenic faults in tectonically active areas around the world and the length\nof the historic record in some regions (i.e., thousands of years in Turkey, Jordan, China,\nand Japan).”\nThis suggests that it is unlikely that we are systematically missing very short inter-event\ntimes and hence biasing our analysis. Note that this is a different issue to the question\nof correlating events between paleoseismic study sites and resolving whether events\nwith overlapping age uncertainties at multiple sites are the same event or different\nevents that occurred closely together in time.\nBerryman, K.R., Cochran, U.A., Clark, K.J., Biasi, G.P., Langridge, R.M. and Villamor, P.,\n2012. Major earthquakes occur regularly on an isolated plate boundary fault. Science,\n336(6089), pp.1690-1693.\nBiasi, G.P. and Weldon, R.J., 2009. San Andreas fault rupture scenarios from multiple\npaleoseismic records: Stringing pearls. Bulletin of the Seismological Society of America,\n99(2A), pp.471-498.\nBiasi G. P., and Wesnousky S. G. 2021. Rupture passing probabilities at fault bends and steps,\nwith application to rupture length probabilities for earthquake early warning, Bull. Seismol.\nSoc. Am. 111, 2235–2247.\nField, E.H., 2015. Computing elastic‐rebound‐motivated earthquake probabilities in\nunsegmented fault models: A new methodology supported by physics‐based simulators.\nBulletin of the Seismological Society of America, 105(2A), pp.544-559.\nField, E.H., Biasi, G.P., Bird, P., Dawson, T.E., Felzer, K.R., Jackson, D.D., Johnson, K.M.,\nJordan, T.H., Madden, C., Michael, A.J. and Milner, K.R., 2015. Long‐term time‐dependent\nprobabilities for the third Uniform California Earthquake Rupture Forecast (UCERF3).\nBulletin of the Seismological Society of America, 105(2A), pp.511-543.\nGoh, K.-I., & Barabási, A.-L. (2008). Burstiness and memory in complex systems. Europhysics\nLetters, 81(4), 48002. https://doi.org/10.1209/0295-5075/81/48002\nGriffin, J.D., Stirling, M.W. and Wang, T., 2020. Periodicity and clustering in the long‐term\nearthquake record. Geophysical Research Letters, 47(22), p.e2020GL089272.\nKearse, J., Little, T.A., Van Dissen, R.J., Barnes, P.M., Langridge, R., Mountjoy, J., Ries, W.,\nVillamor, P., Clark, K.J., Benson, A. and Lamarche, G., 2018. Onshore to offshore ground‐\nsurface and seabed rupture of the Jordan–Kekerengu–Needles fault network during the 2016\nM w 7.8 Kaikōura earthquake, New Zealand. Bulletin of the Seismological Society of America,\n108(3B), pp.1573-1595.\nKempf, P. and Moernaut, J., 2021. Age uncertainty in recurrence analysis of paleoseismic\nrecords. Journal of Geophysical Research: Solid Earth, 126(8), p.e2021JB021996.\nMcCalpin, J.P. and Nelson, A.R., 2009. Introduction to paleoseismology. International\nGeophysics, 95, pp.1-27.\nWalsh, E., Stahl, T., Howell, A. and Robinson, T., 2023. Two‐Dimensional Empirical Rupture\nSimulation: Examples and Applications to Seismic Hazard for the Kaikōura Region, New\nZealand. Seismological Society of America, 94(2A), pp.852-870.\nWilliams, R.T., 2022. Poisson behavior leads to bias when testing for periodicity in the\npaleoseismic record of large earthquakes. Seismological Research Letters, 93(1), pp.118-125.\n- Although I am very positive on the \"averaging\" approach proposed by the authors I think that\nthe authors should modify their terminology and their procedure to assign weights. The\nBayesian Model Averaging (BMA) is something different from what is made by the authors,\nbut, more important, BMA is not an averaging tool, but a model selector. This problem has\nbeen deeply discussed by Monteith et al. (2011; Monteith, K., J. L. Carroll, K. Seppi, and T.\nMartinez (2011). Turning Bayesian model averaging into Bayesian model combination, in Proc.\nof International Joint Conference on Neural Networks, San Jose, California, 31 July–5 August\n2011, 2657–2663.).\nReply: We agree that we should have been more careful about our use of the term Bayesian\nModel Averaging. We have amended the text to make this clearer. In particular\n1. We make it clearer that our approach is not to be confused with what we call classical\nBayesian Model Averaging (CBMA), the limitations of which have been discussed at\nlength by one of the authors of this paper (Fletcher 2018).\n2. We now discuss the connection between WAIC and Bayesian cross validation.\nSee Lines 429-439 which is also included below.\n“To calculate Bayesian model-averaged forecasts, we combine the posterior distributions of\nthe forecast under each model using model weights. We use prediction-based Bayesian\nmodel-averaging (PBMA) (Fletcher, 2018) with the model weights calculated using the\nWatanabe-Akaike Information Criterion (WAIC, Watanabe & Opper, 2010). This is far less\nsensitive to the priors for the parameters than classical Bayesian model-averaging (CBMA),\nwhich uses posterior model probabilities (Fletcher, 2018). PBMA is sometimes referred to as\nBayesian model combination (Monteith et al., 2011). Unlike CBMA, it does not involve the\nassumption that one of the models is true. For large samples PBMA will lead to use of the\nbest linear combination of the posterior distribution of the forecast from each model, whereas\nCBMA will lead to use of the posterior distribution of the forecast from the single best model\n(Monteith et al., 2011; Marzocchi et al., 2012; Yao et al., 2018), which is why some authors\nrefer to CBMA as a tool for model selection (Monteith et al., 2011).”\nFletcher, D. (2018). Model averaging. Springer.\nThis problem has been found and discussed also in seismology (Marzocchi Zechar and Jordan,\n2012, BSSA). The way in which the authors calibrate the weights, using a quantity that is\nproportional to probability (see their equation 8), is similar to BMA and it might have the same\nproblems (see the literature mentioned above).\nReply:We have revised the manuscript to make it clear that our approach will not suffer from\nthe same problems as CBMA. See Lines 429-439 which is also included in the responses to the\nprevious comment.\nIt might be interesting for the authors to know that more recently, Herrmann and Marzocchi\n(GJI, 2023) discussed the optimal weighting scheme for an averaging procedure, which is based\non a fundamentally different logic. Here I suggest a more careful re-thinking about the\nweighting scheme and a solid justification of the method used. Maybe just using a quantity\nproportional to WAIC instead of using the exponential, it might improve the problem (This is\njust a suggestion).\nReply: We thank the reviewer for this suggestion. We agree that the approach in Herrmann\nand Marzocchi (GJI, 2023) is different. We added further discussions around this point in the\nrevised manuscript. See below which is also in the revised text, Lines 448-452.\n“WAIC is a prediction-based criterion, analogous to AIC in the non-Bayesian setting, and use\nof Equation (8) to define the model weights is motivated by the form of AIC weights (Yao et al.,\n2018). Alternative approaches to prediction-based model averaging in seismology have been\nproposed by (Marzocchi et al., 2012) and (Herrmann & Marzocchi, 2023). We prefer to make\nuse of WAIC weights because WAIC has the desirable property of being equivalent to Bayesian\ncross-validation when used for model selection (Vehtari et al., 2017).”\n- Again, despite being positive about the method, I am wondering what is the reason for what\nfaults behave differently (they follow different distributions); is there a different physics? Or\nwhat else? In my understanding, the Weibull distribution is the best performing on faults which\nhave a longer dataset. I am wondering what happens if we use the Weibull model for all faults.\nOf course, in some cases, the Weibull model will perform less well than other models,\nincluding the average model, but what is the information gain per event? One model can\nperform worse than another for a small difference. In this case, the differences do not have any\npractical importance. I suggest the use of some metric, like the information gain per event, to\nshow how much a model performs better than others and discussing the results.\nReply: We did a regression analysis on the single-best model against tectonic regions and\nfaulting styles, but didn’t find significant relationship there. It is difficult to draw conclusion\non the underlying physics with the current data.\nWAIC provides an estimate of the predictive performance of a Bayesian model (Vehtari 2017).\nIf the WAIC-value for model A is smaller than that for model B, we estimate that model A will\nprovide better predictions than model B. WAIC weights provide a numerical comparison of\nthe amount by which model A is better at prediction than B, as they show how much weight\nshould be given to the prediction from each of these models when calculating a model-averaged\nprediction. Based on the WAIC weights in Figure 1, we can see that the predictive performance\nof the Weibull model is not uniformly better than the others. For 33 fault segments, the WAIC\nweight for the Weibull model is close to 0, which suggests that for these fault segments the\nestimated predictive accuracy of the other models is much better than that of the Weibull model.\nWe have added this discussion in the revised manuscript, Lines 160-167.\nWe also calculated the information gain per event for the Weibull renewal process with each\nof the other four models as a reference for each fault, which range between -0.17 and 0.7.\nDue to the word limit, this is only for the reviewer’s interest and was not included in the revision.\nAlthough the information gain per event is positive more often than negative, and the absolute\nvalues of the negative gain are mostly smaller than that of the positive ones except for when\nthe Gamma renewal process is used as a reference, we get the same conclusion as when using\nthe WAIC weights: the predictive performance of the Weibull model is not uniformly better\nthan the others.\n[REDACTED]\n- on a minor side, LINE 30 in the Abstract. I never liked too much the word \"robust\", because\nit is unclear what it means. At least to me.\nReply:We have revised the abstract to shorten it and clarified the relevant sentence.\nReviewer #3 (Remarks to the Author):\nThe manuscript by Wang et al presents a Bayesian statistical analysis of a large number of\npaleoseismic data sets with the apparent goal of developing an improved forecasting method.\nThe central claims of this manuscript, as I see them, are as follows: 1) Bayesian analysis allows\nuncertainties in earthquake-recurrence model selection to be more fully accounted for in\nsubsequent forecasting; 2) No single model can account for all of the observed paleoseismic\nsequences, even those from faults with similar kinematics, strain rates, or tectonic settings; 3)\nRetrospective forecasting shows that a “model averaging” approach to estimating earthquake\noccurrence is more informative and accurate than any single other model, including the Poisson\n(i.e. random earthquakes) model; 4) Various segments of the San Andreas fault system have a\nhigh likelihood of rupture (i.e. up to 70% chance) over the course of the next 50 years; and 5)\nthe large dataset analyzed here allows a variety of generalizations about the nature of\nearthquake recurrence as a function of strain rate, kinematics, and/or tectonic setting.\nOf these claims, it is my opinion that (1) and (2) are generally well supported by the data and\nassociated arguments. Claim (3), however, is a much more complicated question, and to be\nhonest, I feel that the manuscript on some level over intellectualizes it. There are\npractical/societal matters that need to be more fully considered (discussed below). Claim (4) is\nalso well supported by the data and associated discussion, but I do not see where it departs\nappreciably from a variety of other active forecasting methods that have been applied to the\nSan Andreas system. This kind of begs a question (and I do mean this respectfully) of what\nexactly is gained by using this new approach when compared with other methods? Finally,\nclaim (5) is also generally well supported here, although at least some of the subclaims related\nto it have been proposed by other authors who need to be cited. Considering aspects of the\npresentation, the writing is clear and concise. The quality and utility of the illustrations is\nsimilarly top notch. My recommendation is that this manuscript be further considered for\npublication following moderate revisions.\nReply:We thank the reviewer for these comments and we provide detailed responses below.\nI could probably summarize many of my concerns with this manuscript using a single question\nthat I would ask the authors to consider: Why is this approach fundamentally and necessarily\nbetter than what has been done previously? In this way, we need to realize that earthquake\nforecasting is not just an interesting scientific exercise. Rather, it has explicit practical goals,\nwhich vary depending on who is asking the question. As a society we make important decisions\nabout how and where to distribute finite resources based on such forecasts. It may be most\n“correct” for a forecasting approach to more fully consider uncertainties in recurrence-model\nselection, but that only really matters if doing so leads to a demonstrably more accurate forecast.\nThis is the central issue with all earthquake forecasting work - it's going to be a very long time\nbefore we really KNOW if these approaches are of any real utility or which are better than\nothers. Retrospective forecasting is certainly one way to get a tiny glimpse at the answer. It's\noften unsatisfying though, and it is for me in this case as well, in large part because the\nconfidence intervals for the next event are so large (an unfortunate reality of uncertain datasets,\nas this manuscript often points out). If our confidence interval for the occurrence of the most\nrecent event in the data set spans hundreds or thousands of years, and we apply it to a record\nwith a mean recurrence interval of about 100 years, is it really surprising/revealing that we\nsuccessfully capture the event within the predicted interval? So respectfully, why is the\napproach presented here preferred? How would we pitch it to an insurer? A local government?\nA builders association? The factor that all of these entities are interested in is less a matter of\nthe binary occurrence of an earthquake within some window of time, but rather the per-year\nlikelihood of a catastrophic event (often compared with the anticipated useful lifetime of a\nstructure). By decreasing the confidence window of occurrence, the method explicitly raises\nthe per year likelihood of an event, and that has implications. I understand this is a manuscript\nthat is focused on science, so maybe these questions do not need to be explicitly answered, but\nI do think that the subject matter really does demand at least some engagement with the\npractical matters at hand.\nReply:We are grateful to the reviewer for giving us the opportunity to make clarifications. We\nhave revised the manuscript accordingly to further clarify the aim of our research and thus\nanswer the question “Why is this approach fundamentally and necessarily better than what has\nbeen done previously?” Please see Lines 57-64, 88-99, 171-182, 216-218 in the revised\nmanuscript. We also provide more detailed responses below.\nNational earthquake rupture forecasts for local faults (used for insurance purposes, or for\nproviding advice to local government or builders associations) are normally done by building\nin not only a statistical earthquake probability model, but also fault models, deformation\nmodels and earthquake rate models. Taking UCERF-3 for California earthquake rupture\nforecasts for example (Field et al., 2015), the fundamental elastic-rebound theory part of the\nearthquake probability model is a BPT renewal model.\nWe aim to use a global record to find a more robust approach for the elastic-rebound theory\npart of the earthquake probability model, that is, a more robust way to model the inter-arrival\ntimes which have been frequently modelled using either the BPT renewal process, because of\nits physical explanation of the earthquake process (Matthews et al., 2002; Nomura et al., 2011;\nField et al., 2015), or the best model among a few candidate renewal processes, based on an\ninformation criterion (Ogata, 1999; Rhoades et al., 2011). Improving the forecast from this\nlayer of modelling will result in a holistic hazard model that includes fault models, deformation\nmodels, earthquake rate models, and earthquake probability models. We use the global data set\nto demonstrate that there is no universal model and the BPT renewal process is not necessarily\nthe best model. Choosing a different renewal model or choosing to average over several\nrenewal processes may produce a better forecast.\nIt is also noted that our study attempts to present a globally applicable method, and that in some\nmore well-studied regions such as the San Andreas Fault System it may be appropriate to\ndevelop fault-specific approaches such as that presented by Field et al (2015) that build in\nadditional information about fault segmentation and geometry. It is not our intention to repeat\nsuch studies. The method developed in our study could be used as an input to a UCERF3 style\nmodel, in terms of providing more accurate constraints on the mean recurrence intervals and\naperiodicty, but it is beyond the scope of this study to implement a UCERF3 approach.\nToward the point above, it also seems that any manuscript presenting a new approach to\nearthquake forecasting should provide some detailed comparisons with previous forecasting\nwork. For example, the earthquake occurrence probabilities for various segments of the San\nAndreas presented here strike me as largely the same as those presented by the Uniform\nCalifornia Earthquake Rupture Forecast. UCERF also includes a time-dependent forecast in\naddition to a Poisson. I`ll grant that the approach here presents a less uncertain forecast than a\nPoisson alternative whilst also being broadly as accurate (per the retrospective forecast), but\nwhat about other time-dependent models such as UCERF? Again respectfully, why is the\napproach presented here better (this is not to express doubt that it may be, but rather that it is\nnot clear WHY it necessarily is)?\nReply: We hope our response to the previous question also answers this question here.\nUCERF-3 consists of four components, namely the fault model, the deformation model, the\nearthquake rate model, and underlying these three, the earthquake probability model. We aim\nto use a global record to find a better approach for the earthquake probability model. It is not\nour aim to build a full hazard model for all the 93 fault segments worldwide. We use the global\ndata set to demonstrate that there is no universal earthquake probability model and the BPT\nrenewal process is not necessarily the best model. Choosing a different renewal model or\nchoosing to average over several renewal processes may produce a better forecast. The method\ndeveloped in our study could be used as an input to a UCERF3 style model, in terms of\nproviding more accurate constraints on the mean recurrence intervals and aperiodicty.\nFinally, the manuscript consistently does a bit of a hit job on the use of the Poisson model in\nearthquake forecasting. I think that should be dialed back a bit. Yes, I too think the idea that\nearthquakes would occur randomly in time is a bit rubbish from a mechanics perspective, and\nit is not supported by the current state of the paleoseismic literature as the authors seem to agree.\nThat being said, there are practical reasons for the adoption of a Poisson model within\nearthquake forecasting as described in the line referenced comments (i.e. it is a decent choice\nwhen the date of the most recent event is not known).\nReply: We agree that Poisson process is still a valuable model. We have added some\ndiscussions around this in our revised manuscript. See Lines 333-335 which is also provided\nbelow.\n“Having said this, the Poisson process may still be a valuable model when limited data are\navailable, which is the case for many fault segments that are not included in this study because\nthey have less than five events in the record.”\nIn summary, this is interesting and high quality work. Having communicated with some of the\nauthors on the subject of earthquake statistics and forecasting previously, I know that their\ndesire to do good science is genuine. I think this work too will soon be worthy of publication,\nbut given its nature, some considerably more nuanced discussion of why this approach is better\nthan the alternatives is required.\nReply:Thank you. We really appreciate all your comments on our work.\nRandy Williams\nLine Referenced Comments\n18-20: Interesting question really - does better incorporation of uncertainty lead to a more\nrobust forecast? Depends on what a “robust” forecast is. Seems at least possible that we end up\nwith a forecast that is no more or less accurate than was the case prior to using this Bayesian\napproach.\nReply:We agree with the reviewer that the use of “robust” is vague here. We have revised the\nabstract to shorten it and clarified the relevant sentence.\n21-22: Again, this is cool, but minimizing the forecast uncertainty is only a useful thing if the\nforecast itself is accurate. The central dilemma of earthquake forecasting! We’ll need decades\n(if not centuries) to find out.\nReply: Here we are comparing the model-averaging approach against the single best-model\napproach. When there is some uncertainty as to the best model, model averaging (MA)\noutperforms a single best-model (BM) primarily in terms of better representing all the\nuncertainty. This typically means the coverage of intervals is much better with MA, and can be\npoor with BM.\n40-47: Could simplify this section for the novice with some more plaintext terminology. The\nPoisson process describes events that effectively occur randomly in time. And the antithesis of\nthat for periodic sequences of course.\nReply:We have made changes to the text. See Lines 49-54 which is also included below.\n“A renewal process is a statistical model that describes event occurrences in time, treating\neach new occurrence as a renewal in which the system is reset. It assumes that the time between\nevents (the inter-event time, such as the time between two consecutive earthquake occurrences)\nare independent and identically distributed. The Poisson point process is a special type of\nrenewal process that describes events effectively occurring randomly in time, and thus does\nnot have any memory of the history.”\n57-60: Maybe rephrase to better emphasize the severity of the problem. I’d wager that relatively\nfew of those 90 records actually have anywhere near 20 distinct events. In my experience\nrecords with more than 10 events are quite rare. Many have fewer than 5.\nReply: Thank you for this suggestion. We have added more description on the data here to\nemphasize the problem. See Lines 79-80 which is also included below.\n“(90 out of a total of 93 earthquake records in this study have evidence for fewer than 20\nearthquakes; 79 records have evidence for no more than 10 earthquakes; and 30 have evidence\nfor only 5 earthquakes)”\n60-61: It might, certainly, but does it? Again, we really don’t know at this point.\nReply: We agree, which is why we used “may produce erroneous or unrealistically precise\nforecasts”.\n77: First occurrence of (what I assume is) Brownian passage time. Please define.\nReply:Thank you. We have defined it when it first appeared in the revised manuscript.\n98-99: See comment above about record length.\nReply: We have now used “only 3 fault segments having more than 20 events, and 14 having\nmore than 10 events” to emphasize the severity of the problem.\n108: Perhaps this comment is too “in the weeds” but I feel like we have known for a while that\na lognormal distribution of interevent times is not an appropriate model for earthquakes, largely\nbecause its hazard function actually decreases through time from some maximum value.\nReply:We thank the reviewer for this comment. Given that lognormal has been considered by\nmany studies in literature for paleoearthquakes, we think it’s useful to include it for this global\nstudy. After reading this comment, we doubled checked our code and found an error when\ncalculating the WAIC weight for the lognormal model: it was a mistake when specifying the\nstandard deviation of the lognormal distribution for the calculation of the WAIC weight. We\nhave thus corrected the error and recalculated everything. Although some numbers reported in\nthe results changed, our conclusions do not change, and are better supported by the new\ncalculations.\n147-150: Quite interesting though.\nReply:We have now revised this part to add more discussion.See Lines 171-182 which is also\nincluded below.\n“It is unclear if this is simply an outcome of sampling or due to real differences in recurrence\nbehaviour between the different segments, including how neighbouring fault segments interact.\nVariability in observed earthquake recurrence behaviour at paleoseismic sites on the San\nAndreas Fault has been proposed to be (at least partially) due to overlap of ruptures occurring\non neighbouring segments (Weldon et al., 2004) and this has been supported by studies using\nearthquake simulators (Field, 2015). In contrast, relatively strongly quasi-periodic recurrence\non the Alpine Fault has been attributed to its geometric simplicity and relative isolation from\nother faults (Berryman et al., 2012). The persistence (or otherwise) of rupture barriers between\nsegments (Philibosian & Meltzner, 2020) may also be a significant factor controlling the\ndistribution of inter-event times observed on a fault segment. Therefore, because of the\nlimitations of the available data, even within a single well-studied fault system, we cannot use\na universal single best model, and it is not clear that one exists.”\n169: Could just call this “27” time as likely?\nReply:Thank you. We have changed this.\n207-222: Some of this discussion makes one speculate as to whether variations in strain rate or\nkinematics have an appreciable impact on record length and/or event preservation. Maybe\nworth a teaser sentence?\nReply: Thank you. We have added in the revised manuscript the following sentence, Lines\n282-284.\n“With more data becoming available in the future, one could investigate whether variations in\nstrain rate or kinematics have a significant impact on record length and event preservation.”\n211-215: Pretty sure that Berryman made a similar argument in the 2011 Alpine fault\npaleoseismology paper. Should probably give them a citation.\nReply:Thank you. We have now cited Berryman et al., 2012.\n251-252: I think this statement needs a bit of justification. Why is having smaller errors\nnecessarily preferable from a hazards forecasting perspective? As stated above, it’s only really\na good thing if the forecast is accurate.\nReply: Thank you. We have revised this part to reflect this suggestion. See Lines 326-335\nwhich is also included below.\n“Even though in some situations with a small number of events in the record, the less\ninformative, more uncertain Poisson-based forecasts seem to cover the true value, the majority\nof faults with small numbers of events are still better represented by the model-averaging\napproach (e.g., 23 out of the 30 fault segments which had 4 events in the retrospective forecasts\nhave model-averaged forecasts with much smaller MSE than the Poisson forecasts). It is\nanticipated that for most hazard modelling purposes the smaller errors associated with the\nmodel-averaged forecasts favours their use. Having said this, the Poisson process may still be\na valuable model when limited data are available, which is the case for many fault segments\nthat are not included in this study because they have less than five events in the record.”\n253: Although I am also in favor of time-dependent approaches, we can probably cut the\nPoisson model a little bit of slack here. It does have its uses from a hazards perspective,\nparticularly if we have some sense of the average interseismic period length but not a date on\nthe most recent slip event. We see in the retrospective forecast that in a select few cases smaller\nerrors caused the analysis to miss the most recent event.\nReply:Please see the reply above.\n235-240: Yes, the Poisson model does yield larger uncertainties and thus more frequent\n“overlap” with the actual events, but to be fair, the model average approach also yields fairly\nlarge uncertainties in an absolute sense. In the example of the San Andreas, if the forecast\nproduces a 95% confidence interval for the next event that is a period of several hundred years,\nwhilst the actual data suggest that a typical recurrence interval is on the order of about 100\nyears, it is a bit unsurprising that the forecast would capture the event within its uncertainty.\nRespectfully, I feel that I have to ask what exactly we learn from this? Certainly it could be\ninterpreted to show that your approach isn’t fundamentally off base, but why is this necessarily\nbetter than all of the other hazards forecasts that are available for the San Andreas (or other\nfault systems)?\nReply:As mentioned earlier in one of our replies, we aim to use a global record to find a more\nrobust approach for the elastic-rebound theory part of the earthquake probability model, which\ncan then be used to aid development of a better full hazard model like UCERF3 for local faults\nthat includes fault models, deformation models, earthquake rate models, and earthquake\nprobability models.\nAlso, the bit about the model averaging approach “outperforming” the Poisson is a bit\ncomplicated as written. You state that the model averaging approach successfully covered 79\nof the true occurrence times, but then state that the model averaging approach outperforms\nPoisson because 81 of the model averaging results had lower mean squared errors. Does that\nnot then by definition include 2 records for which your retrospective forecast failed (i.e. the CI\ndid not include the true occurrence time)? Strange to count those as successes of the approach.\nReply:We have now added an figure in the supplement and added the following discussion to\nclarify this. See Lines 337-339 which is also included below.\n“For about half of the fault segments (45 of the 93), the MSEs of the forecasts from the Poisson\nprocess are more than twice of that from the model-averaged forecasts (Supplementary Fig.\n7).”\nFinally, please consider comments above regarding the utility of the Poisson process in cases\nwhere the date of the most recent event is unknown. There is practical utility in the time-\nindependent approach in some cases, even if we believe that fault behavior is, in fact, time\ndependent.\nReply:Please see replies above.\n251-252: Again, this is a really tricky area of application. Please justify why this is favored,\nremembering WHY we want to forecast earthquakes accurately in the first place. We could of\ncourse be conservative/cautious and assume that any place that has a large fault that moved\nwithin the last 50 ka has a very high earthquake risk, but that would not do good things for\nlocal economies / insurance / building budgets. Although I do like this manuscript, I feel it has\na tendency to treat seismic forecasting as little more than an intellectual exercise by focusing\non what is most “correct” about fault behavior and uncertainty minimization. Neither of these\nfactors are readily translated into accuracy/reliability/usefulness of the forecast, although they\nmay be components of them. In a practical sense, its less about correctly predicting an\nearthquake within some interval of time. Rather, it's about calculating a per year likelihood of\na catastrophic event for a particular area. Your approach by definition maximizes that\nlikelihood by minimizing the predicted interval.\nReply: Please see replies above. To reiterate, we aim to use a global record to find a more\nrobust approach for the elastic-rebound theory part of the earthquake probability model, which\ncan then be used to aid development of a better full hazard model like UCERF3 for local faults\nthat includes fault models, deformation models, earthquake rate models, and earthquake\nprobability models. Currently, many studies, including UCERF3, use the BPT for the elastic-\nrebound theory part of the probability model. Our research found that the BPT renewal process\nis not necessarily the best model. Choosing a different renewal model or choosing to average\nover several renewal processes may produce a better forecast\n265-266: I did not see a comparison between the model averaging approach and the “single\nbest model” approach in the retrospective forecasting data?\nReply: We have added Supplementary Fig. 8 in the revised supplement to show the relative\ndifferences in MSEs between the single best model and the model-averaging approach, and\nadded this discussion in the revised manuscript. See Lines 336-341, which is also included\nbelow.\n“The MSEs of the retrospective forecasts from the single-best model approach are very close\nto those from the model-averaging approach for the majority of the fault segments\n(Supplementary Fig. 8), all within two times relative difference. Model-averaging with WAIC\nweights is not usually designed to achieve a better MSE than the single best-model approach.\nHowever, when there is some uncertainty as to the best model, model-averaging outperforms\na single best-model primarily in terms of better representing all the uncertainties.”\n267-269: This might be a stronger point to the work that is not being fully appreciated by the\nline in question. If we are ever to validate whether our hazards forecasts are meaningful (in the\nabsence of waiting a few hundred years to find out), it may be best achieved by looking at the\nglobal distribution of forecasts and how often they ring true.\nReply:Thank you for this suggestion. We agree with you that this may be the way to validate\nhazard forecasts. However, we decide to leave the sentence as it is as we think that statement\nshould be sufficient.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have well addressed the reviewers' concerns and comments. The changes and the\naddifional discussion have significantly improved the quality and clarity of the manuscript. However,\nthere sfill remains one point that should be further addressed to enhance its contribufion to the field.\nAs a response to my comment, the authors have added Supplementary Fig. 4, which compares the\nnormalized inter-event fimes between the single-best models and concluded that there does not seem\nto be any characterisfic recurrence behaviour for the fault segments that have the same single-best\nmodel.\nAlthough the conclusion is not wrong, the differences between the the single-best models should be\ndiscussed.\nFor example, though not menfioned in the revised manuscript, it is clear that the SDs of scaled inter-\nevent fimes in the two major models, Weibull and lognormal distribufions (Supplementary Fig. 4 (a,d)),\nare much different in their averages.\nWhile the SDs in the Weibull cases are smaller than 0.8 except for one case, more than half of the SDs in\nthe lognormal cases are larger than 0.8.\nThe authors should discuss this point.\nIn addifion, I guess the reason for this difference might be related to the long tail of the lognormal\ndistribufion.\nIf there is an excepfionally long inter-event fime in the records, it solely increases the coefficient of\nvariafion and also provides much higher likelihood to the lognormal model than the other models\nbecause of its long tail.\nI think this might be the case for the lognormal cases with large SDs, and that is why I suggested in the\nprevious round of review that a histogram of the normalised inter-event fimes be drawn for each single\nbest model.\nReviewer #2 (Remarks to the Author):\nI think that the authors address well my main concerns, maybe except one.\nI do appreciate the fact that the authors now make a disfincfion between CBMA and PBMA.\nHowever, I noficed one possible inconsistency; or, at least something that is not clear to me. The authors\nclaim that \"For large samples PBMA will lead to use of the best linear combinafion of the posterior\ndistribufion of the forecast from\neach model, ...\"; however, according to the authors \"WAIC weights provide a numerical comparison of\nthe amount by which a model is befter at predicfion than another, as they show how much weight\nshould be given to the predicfion from each of these models when calculafing a model-averaged\npredicfion.\"\nIt is not clear to me how using the weights that reflect the forecasfing performance of each single model\nwill guarantee, with a large sample, to get the best performance of the linear combinafion. Of course,\nthere must be a connecfion and the ranking of the models will be the same (the best performing model\nwill likely have also the largest coefficient in the linear combinafion), but I think that using the weights\nthat account for the individual performance of each single model does not guarantee that the linear\ncombinafion would be the best one. If I am wrong (and I could be!), it would be good to add more\nexplanafions and/or references that address this issue in depth.\nReviewer #3 (Remarks to the Author):\nThank you for taking the fime to address my previous comments. I think in general this manuscript has\nbeen improved relafive to the original version. I would suggest that a bit of revising be done in the\nintroducfion, largely because although useful addifional details and potenfial for clarity have been\nadded, my opinion is that the quality of the wrifing went down slightly. It reads a bit as though the\naddifional pieces of informafion requested by the reviewers were added without considerafion of overall\nflow. This is not a huge problem though, and I do not feel it should hold up publicafion.\nThe only suggesfion I made that was not implemented was to provide some comparisons between the\nprobabilisfic forecasfing provided by this model (beginning on line 204) and that provided by other\napproaches (UCERF, etc) where possible. This suggesfion was dismissed on the basis that other efforts\nlike UCERF are mulfi component forecast models, whereas this manuscript seeks to model the renewal /\nelasfic rebound process only. Respecffully, I don’t understand this argument. Other models like UCERF\nprovide, amongst other things, probabilisfic forecasts of the likelihood of large earthquakes over some\nforward period of fime. This manuscript also provides such forecasts. As such, I do not see why they\ncould not or should not be compared. It is really a prefty basic and fundamental quesfion. I do not\nnecessarily think this manuscript has to include probabilisfic forecasts as a result/product, but it does in\ncurrent form, and as such you have to expect that people will consider them on a case by case basis for\nan area of interest/concern. A comparison with other forecasts (where available) is warranted. I do not\nthink this would be complicated to do in a qualitafive way, and as such my recommendafion is that this\nmanuscript be further considered for publicafion following minor revisions. I think at a minimum, an\nexplicit statement about why that comparison is not included should be provided.\nOverall though, I`ll look forward to seeing this in print. Best of luck with the remainder of the process.\nRandy Williams\nRESPONSES TO REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have well addressed the reviewers' concerns and comments. The changes and the\nadditional discussion have significantly improved the quality and clarity of the manuscript.\nHowever, there still remains one point that should be further addressed to enhance its\ncontribution to the field.\nAs a response to my comment, the authors have added Supplementary Fig. 4, which compares\nthe normalized inter-event times between the single-best models and concluded that there does\nnot seem to be any characteristic recurrence behaviour for the fault segments that have the same\nsingle-best model. Although the conclusion is not wrong, the differences between the single-\nbest models should be discussed. For example, though not mentioned in the revised manuscript,\nit is clear that the SDs of scaled inter-event times in the two major models, Weibull and\nlognormal distributions (Supplementary Fig. 4 (a,d)), are much different in their averages.\nWhile the SDs in the Weibull cases are smaller than 0.8 except for one case, more than half of\nthe SDs in the lognormal cases are larger than 0.8. The authors should discuss this point.\nIn addition, I guess the reason for this difference might be related to the long tail of the\nlognormal distribution. If there is an exceptionally long inter-event time in the records, it solely\nincreases the coefficient of variation and also provides much higher likelihood to the lognormal\nmodel than the other models because of its long tail. I think this might be the case for the\nlognormal cases with large SDs, and that is why I suggested in the previous round of review\nthat a histogram of the normalised inter-event times be drawn for each single best model.\nReply:We thank the reviewer for this observation and for the thoughtful suggestions. We have\nadded the suggested discussion about the differences in the median standard deviations of the\nscaled inter-event times between the fault segments that were best fit by a Weibull model and\nthose that were best fit by a lognormal model in the revised manuscript. We also added some\ndiscussion around the reason for this difference. The reviewer is right that this is related to fault\nsegments that have inter-event times showing a long-tailed distribution which is best captured\nby a lognormal renewal process. The long inter-event times at the tail end increase the spread\nof the data and thus we see larger standard deviations of the scaled inter-event times for the\nfault segments that were best fit by a lognormal model. We added this discussion in the revised\nmanuscript as well. See lines 182-194 which is also included below.\n“For each of the 93 fault segments, we calculated the scaled inter-event times for each MC\nsample, and then reported the median and the 2.5% and 97.5% quantiles of the standard\ndeviations of the scaled inter-event times calculated for the 100 MC samples. Fault segments\nwith higher rate of earthquake occurrences appear to have smaller standard deviations of the\nscaled inter-event times (Supplementary Fig. 1 (e)). The median standard deviations of the\nscaled inter-event times for the 41 fault segments that were best fit by a Weibull model are all\nsmaller than 0.8 except for one case. In contrast, more than half of the median standard\ndeviations of the scaled inter-event times for the 15 fault segments that were best fit by a\nlognormal model are larger than 0.8 (Supplementary Fig. 1 (a,d)). The large median standard\ndeviations are related to fault segments that have inter-event times with long-tailed\ndistributions which result in a large spread of data values. These long-tailed distributions are\nbest fit by a lognormal renewal process.”\nReviewer #2 (Remarks to the Author):\nI think that the authors address well my main concerns, maybe except one.\nI do appreciate the fact that the authors now make a distinction between CBMA and PBMA.\nHowever, I noticed one possible inconsistency; or, at least something that is not clear to me.\nThe authors claim that \"For large samples PBMA will lead to use of the best linear combination\nof the posterior distribution of the forecast from each model, ...\"; however, according to the\nauthors \"WAIC weights provide a numerical comparison of the amount by which a model is\nbetter at prediction than another, as they show how much weight should be given to the\nprediction from each of these models when calculating a model-averaged prediction.\"\nIt is not clear to me how using the weights that reflect the forecasting performance of each\nsingle model will guarantee, with a large sample, to get the best performance of the linear\ncombination. Of course, there must be a connection and the ranking of the models will be the\nsame (the best performing model will likely have also the largest coefficient in the linear\ncombination), but I think that using the weights that account for the individual performance of\neach single model does not guarantee that the linear combination would be the best one. If I\nam wrong (and I could be!), it would be good to add more explanations and/or references that\naddress this issue in depth.\nReply: We are grateful for this comment, as it made us realise that we could have made this\npoint more clearly. In particular, we now state the following on lines 432-442.\n\"We prefer to make use of WAIC weights for the following reasons. Model selection using\nWAIC has the desirable property, in large samples, of being equivalent to Bayesian leave-one-\nout cross-validation (B-LOO)43. When B-LOO is used in model averaging it is known as\nBayesian stacking, and has the useful property that, for large samples, it leads to the best linear\ncombination of the posterior distributions of the forecasts from each model, whereas CBMA\nwill lead to use of the posterior distribution of the forecast from the single best model39-41\n(which is why some authors refer to CBMA as a tool for model selection39). We would\ntherefore expect WAIC weights to provide a close-to-optimal linear combination of the\nposterior distributions of the forecasts from each model, whilst being much less\ncomputationally-intensive than Bayesian stacking.\"\nReviewer #3 (Remarks to the Author):\nThank you for taking the time to address my previous comments. I think in general this\nmanuscript has been improved relative to the original version. I would suggest that a bit of\nrevising be done in the introduction, largely because although useful additional details and\npotential for clarity have been added, my opinion is that the quality of the writing went down\nslightly. It reads a bit as though the additional pieces of information requested by the reviewers\nwere added without consideration of overall flow. This is not a huge problem though, and I do\nnot feel it should hold up publication.\nReply:Thank you for this suggestion. We have revised the introduction to make it flow better.\nPlease see lines 43-54 and 66-91.\nThe only suggestion I made that was not implemented was to provide some comparisons\nbetween the probabilistic forecasting provided by this model (beginning on line 204) and that\nprovided by other approaches (UCERF, etc) where possible. This suggestion was dismissed on\nthe basis that other efforts like UCERF are multi component forecast models, whereas this\nmanuscript seeks to model the renewal / elastic rebound process only. Respectfully, I don’t\nunderstand this argument. Other models like UCERF provide, amongst other things,\nprobabilistic forecasts of the likelihood of large earthquakes over some forward period of time.\nThis manuscript also provides such forecasts. As such, I do not see why they could not or\nshould not be compared. It is really a pretty basic and fundamental question. I do not necessarily\nthink this manuscript has to include probabilistic forecasts as a result/product, but it does in\ncurrent form, and as such you have to expect that people will consider them on a case by case\nbasis for an area of interest/concern. A comparison with other forecasts (where available) is\nwarranted. I do not think this would be complicated to do in a qualitative way, and as such my\nrecommendation is that this manuscript be further considered for publication following minor\nrevisions. I think at a minimum, an explicit statement about why that comparison is not\nincluded should be provided.\nReply: Thank you for this suggestion. We agree with the reviewer that a comparison with the\nforecasts from other studies is important. We did try to compare but found that many\npublications of the paleoearthquake records only provide an estimate of the mean recurrence\ninterval which is equivalent to a forecast of the next occurrence time using a Poisson process\n(this comparison is shown in Supplementary Fig. 2-4). For studies that also provide a forecast\nof the probability of the next earthquake occurrence within a time period, they use different\nforecast intervals, which make it difficult to provide a meaningful comparison. For example,\nUCERF3 provides 30 year probabilities from 2014 aggregated by parent fault section, whereas\nwe provide 50 year probabilities from 2022 on each fault segment. For the San Andreas fault\nsegments, UCERF3 gives 30 year M>=6.7 mean probabilities (min, max):\nBig Bend: 18% (0.2%, 69%)", "reviewer_3": "The manuscript by Wang et al presents a Bayesian stafisfical analysis of a large number of paleoseismic data sets with the apparent goal of developing an improved forecasfing method. The central claims of this manuscript, as I see them, are as follows: 1) Bayesian analysis allows uncertainfies in earthquake-recurrence model selecfion to be more fully accounted for in subsequent forecasfing; 2) No single model can account for all of the observed paleoseismic sequences, even those from faults with similar kinemafics, strain rates, or tectonic seftings; 3) Retrospecfive forecasfing shows that a “model averaging” approach to esfimafing earthquake occurrence is more informafive and accurate than any single other model, including the Poisson (i.e. random earthquakes) model; 4) Various segments of the San Andreas fault system have a high likelihood of rupture (i.e. up to 70% chance) over the course of the next 50 years; and 5) the large dataset analyzed here allows a variety of generalizafions about the nature of earthquake recurrence as a funcfion of strain rate, kinemafics, and/or tectonic sefting.\nOf these claims, it is my opinion that (1) and (2) are generally well supported by the data and associated arguments. Claim (3), however, is a much more complicated quesfion, and to be honest, I feel that the manuscript on some level over intellectualizes it. There are pracfical/societal mafters that need to be more fully considered (discussed below). Claim (4) is also well supported by the data and associated discussion, but I do not see where it departs appreciably from a variety of other acfive forecasfing methods that have been applied to the San Andreas system. This kind of begs a quesfion (and I do mean this respecffully) of what exactly is gained by using this new approach when compared with other methods? Finally, claim (5) is also generally well supported here, although at least some of the subclaims related to it have been proposed by other authors who need to be cited. Considering aspects of the presentafion, the wrifing is clear and concise. The quality and ufility of the illustrafions is similarly top notch. My recommendafion is that this manuscript be further considered for publicafion following moderate revisions.\nI could probably summarize many of my concerns with this manuscript using a single quesfion that I would ask the authors to consider: Why is this approach fundamentally and necessarily befter than what has been done previously? In this way, we need to realize that earthquake forecasfing is not just an interesfing scienfific exercise. Rather, it has explicit pracfical goals, which vary depending on who is asking the quesfion. As a society we make important decisions about how and where to distribute finite resources based on such forecasts. It may be most “correct” for a forecasfing approach to more fully consider uncertainfies in recurrence-model selecfion, but that only really mafters if doing so leads to a demonstrably more accurate forecast. This is the central issue with all earthquake forecasfing work - it's going to be a very long fime before we really KNOW if these approaches are of any real ufility or which are befter than others. Retrospecfive forecasfing is certainly one way to get a finy glimpse at the answer. It's often unsafisfying though, and it is for me in this case as well, in large part because the confidence intervals for the next event are so large (an unfortunate reality of uncertain datasets, as this manuscript often points out). If our confidence interval for the occurrence of the most recent event in the data set spans hundreds or thousands of years, and we apply it to a record with a mean recurrence interval of about 100 years, is it really surprising/revealing that we successfully capture the event within the predicted interval? So respecffully, why is the approach presented here preferred? How would we pitch it to an insurer? A local government? A builders associafion? The factor that all of these enfifies are interested in is less a mafter of the binary occurrence of an earthquake within some window of fime, but rather the per-year likelihood of a catastrophic event (often compared with the anficipated useful lifefime of a structure). By decreasing the confidence window of occurrence, the method explicitly raises the per year likelihood of an event, and that has implicafions. I understand this is a manuscript that is focused on science, so maybe these quesfions do not need to be explicitly answered, but I do think that the subject mafter really does demand at least some engagement with the pracfical mafters at hand.\nToward the point above, it also seems that any manuscript presenfing a new approach to earthquake forecasfing should provide some detailed comparisons with previous forecasfing work. For example, the earthquake occurrence probabilifies for various segments of the San Andreas presented here strike me as largely the same as those presented by the Uniform California Earthquake Rupture Forecast. UCERF also includes a fime-dependent forecast in addifion to a Poisson. I`ll grant that the approach here presents a less uncertain forecast than a Poisson alternafive whilst also being broadly as accurate (per the retrospecfive forecast), but what about other fime-dependent models such as UCERF? Again respecffully, why is the approach presented here befter (this is not to express doubt that it may be, but rather that it is not clear WHY it necessarily is)?\nFinally, the manuscript consistently does a bit of a hit job on the use of the Poisson model in earthquake forecasfing. I think that should be dialed back a bit. Yes, I too think the idea that earthquakes would occur randomly in time is a bit rubbish from a mechanics perspecfive, and it is not supported by the current state of the paleoseismic literature as the authors seem to agree. That being said, there are pracfical reasons for the adopfion of a Poisson model within earthquake forecasfing as described in the line referenced comments (i.e. it is a decent choice when the date of the most recent event is not known).\nIn summary, this is interesfing and high quality work. Having communicated with some of the authors on the subject of earthquake stafisfics and forecasfing previously, I know that their desire to do good science is genuine. I think this work too will soon be worthy of publicafion, but given its nature, some considerably more nuanced discussion of why this approach is befter than the alternafives is required.\nRandy Williams\nLine Referenced Comments\n18-20: Interesfing quesfion really - does befter incorporafion of uncertainty lead to a more robust forecast? Depends on what a “robust” forecast is. Seems at least possible that we end up with a forecast that is no more or less accurate than was the case prior to using this Bayesian approach.\n21-22: Again, this is cool, but minimizing the forecast uncertainty is only a useful thing if the forecast itself is accurate. The central dilemma of earthquake forecasfing! We’ll need decades (if not centuries) to find out.\n40-47: Could simplify this secfion for the novice with some more plaintext terminology. The Poisson process describes events that effecfively occur randomly in fime. And the anfithesis of that for periodic sequences of course.\n57-60: Maybe rephrase to befter emphasize the severity of the problem. I’d wager that relafively few of those 90 records actually have anywhere near 20 disfinct events. In my experience records with more than 10 events are quite rare. Many have fewer than 5.\n60-61: It might, certainly, but does it? Again, we really don’t know at this point.\n77: First occurrence of (what I assume is) Brownian passage fime. Please define.\n98-99: See comment above about record length.\n108: Perhaps this comment is too “in the weeds” but I feel like we have known for a while that a lognormal distribufion of interevent fimes is not an appropriate model for earthquakes, largely because its hazard funcfion actually decreases through fime from some maximum value.\n147-150: Quite interesfing though.\n169: Could just call this “27” fime as likely?\n207-222: Some of this discussion makes one speculate as to whether variafions in strain rate or kinemafics have an appreciable impact on record length and/or event preservafion. Maybe worth a teaser sentence?\n211-215: Prefty sure that Berryman made a similar argument in the 2011 Alpine fault paleoseismology paper. Should probably give them a citafion.\n251-252: I think this statement needs a bit of jusfificafion. Why is having smaller errors necessarily preferable from a hazards forecasfing perspecfive? As stated above, it’s only really a good thing if the forecast is accurate.\n253: Although I am also in favor of fime-dependent approaches, we can probably cut the Poisson model a liftle bit of slack here. It does have its uses from a hazards perspecfive, parficularly if we have some sense of the average interseismic period length but not a date on the most recent slip event. We see in the retrospecfive forecast that in a select few cases smaller errors caused the analysis to miss the most recent event.\n235-240: Yes, the Poisson model does yield larger uncertainfies and thus more frequent “overlap” with the actual events, but to be fair, the model average approach also yields fairly large uncertainfies in an absolute sense. In the example of the San Andreas, if the forecast produces a 95% confidence interval for the next event that is a period of several hundred years, whilst the actual data suggest that a typical recurrence interval is on the order of about 100 years, it is a bit unsurprising that the forecast would capture the event within its uncertainty. Respecffully, I feel that I have to ask what exactly we learn from this? Certainly it could be interpreted to show that your approach isn’t fundamentally off base, but why is this necessarily befter than all of the other hazards forecasts that are available for the San Andreas (or other fault systems)?\nAlso, the bit about the model averaging approach “outperforming” the Poisson is a bit complicated as wriften. You state that the model averaging approach successfully covered 79 of the true occurrence fimes, but then state that the model averaging approach outperforms Poisson because 81 of the model averaging results had lower mean squared errors. Does that not then by definifion include 2 records for which your retrospecfive forecast failed (i.e. the CI did not include the true occurrence fime)? Strange to count those as successes of the approach.\nFinally, please consider comments above regarding the ufility of the Poisson process in cases where the date of the most recent event is unknown. There is pracfical ufility in the fime-independent approach in some cases, even if we believe that fault behavior is, in fact, fime dependent.\n251-252: Again, this is a really tricky area of applicafion. Please jusfify why this is favored, remembering WHY we want to forecast earthquakes accurately in the first place. We could of course be conservafive/caufious and assume that any place that has a large fault that moved within the last 50 ka has a very high earthquake risk, but that would not do good things for local economies / insurance / building budgets. Although I do like this manuscript, I feel it has a tendency to treat seismic forecasfing as liftle more than an intellectual exercise by focusing on what is most “correct” about fault behavior and uncertainty minimizafion. Neither of these factors are readily translated into accuracy/reliability/usefulness of the forecast, although they may be components of them. In a pracfical sense, its less about correctly predicfing an earthquake within some interval of fime. Rather, it's about calculafing a per year likelihood of a catastrophic event for a parficular area. Your approach by definifion maximizes that likelihood by minimizing the predicted interval.\n265-266: I did not see a comparison between the model averaging approach and the “single best model” approach in the retrospecfive forecasfing data?\n267-269: This might be a stronger point to the work that is not being fully appreciated by the line in quesfion. If we are ever to validate whether our hazards forecasts are meaningful (in the absence of waifing a few hundred years to find out), it may be best achieved by looking at the global distribufion of forecasts and how often they ring true.", "reviewer_author_response": "RESPONSES TO REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have proposed a Bayesian model-averaging approach for forecasting large\nearthquakes along active faults. They emphasise the importance of model uncertainty in the\nconventional single-model approaches and ensemble all the candidate models with their\nweights derived from the model selection criterion WAIC to mitigate the effect of model\nuncertainty on forecasts.\nSelection of recurrent-time distribution has long been controversial for long-term earthquake\nforecasts. The authors argue that no single model consistently provides the best fit across the\nglobal catalogs, or even within the individual fault systems. Their approach presents a plausible\nsolution in this situation. Additionally, the authors discussed the periodicity of earthquake\nrecurrence and revealed through regression analysis that the degree of periodicity depends not\nonly on the occurrence rate (as discussed in the previous studies) but also on the faulting type\nand the tectonic region. The authors provide the probability forecasts in the next 50 years for\nthe global 93 fault segments. Their retrospective analysis shows good accuracy of the proposed\napproach for most of the fault segments.\nThe paper is well organised and insightful. However, I'd like to provide feedback on several\nspecific points:\nL86: The authors describe that the AIC can’t provide guidance as to how much better a\nparticular model is compared to other models. However, since the AIC is an estimator of the\nprediction error of a fitted model in terms of the log-likelihood of a new observation, it can\nindeed assess the relative predictive competence of different candidate models.\nReply:We agree that this statement is misleading, and we have therefore removed it.\nL148: Though the authors have revealed relationships between the number of events and the\nsingle-best model in their catalogs, they do not provide discussion on characteristic behavior\nof earthquake recurrence in each single model. For example, calculating the normalised\nrecurrence intervals, which are divided by the mean recurrence intervals of the respective fault\nsegments, grouping them by the selected single-best models, and comparing the histograms\nbetween the single-best models might help us to observe some characteristic recurrence\nbehavior in each model. In this case, I suggest referring the following paper:\nNishenko, S. P., & Buland, R. (1987). A generic recurrence interval distribution for earthquake\nforecasting. Bulletin of the Seismological Society of America, 77(4), 1382-1399.\nReply:We thank the reviewer for this suggestion. We carried out this analysis as suggested by\nthe reviewer. For each fault segment, wecalculated the normalised recurrence intervals for each\nMC sample, and then reported the median and the 2.5% and 97.5% quantiles of the standard\ndeviations of the normalised recurrence intervals calculated for the 100 MC samples. However,\nwe didn’t find any characteristic recurrence behaviour for the fault segments that have the same\nsingle-best model. See the figure below. We did observe that fault segments with higher rate\nof earthquake occurrences appear to have smaller standard deviations of the normalise\nrecurrence intervals. We have added this discussion in the revised manuscript, Lines 196-203,\nwhich is also included below.\n“Previous studies (Sykes and Nishenko, 1984; Nishenko, 1985; Nishenko and Buland, 1987)\nhave shown that the standard deviations of the scaled (divided by the mean value) inter-event\ntimes along several fault segments appear to be constant. For each of the 93 fault segments,\nwe calculated the scaled inter-event times for each MC sample, and then reported the median\nand the 2.5% and 97.5% quantiles of the standard deviations of the scaled inter-event times\ncalculated for the 100 MC samples. There does not seem to be any characteristic recurrence\nbehaviour for the fault segments that have the same single-best model. Fault segments with\nhigher rate of earthquake occurrences appear to have smaller standard deviations of the scaled\ninter-event times (Supplementary Fig. 4).”\nL207-222: The ranges of the CIs in this paragraph appear too short compared to the sample\nsize of 93 and the dispersion of the data points in Figure 4. I suspect that this might be because\nall the 100 MC samples and/or 5000 MCMC samples for each segment are used as the data\npoints in evaluating the CIs. If this is the case, it is not appropriate to discuss the statistical\nsignificance using the CIs from these leveraged samples.\nReply: We are grateful for the reviewer's comment related to this part of the analysis. This\nmade us realise that we had not allowed for both sources of errors in the regression model:\nprocess error (which would be present even if we knew the true values of the shape parameter)\nand sampling error, i.e. the uncertainty in our estimation of the true value of the shape\nparameter. In the revised manuscript, we added both sources of errors in the regression model\nand have updated the CIs which are wider than in the original manuscript. See Lines 273-296\nin the revised manuscript, which is also included below. R code related to all the data analysis\nis provided in the supplement.\n“A regression model for the relationship between the shape parameter of a Weibull renewal\nprocess and the earthquake rate, tectonic setting, faulting type, and the number of earthquakes\nof each of the 93 fault segments suggests that when the earthquake rate increases 10 fold (e.g.,\nfrom 0.0001 to 0.001), the shape parameter of the Weibull renewal process increases by about\n24% (95% CI: 5%-45%) (Figure 4). This suggests that fault segments with higher earthquake\nrates tend to have more periodic behavior. The shape parameter of the Weibull renewal process\nfor fault segments located in an intraplate setting is about 87% (95% CI: 51%-148%) of that\nfor fault segments located at or near plate boundary. Although this the 95% CI is wide and\ncovers 100%, the posterior density plot (Supplementary Fig. 5) suggests that there is a high\nprobability (0.7) that the latter may be more periodic than the former. There are only five fault\nsegments from an intraplate setting, so to reach a more robust conclusion, more data from\nintraplate fault segments are needed. The shape parameter of the Weibull renewal process for\nfault segments located in an active intraplate setting (predominantly faults in central China) is\nabout 1.4 (95% CI: 1.1-1.9) times of that for fault segments located at or near plate boundary,\nsuggesting that the former appear to be more periodic than the latter. The shape parameter of\nthe Weibull renewal process for fault segments located in a subduction region is about 1.9 (95%\nCI: 1.1-3.3) times of that for fault segments located at or near plate boundary, suggesting that\nthe former appear to be more periodic than the latter. On average, the Weibull shape\nparameter for reverse faults is about 69% (95% CI: 46%-104%) of that for normal faults,\nsuggesting that large earthquakes recur more periodically on normal faults than on reverse\nfaults.”\nL244: From Figure S4, the authors have shown that the forecast values from the model-\naveraging approach have much smaller mean square errors (MSEs) than those from the Poisson\nprocess. However, it is difficult to see the relative improvement in accuracy from this figure\nbecause the fault segments have a wide range of scales in their recurrence times. Therefore, I\nsuggest adding a figure and comments on the relative differences in MSEs between the model-\naveraging approach and the Poisson process.\nReply:We have included Supplementary Fig. 7 in the revised supplement to show the relative\ndifferences in MSEs between the Poisson process and the model-averaging approach. We have\nalso added the following comments in the revised manuscript, Lines 319-321.\n“For about half of the fault segments (45 of the 93), the MSEs of the forecasts from the Poisson\nprocess are more than twice of that from the model-averaged forecasts (Supplementary Fig.\n7).”\nL252: favours -> favour ?\nReply:Thank you. We have changed this.\nL265: The accuracy of forecasts for the single-best model approach does not seem to be\ndiscussed in the retrospective tests in Section 5.\nReply: Thank you. We have added Supplementary Fig. 8 in the revised supplement to show\nthe relative differences in MSEs between the single best model and the model-averaging\napproach, and added this discussion in the revised manuscript. See Lines 336-341, which is\nalso included below.\n“The MSEs of the retrospective forecasts from the single-best model approach are very close\nto those from the model-averaging approach for the majority of the fault segments\n(Supplementary Fig. 8), all within two times relative difference. Model-averaging with WAIC\nweights is not usually designed to achieve a better MSE than the single best-model approach.\nHowever, when there is some uncertainty as to the best model, model-averaging outperforms\na single best-model primarily in terms of better representing all the uncertainties.”\nL292: The sentence ends with ‘..’.\nReply:Thank you. We have fixed this.\nReviewer #2 (Remarks to the Author):\nThe paper aims at introducing a new methodology to estimate large earthquake occurrence on\nfaults. Although the topic is of great interest and the paper contains a very interesting approach,\nI have some concerns that should be carefully addressed; some of them might be critical.\nBelow I list my main concerns.\n- This is likely the most critical point. In the \"Introduction\" section, the authors did not quote\nrecent bibliography on the earthquake recurrence models on faults, such as the recent Neely et\nal (SRL, 2023) and the work made in UCERF3 (Field et al.", "reviewer_ci)": "nBig Bend: 41% (36%, 46%)", "reviewer_the\nprocess\nreply": "Thank you for all your constructive suggestions and comments, which have improved\nour manuscript significantly.\nRandy Williams\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have addressed my concerns well and revised the manuscript accordingly. I believe that the\nreviewers' concerns have been sufficiently addressed and that the revised manuscript is worthy of\npublicafion as it stands.\nReviewer #2 (Remarks to the Author):\nThank you to the authors to have addressed safisfactorily all my main concerns.\nReviewer #3 (Remarks to the Author):\nThank you for considering my comments and addressing my last remaining concern with the manuscript.\nI agree that the forecast comparisons that you are able to provide are not \"apples to apples\" in terms of\ntheir temporal formulafion, but I do think that the comparison is sfill useful. The revised text does a good\njob of discussing those issues succinctly. I appreciate your efforts to provide these data, and I have\nnothing further to add in terms of comments. Good luck with the remainder of the publicafion process!\nRandy Williams"}, "author_response": ""}]

Decision letter