Thank you for your insightful comments and suggestions. We appreciate your positive feedback and your recommendations for improving the manuscript. We have rewritten the governing equations using arrows for vectors as suggested in the revised manuscript. We have expanded the explanation of the boundary conditions and their role in the numerical simulations. Meanwhile, we have added the corresponding interpretation to the caption of Figure 3 as suggested.  The English errors have also been modified as suggested. We appreciate your time and effort in reviewing our work and providing these valuable suggestions. We believe these modifications have improved the clarity and technical accuracy of the manuscript. Thank you once again for your constructive feedback.

Thanks for giving me the chance to review the MS by Su et al. I think this MS has real merit and contribution to a certain extent. However, some concerns need to be addressed in the current version.

Reply: Thank you very much. I appreciate the thoughtful and constructive comments. I have carefully addressed each comment and revised the manuscript accordingly.

Major concerns:

1. Section Abstract. The current abstract is concise but could be more informative. I think the key quantitative results or/and specific findings should be given.

Reply: Thank you for bringing this to our attention. The magma beneath the fault ascended rapidly, reaching the lower crust within 20 million years, with a cooling rate of approximately ∼35°C/Myr. Conversely, the thickened magma took 40-50 million years to ascend to the lower crust, cooling at a rate of ∼10°C/Myr. In contrast, magma without thickening and fault would take considerably longer time to reach the lower crust. Lines 19-24.

1. Section Numerical simulation and model setup. First, authors should introduce the specific location of the cross section (Fig. 2). Second, I suggest that a single section is not enough to reflect the tectonic activity of the whole study area, and additional sections are necessary. Third, Section model setup should include a justification for the chosen parameters and boundary conditions. Explain why these choices are appropriate for the study area and the geological processes being modeled.

Reply: Thank you for bringing this to our attention. We have added an approximate location of the cross-section Fig. 2b within Figure 1. As depicted in Figure 1, the Changle-Nan'ao Fault Zone significantly influences the magma distribution along the coast, showing a similar pattern in magma planes. Hence, selecting a profile for simulation can effectively illustrate the magma evolution process along the coast. The coastal fault zone of the SCB is approximately 40-60 km wide (Cui et al., 2013). To better reflect this, we increased the fault width of models from 4 km to 45 km (Fig. 2). This adjustment slightly altered the ascent style of magma, but it did not significantly affect their uplift and cooling rates. We have also included relevant references to justify our boundary conditions and parameter selections, enhancing the realism of the model. See Lines: 144,151, 168, 186-187, 219, 285-293.

1. Section Result. Due to the absence of time and fault information in Fig. 3, I cannot accurately judge the accuracy of the interpretation of the results.

Reply: Thank you for your feedback. We have updated each figure to include the evolutionary time. In addition, we added the evolutionary results of the bottom magma. “In Model 1, after 10 Myr, the bottom magma was divided into five upwelling magma bodies: 1, 2, 3, 5, and 6. Magma body 1 corresponds to the fault location, bodies1, 2 and 5 correspond to pre-existing thickened magma bodies, body 4 corresponds to relatively shallow pre-existing magma, and body 3 is located below body 4. At this stage, bodies 1, 2, and 5 are uplifted higher than bodies 3 and 6. As evolution progresses, body 3 reaches the position of body 4 and rapidly achieves heights comparable to those of bodies 2 and 5 (Fig. 3e), whereas body 6 is significantly less uplifted. In Model 2, under compressive stress, four rising magma bodies formed at the bottom. No rising magma bodies formed at the bottom right. The formed magma bodies tilted and ascended to the left (Figs. 3f-j).” Lines 224-236.

1. Section Discussion. I suggest that authors should add some text to explicitly address potential limitations of the study or models and suggest areas for future research.

