We would like to thank the anonymous referee for the thoughtful review of our manuscript and for giving these constructive comments and suggestions, which substantially helped us improve the quality of the paper. All the points that were raised have been adopted in the revised manuscript. We believe the new version of the manuscript has been significantly improved. Below is a point-by-point answer to the comments and suggestions raised by the reviewer.

Major comments：

1. Aftershock selection and illustration.

Number and locations of early aftershocks are critical in this algorithm. I am not sure how the early (within 2 h) aftershocks could accurately reveal the rupture pattern of earthquakes. I do recommend the authors make plots of the early aftershocks for the events shown in this manuscript. By doing this, readers can easily judge how the early aftershocks reflected the source dimension. Statistical analysis might be required to demonstrate this question.

Reply: In the new version of the manuscript, we have added temporal and spatial distribution plots of the early aftershock sequences of the Wenchuan Mw7.9 and Kaikōura Mw7.8 earthquakes in Section 2.2.2, as well as interpreted the insets. The added content is as follows:

“The early aftershocks of these two earthquakes were mainly distributed along the direction of the surface rupture zone and within a certain range on both sides of the surface rupture zone, based on the spatial and temporal distribution of the aftershock sequences (Fig. 3). The Wenchuan earthquake's early aftershocks mostly occurred within 300 kilometers of the epicenter and were concentrated along the main rupture direction. And the early aftershocks of the Kaikōura earthquake were distributed within 200 kilometers of the epicenter and were relatively dispersed along the main rupture strike, which was primarily caused by the earthquake's complex fault system (Wallace et al., 2018). Early aftershocks in both cases exhibit a pattern of rupture in a single direction. Based on this, we believe that the fitting results in figure 2 can roughly depict the length and direction information of the earthquake rupture.”

Figure 3: Spatial and temporal distribution of early aftershock sequences after the mainshock. (a) and (c) depict the spatial distribution of aftershocks of the 2008 Wenchuan Mw 7.9 earthquake and the 2016 Kaikōura Mw 7.8 earthquake, respectively. (b) and (d) depict the temporal distribution of aftershocks of the Wenchuan and Kaikōura earthquakes, respectively. The s in Log(s) denotes the time in seconds between the aftershock and the mainshock.

2. The authors have proven that the accuracy of the estimated intensity map by this method by comparing it with other results. That is good. We can estimate the damage levels in space. But the time efficiency is less discussed or demonstrated. How fast you could deliver this result? And comparing it with other approaches would greatly enhance the importance of this work.

Reply: Some discussion of time efficiency and comparisons with other methods are added to Section 4.2, and the corresponding references were added. The added contents are as follows:

“The primary goal of developing this method is to provide information services to response workers during the black box period of an earthquake emergency. Lessons learned from the calculations of all the cases in this study and from actual earthquakes emergency work (Zhao et al., 2022b; Zhao et al., 2023). If seismic stations in the seismogenic region are as sparse and uneven as those in western China, once the earthquake is determined to be suitable for use with AL-SM99, a reliable seismic intensity assessment map can be produced within 1-1.5 h of the mainshock. The majority of the time required to produce the results is spent acquiring aftershock data, while processing the aftershock data takes only a few seconds and the calculation and output of the maps takes about five minutes.”

“Chen et al (2022a) proposed a rapid assessment method that can generate intensity assessment maps within 30 minutes. The spatial location of the rupture trajectories obtained from the inversion of rupture processes for earthquakes of small magnitudes may be less satisfactory (Honda et al., 2011; Yao et al., 2019). In small-magnitude earthquakes, however, seismic intensities assessed using aftershock data may be more accurate (Kang et al, 2023). The AL-SM99-fitted curves of the spatial distribution of aftershocks can be used as a cross-reference for the correction of the above inversion results, speeding up the operation using both methods. For global earthquakes, when the magnitude reaches the trigger threshold of the ShakeMap system, it will generate the first version of the assessment results through the original solution built into the system within minutes after the earthquake, and will be continuously updated as data is aggregated and accumulated (Worden et al, 2020). It inspires us to combine AL-SM99 with aftershock monitoring to dynamically present intensity assessment results, since for earthquakes with small rupture scales, relying on the epicenter coordinate or a small number of aftershocks can provide very useful shaking distribution estimates.”

Reference:

“Honda, R., Yukutake, Y., Ito, H., Harada, M., Aketagawa, T., Yoshida, A., Sakai, S. I., Nakagawa, S., Hirata, N., and Obara, K.: A complex rupture image of the 2011 off the pacific coast of Tohoku earthquake revealed by the meso-net, Earth Planet Sp., 63(7), 583–588, https://doi.org/10.5047/eps.2011.05.034, 2011.”

