We thank the reviewer for the comments. We provide our point-by-point response below:

• Point 1: Some of the reasoning given for the exclusion of three faults of the Central Apennines is unconvincing and speculative, particularly the argument made about the minor contribution of these faults to the regional seismic hazard (lines 143-144). If the purpose of the study is to examine the influence of fault network geometry, the slow slip rates should not be a condition for exclusion. First, the impact on the hazard will depend on the target probability of exceedance (i.e., on the purpose of the hazard assessment). Second, their contribution to the hazard might not be straightforward, especially in a complex fault network context. That is, the activity of these faults might still impact earthquake cycles of nearby faults independently of their low activity rates, thereby affecting the rates of earthquakes and behavior within the system. Consider adapting the tone of this justification.
  
  We have modified the Methods section to focus only on the numerical reasons behind the exclusion of these faults (L140-L144):
  
  “In preliminary trials, we identified that the inclusion of three small, slowly slipping active faults in the Central Apennines (Parasano-Pescina, Tremonti and San Sebastiano faults) led to unrealistically high stress concentrations in the Fucino fault, which caused numerical instability and forcing the simulation to stop prematurely. For these computational reasons, we exclude the three active faults from the Central Apennines model. The potential implications of this exclusion are addressed in the Discussion section”
  
  We have additionally included in the Discussion section a paragraph discussing in a qualitative way the implications of including these faults (L518-L525):
  
  “Although not included in the simulations, the Parasano-Pescina, Tremonti and San Sebastiano may influence the seismic sequences of the Central Apennines fault network. Due to their limited area, the Parasano-Pescina and Tremonti faults would likely produce full-rupture events, similar to the Irpinia Antithetic fault in the Southern Apennines (Figure 4a-c). The inclusion of these tree faults would increase the number of across-strike interacting fault segments, thus promoting more stress heterogeneities in neighbouring faults such as Fucino, Scanno and Trasacco, and potentially leading to more complex seismic cycles or multi-fault ruptures. Given the limitations of the modelling framework, it is currently not feasible to investigate these interactions within a network-wide approach.”
•  Point 2: The Alburni fault depth is limited at 3km to avoid intersection with the Vallo di Diano fault. This is a very short depth to effectively generate earthquakes. In fact, in your study this fault does not generate earthquakes.  Have you considered resolving this structure better at depth? For instance, propagating it up to the intersection with the Vallo di Diano fault. Otherwise, I do not see the point of including this fault.
  
  In previous trials, we attempted to model the Alburni fault up to its intersection with the Vallo di Diano fault, where the fault depth remains shallower than ~4 km. This set up led to unrealistically high stress concentrations along the Vallo di Diano fault, which caused the numerical solution to become unstable, forcing the simulation to stop prematurely. To ensure numerical stability, we shortened the fault without necessarily excluding it to determine whether the truncated segment would participate in multifault ruptures, which was not observed in the simulations. As mentioned by the reviewer, one alternative would be to exclude this fault, but this requires rerunning the full models, which is computationally demanding given the two-month runtime of the simulation.
• Point 3: Line 247: “fault sin” should be “faults in”
  
  The typo was fixed.