Reply: Thank you for highlighting this issue. The deep structure of the late Mesozoic SCB is poorly constrained, resulting in speculative assumptions about key parameters such as fault depth and magma thickness. The model might oversimplify the complex geological features, potentially leading to inaccurate results. Additionally, assuming a uniform crustal thickness may not capture the true variability of the crust. The geometry of the lithospheric faults in the model is simplified, and important details that could affect magmatic processes are neglected. Variations in interpretations of fault characteristics and magma properties contribute to the idealized nature of our model results, which may differ significantly from actual geological conditions. This study aims to use these idealized scenarios to illustrate the complexity and diverse interpretations of magma evolution processes along the South China coast. Moreover, given the model's two-dimensional nature, it simplifies the intricate three-dimensional processes that likely influence magmatic evolution. Future research will focus on addressing these complexities to provide a more comprehensive understanding. Lines 323-336.

Specific concerns:

1. The figure captions are simple and need more detail information, such as Fig. 3.

Reply: Thank you very much for pointing out this issue. The captions have been modified as follows:

Fig. 2 (A) Schematic cross section illustrating the subduction of the Paleo-Pacific plate. (B) Reference model geometry depicting temperature, density and viscosity variations with depth (see location in Fig. 1). The boundary conditions are the same between Model 1 and Model 2, except that Model 2 is assigned a horizontal force on the right side.

Fig. 3 Results of Models 1 and 2, illustrating magma upwelling for 10-50 myr, respectively. (a)-(e): Underplating magma rising, forming five magma bodies of varying heights in Model 1; (f)-(j): Underplating magma tilting under right-sided compression in Model 2. The crust, mantle and magma materials are modelled as phases (fluids) on a dimensionless scale, with values of 0, 0.5 and 1, respectively. The contours denote the flow distribution of mantle fluid.

Fig. 4 (a) Comparison of the observed cooling histories of the CNB magmatic plutons (data from Chen et al. 2020) with time‒temperature paths generated by models of rising magma; the green shadow represents the actual cooling age of the CNB magmatic plutons. (b) Zircon Hf isotopes and ages of coastal magmatic rocks in the SCB (data from Li et al., 2023). (c) Sketch illustrating the formation stage of the underplating magma and tectonic background during 80-110 Ma, 110-130 Ma and 130-160 Ma. The olive shadow crossing sections (a) and (b) denote the age range of 110-130 Ma..

1. Fig. 1 is not cited in the article.

Reply: Thank you very much for pointing out this issue. Fig.1 has been cited on lines 95 and 137.

Reply: Thank you very much. I appreciate the thoughtful and constructive comments. I have carefully addressed each comment and revised the manuscript accordingly.

Here are my detailed comments:

1. The settings for heat capacity and thermal conductivity are declared twice in Lines 159-161 and Lines 179-181, respectively. I believe the thermal conductivity for the mantle should be larger than that for the crust, although the authors cited Chapman, 2021.

Reply: Thank you for highlighting this issue. We have removed the redundant content from Lines 179 to 181. While heat conduction in the mantle is indeed greater than that in the crust, the difference is minimal and does not significantly affect the model results. This study focuses on the impact of initial magma position differences on magma intrusion dynamics. Consequently, we simplified the parameter variations and adopted parameters similar to those in Chapman (2021).

1. It would be better to provide more details about the initial model, such as the boundary conditions, mesh resolution, etc. The contrast between the two models should be reflected in Figure 2b.

Reply: Thank you very much for this suggestion. We have incorporated the mesh resolution details as suggested: " The resolution of the model is physically controlled but generally represented by a triangular mesh, with side lengths of 2-5 km, a minimum area of 1.8 km² and a maximum area of 3.2 km². Additionally, we have added the corresponding boundary conditions: " Model 2 has similar boundary conditions and is assigned a 2e9 Pa horizontal force on the right side, representing the compression effect induced by subduction." in Lines 202-203, 211-212. We have also added the fault condition in the model: " The model fault was assigned a low viscosity of 1e19 Pa·s to represent the CNB (Vissers et al., 1995; Columbu et al., 2015)" Lines 186-187. Furthermore, the difference between the two models has been added to Figure 2b and its caption.

The coastal fault zone of the SCB is approximately 40-60 km wide (Cui et al., 2013). To better reflect this, we increased the fault width of models from 4 km to 45 km (Fig. 2). This adjustment slightly altered the ascent style of magma, but it did not significantly affect their uplift and cooling rates.