“Yao, Q., Wang, D., Fang, L. H., and Mori, J.: Rapid estimation of magnitudes of large damaging earthquakes in and around Japan using dense seismic stations in China, Bull. Seismol. Soc. Am. 109, 2545–2555. https://doi.org/10.1785/0120190107, 2019.”

3. To better validate the accuracy of the source dimension estimated from the early aftershocks, the authors could compare your results with source ruptures, at least for large earthquakes. I believe there are many cases that can be utilized for such comparison.

Reply: Nine earthquakes with Mw≥7.0 were used as an example in Section 3.2.Our results were compared to surface rupture lengths calculated using an empirical formula for wells and those documented in the literature, and the average linear direction of surface rupture was calculated using ArcGIS software (Table 2). Subject to the conditions of use, the results of our method's fitting can provide reasonably accurate information on the length and direction of surface rupture. In addition, we include a comparison with the back-projection results of Chen et al (2022a) in this section, which supports our conclusions. However, since Lowess is essentially a nonparametric regression method that ignores the complex physical relationships contained in the aftershock sequence, we believe its results cannot fully replace those obtained through physical means (e.g., back-projection techniques). But the different methods can be cross-referenced to make further corrections to the results. The added content is as follows:

“We gathered the source rupture data from Chen et al (2022a) using the back-projection technique. And using the same technique, a set of results reflecting the surface rupture of the 2016 Kaikura Mw7.8 earthquake was calculated, using waveform data from high sensitivity seismograph network in Japan. Both the Lowess and back-projection results show rupture directions similar to those indicated by the long axis of the isoseismic line in the area with intensity VIII of the Wenchuan earthquake, but the former estimates a longer rupture length (Fig.11(a)). Furthermore, the back-projection results reveal more details about the rupture. For example, the back-projection results point to a possible fracture near the IX-degree intensity anomaly in the long-axis direction. This method has also demonstrated benefits in determining the intensity anomaly area in the application of the 2022 Maduo Mw7.3 earthquake (Chen et al, 2022b). As a nonparametric method, the points fitted by Lowess are clearly distributed along a curve. However, when the fault system in the seismogenic region is complex, the dominant orientation of the rupture traced using the back-projection method may be problematic(fig.11(b)). A clear guide to array data selection may be required when using the back-projection method, and we recognize that the results of array data calculations will be more accurate if the appropriate region is chosen. We believe that the two methods can be cross-referenced in their application to obtain more accurate intensity assessment results.

Figure 11：Comparison of surface rupture results obtained using the lowess and inverse projection methods for the (a) 2008 Wenchuan Mw 7.9 and (b) 2016 Kaikōura Mw 7.8 earthquake.”

4. Comparison of your results with Chen et al. (2022a, b) that you already cited in this work is also beneficial.

Reply: We have added a comparison with Chen et al's (2022a,b) work to both the examination of the source rupture results and the discussion of time efficiency, which adds to the richness of our manuscript.The additions are mentioned above in the responses to major comment 2 and 3.

Minor comments:

1. The English needs improvements.

Line 10, mainshocks

Line 13, of 59 M XXX~XXX earthaukes that occurred from 2000-2022

Line 21, Our study suggest that with early accessible aftershocks, we are able to rapidly

determine the rupture fault plane (s), thus have better estimae of the seismic intensities.

Line 44, of an earthquake is limited,

Line 47, after earthquakes

Line 94, We selected Mw ≥ 6.6 shallow earthquakes that occurred during 2000-2022 in this study.

...

Reply: We have checked the language errors in the manuscript and polished it.

The complete response is available in the supplementary file, which also contains the accompanying illustrations.

We are grateful to the reviewer for the insightful comments and suggestions, which substantially helped us improve the quality of the paper. Almost all the points that were raised have been adopted in the revised manuscript. Almost the details in the supplementary documents have also been adopted and revised. We believe the new version of the manuscript has been significantly improved. Detailed responses to the comments are presented below.

Major comments:

1. I believe the introductory part of the work should be greately improved before being published.

Reply: We appreciate your comments regarding the introduction of our manuscript. In accordance with the comments in the supplement, we made considerable revisions to the introduction section. Almost all of the sentences were rewritten without altering the original text's intended meaning. In the new version of the introduction, we have paid special attention to sentence structure, grammar, and the transitions between texts. Several details have been elaborated in light of the additional document's comments. The updated introduction was expanded from four to five paragraphs, making the content of each paragraph more distinct. The obvious modifications are shown below.