•  Point 4: For the analyses in figure 4, I suggest to compare the MFDs obtained with the simulator with analytical MFDs based on the seismic moment rate balancing. Given that your models do not generate multi-fault ruptures, such a comparison would be fair and would be useful to clearly determine what are the differences/improvements of your physics-based approach against more simpler, but widely used, approaches in seismic hazard. Physics-based models are expensive and not always easy to implement, so this exercise might help you further justify why we need them in favor of less expensive approaches that might provide similar occurrence statistics.
  
  Following the reviewer suggestion we have compared the earthquakes rates based on the seismic moment rate balancing to the ones obtained from the SEAS simulations (Appendix H, Fig. 4c,f). We address this comparison in the discussion section 5.4 (“Implications of fault network geometry on seismic hazard assessment”) (L552-L570):
  
  “Physics-based earthquake simulators based on SEAS models provide an alternative framework to constrain earthquake rates and magnitude-frequency distributions for seismic hazard assessment (SHA, e.g., Milner et al., 2021; Shaw et al., 2018, 2022). To explore this potential, we compare magnitude-frequency distributions derived from our SEAS simulations with those obtained from an approach (Pace et al., 2016) based on seismic moment conservation and a truncated Gutenberg-Richter distribution (“fault-based method”, Appendix H). 
  
  For both fault networks, the two approaches yield systematically different magnitude-frequency distributions (Fig. 4c,f). The fault-based method predicts lower recurrence rates than the SEAS simulations for Mw<~6.3, but higher rates at larger magnitudes. One likely explanation is that, unlike the fault-based method which assumes that the entire long-term slip rate is released seismically, the SEAS models allow part of the moment budget to be accommodated aseismically through creep or slow slip transients (e.g. Fig. 5c), thereby reducing the seismic budget for small earthquakes (Rodriguez Piceda et al., 2025b). Additionally, the tested frictional properties in our SEAS models (Table A1) yield a W_s/L_∞   ratio (Eq. A5) that may favour more characteristic rupture behaviour, limiting the occurrence of smaller events relative to Gutenberg-Richter assumptions (Barbot, 2019; Cattania, 2019). 
  
  At the high-magnitude end, SEAS simulations produce larger magnitudes and higher rate of events than the fault-based method. This difference is likely due to the fault-based method assuming stationary recurrence and a priori prescribed maximum magnitude (Pace et al., 2016), while in the SEAS models rupture sizes emerge dynamically from stress evolution over a finite catalog duration. Therefore, this comparison, although limited, shows how physics-based models can provide independent constraints both on the shape of the magnitude-frequency distribution and the maximum magnitude for SHA.”
•  Point 5: In the paragraph starting in line 317, you mention that fault networks with fewer partially overlapping across-strike nearby faults tend to nucleate earthquakes near the fault tips. In the appendix, you also mention that QDYN uses the back-slip approach to apply loading on faults. The back-slip approach is known to generate stress singularities near the fault tips resulting in unrealistic earthquake nucleations in these fault regions, especially when slip rates are not tapered (as your case). I was wondering to what extent these fault-tip nucleations in more isolated fault systems like the Southern Apennines can be linked to effects of the back-slip loading? Does QDYN or your post-simulation analysis consider any internal correction to mitigate this effect?
  
  Our models include a velocity-strengthening region around fault edges to prevent the infinite stress rates that would otherwise arise there from the backslip method (Rodriguez Piceda et al., 2025a). Despite the smoothening effect of these velocity-strengthening regions, the resulting stressing rates are still large near their edges, which promotes earthquake nucleation there.
  
  We added in the Appendix A section a line explaining the model characteristics that prevents the unrealistic stressing rates from the backslip method (L560-L562):
  
  “Additionally, we set a 1 km transition zone with velocity-strengthening friction properties along the edges of the velocity-weakening regions to prevent infinite stress rates at the fault edges that could arise from the backslip method (Rodriguez Piceda et al., 2025a).”
  
  We included in the Results section a description attributing the nucleation location towards the fault tips to the high stressing rates due to the backslip loading method (L322-L327):
  
  “Faults such as Volturara, which have few nearby, partially overlapping, across-strike faults, tend to show coseismic stress changes that remains broadly similar through successive earthquake cycles. Seismicity comprises mainly full ruptures that nucleate near fault tips (Fig. 7a, c). These nucleation locations correspond to regions of elevated stressing rates produced by the backslip loading method, which remain highest near fault tips despite the inclusion of velocity-strengthening buffer regions (Appendix A).
•  Point 6: In figure 9, the mean recurrence time dashed line is grey in the figure but black in the legend. Consider homogenization.
  
  The colour black in the legend was intended to convene a general colour for the three types of recurrence time (for all ruptures, full and partial ruptures). To avoid future confusion, we modified the legend to include the line stroke of the three types of recurrence times (Figure 9).
• Point 7: Across-strike interaction index is referred to as both AI and GI in the manuscript (e.g., line 432). Consider homogenization.
  
  All the references to GI were updated as AI (across-interaction index)