1. The mesh has side lengths of 2-5 km. However, the square representing the previously intruded magma is only 10 km thick. The grid seems too coarse for the pre-existing magma. The authors should conduct a test on the grid size.

Reply: Thank you very much for pointing out this issue. The grid we used is approximately 2 km in length and 1 km in height for the intruded magma. After testing, we confirmed that this resolution meets the requirements for accurate simulation.

1. I only see two thick bodies located at the bottom of the model, which develop into the 2 and 5 magma intrusions. Why is there no thick body under the location of magmas 1 and 3 in the initial model? One conclusion of this manuscript is that pre-existing magma can accelerate underplating magma emplacement. The authors do not provide details on how this conclusion was reached. I assume it is through the comparison between magma 3 and magmas 2/5. However, the initial conditions are different between magma 3 and magmas 2/5.

Reply: Thank you very much for this question. The thickened magma bodies are located at positions 1, 2 and 5, with an additional magma body assigned at position 4 above position 3 (see new Fig. 2). A fault is set up at position 1, causing the magma to rise at position 1. At position 4, the magma also rose, leading to the accumulation and ascent of magma at position 3. The fault and pre-existing magma are responsible for the gathering and upflow of magma at positions 1 and 3. In Model 1, the initial thickness of the magma at position 3 is less than that at positions 2 and 5, resulting in a slower initial uplift compared to positions 2 and 5 (Figure 3a). However, in the later stages of the uplift process, the height of the magma at position 3 exceeded that of the magma at positions 2 and 5 (Figure 3e). Additionally, the velocity of magma at position 3 is greater than that of the magma at the right boundary, position 6. Therefore, the presence of pre-existing magma significantly promoted and accelerated the upwelling process. We have added the evolutionary result of these magmas in Lines 229-236.

Line 199: The presence of Model 2 before Model 1 can be confusing. Since the next paragraph describes the differences between the two models, it is recommended to place the differing content of the boundary conditions in the next paragraph.

Reply: Thank you for your suggestion. We have added the initial boundary condition to Model 1 and moved the boundary condition for Model 2 to the next paragraph as recommended. Please see Lines 204-205 and 211-212 for the updates.

1. The description of the simulation results is very important, but the article lacks detail in this area. It is suggested to add a paragraph describing in detail the evolution of magma from 1 to 5 and the different results between the two models.

Reply: Thank you very much for the suggestion. We have added a new paragraph describing the evolution results of magma bodies 1-6. In Model 1, after 10 Myr, the bottom magma was divided into five upwelling magma bodies: 1, 2, 3, 5, and 6. Magma body 1 corresponds to the fault location, bodies 1, 2 and 5 correspond to pre-existing thickened magma bodies, body 4 corresponds to relatively shallow pre-existing magma, and body 3 is located below body 4. At this stage, bodies 1, 2, and 5 are uplifted higher than bodies 3 and 6. As evolution progresses, body 3 reaches the position of body 4 and rapidly achieves heights comparable to those of bodies 2 and 5 (Fig. 3e), whereas body 6 is significantly less uplifted. In Model 2, under compressive stress, four rising magma bodies formed at the bottom. No rising magma bodies formed at the bottom right. The formed magma bodies tilted and ascended to the left (Figs. 3f-j). Lines 224-236.

1. Lines 221-223: Figure 3 only shows the evolution from 1 to 5 Myr. How do the authors know the evolution within 20 Myr? The captions for Figure 3 are too simple. The authors should add descriptions for the different colors and the meanings of the values in the color bar. It would be better to add the time of each snapshot. What do the contours and white lines in Figure 3 represent? These lines are not clearly visible.

Reply: Thank you for pointing out this issue. The actual age of evolution is 10-50 Myr. We have revised the caption of Figure 3 as suggested. The colour bar has been described, and the time of each snapshot has been added. The contours denote the flow distribution of mantle fluid. Initially, the model was divided into many modules. Although these modules were not used in the final model, their boundary lines were recorded during processing. We have now deleted these unnecessary white lines.