(1) Based on the review, we rewrote the first paragraph of the introduction and included references.

Original text: “Seismic intensity reflects the strength of ground motion and its influence. Rapid seismic intensity assessment helps in formulating an early emergency response after a destructive earthquake. The rapid and accurate output of seismic intensity assessment could notably reduce the loss of life and property in disaster areas. Therefore, it is necessary to develop methods for the faster assessment of seismic intensity and the efficient use of disaster data in the early post-earthquake period.”

After modification: “Seismic intensity reflects the strength of ground motion caused by an earthquake at certain location and its influence. Rapid and accurate assessment of seismic intensity facilitates the development of emergency measures in the aftermath of a destructive earthquake, thereby reducing the number of fatalities and property damage (Erdik et al., 2011; Poggi et al., 2021). Therefore, it is necessary to develop methods for the faster assessment of seismic intensity and the effective utilization of disaster data in the early post-earthquake period.”

References:

Erdik, M., Şeşetyan, K., Demircioğlu, M. B., Hancılar, U., and Zülfikar, C.: Rapid earthquake loss assessment after damaging earthquakes. Soil Dyn. Earthq. Eng., 31(2), 247-266. ,2011

Poggi, V., Scaini, C., Moratto, L., Peressi, G., Comelli, P., Bragato, P. L., and Parolai, S.: Rapid damage scenario assessment for earthquake emergency management. Seismol Res Lett., 92(4), 2513-2530. https://doi.org/10.1785/0220200245, 2021.

(2) Lines 41-44 of the original text were rewritten, and the paragraph was split into two paragraphs from here.

Original text: “From the perspective of data acquisition, the time from the occurrence of the earthquake to the first acquisition of disaster data from the disaster area (generally within a few hours after the mainshock) is considered as the black box period of earthquake emergency disaster service.”

After modification: “The time between the occurrence of an earthquake and the first acquisition of disaster data from the seismogenic region, usually within 2-3 hours of the mainshock, is defined as the black box period for earthquake emergency response (Nie and An, 2013).”

Reference:

Nie G, and An J.: Basic theoretical model of earthquake emergency response (in Chinese). Urban Disaster Reduct. 3:25–29. 2013.

(3) Lines 50-58 of the original text have had their content optimized.

Original text: “To expand the method of rapid seismic intensity assessment and improve its timeliness and accuracy, Chen et al. (2022a) proposed a method to predict the source rupture process by using the far-field seismic array data back-projection technique, and combining it with the ground motion prediction equation (GMPE) for rapid assessment of seismic intensity. This method was validated in the 2021 Maduo Mw 7.3 earthquake in Qinghai province and the Yangbi Mw 6.1 earthquake in Yunnan province (Chen et al., 2022b). However, accurate inversion of the source rupture process for earthquakes that occur in different regions and selection of more applicable regional GMPEs are the key points that still need to be addressed and improved in the Chen et al. method.”

After modification: “Back projection could image the fault geometry of large earthquakes at high resolution and is frequently used to trace surface rupture processes and source durations (Ishii et al., 2005; Wan et al., 2022). Using the back-projection technique and ground motion prediction equation (GMPE), Chen et al. (2022a) developed a new algorithm for quickly obtaining intensity maps of destructive earthquakes. The algorithm was validated during the emergency response phase of the 2021 Maduo Mw 7.3 earthquake and the 2021 Yangbi Mw 6.1 earthquake in China, and it was confirmed to be suitable for intensity assessment in regions with sparse observation networks (Chen et al., 2022b). However, accurate inversion of the source rupture process for earthquakes that occur in different regions and the selection of more applicable regional GMPEs are the key aspects of the algorithm that still require improvement.”

References:

Ishii, M., Shearer, P. M., Houston, H., and Vidale, J. E.: Extent, duration and speed of the 2004 Sumatra–Andaman earthquake imaged by the Hi-Net array. Nature, 435(7044), 933-936. https://doi.org/10.1038/nature03675, 2005.

Wan, Z. F., Wang, D., Zhang, J. F., Li, Q., Zhao, L. F., Cheng, Y. F., Mori, J., Chen, F., and Peng, Y. Y.: Two‐Staged Rupture of the 19 October 2020 M w 7.6 Strike‐Slip Earthquake Illuminated the Boundary of Coupling Variation in the Shumagin Islands, Alaska. Seismol. Res.Lett., 94(1), 52-65. https://doi.org/10.1785/0220220203, 2022.

(4) Lines 87-91 of the original text have been rewritten.