• Point 8: Line 435 should say Central Apennines? Spearman correlation indices in Table F1 indicate poorer correlation for the Central Apennines than the Southern.
  
  We changed L442 to “Central Apennines”
• Point 9: I am curious to know a bit more on why the simulator does not produce multi-fault ruptures in the system. There are other quasi-dynamic simulators (RSQSim or MCQsim) that are demonstrating to produce such ruptures in across-strike fault systems, so I am not sure how the quasi-dynamic approximation in QDYN is inhibiting this behavior. Does it have to do with the specific friction/stressing parametrization of your models? Have you verified that other friction parameters also inhibit multi-fault ruptures? I think the paper should discuss more details on this topic
  
  Quasi-dynamic simulators produce fewer multi-fault ruptures than fully-dynamic models because they do not account for dynamic triggering by seismic waves, a mechanism that facilitates rupture jumping at larger fault separations than static triggering. Moreover, our models do not account for dynamic weakening mechanisms (e.g. flash heating), which are also known to promote rupture across multiple fault segments (Thomas et al., 2014).
  
  To our knowledge no studies using MCQSIM or QDYN have reported multi-fault rupture events under standard quasi-dynamic assumptions. However, as pointed out by the reviewer, multi-fault ruptures are commonly reported in RSQSim simulations (Herrero-Barbero et al., 2021; Milner et al., 2021; Shaw et al., 2025). We identify multiple reasons of why this could be the case and included them in the Discussion (L471-L489):
  
  “While our models reproduce the magnitude range of historical single-fault events, they do not produce multi-fault ruptures that have been documented in the region (Benedetti et al., 1998; Cello et al., 2003; Galli et al., 2006). This likely stems from the quasidynamic approximation, which does not simulate dynamic triggering by seismic waves and limits stress drops, slip velocities and rupture jumping compared to fully-dynamic simulations (Thomas et al., 2014). However, studies using other quasidynamic simulators (e.g. RSQSim, Dieterich and Richards-Dinger, 2010) have reported multi-fault rupture scenarios (Herrero-Barbero et al., 2021; Milner et al., 2021; Shaw et al., 2025). In these models, faults are more closely spaced and subject to higher effective normal stresses than our QDYN simulations (~100 MPa), which increases rupture energy favouring multi-fault ruptures. In addition, model approximations of RSQSim, such as discrete state transitions and prescribed slip rates during the sliding stage might also promote larger ruptures. In comparison, our simulations resolve continuous rate-and-state friction, include velocity strengthening barriers at fault tips and prescribe a lower effective normal stress, all of which would favour rupture arrest at the fault boundaries. This limitation likely leads to an underestimation of the largest seismic events and biases the upper tail of the magnitude-frequency distribution toward single-fault ruptures (Fig. 4). While non-periodic multi-fault ruptures would likely increase the recurrence variability and cycle complexity on individual faults, it is not expected to change the first-order role of fault network geometry identified here. Multi-fault ruptures could be promoted by stronger velocity weakening asperities or higher normal stresses, but such parameter choices would increase the computational cost, limiting their feasibility in these regional-scale simulations. A second alternative to address this would be using simulated cycles as initial conditions for fully dynamic rupture models (Galvez et al., 2020).”
• Point 10: Concerning the preferred nucleation locations discussed in section 5.4, I suggest discussing the role of including fault-specific geometric complexities in future models as well. Complex, non-planar, fault geometries (e.g., variable dip, fault roughness) would most likely impact earthquake nucleation and rupture behavior significantly in the studied fault system networks.

The role of fault geometric complexities (including fault roughness, variable fault geometry) on the rupture growth and nucleation of earthquakes was mentioned in the Method section (L180-L185):

“The set up does not contain further complexities such as variable fault geometry, heterogeneous slip-rate distribution along the fault strike, and or complex heterogeneous distribution of frictional properties. Although these elements may generate heterogeneous stress concentrations that might act as barriers to rupture propagation or as regions of earthquake nucleation (Delogkos et al., 2023; Hillers et al., 2007; Luo and Ampuero, 2018; Mildon et al., 2019; Rodriguez Piceda et al., 2025b), we chose to omit these complexities to primarily focus on the impact of fault network geometry and fault stress interactions on rupture dynamics.”

