We thank the referee for the careful review and the helpful remarks. In the following, we address all comments and outline the foreseen amendments of the manuscript. After having received comments from additional reviewers, we will thoroughly state all the amendments in our final rebuttal and upload the updated manuscript.  

Comment 1: Line 46: there is not the reference year for Baker and Chen (it should be 2020).

Response: Thank you. We will add the reference year.

Comment 2: Line 72: do authors refer to the case of the same IM at different sites? Please specify (if different IMs at multiple sites are considered, then a spatial-cross correlation model is needed to model correlation among GMM residuals).

Response: Yes, the manuscript focuses on correlations of the same IM at different sites. We will carefully specify this in the revised manuscript. For further information, please see the response to the following comment.

Comment 3: In fact, following previous comment, it is not clear whether the correlation models described in Section 2 can apply when more than one IM at different sites are of concern. If it is not the case, can the models be extended? The “model E” of course could consider the case of more IMs (literature gives spatial-cross correlation models, one of which published by the second author of this manuscript; i.e., Loth and Baker, 2013), but what about “model EA” and “model EAS”? Please discuss.

Response: We thank the reviewer for raising this point. As mentioned in our response to the previous comment, the present study focuses on correlation of the same IM at different sites. However, all three models can be extended to include spatial cross-correlations between multiple IMs. Specifically, the linear co-regionalization technique, used by Loth and Baker (2013) for isotropic models, can also be employed for models EA and EAS. In response to the reviewer’s comment, we will include additional clarifications in the revised version of the manuscript.

References

• Loth C, Baker JW. (2013) A spatial cross-correlation model of spectral accelerations at multiple periods. Earthquake Engineering and Structural Dynamics.

Comment 4: Lines 139-140: “We also tried a model with two separate isotropic components E1 and E2 […]” Can the authors better explain the meaning of E1 and E2 components?

Response: We will precise this statement in the revised version. Equation (7) states the model EAS as it is used in the manuscript. This model has one isotropic component ρ_E that is multiplied with both, ρ_A and ρ_S respectively. Additionally, we experimented with a model that has two separate isotropic components,  ρ_E1 and ρ_E2 , each with two parameters. One component is multiplied with ρ_A, while the other is multiplied with ρ_S. The increased flexibility of this model comes at the cost of two additional parameters that need to be estimated. The results indicated that this additional complexity does not improve the performance.    

Comment 5: Line 155: A reference about the mixed-effects regression should be provided.

Response: We will add a carefully selected reference to the revised manuscript.

Comment 6: Line 156: “a maximum usable period within the period range of interest.” What is the period range of interest?

Response: In Section 4.2 we estimate correlation models for spectral accelerations at nine vibration periods in the range [0.01s, 6s], which is the period range of interest. While we explain this in line 274, we agree with the reviewer that readers might benefit from a similar statement in line 156. We will modify the manuscript accordingly.  

Comment 7: Pag 7, line 162: “For Sa(1s), the pooled data set consists of 13,342 records from 128 events.” Do the number of data vary with the considered IM? Can the authors provide some additional explanation about that?

Response: The size of the data set varies depending on the period. As stated in line 156, we restrict the entire set of records based on the maximum admissible period. Specifically, this results in lower numbers of available records for longer periods. In response to the reviewer’s comment, we will specify this in Section 3.1 of the revised manuscript.

Comment 8: Figure 2: Although some characteristics of the figure are discussed in the paper, it is not clearly motivated that fact that the two panels of the figure have a similar trend. For instance, is it somehow expected that the number of stations pairs at Euclidean distance larger than 35km and angular distance lower than 20° are significantly larger than those with the same Euclidean distance and angular distance larger than 160? A similar reasoning can be done for station pairs with low and large soil dissimilarity. Can the authors add some comments about that?

Response: We thank the reviewer for this question. The following explanation focuses first on joint bins of Euclidean distance and angular distance (Figure 2a), followed by joint bins of Euclidean distance and soil dissimilarity (Figure 2b).