Original text: “We used of the interquartile range (IQR) to exclude outliers from the aftershock sequence that occurs within 2 h of the mainshock, utilise Lowess to fit the spatial distribution trend of aftershocks, combine the GMPE and seismic intensity scale to assess the seismic intensities, demonstrate the implementation process and intensity assessment results through specific earthquake cases, and finally discuss the applicability of this method.”

After modification: “We utilize the interquartile range (IQR) to exclude aftershocks with abnormal geographic coordinates in the aftershock sequence that occurred within 2 hours of the mainshock. The geographic coordinates of the processed aftershocks are then fitted with Lowess, and the results of the fitting are used in the GMPE calculation. Finally, the ground motion calculation results are converted to seismic intensity using the seismic intensity scale. The new method's implementation process and the effect of intensity assessment will be demonstrated by specific earthquake cases, and the method's applicability will be discussed.”

Furthermore, language issues present in other sections of the manuscript will be reviewed and corrected.

2. In addition, arrangement and presentation of tables and figures for chosen earthquakes needs to be enhanced.

Reply: The figures in the manuscript have been rearranged and two illustrations have been added, including a map of the spatial and temporal distribution of the early aftershock sequence and a comparison of the fitted results with the surface rupture results inverted by physical means. We are also considering adopting the supplementary comments. The revised manuscript will include a new table that tallies the number of aftershocks eliminated from each case and provides additional information about these cases. And we will attempt to improve the presentation quality of Figure 11 by utilizing a new format. The newly added illustrations can be found in the response's supplementary file.

3. I was not entirely convinced if the propsed technique was efficient but I hope re-writing the results could clear up the benefits.

Reply: In the results section, we have added a comparison between the rupture results of the physical means inversion and the Lowess fit result. The comparison also demonstrates that the seismic intensities estimated by the method are reasonable under the conditions of use.We will accept your comments for further optimization of the results section.

There are conditions for the use of this method, and the Lowess is essentially a nonparametric regression method that ignores the complex physical relationships contained in the aftershock sequence, we believe its results cannot fully replace those obtained through physical means (e.g., back-projection techniques). But the different methods can be cross-referenced to make further corrections to the results.

We put our method to the test in the 2022 Luding Mw 6.6 earthquake in China, as well as two great earthquakes (Mw 7.8, Mw 7.5) in Turkey in 2023. The results show that this method is feasible in the given conditions. Furthermore, we improved the method proposed in this study so that it could be used to assess the seismic intensity of small magnitude earthquakes as well. The effects of the intensity assessment of the three aforementioned earthquakes can be viewed in the supplementary file. Our recent publications that contain related works are listed below.

Kang, D.J., Chen, W.K., Zhao, H.Q., and Wang, D.: Rapid assessment of the September 5, 2022 Ms 6.8 Luding earthquake in Sichuan, China. Earthquake Research Advances, https://doi.org/10.1016/j.eqrea.2023.100214, 2023.

Zhao, H., Jia, Y., Chen, W., Kang, D., and Zhang, C.: Rapid mapping of seismic intensity assessment using ground motion data calculated from early aftershocks selected by GIS spatial analysis, Geomatics, Natural Hazards and Risk, 14:1, 1-21, https://doi.org/10.1080/19475705.2022.2160663, 2023.

4. Scientifically, I think a physics-based simulation would be an approperaite way of proving the point and it would not be a entirly hard task to do.

Reply: We were inspired by the 2021 Journal of Geophysical Research: Solid Earth article on the simulation of the relationship between earthquake sequences and geometrically complex faults, and the 2017 Earth, Planets, and Space article on the relationship research between mainshock ruptures and aftershock sequences based on dense seismic observations. The goal of this study is to broaden the use of early aftershock sequences and to serve as a reference for testing the accuracy of the method proposed earlier by our team (Chen et al., 2022). We will focus on using physical simulations to validate our approach in the next step of our work.

References:

Ozawa, S. and Ando, R.: Mainshock and aftershock sequence simulation in geometrically complex fault zones, J. Geophys. Res. Solid Earth, 126, e2020. https://doi.org/10.1029/2020JB020865, 2021.

Yukutake, Y. and Iio, Y.: Why do aftershocks occur? Relationship between mainshock rupture and aftershock sequence based on highly resolved hypocenter and focal mechanism distributions, Earth Planets Space, 69, 1–15. https://doi.org/10.1186/s40623-017-0650-2, 2017.

Chen, W., Wang, D., Si, H., and Zhang, C.: Rapid estimation of seismic intensities using a new algorithm that incorporates array technologies and ground‐motion prediction equations (GMPEs), Bull. Seismol. Soc. Am., 112, 1647–1661. https://doi.org/10.1785/0120210207, 2022.