References

Barbot, S.: Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault, Tectonophysics, 768, 228171, https://doi.org/10.1016/j.tecto.2019.228171, 2019.

Cattania, C.: Complex Earthquake Sequences On Simple Faults, Geophysical Research Letters, 46, 10384–10393, https://doi.org/10.1029/2019GL083628, 2019.

Chan, C.-H., Sørensen, M. B., Stromeyer, D., Grünthal, G., Heidbach, O., Hakimhashemi, A., and Catalli, F.: Forecasting Italian seismicity through a spatio-temporal physical model: importance of considering time-dependency and reliability of the forecast, Annals of Geophysics, 53, 129–140, https://doi.org/10.4401/ag-4761, 2010.

Cinque, A., Ascione, A., and Caiazzo, C.: Distribuzione spazio-temporale e caratterizzazione della fagliazione quaternaria in Appennino meridionale., ITA, 2000.

Delogkos, E., Howell, A., Seebeck, H., Shaw, B. E., Nicol, A., Mika Liao, Y.-W., and Walsh, J. J.: Impact of Variable Fault Geometries and Slip Rates on Earthquake Catalogs From Physics-Based Simulations of a Normal Fault, Journal of Geophysical Research: Solid Earth, 128, e2023JB026746, https://doi.org/10.1029/2023JB026746, 2023.

Faure Walker, J. P.: Mechanics of continental extension from Quaternary strain fields in the Italian Apennines, Doctoral, UCL (University College London), 405 pp., 2010.

Faure Walker, J. P., Visini, F., Roberts, G., Galasso, C., McCaffrey, K., and Mildon, Z.: Variable Fault Geometry Suggests Detailed Fault‐Slip‐Rate Profiles and Geometries Are Needed for Fault‐Based Probabilistic Seismic Hazard Assessment (PSHA), Bulletin of the Seismological Society of America, 109, 110–123, https://doi.org/10.1785/0120180137, 2019.

Faure Walker, J. P., Boncio, P., Pace, B., Roberts, G., Benedetti, L., Scotti, O., Visini, F., and Peruzza, L.: Fault2SHA Central Apennines database and structuring active fault data for seismic hazard assessment, Sci Data, 8, 87, https://doi.org/10.1038/s41597-021-00868-0, 2021.

Galli, P. A. C., Peronace, E., Quadrio, B., and Esposito, G.: Earthquake fingerprints along fault scarps: A case study of the Irpinia 1980 earthquake fault (southern Apennines), Geomorphology, 206, 97–106, https://doi.org/10.1016/j.geomorph.2013.09.023, 2014.

Herrero-Barbero, P., Álvarez-Gómez, J. A., Williams, C., Villamor, P., Insua-Arévalo, J. M., Alonso-Henar, J., and Martínez-Díaz, J. J.: Physics-Based Earthquake Simulations in Slow-Moving Faults: A Case Study From the Eastern Betic Shear Zone (SE Iberian Peninsula), Journal of Geophysical Research: Solid Earth, 126, e2020JB021133, https://doi.org/10.1029/2020JB021133, 2021.

Hillers, G., Mai, P. M., Ben-Zion, Y., and Ampuero, J.-P.: Statistical properties of seismicity of fault zones at different evolutionary stages, Geophysical Journal International, 169, 515–533, https://doi.org/10.1111/j.1365-246X.2006.03275.x, 2007.

Iacoletti, S., Cremen, G., and Galasso, C.: Advancements in multi-rupture time-dependent seismic hazard modeling, including fault interaction, Earth-Science Reviews, 220, 103650, https://doi.org/10.1016/j.earscirev.2021.103650, 2021.

Luo, Y. and Ampuero, J.-P.: Stability of faults with heterogeneous friction properties and effective normal stress, Tectonophysics, 733, 257–272, https://doi.org/10.1016/j.tecto.2017.11.006, 2018.