For geometrical reasons, station pairs with small angular distances are more likely than station pairs with large angular distances. Figure R1 illustrates this for joint bins of Euclidean distance from 35 km to 40 km, and angular distances from 0° to 20° (a), and 160° to 180° (b), respectively. The red shaded area indicates the joint bin and we observe that this area is larger for smaller angular distances. Thus, the chance of a second station being in the joint bin is larger for smaller angular distances. We agree with the reviewer that this fact was missing in the submitted version of the manuscript. We consider adding Figure R1 as an Appendix and referring to it in Section 3.1.

Figure R1. Joint bins of Euclidean distance, d_E , and angular distance, d_A , for a reference site: Euclidean distances from 35 km to 40 km, and angular distances from 0° to 20° (a), and from 160° to 180° (b).

For bins of soil dissimilarity and Euclidean distance, our explanations for the trend in Figure 2b are twofold: First, we note that it is more likely to observe similar vs30 values for pairs that are close to each other. Note that [35, 40] km is the largest Euclidean distance bin plotted in Figure 2. Second, low soil dissimilarities are more frequent in the training data set described in Section 3.1 with soil dissimilarities larger than 400 m/s accounting for less than 2% of all records. In the revised manuscript, we will add further details to our description of Figure 2b.

Comment 9: Lines 191-195: To increase the readability of the work, authors may consider providing some additional details describing the simulation procedure apart from the cited references.

Response: Parameter inference was performed with Markov-chain Monte Carlo (MCMC). MCMC is a sequential sampling algorithm that is often used in Bayesian statistics to draw samples from a certain target posterior (Neal, 1993). Specifically, MCMC algorithms generate correlated series of samples that will converge in distribution to the target posterior. While there exist many different MCMC algorithms, we implemented the No-U-Turn Sampler algorithm that allows for an efficient exploration of the parameter space without requiring costly tuning runs (Hoffman and Gelman, 2014). Implementation details can be found in the cited Zenodo repository (see "Assets" on NHESS). In response to the reviewer’s comment, we will add further details to Section 3.4.

References

• Hoffman, MD. and Gelman, A. (2014) The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo, Journal of Machine Learning Research
• Neal, RM. (1993) Probabilistic Inference Using Markov Chain Monte Carlo Methods, Tech. rep., University of Toronto, Canada

Comment 10: Line 196: “The results shown in Fig. 3 illustrate slight differences in the three event-specific models”. Are the differences so slight? They seem significant, at least from a graphical point of view. Please explain.

Response: We agree with the reviewer. Indeed, we think that the term “slight” is not an accurate description of the differences in the three event-specific models. We will either delete this term or explain the differences in more detail.

Comment 11: Line 190 states that both paragraphs 3.4.1 and 3.4.2 refer to Sa(1s), so the parameters given in Table 2 should refer to such an IM, shouldn’t they? If it is the case, it should be recalled in the table’s caption. This comment also applies to Figure 4, which is supposed to refer to Sa(1s).

Response: We will modify the captions of Table 2 and Figure 4 accordingly.

Comment 12: Lines 210-220 give a brief discussion about the trend of the correlation models represented in Figure 4, which should refer to Sa(1s) (see also previous comment). Have authors investigated such trends for other spectral ordinates? Authors should briefly mention if and how the trend of the correlation models vary with the considered spectral ordinate. For example, is the discussed similarity between “model E” and “model EAS” at small angular distances (vibration) period-invariant?

Response: The authors investigated such trends for other spectral ordinates. In general, the similarity between model E and model EAS at small angular distances holds true for all examined vibration periods. As we state on line 220, this is because most station pairs with Euclidean distances smaller than 40 km (the interval shown in Figure 4) have small angular distances.

For the other discussed trends, we would like to refer to Section 4.2, where we investigate the role of site and path effects. Compared to Sa(1s), shown in Figure 4, the importance of site effects increases with periods longer than one second, and decreases with periods shorter than one second. This is also reflected by the estimated parameters for model EAS. Specifically, the weight parameter (w) decreases for longer periods, which means the model assigns more weight to the site-effect term (i.e., the term that accounts for soil dissimilarities), compared to the path-effect term (i.e., the term that accounts for angular distances). In response to the reviewer’s comment, we will add a notebook to the cited Zenodo repository (see "Assets" on NHESS) that shows Figure 4 for other spectral ordinates.