1. Lines 230-232: Whether or not there was magma, the cooling rate for the underplating magma looks similar in Figure 4a.

Reply: Thank you for pointing out this issue. The cooling rate of thickened magma is faster than that of magma without thickening (e.g., magmas 3 and 5 in Fig. 4). Additionally, magma at the fault cools more rapidly (e.g., magma 1 in Fig. 4). However, under compression, the cooling rate of the magma body is relatively slower. Therefore, we have revised the sentence to: "The results indicate a more rapid cooling rate for the ascending magma through the fault (∼35°C/Myr) and a slower cooling rate for underplating magma (∼10°C/Myr)." Lines 247-249.

Lines 234-235: What does this mean? Is there a fault under the pre-existing magma? I did not see it.

Reply: Thank you for bringing this to our attention. The sentence has been revised to: " Moreover, the ascent of magma to the location of pre-existing magma is also faster than that without pre-existing magma (e.g., magma 3 and magma 6 in Fig. 3a-e)."

The captions for Figures 4a and 4b are in the wrong order. The captions need more details. What is the meaning of the green shadow in Figure 4a? What is the meaning of the olive shadow crossing Figures 4a and 4b?

Reply: Thank you for bringing this to our attention. The caption for Figure 4 has been revised accordingly. Fig. 4 (a) Comparison of the observed cooling histories of the CNB magmatic plutons (data from Chen et al. 2020) with time‒temperature paths generated by models of rising magma; the green shadow represents the actual cooling age of the CNB magmatic plutons. (b) Zircon Hf isotopes and ages of coastal magmatic rocks in the SCB (data from Li et al., 2023). (c) Sketch illustrating the formation stage of the underplating magma and tectonic background during 80-110 Ma, 110-130 Ma and 130-160 Ma. The olive shadow crossing sections (a) and (b) denote the age range of 110-130 Ma.

1. Lines 242-243: What does “This period” refer to? The transition from negative to positive values? It is during 110-125 Ma. How does it correspond to 130-105 Ma?

Reply: Thank you for your feedback. We have revised the sentences as follows: " The εHf(t) values of the late Mesozoic igneous rocks along the coastal South China Block (SCB) tend to increase, peaking at positive values during 110-130 Ma (Fig. 4b)." Additionally, " This period corresponds to the magmatic lull in the SCB, which coincided with a compression phase in the CNB during 130-105 Ma (Wei et al., 2023)." Lines 260 and 270.

1. Lines 257-258: “The lower viscosity of the CNB…” How low is the viscosity? What constraints are there for the lower viscosity? How is this reflected in the model?

Reply: Thank you very much for pointing out this issue. It was suggested that once the shear zone network went through, the shear strength of the lithospheric mantle drastically decreased (Vissers et al., 1995). The CNB indicates long-term tectonic shearing activity, which would result in a lower viscosity than that in the interior of the SCB. Therefore, we assigned a low viscosity of 1e19 Pa·s to represent the active fault in the model. If the fault only changes in viscosity and the bottom magma is not thickened, then the magma at the bottom will not upflow. Shear zones control the ascent and emplacement of magmas (Weinberg et al., 2004). This implies that the shear zone should be a pathway for thermal fluid and have higher thermal expansivity (Afonso et al., 2005). Therefore, the weak CNB compared with that of the interior SCB facilitated the emplacement of mantle magma during the magmatic lull. We have added the fault condition in the model: "The model fault was assigned a low viscosity of 1e19 Pa·s to represent the CNB" Lines 186-187. In addition, we have modified the sentences and added the above discussion into this paragraph, Lines 285-295.

1. Lines 259-261: How was this conclusion reached? If it comes from the modeling results, Figure 3 only shows the results during 1-5 Myr.

Reply: Thank you for bringing this to our attention. Figure 3 has indeed been updated to accurately reflect the results spanning 10-50 Myr.

1. It would be better to separate the discussion from the conclusion.

Reply: Thank you for your suggestion. We have separated the discussion and conclusion sections accordingly. The final paragraph has been revised to serve as the Conclusion.