Mignan, A., Danciu, L., and Giardini, D.: Considering large earthquake clustering in seismic risk analysis, Nat Hazards, 91, 149–172, https://doi.org/10.1007/s11069-016-2549-9, 2018.

Mildon, Z. K., Roberts, G. P., Faure Walker, J. P., and Toda, S.: Coulomb pre-stress and fault bends are ignored yet vital factors for earthquake triggering, Nat. Commun., 1–9, https://doi.org/10.31223/osf.io/pt829, 2019.

Milner, K. R., Shaw, B. E., Goulet, C. A., Richards‐Dinger, K. B., Callaghan, S., Jordan, T. H., Dieterich, J. H., and Field, E. H.: Toward Physics‐Based Nonergodic PSHA: A Prototype Fully Deterministic Seismic Hazard Model for Southern California, Bulletin of the Seismological Society of America, 111, 898–915, https://doi.org/10.1785/0120200216, 2021.

Morewood, N. C. and Roberts, G. P.: The geometry, kinematics and rates of deformation within an en échelon normal fault segment boundary, central Italy, Journal of Structural Geology, 22, 1027–1047, https://doi.org/10.1016/S0191-8141(00)00030-4, 2000.

Pace, B., Visini, F., and Peruzza, L.: FiSH : MATLAB Tools to Turn Fault Data into Seismic‐Hazard Models, Seismological Research Letters, 87, 374–386, https://doi.org/10.1785/0220150189, 2016.

Papanikolaou, I. D., Roberts, G. P., and Michetti, A. M.: Fault scarps and deformation rates in Lazio–Abruzzo, Central Italy: Comparison between geological fault slip-rate and GPS data, Tectonophysics, 408, 147–176, https://doi.org/10.1016/j.tecto.2005.05.043, 2005.

Roberts, G. P. and Michetti, A. M.: Spatial and temporal variations in growth rates along active normal fault systems: an example from The Lazio–Abruzzo Apennines, central Italy, Journal of Structural Geology, 26, 339–376, https://doi.org/10.1016/S0191-8141(03)00103-2, 2004.

Rodriguez Piceda, C., Mildon, Z. K., van den Ende, M., Ampuero, J.-P., and Andrews, B. J.: Normal Fault Interactions in Seismic Cycles and the Impact of Fault Network Geometry, Journal of Geophysical Research: Solid Earth, 130, e2024JB030382, https://doi.org/10.1029/2024JB030382, 2025a.

Rodriguez Piceda, C., Mildon, Z., Andrews, B. J., Visini, F., Ampuero, J. P., and van den Ende, M.: Spatially heterogenous Holocene slip rates drive seismic sequence variability on normal faults, Seismica, 4, https://doi.org/10.26443/seismica.v4i2.1682, 2025b.

Shaw, B. E., Milner, K. R., Field, E. H., Richards-Dinger, K., Gilchrist, J. J., Dieterich, J. H., and Jordan, T. H.: A physics-based earthquake simulator replicates seismic hazard statistics across California, Science Advances, 4, eaau0688, https://doi.org/10.1126/sciadv.aau0688, 2018.

Shaw, B. E., Fry, B., Nicol, A., Howell, A., and Gerstenberger, M.: An Earthquake Simulator for New Zealand, Bulletin of the Seismological Society of America, 112, 763–778, https://doi.org/10.1785/0120210087, 2022.

Shaw, B. E., Milner, K. R., and Goulet, C. A.: Deterministic Physics‐Based Earthquake Sequence Simulators Match Empirical Ground‐Motion Models and Enable Extrapolation to Data‐Poor Regimes: Application to Multifault Multimechanism Ruptures, Seismological Research Letters, 96, 2431–2444, https://doi.org/10.1785/0220240141, 2025.

Spagnuolo, E., Herrero, A., and Cultrera, G.: The effect of directivity in a PSHA framework, Geophysical Journal International, 191, 616–626, https://doi.org/10.1111/j.1365-246X.2012.05630.x, 2012.