Comment 13: Section 4.2 discusses the pooled model performance for other spectral ordinates. Why don’t authors provide the parameters of the correlation models for all considered spectral ordinates? This would be helpful to the reader interested in implementing the EAS correlation model.

Response: Thank you for pointing out this issue. We now provide the parameters for all spectral ordinates in the GitHub repository, and we will add a new copy to the Zenodo repository (see "Assets" on NHESS) when submitting the revised version. We intended to provide this file already at submission and apologize for any inconveniences caused. In addition, we consider providing the parameters as a Table in the appendix of the revised version.

Comment 14: Section 6: Choosing the same threshold value for all the sites of an area means that the return periods of the chosen thresholds are different from site to site and are all dependent on the rate of earthquakes of the seismic source. In fact, it was demonstrated that the effects of spatial correlation models may have different significance varying the return period of the considered thresholds. Thus, it may be relevant to define threshold values avoiding the dependency of the source rate. To this aim, a possible strategy is discussed in Cito et al. (2023).

Response: We thank the reviewer for pointing us to the recently published study of Cito et al. (2023). We agree with the reviewer that the choice of appropriate thresholds can be important. For that reason, we performed an extensive sensitivity study with different threshold values prior to submission. As we state on line 337, varying threshold values did not affect the major trends in the different correlation models that we reported in Section 6.

The reviewer suspects that our results depend on the rate of earthquakes of the considered seismic source. To reply to this comment, we first briefly summarize the study of Cito et al., before discussing some key differences to the case study presented in our manuscript.

Cito et al. evaluated the accuracy of an approximate method to generate spatially cross-correlated samples of different ground-motion intensity measures at multiple sites. For that purpose, they compare the joint exceedance of a threshold value conditional on all possible rupture scenarios in a specific region over a certain time horizon. In that context, they show the importance of choosing site-specific threshold values that are independent of the rates of earthquake events on a particular seismic source.

In our case study, we generate spatially correlated samples of the same intensity measure at multiple sites. For the different proposed correlation models, we compare the joint exceedance of a threshold value conditional on one specific rupture scenario (Figure 9), instead of all possible rupture scenarios that are generated by a particular source over a certain time horizon. Thus, the rates of earthquake events do not affect the presented results.

As stated above, we experimented with different threshold values. These experiments also included cases with site-specific threshold values. Those thresholds were chosen such that they have identical site-specific exceedance probabilities for the considered rupture scenario. We found that the reported trends for the different spatial correlation models do not depend on the threshold values. In response to the reviewer’s comment, we will include additional explanations in Section 6 of the revised manuscript.

References

• Cito P., Chioccarelli E., Iervolino I. (2023) Conditional hazard for simplified multi-site seismic hazard and risk analyses. Earthquake Engineering and Structural Dynamics

Comment 15: Section 6 is entitled "Implications for regional seismic risk", On the other hand, the presented case-study deals with seismic hazard rather then risk (Figure 10 factually shows hazard curves). Although hazard is necessary for risk assessment, authors should specify that the results they find pertain specifically to the hazard side. In fact, the way the covariance matrix of IMs is modelled does not necessarily affect hazard and risk results in the same way: for example, the difference between the curves obtained via the “E” and “EAS” correlation models discussed in the manuscript can be different when (a metric of) risk is considered in lieu of hazard.

Response: We agree with the reviewer that Section 6 may benefit from additional clarifications on the difference between hazard and risk, which we will add in the revised version of the manuscript.

We thank the reviewer for the careful review and appreciate the helpful remarks for further improvement of our manuscript. In the following, we provide preliminary replies to all comments and outline the foreseen amendments of the manuscript. We will thoroughly state all the amendments in our final rebuttal once the open discussion is over.