Stein, R. S., Barka, A. A., and Dieterich, J. H.: Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering, Geophysical Journal International, 128, 594–604, https://doi.org/10.1111/j.1365-246X.1997.tb05321.x, 1997.

Thomas, M. Y., Lapusta, N., Noda, H., and Avouac, J.-P.: Quasi-dynamic versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced coseismic weakening, Journal of Geophysical Research: Solid Earth, 119, 1986–2004, https://doi.org/10.1002/2013JB010615, 2014.

Thompson, E. M. and Worden, C. B.: Estimating Rupture Distances without a Rupture, Bulletin of the Seismological Society of America, 108, 371–379, https://doi.org/10.1785/0120170174, 2017.

Toda, S., Stein, R. S., Reasenberg, P. A., Dieterich, J. H., and Yoshida, A.: Stress transferred by the 1995 M = 6.9 Kobe, Japan, shock: Effect on aftershocks and future earthquake probabilities, Journal of Geophysical Research: Solid Earth, 103, 24543–24565, https://doi.org/10.1029/98JB00765, 1998.

Valentini, A., Visini, F., and Pace, B.: Integrating faults and past earthquakes into a probabilistic seismic hazard model for peninsular Italy, Natural Hazards and Earth System Sciences, 17, 2017–2039, https://doi.org/10.5194/nhess-17-2017-2017, 2017.

Yin, Y.: Seismicity and fault interaction through observation and simulation, ETH Zurich, https://doi.org/10.3929/ETHZ-B-000578638, 2022.

We thank the reviewer for the comments. We provide our point-by-point response below:

•    Point 1: First, regarding the scaling factor of ×100 adopted for numerical reasons: although the explanation is present, the paragraph is quite dense, and not all readers familiar with QDYN modelling will immediately recognise the implications of this choice. It would be helpful if the authors could elaborate on why this specific scaling factor was chosen rather than another value, and, if available, include a brief robustness check showing results obtained with a different scaling factor (e.g. ×50 or ×150). This would provide readers with a clearer sense of the sensitivity of the results to this numerical adaptation.

We modify the paragraph in the Methods section to justify our choice of scaling factor (L166-L179):

“We assign a constant loading rate value to each fault, corresponding to the maximum Holocene slip-rate (Figure 2; Cinque et al., 2000; Faure Walker, 2010; Faure Walker et al., 2019, 2021; Galli et al., 2014; Morewood and Roberts, 2000; Papanikolaou et al., 2005; Roberts and Michetti, 2004; Valentini et al., 2017). Preliminary trials conducted in this study and from previous work (e.g., Yin, 2022), show that in the (implicit) LSODA solver used here, slip rates below 1 mm/yr can cause velocity values to drop below numerical precision, especially when neighbouring faults also exhibit low slip rates. We suspect that this is due to the difference between the minimum and maximum slip rates present in the system, which span many orders of magnitude, causing precision underflow when jointly solving the system of equations. Therefore, we scale up the loading rate by a factor of 100 which represents the minimum value required to ensure numerical stability across the full fault network, while remaining as close as possible to the geologically inferred loading rates. Based on these trials and on previous studies (Yin, 2022) we find that this scaling has minimal impact on seismicity statistics, once we correct simulations outputs (such as recurrence times) for this upscaling factor (see Appendix B). All time scales reported hereafter account for this correction factor.”

•    Point 2: A second point concerns the exclusion of the three smaller faults (Parasano–Pescina, Tremonti, and San Sebastiano) from the Central Apennines network. While I understand the need to reduce intersections and avoid numerical instabilities, the justification feels somewhat incomplete given the stated aim of analysing system-scale behaviour. Even relatively short or low-slip-rate segments may influence the connectivity of the network and the topology of rupture paths. Previous studies have shown that such linking structures can significantly affect recurrence intervals, segmentation patterns, and fault-to-fault interactions. For these reasons, it would be beneficial for the authors to provide a clearer rationale for their exclusion and, if any preliminary sensitivity tests were performed (even qualitatively), to briefly comment on whether the presence or absence of these segments modifies the system behaviour. A short assessment of this point would greatly strengthen confidence in the selected network geometry.