Comment 1: Introduction: As the topic focuses on site and path effects in the residuals of a GMM and consequently their spatial correlation, a discussion could be held on the different approaches in the literature to capture such anisotropic and non-stationary effects in spatial ground motion modelling. The issue of spatial correlation is also critical for non-ergodic analyses that are based on observations, however, it has recently been shown that the exponential growth of seismological records in the vicinity of earthquakes provides an opportunity to separate the source, regional propagation and site factors that control the variability of ground motion, thus not only from physics-based simulations. For example, in the non-ergodic framework there are several approaches to capture repeatable site and path effects (useful references can be found in Sgobba et al. Earthquake Engng Struct Dyn. (2021) - https://doi.org/10.1002/eqe.3362; Parker and BaltayBSSA (2022) - https://doi.org/10.1785/0120210175; Lavrentiadis, et al. Bull Earthquake Eng (2022) - https://doi.org/10.1007/s10518-022-01485-x, among others). The authors make some reference to this point in section 2.2, but without a comprehensive distinction between the spatial modelling approaches in the two main contexts (ergodic and non-ergodic) and their implications. This is a crucial aspect for ongoing research, so the manuscript would benefit from a fuller discussion on this point.

Response: We thank the reviewer for raising this point. While we state some differences between ergodic and non-ergodic models in Section 2.2, we agree that this aspect deserves further discussion. Therefore, we plan to add a paragraph to the introduction of the revised manuscript.

Comment 2: Section 2.1: The authors introduce the gamma-exponential function as an isotropic model based on Euclidean distance. Why did they choose this type of model? This assumption should be better justified. Furthermore, I recommend that the authors refer to and consider the notations of previous studies or define which terms of their model are the same as the notations of, for example, Schiappapietra and Douglas (2020), which define the most typical parameters of the exponential model, i.e. sill and the range of the correlation model;

Response: The gamma-exponential model was used in the present study because of its generality. It includes the exponential model (as referred to by the reviewer) as a special case. We also note that a similar model was employed in the ground-motion correlation studies of  Heresi and Miranda (2019), and of Goda and Hong (2008). The gamma-exponential model, defined in Eq. (3) of the manuscript, has two parameters: The lengthscale, l, and the exponent, γ. The lengthscale is the distance at which the correlation is 0.368.

The exponential model, used for example in Jayaram and Baker (2009), is commonly defined as exp(-3d_E / b), where b is the distance at which the correlation is 0.05 and commonly referred to as the correlation range. We note that this model is equivalent to the gamma-exponential model with an exponent of one and a lengthscale of b/3.

In response to the reviewer’s comment, we will specify the above-mentioned relation to the exponential model and to the range parameter in Section 2.1 of the revised manuscript.  

References

• Heresi P., Miranda E. (2019): Uncertainty in intraevent spatial correlation of elastic pseudo-acceleration spectral ordinates. Bulletin of Earthquake Engineering.
• Goda K., Hong HP. (2008): Spatial correlation of peak ground motions and response spectra. Bulletin of the Seismological Society of America.
• Jayaram N., Baker JW. (2009): Correlation model for spatially distributed ground-motion intensities. Earthquake Engineering and Structural Dynamics.

Comment 3: Section 3.1: please, specify which type of Vs30 is considered in the study to compute the soil dissimilarity. are they homogeneous estimates from direct measurements or inferred from other quantities?

Response: We thank the reviewer for this question. We used vs30 values as provided in the NGA-West2 database (Ancheta et al., 2014) and did not separate between stations with measured or inferred vs30 values. We will clarify this in the revised version of the manuscript.

References

• Ancheta TD., Darragh RB., Stewart JP. et al. (2014): NGA-West2 Database. Earthquake Spectra.

Comment 4: Section 3.4: explain why Sa(1s) was considered as an example parameter in the study. Furthermore, it would be interesting to add some comments about the variability of the results with the vibration period;

Response: The focus on one spectral ordinate as an example parameter allowed us to explain the different steps in sufficient detail. We chose Sa(1s) because the analysis presented in Section 4.1 showed that it offered a balanced mix in terms of the relative importance of site and path effects. We agree with the reviewer that this aspect deserves further clarification in the revised version of the manuscript.

Based on this comment, and comment #12 from the first reviewer, we will add supplementary material that illustrates the estimated models E and EAS for other spectral ordinates (similar to Figure 4). In the revised manuscript, we also plan to add a brief summary of this supplementary material in Section 3.4.2.

Comment 5: Section 4: The validation procedure is not straightforward. Some concepts are briefly introduced in Section 4 but are then clarified in Section 4.1 (e.g. it is not immediately clear what the basic BL model is); why did the authors choose to apply an in-sample and an independent model performance later? Why did they not choose a cross-validation technique? Please give reasons for these choices;

Response: We thank the reviewer for raising this point. Based on the reviewer's comment, we will improve our explanation of the different comparisons that we performed. As stated by the reviewer, we compare the models in terms of their in-sample performance (Section 4), and out-of-sample performance (Section 5). The former allowed us to compare the performance of event-specific and pooled models in Section 4.1. A crucial step to address the question whether the previously found event-to-event variability in isotropic model parameters may be explained by the lack of accounting for site and path effects.

In Section 5, we compare the model performance on unseen, out-of-sample data from the 2019 Ridgecrest earthquake sequence. We chose this approach (rather than a cross-validation technique) for following reasons: First, the ergodic GMM, used to compute the within-event residuals, differs for both data sets. As stated in Section 3.1, we used the GMM of Chiou and Youngs (2014) for the NGA-West2 training data set, while the GMM of Liu et al. (2022) was used for the Ridgecrest data set (explained in Section 5). Second, choosing a suitable strategy to perform the train/test splits for cross-validation is not trivial. Special care would be required to ensure that the resulting splits have a balanced number of events and records, while remaining representative of the overall data set. Finally, we note that the choice of using a fixed test set was also motivated by the similar procedure used in the study of Liu et al. (2022). 

References

• Chiou BSJ., Youngs RR. (2014): Update of the Chiou and Youngs NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra. Earthquake Spectra.
• Liu C., Macedo J., Kottke AR. (2022) Evaluating the performance of nonergodic ground motion models in the ridgecrest area. Bulletin of Earthquake Engineering.

Comment 6: Section 4.2: results in Figure 7 are interesting as they show the improved performance of EAS model compared to EA in capturing the site effects at long-periods. On the other side, it appears that the path effects are better captured by EA model at short-periods compared to the isotropic model E (it can be noted a peak at T=0.1s). Please, provide a comment on this point;

Response: We thank the reviewer for this remark. Indeed, the path effects (as captured by models EA) appear to be more relevant for periods shorter than T=0.3s, and we will include a comment on this aspect in the revised version of the manuscript.

Comment 7: Section 6: why is there talk of regional risk implications, when the assessments being made are about hazard? certainly the analysis shown can be seen in a broader context of risk assessment, but i think the terminology adopted here may be misleading;

Response: We agree with the reviewer that Section 6 would benefit from additional clarifications on the difference between hazard and risk, which we will add in the revised version of the manuscript.

Comment 8: Section 6: in this case-study, there are no plots to show maps of the correlation fields. I recommend the authors to plot the EAS model to highlight the resulting spatial correlation structure in the region and extend the considerations at least at two IMs (at short and long-periods, respectively).

Response: In response to the reviewer’s comment, we will add maps of the correlation coefficient for models E and EAS to Figure 9 of the revised manuscript (see Figure R1 below for a first draft). We will also provide a supplementary material which compares the results presented in the manuscript (Figures 10 and 11 for Sa(1s)) to at least two other IMs, and we plan to add a brief summary of these results to Section 6.

Figure R1: Map of the case study area used for regional risk assessment and the considered M6.25 rupture on the Hayward fault: (a) Soil conditions, quantified via vs30, and correlation coefficient from the indicated reference site to all other sites obtained with posterior mean parameters of model E (b) and model EAS (c). Numbered boxes indicate the subregions of interest (blue areas are water).

We thank the reviewer for the careful review and appreciate the helpful remarks for further improvement of our manuscript. In the following, we provide preliminary replies to all comments and outline the foreseen amendments of the manuscript. We will thoroughly state all the amendments in our final rebuttal once the open discussion is over.

• Comment 1: The authors have primarily compared the three models against each other; however, the absolute goodness-of-fit and prediction power are not clear. The authors could use conditional R²or MAP to obtain R² for presenting the statistical goodness of fit. Relative improvements are not helpful if the models do not possess high absolute predictive power on unseen/test datasets.

Response: We thank the reviewer for raising this point. We agree that there are many metrics to assess the goodness-of-fit. The R-squared metric quantifies the proportion of total variation about the mean explained by a regression model (Draper and Smith, 1998). As such, it would be a useful metric for validating a novel ground-motion model (GMM), where one estimates systematic and repeatable effects. In the submitted study, however, we focus on the spatial correlation of scaled within-event residuals. As explained in Section 2 (lines 72-81), the proposed correlation models estimate the off-diagonal entries of the correlation and covariance matrices. In order to assess and compare different models, we need a metric that accounts for the full joint distribution, including these off-diagonal terms. The logarithmic posterior predictive density (LPPD), as defined in Eq. (9), is such a metric. The R-squared metric, on the other hand, focuses solely on the variance at each data point, i.e., the diagonal entries of the covariance matrix. In response to the reviewer’s comment, we plan to add a summary of the above stated justifications to the revised manuscript.

We also note that all absolute values of LPPD can be found in the cited online repository (see "Assets" on NHESS). In the manuscript, however, we prefer to compare the LPPDs of different models rather than stating their absolute values. Absolute values may be helpful and important if they are applied to tasks and data sets for which previous studies reported the same metric. In that case, there may be some common perception on absolute values that indicate sufficient or insufficient performance, such as for the MNIST data set for visual recognition of handwritten digits (LeCun et al., 2010). Because our study is the first one that applies the LPPD metric to spatial correlation models, absolute values would lack context and reference points. 

• Comment 2: In section 2, the authors mention that their spatial model is based on GMM’s functional form with between-event and within-event residuals, then they have developed event-specific and pooled models too. This isn’t very clear:
  • Comment 2a: If the models are based on mixed effects, then given the Bayesian MCMC is a sampling approach, the authors can try to test the assumption of independence between the two residuals, as mentioned in lines 65-70. It would be much more interesting to see the results if the correlation between the two types of residuals is modeled, unlike the frequentist mixed-effects modeling approaches used conventionally.

Response: We thank the reviewer for this comment. Our study focuses on spatial correlation of within-event residuals. The proposed spatial correlation models can be combined with an existing GMM to perform regional risk assessments. The separation of between- and within-event residuals is common in existing GMMs, see for instance Chiou and Youngs (2014). The development of a novel GMM that may feature correlated within-event and between-event residuals is considered as out of scope of the present study.

• Comment 2b: If the models are only developed as event-specific and pooled, then the authors should justify the reason. Hierarchical modeling using Bayesian approaches is relatively easy as compared to frequentist approaches. Why haven’t the authors separated the two types of residuals since the multiple records from multiple events need a hierarchical kind of modeling?

Response: As mentioned in our reply to comment 2a above, the present study focuses on spatial correlation of within-event residuals. In this context, we refer to event-specific models as models that are estimated with the within-event residuals from one event only (such as in the study of Heresi and Miranda, 2019). A pooled model, on the other hand, is estimated with the within-event residuals from multiple earthquake events (such as in the study of Esposito and Iervolino, 2012). In both cases, however, we did separate the two types of residuals based on the employed GMM.

• Comment 3: The authors have used the three parameters to model the correlations. It is clear from Figure 5 that the inclusion of three parameters model is better; however, it would be interesting to see a sensitivity analysis of which parameter has the most impact on the permutations. Also, using a different weighted objective function for the three factors can be checked.

Response: We thank the reviewer for raising this point. In Section 4.2, we compared three different models: A model that only accounts for spatial proximity (model E), one that additionally accounts for path effects (model EA), and finally the model that additionally accounts for site effects (model EAS). The objective of Section 4.2 is to assess the relative importance of the three effects: spatial proximity, path and site effects (see line 238). In lines 280-281, for example, we state: “The major part of the benefit in LPPD stems from the isotropic model E, although accounting for path and site effects … leads to further increases in LPPD”.

As explained in Section 3, we compute the posterior distribution of the model parameters based on Bayes’ theorem. As such, we do not prescribe specific weights for an objective function. If the reviewer refers to how we combine the above mentioned effects in the functional form of model EAS: In Section 2.3, we describe the different functional forms that we tested for model EAS. Please see also our reply to comment #6 below. 

• Comment 4: Apart from the considered parameters, including wall effects could have made the models more physics-oriented. E.g., footwall vs. hanging wall effects. Also, besides Vs30, basin effects could be valuable to improve the model, e.g., parameters like Z1 or Z2.5. These parameters are accessible in the NGAWest2 database and should have been included in the models or tested.

Response: We agree with the reviewer that there are many variables, which could be additionally included in models similar to the ones presented in this study. We reflect on our choice of vs30 as a proxy for soil conditions in Section 2.3 (lines 129-131), which was not only motivated by its availability in the training data set (NGA-West2), but also in the case study regions employed in this study. The proposed model EAS is based on three dependent variables (Euclidean distance, angular distance and soil dissimilarity), whereas most previous spatial correlation models were exclusively based on Euclidean distance. We discuss the limitations of these dependent variables in the conclusion (lines 378-379), but consider models with additional variables as out-of-scope of the present study. Instead, we hope that our work motivates future studies that explore refined model parametrizations and additional dependent variables.

• Comment 5: The dataset is not clearly explained and explored in the paper. Figure 2 and some results indicate biases in the datasets. The authors should consider sampling the dataset in a more balanced manner to minimize any biases. Furthermore, information about the distributions of the dataset parameters can help understand the study.

Response: The revised version of the manuscript will include further details on the trends seen in Figure 2. As we argue in our preliminary reply to comment #8 from the first reviewer, these trends, or apparent “biases”, have geometrical and geological reasons. Note that we use Bayesian inference to estimate model parameters and the uncertainty in the posterior parameters reflect the above-mentioned trends. We discuss this in lines 215-218 for the case of large angular distances and small Euclidean distances. For these reasons, we do not think that an over- or under-sampling technique is required for estimating the model parameters.

While we will consider adding more details to our description of the data set, we note that the same data set was used in many previous studies for GMMs (e.g., Chiou and Youngs, 2014), as well as for spatial correlation models (e.g., Baker and Chen, 2019). A very detailed description of the data set is available in Ancheta et al. (2013), and we do not think that repeating the information in the present manuscript adds further benefit. Instead, we will focus on the aspects that are important for the correlation models presented in the study (such as the ones mentioned in the paragraph above).

• Comment 6: In line 140, the authors claim the isometric parameters increased the predictive power. Can some details be provided on how much improvement was observed?

Response: We will precise this statement in the revised version. Please see also our preliminary reply to comment #4 from the first reviewer, where we give a more detailed explanation.

• Comment 7: Details on the MCMC NUTS process are not quite thorough such as the number of samples, thread, etc., are not clear. Typically, in Bayesian approaches, the hyperparameters need to be stabilized through sensitivity analysis. Such details are missing in the paper.

Response: In response to the reviewer’s comment, we will add a more detailed description of the employed Markov-chain Monte Carlo (MCMC) procedure. Please see also our preliminary reply to comment #9 from the first reviewer. Note that the full implementation and detailed model specifications are available in the cited online repository that supplements the manuscript (see “Assets” on NHESS), such that interested readers can reproduce the presented results.  

• Comment 8: The train-test split of the data is not clear. Can the authors please provide more details on the datasets after the split and also provide the statistical analysis using the events of test sets?

Response: Based on the reviewer’s comment, we will improve our explanation of the different comparisons that we performed. We also refer to our preliminary reply to comment #5 from the second reviewer. As explained in Section 3.1, we used the NGA-West2 data set as the training data (see also our reply to comment #5 above). We did not perform any train-test splits on the NGA-West2 data set. Instead, we used data from the 2019 Ridgecrest earthquake sequence (Rekoske et al., 2020) as test data. This data set was not used to train the models. In lines 292-293, we refer to Rekoske et al. (2020) and Liu et al. (2022), which both present thorough statistical analyses of the test data set.

• Comment 9: It’s unclear, but it seems that the authors have not considered cross-correlation between the Sas, so having independent models can have issues like those raised by the article’s co-authors?

Response: The manuscript proposes a novel correlation model that accounts for path and site effects. As an important first step, the study focuses on correlations of the same IM at different sites, while the model can be extended to include spatial cross-correlations between multiple IMs. Specifically, the linear co-regionalization technique, used by Loth and Baker (2013) for isotropic models, can also be employed for the newly proposed model EAS. We will carefully specify this in the revised manuscript. Please see also our preliminary reply to comments #2 and #3 from the first reviewer.

• Comment 10: On the same note, the authors use RotD50Sa as the IM for the spatial correlation. However, it is a bit inconsistent with the objective since spectral ordinates of the RotD50Sa spectrum do not belong to the same rotation angle and hence arguably not the same ground motion. The authors should reflect on this in terms of spatial correlation.

Response: We mainly chose RotD50 as an intensity measure (IM) for following two reasons: First, the underlying GMMs are based on this IM. Second, the same IM was used in previous spatial correlation studies (e.g., Baker and Chen, 2019), as well as in non-ergodic GMMs that quantify systematic path effects (e.g., Liu et al., 2022). While we think that the derivation of additional models for other IMs (such as RotD100) is interesting, we consider it not to be within the scope of the present study.

• Comment 11: While it doesn’t make a significant difference for small distances, the authors should check more realistic Haversine distance measures compared to the Euclidean distance in the correlation model.

Response: We thank the reviewer for this comment. As mentioned in Section 2.1, we compute the Euclidean distances between projected Cartesian coordinates. Before submission, we validated the accuracy of this approach by computing the geodesic distance between every station pair on the WGS 84 ellipsoid.  For the latter computations, we used the PyProj software library (PROJ contributors, 2023). The maximum relative difference between both approaches amounts to less than 0.02%, and is deemed as negligible for the present purpose.   

• Comment 12: Lines 235-239: Why a subset from the training dataset is used as test data? Why aren’t they selected from a separate database?

Response: The statement on lines 235-239 pertains to Section 4 where we assess the in-sample model performance. In this section, the testing data set is a subset of the training data. In Section 5, however, we compare the different models on out-of-sample data that was not used for model training. Based on the reviewer’s comment, we will improve our explanation of the different comparisons in the revised manuscript. See also our reply to comment #8 above.

References cited in this preliminary reply

Ancheta TD., Darragh RB., Stewart JP. et al. (2014): PEER NGA-West2 Database. PEER Report 2013/03, Pacific Earthquake Engineering Research Center.

Baker JW., Chen Y. (2019): Ground motion spatial correlation fitting methods and estimation uncertainty, Earthquake Engineering and Structural Dynamics.

Chiou BSJ., Youngs RR. (2014): Update of the Chiou and Youngs NGA Model for the Average Horizontal Component of Peak Ground Motion and Response Spectra. Earthquake Spectra.

Draper NR., Smith H. (1998): Applied Regression Analysis (3rd ed.). Wiley.

Esposito S., Iervolino I. (2012): Spatial Correlation of Spectral Acceleration in European Data. Bulletin of the Seismological Society of America.

Heresi P., Miranda E. (2019): Uncertainty in intraevent spatial correlation of elastic pseudo-acceleration spectral ordinates. Bulletin of Earthquake Engineering.

LeCun Y., Cortes C., Burges CJ. (2010): MNIST handwritten digit database. http://yann.lecun.com/exdb/mnist

Liu C., Macedo J., Kottke AR. (2022) Evaluating the performance of non-ergodic ground motion models in the Ridgecrest area. Bulletin of Earthquake Engineering.

Loth C, Baker JW. (2013) A spatial cross-correlation model of spectral accelerations at multiple periods. Earthquake Engineering and Structural Dynamics.

PROJ contributors. (2023): PROJ coordinate transformation software library. Open Source Geospatial Foundation. URL https://proj.org/.

Rekoske JM., Thompson EM., Moschetti MP. et al. (2020): The 2019 Ridgecrest, California, Earthquake Sequence Ground Motions: Processed Records and Derived Intensity Metrics. Seismological Research Letters.