We refer to answer to Point 1 of Reviewer#1 for a detailed analysis of this modelling choice and the corresponding modifications made to the manuscript.

•    Point 3: A related and important limitation is that the current numerical framework does not allow ruptures to propagate across fault intersections (namely multi-fault ruptures). In the tectonic context of the Central Apennines, multi-fault ruptures for example have been proposed to explain several aspects of the paleoseismic record. Recent fault-based PSHA studies have also shown that relaxing segmentation rules improves consistency with observed rupture histories in the region. A more explicit acknowledgement of how this limitation affects the interpretation of the results would be valuable. In particular, discussing whether this might influence the reported patterns of recurrence variability and complexity would help readers understand the scope of the conclusions. Even a brief paragraph in the Discussion would make this aspect clearer.

This comment has substantial overlap with the point 9 of Reviewer #1. Our response to Reviewer #1 was quite lengthy, and so we will not repeat it in full here, but we have also included in the paragraph cited in that answer our qualitative interpretation of how multifault ruptures would affect the results (L484-L486):

“In summary, while non-periodic multi-fault ruptures would likely increase the recurrence variability and cycle complexity, it is not expected to change the first-order role of fault network geometry identified here.”  

•    Point 4: Finally, in the last part of the manuscript, the implications for fault-based seismic hazard are interesting but could be better articulated. While the scientific insights into geometric control on network behaviour are compelling, the practical consequences for hazard modelling remain somewhat unclear. Probabilistic fault-based models typically require magnitudes, recurrence intervals, rupture dimensions, and structured epistemic uncertainties. It is not yet evident how the results of the present SEAS models can be translated into these quantities. For example, can the synthetic recurrence behaviour be used to inform logic-tree branches? Do the simulated ruptures provide magnitude–frequency constraints? Or should the results be interpreted mainly as conceptual insights into fault-network dynamics? Clarifying what aspects of the modelling are directly actionable for PSHA, and what aspects remain at the conceptual or physics-modelling level, would significantly improve the applied relevance of the work.

Following this point and Point 4 of Reviewer 1, we have restructured the Discussion section 5.4 (Implications of fault network geometry on seismic hazard assessment) which includes the potential use of physics-based earthquake simulators to constrain earthquake rates and magnitude frequency distributions (see response to Point 4 of Reviewer 1). We have additionally rephrased the discussion section to better convey the possible applications to seismic hazard assessments (L571-L587):

“Beyond these general applications, our results also show how fault-network geometry modulates earthquake recurrence and magnitude variability challenging the common practice of applying uniform recurrence parameters (e.g. a single mean recurrence time or coefficient of variation) across an entire fault system (e.g., Nishenko and Buland, 1987, Ellsworth et al., 1999; Matthews et al., 2002). In networks with numerous across-strike faults and high long-term slip rates, recurrence variability can vary greatly between faults (e.g., Central Apennines, Fig. 9b), suggesting that hazard models should allow for broader epistemic uncertainty in recurrence and magnitude distributions. Therefore, probabilistic SHA could further benefit from integrating network-derived metrics, such as the across-strike interaction index, as quantitative guides for weighting logic-tree branches. While such metrics do not directly prescribe recurrence or magnitude parameters, they may provide physically grounded constraints to better reflect the complexity of fault interactions. Physics-based simulators are well suited to quantify the recurrence and magnitude variability (Milner et al., 2021; Shaw et al., 2018, 2022), and to inform the weighting of alternative fault-based SHA models. In addition, time-dependent SHA approaches, such as those based on Coulomb stress transfer histories (Chan et al., 2010; Iacoletti et al., 2021; Mignan et al., 2018; Stein et al., 1997; Toda et al., 1998), should also account for spatial variability in stress changes. Finally, our simulations indicate that overlap zones between faults can act as preferred nucleation sites and constrain rupture extent (Fig. 7), highlighting their importance for assessing rupture scenarios and directivity effects (Spagnuolo et al., 2012; Thompson and Worden, 2017).”