OVERVIEW

This manuscript presents a simple 1D thermokinematic modeling tool that predicts the thermochronologic cooling ages of crust that has experienced lithosphere-scale processes. These processes, such as burial, fault, and delamination, produce prescribed km-scale burial and/or exhumation that perturbs the crustal thermal field. The thermal and thermochronologic consequences of these different processes are quantified by forward modeling particle time, temperature, and depth paths relative to Earth’s surface and varying the lithospheric thermal field with changes in exhumation rate, heat production, Moho depth, and/or fault motion. The goal of this tool is to provide a relatively simple thermokinematic modeling tool that permits exploration of how specific geologic processes are predicted to impact low-temperature cooling ages. The authors argue that although simple, the results of their demonstration models suggest that Moho depth evolution over tens of millions of years can impact upper crustal thermal fields in ways that are resolvable by low-T thermochronology. These models also demonstrate in what ways the thermal history of rocks, and the cooling ages that document these histories, can be decoupled in time from the processes that (eventually) produced rock cooling.

This work demonstrates a new tool for investigating one of the central problems in low-temperature thermochronology: the challenge of inferring the timing and magnitude of lithosphere-scale geological processes through their impact on the thermal history of upper-crustal rocks, as documented by He and FT cooling ages. I think that simple 1D approaches, like the one used here, can provide an outstanding first-order tool for building intuition about how cooling ages should be expected to document these processes. This is a welcome addition. However, I find the manuscript as currently written and illustrated has several substantial limitations. First, I do not think the model results support the conclusion that Moho depth evolution over tens of millions of years can impact upper crustal thermal fields in ways that are resolvable by low-T thermochronology in real rocks, because the predicted variability among “perfect” model ages is small (< 15% of ages in most cases) compared to the typical uncertainty in real thermochronologic datasets. Second, 1D approaches are only useful if their simplifications, and the corresponding limitations, are very clearly articulated at the outset. I find that this paper does not do this effectively as currently written in several key ways, which would substantially limit reader engagement with the paper and the tool it is promoting. Finally, the figures, tables, and text that describe the model design and results are challenging to follow. I expand upon these comments below, and hope the authors find it useful as they revise.

We greatly appreciate the thorough review and constructive comments. These raised issues are covered in more detail below, so we respond to them individually below.

MAJOR COMMENTS.

1 Results do not support one of the conclusions

I disagree that a main takeaway from the results that the “moving Moho” produces “significant” differences cooling ages. Yes, the MM produces systematic differences in the predicted ages for especially the two highest-T systems, but the cooling age differences are, in the vast majority of cases, less than what is considered reproducible in most real low-T thermochronometer datasets. This signal would simply get swamped out by other things.

For the example highlighted in the discussion (line 416), EE1: “the resulting ZHe and ZFT ages are 3-7 Myr younger for the MM variant”. However, these ages appear to be (reading off the figure) ZHe: 21 Ma vs 18 Ma, and ZFT: 42 Ma vs. 37 Ma. In real ZHe data, for example, one would be hard pressed to convince me that such a 3 Myr age difference has geological significance, because the intra-sample reproducibility of “well-behaved” samples is generally ~10-20%. Add in the fact that ZHe ages can vary much more significantly due to U-Th compositional variability, and that makes this even more tenuous. 

Thanks for this comment. We argue that we have clearly demonstrated, with quantification, that the evolution of the Moho depth has an effect on rock thermal histories, and we have outlined scenarios for which an exhuming Moho showed a greater or less significant effect on thermal histories and low-temperature thermochronometer (LTT) ages compared to a fixed Moho. We agree that intra-sample reproducibility can be large in low T thermochronology datasets, however these forward models have the advantage of not being limited by the variability inherent in empirical datasets. Still, in the results we present, we observe percent differences in age of up to 21%, 26%, 34%, and 66% for the predicted AHe, AFT, ZHe, and ZFT ages, respectively. Though it is true in numerous cases that the uncertainties of measured ages might mask the effect of the Moho depth varying, it is notable that the average percent difference in ages for all chronometers is ~9% for all model variants, which would likely be resolvable in less dispersed measured bedrock ages. We thus argue that the high inherent uncertainties in empirical datasets are not in conflict with our reported study outcome that making an informed decision about the Moho evolution during a particular geological history of interest can increase the accuracy of the thermal history modeling exercise. We will highlight the magnitude of the differences in predicted age for the moving and fixed Moho variants more clearly in the revised text.

In addition, based on a comment from reviewer 2, we plan to calculate the Peclet number for each model, which will further emphasize the difference in heat transfer for the fixed and moving Moho model variants. We hope that this will aid readers in applying these concepts to future modeling exercises.

2 Several key assumptions and limitations are not clearly explained

Section 2.1 of the methods falls significantly short of explaining key assumptions and limitations of this 1D approach. Instead, the explanations that are given for those I highlight below (and others) are scattered throughout the methods and results sections is a way that is difficult to keep track of; it is confusing to know what applies to a subset of scenarios and what applies more broadly. This is even more important to lay out clearly for those who want to use this tool to build and interpret their own models. As a reader who has though a lot about some of these simplifications, I’m eager to know what these authors think about them! But I was frustrated and distracted by the need to scour the manuscript for this information.

Thanks for this feedback. We are prepared to introduce a new section on model design (including assumptions) that will consolidate the existing text that is distributed across other sections and address the questions outlined below. We also include specific responses below.

We would also like to note that part of the limited description of some features in Tc1D results from us concurrently working on a technical paper providing a complete description of the code features, design choices, assumptions, etc. Unfortunately, we do not anticipate the technical paper will be accepted and citable prior to the need to address the comments raised in review of this manuscript. The points raised by the reviewer are clearly important, and for those that are relevant to this text we will be adding further description to ensure the text stands alone and readers fully understand how the software works and how the results are produced.

(2a) surface boundary condition and relationship to rock uplift

This is currently described as a fixed boundary condition in the FD solution to the heat transfer equation (line 59) “Temperatures are fixed at the surface…” and in a single sentence (lines 63-64): “A variety of erosion models are built into Tc1D allowing users to explore many different burial (i.e., negative erosion) and exhumation histories (Whipp, 2022).”

But…what is a built-in erosion model? This is very vague, and it sounds like it could be a geomorphic transport law. The reference in Whipp, 2022 is a link to the zenodo page that hosts the code, but I'm not clear on what the purpose of this citation is—what is the reader being referred to exactly? This Zenodo page doesn't provide any obvious narrative details about the "many different burial and exhumation histories" and how they have been conceptually and/or numerically constructed. Perhaps this is embedded in the documentation, but I think it is asking too much of the reader to dig into that.

We agree that this information was unnecessarily hard to find, and in the revised version, will include a clear description of what constitutes an erosion model in Tc1D as well as better description of the relevant erosion models used in this manuscript. In addition, a copy of the erosion model descriptions from the documentation for this version of Tc1D will be included in the supplementary data.

I think the essential thing that needs to be described is how the simple surface boundary condition defines the relationship between rock uplift, surface uplift, and rock exhumation (the latter being only one of these three things that actually cools off rocks).

Conceptually, is the model space “hung” (and fixed) at depth = 0 km, without a sea level datum? And therefore erosion “scenarios” are simply changes in rock uplift rate through this fixed surface? Or, are “erosion scenarios” superimposed changes on independent rock uplift rates?

This is an important point raised by the reviewer, which was not clearly described in the earlier text. We will expand the description of the boundary conditions to note that the surface elevation is fixed within the model at 0 km, and thus any velocity toward the surface reference elevation is an exhumation velocity. Generally, this will also be the rock uplift velocity, however, if the surface elevation were to change in some geological scenario (e.g., isostatic uplift), the rock uplift velocity would be a combination of the exhumation velocity and the velocity due to surface elevation changes. Internally in Tc1D, this difference would not be seen, but could be imposed via changes in the surface temperature to simulate surface uplift or lowering as a proxy. We will clarify this concisely in the text.

A depth = 0 km datum certainly makes sense from a thermochronology perspective (since cooling ages are agnostic about surface elevation relative to sea level), and I think it would make practical sense to anyone who has designed a similar 1D model themselves. But readers who build their own thermokinematic models are not the audience for this work (e.g., line 35). Critically, if I understand what is being done here (surface is fixed at 0 km depth, no surface uplift?), it means that rock uplift is being imposed with no surface uplift (i.e., exhumation is keeping perfect pace with rock uplift). Imposing rock uplift with no surface uplift is a huge simplifying assumption with cascading thermal and geodynamic implications. Moreover, it’s a significant deviation from how the Earth works in many of the scenarios of interest here (or else we’d have no relief at Earth’s surface). This simplification and its consequences are not clear right now.

We thank the reviewer for raising this point, and agree that we were insufficiently clear about the surface boundary condition and its implications. We agree this is an important assumption that should be very clear to the target audience, and will be clarifying these aspects in our revised version as described above.

The Thrust Fault models are the extreme example, because motion on a thrust fault alone cannot cool off rocks; thrusts cannot bring particles closer to Earth’s surface (unless accompanied by erosion, a separate process with its own drivers). But, in these models, when a particle in the HW of a thrust fault moves up a ramp, it is assumed that erosional exhumation is perfectly keeping pace with that motion. Indeed, this assumption is mentioned in passing (line 127), but its implications are not discussed.

In sum, if rock uplift always strictly equal to rock exhumation in these models, in both magnitude and timing, this needs to be explicitly stated and discussed up front. This simplification also needs to be related to the key takeaways of this paper, for example: “our results show that the recorded response of thermal histories/thermochronometers in the upper crust and geological processes that disrupt the crustal thermal field may be disassociated in time, because of the time and length scales of different heat transfer mechanisms.” (abstract, lines 9-11).

We thank the reviewer for the opportunity to clarify how Tc1D deals with surface uplift, rock uplift and exhumation and will clarify these important points as noted in the comments above.

(2b) thermochronometer behavior

The paragraph at line 65 in section 2.1 that describes the basics of the He and FT age calculations is missing some key information. For example, does the model only predict one age per system? Can the user control key decisions about how these ages are calculated, such as grain size and U-Th composition in the He systems? Of course, a lot of the resolving power of these chronometers comes from age-eU trends (He) and track length distributions (FT); which are not not mentioned and I assume not being predicted here; but this is a critical simplification that should be emphasized more. The radiation damage models for the He systems are being implemented; this is great, but it means grain composition is another potential ‘knob’ that will change predicted He ages in some scenarios, not being explored in the examples here. But, a hypothetical future user of this tool, upon sitting down to design a model run, would immediately be presented with the challenge of choosing particular grains to model. I think more clearly describing in this section what the consequences of only predicting single ages are would significantly improve the communication of the method.

Some of this information is currently relegated to section 2.3, after the extensive overviews of each of the geodynamic scenarios in Section 2.2. For example: “Hence, most of the scenarios explored below involve a starting depth for which open system behavior is expected, and zircon or apatite should not have accumulated any alpha damage (∼ 300C). Where this is not the case, the implications are discussed”. However, I cannot find any discussion of these implications in the rest of the text. This is another example of relevant assumptions and their implications being challenging to find and keep track of in the current manuscript.

Thanks for this feedback. We agree that eU and grain size,  and fission track length distributions are key information in LTT studies, particularly in inverse modeling of thermal histories. Since our focus in this contribution is to explore the thermal histories resulting from forward modeling evolving temperature conditions of a rock particle moving through a dynamic geothermal gradient over geological timescales, we felt that it would distract from the focus of the manuscript to also explore the effects on resulting ZHe and AHe ages of different eU and grain size conditions. Rather, it was important to keep those parameters constant between scenarios to be able to compare LTT ages between models (values used are listed in Table 2 of the manuscript). The thermal histories of the forward models are not influenced by those variables, only the calculated LTT ages. However, the Tc1D software allows the user to specify any desired zircon and apatite grain radii, U concentrations and Th concentrations, which we will clarify  in the text for potential users. Furthermore, any modeled thermal history output from Tc1D can be imported into the TcPlotter software (https://doi.org/10.5281/zenodo.6341671) that we presented in Whipp et al. (2022, GChron; https://doi.org/10.5194/gchron-4-143-2022), to explore a full range of eU and grain size relevant to a particular study area. We agree with the reviewer that these specifics were not adequately explained in the manuscript, and we propose to clearly articulate this information in section 2.1 in the revised version.

Since for the purposes of this manuscript we are using Tc1D in a forward modeling capacity, track length distributions of the AFT system are an output, not an input. These are produced for every forward model but have not been presented here. We are happy to include these in the supplementary data and to note that they are standard outputs from Tc1D so that readers are aware.

(2c) starting conditions

The starting thermal and age conditions may have a strong control of the model results, but the starting conditions are not completely explained or justified here. As someone with modeling experience, in most cases I can figure out why these models are designed the way that they are; but, I shouldn’t have to spend time doing this as a reader, and again, experienced modelers are not the primary audience here. Some of my questions are:

Why do some scenarios start with 5 Myr of thermal steady state, whereas other do not? Why 5 Myr? Is that sufficient in duration to equilibrate both thermal and age structure across the model domain? (There are key length-scaling relationships here). And what actually are the “steady state” conditions in each model that uses them to set the thermal and age structure of the lithosphere? Are they different between the various scenarios?

What are the geological implications of that 5 Myr of steady state? For someone designing their own scenario, how should they determine whether they need a period of steady state prior to implementing a geological process of interest?

Thank you for the opportunity to provide clarity on our modeling design. First off, we will add text to clearly state that all models start from an initial conductive steady state. We realized during modeling that initiating faulting or delamination at the start of the simulations made it difficult to see their thermal effects. Thus, we start with 5 Myr of additional time at a thermal steady state prior to faulting and/or delamination to clearly show how the thermal field changes. For comparison to the other models presented for other scenarios, we then simulate 50 Myr of thermal evolution thereafter. This will be made clear in the revised text.

3 The text, figures, and tables that communicate the model designs, results, and implications are difficult to follow and comprehend as written

Starting in Section 2.2, and for the rest of the paper, I found I had to simultaneously be looking at Table 1, the figures presenting the tT and Tdepth info (ex, figure 3), and the relevant text in order to follow any of the main points. Few readers will have the time and patience to constantly flip between Table 1, the text, and the figures.

I think revision of the tables and figures would address this problem. The biggest barrier from my perspective is the figures; especially figures 3, 4, and 7 do not stand on their own at all. The figure captions don’t say anything useful, there is no annotation, and the only way to understand how the plots relate to geologic process, model design, or the variables of interest—let alone the implications and limitations of the results, e.g., don’t pay attention to the ZFT ages (line 220)—is to have Table 1 and multiple sections of the text also visible. There doesn’t seem to be a practical (figure design limitation) reason for this. In figure 3, for example, identical legends take up space in all 12 panels. Scenario IDs, which contain no specific information, are the only other label. The panels are not visually organized by scenario type or question, and none of the model design information is annotated. There are many opportunities for improvement here; for example, the information in Table 1 is laid out on a time axis, why not just annotate the essential information from that table onto these panels?

We appreciate this feedback and the suggestions for improvement and propose to incorporate them into the figures in the revised manuscript. In particular, we plan to remove Table 1 and integrate its information directly into Figures 3-7 as inset tables with relevant model information above the existing plots.

A few other related comments regarding the figures and tables:

Figure 2, All panels are all given the same visual weight in this figure, but in fact panels (a) and (b) apply to c-g.

We propose to note this more clearly in the figure caption.

Figure 3: The x-axis is both time and predicted cooling age (right?). It would be useful to describe what the significance of plotting the predicted cooling ages here means. I assume the ages are simply placed at the depth that corresponds to that time in the forward model? What is the best way to think about that information geologically, if one wants to export this information into a geological framework? This is a not typical way of visualizing thermochron data (though I think it is an effective way to present the model results).

Thanks for the question, and the opportunity to provide additional clarity. The left axis shows the modeled particle’s thermal histories (typically there are two thermal histories for each model, one for a moving Moho scenario and one for a fixed Moho scenario), and the right axis shows the modeled particle’s depth history (same, regardless of the Moho position). The depth history is a function of the exhumation/burial history of the modeled particle, and can vary from the thermal history in some cases. The thermochronometer ages are plotted on the thermal history at the time corresponding to their age. Thus, the reviewer is correct in noting that time on the x-axis is both past time in the thermal model and ages for the predicted thermochronometer systems. We will clarify that the ages are plotted according to where they are in the thermal history at that time.

Where ZFT data should be ignored (in models where total exhumation is 10 km? line 220), modify the symbol in the figures to make this clear. Does this limitation apply to all models, or just those in the EE group?

We appreciate this suggestion and will include a modified symbol for ZFT data in cases when it should be ignored.

Table 1. Eliminate some unnecessary abbreviations. SS = steady state (but there’s enough space to just spell it out). IHP, the same; consider also using numbers for IHP instead of high vs. avg.  IHP is called “volumetric heat production” in table 2; consider using same terminology for clarity, unless a (conceptual) difference is reflected in the nomenclature difference. “CE 5 km” and “Erode 5 km, constant rate” is identical information, but the different phrasing makes it seem different; erosion rates would be a welcome addition to this table.

Thanks for these – we will make the suggested changes in the revised manuscript.

Reviewer 2: Review of manuscript entitled “Simplified modeling of the impact of lithospheric-scale geological processes on thermal histories and low-temperature thermochronometers” by Dawn Kellett and David Whipp.

This manuscript reports the results of a series of numerical experiments in which the heat equation is solved to predict temperature-time paths of rock particles exhumed to the surface from which synthetic cooling ages are derived. Various scenarios are envisaged representing different tectonic and erosional events. The different age patterns are then compared to each other. The manuscript concludes that a factor that has been neglected in previous studies, i.e., whether the Moho depth is kept constant or not during the various events, has a significant impact on the predicted ages. The authors also compare the ages distributions associated with different geological events.

I enjoyed reading this manuscript. It is relatively well written and presented, apart from some figures being at the limit of readability (too small font size), in particular those illustrating the main results (Figures 3 to 7).

We appreciate the thorough review and positive comments.

1. A) I have some important concerns about the content of the manuscript that I will now express, with some recommendations on how to address them:
   
   
   
   Most models include a long phase of imposed exhumation that ends today; this implies that most of the tectonic events considered are very ancient; this is confirmed by the predictions that almost all ages are set during the “exhumation” phase, not the “tectonic” phase. This is not a problem per se, but it needs to be made clearer in the abstract and introduction of the manuscript, i.e., that most scenarios focus on “past” tectonic events. Personally, I think the manuscript would have a much greater impact if a broader range of scenarios that include many cases of recent and/or on-going tectonic events were considered.

This is a good point raised by the reviewer and we appreciate it. For the sake of plotting the results and inter-model comparisons we have opted to use simulations with a duration of 50–55 Myr. This provides time for heat to be conducted from the Moho to shallower depths in delamination models, as well as time for both tectonic and post-tectonic exhumation. We also felt this longer duration was a “Goldilocks” choice, as later stages of the model evolutions could be linked to regions of recent tectonic activity, while the longer duration could be relevant to more slowly eroding areas or those with older tectonic events.

Further, the modeled thermal history should begin at a point of no daughter/alpha damage accumulation for the thermochronometers of interest. We aimed to explore both lower (AFT, AHe) and higher temperature LTTs (ZFT, ZHe), so we aimed to have the tracked particle start or reach ~300 C during the model run and also reach the surface by the end of the model. We agree this is a limitation of our approach.

Regardless, the point raised is an important one, and we intend to explore a few additional model scenarios that might be more relevant for "younger” areas, and would supplement our results with these scenarios to have more examples of cooling ages corresponding to the time of tectonic (and associated erosional) activity.

1. 
   
   The amplitude of the exhumation phase that ends almost all scenarios considered in this analysis appears quite arbitrary and thus, if I understand properly, unrelated to the prior tectonic event. The various tectonic events therefore appear as “thermal preconditioners” for the late exhumation phase rather than being directly responsible for the shape of the cooling curve experienced by the rock particles when their ages are set. This, in my opinion, lessens the impact of the manuscript as the effect of the various tectonic styles/events is much smaller due to this strong “filtering” by the exhumation phase.

We agree that the tectonic scenarios we explored act as thermal preconditioners, and that late exhumation dominates in most scenarios. However, we think this is an important outcome for the growing number of researchers working on orogens/geological settings that are 50 Ma and older. As noted above, we plan to include a few additional results in the revised manuscript with ages that better correlate with the timing of tectonic activity. In such cases, we also intend to better balance the estimated surface uplift with the post-tectonic erosional exhumation, for example.

1. 
   
   This disconnect between the amount of exhumation and the nature/amplitude of the assumed tectonic event may also cause inconsistencies in the proposed scenarios. Indeed, if the exhumation phase is meant to represent a post-orogenic eroding phase, its amplitude should be in proportion to the relief/topography generated during the tectonic event. I take as an example, the delamination scenarios. Delamination causes surface uplift but the recovering of the original lithospheric temperature structure (regrowth of the mantle root) that is modeled during the exhumation phase should lead to a lowering of the surface back to its original value. It is not clear to me how the amplitude of the exhumation phase is set, i.e., is it arbitrary or in connexion to the amplitude of the tectonic event? This needs to be stated more clearly and justified. Ideally I would like to see the 1D thermal model coupled to a simple 1D isostatic model to make sure that the two parts of most scenarios are physically connected.

We appreciate this comment and see the reviewer’s point about a possible disconnect between the erosion history and event giving rise to uplift. Our approach in this manuscript is to present both a suite of generic model predictions for various tectonic and erosional scenarios and to provide similar magnitudes of exhumation and model durations to enable them to be compared to one another easily. That said, these generic scenarios may not be representative of specific geological events in a self-consistent way, as a result. Thus, as noted above, we intend to provide supplementary results to demonstrate scenarios that avoid the perceived disconnect between the tectonic event and subsequent exhumation.to better enable comparisons between models. For example, Tc1D already includes calculation of the Airy isostatic response to changes in lithospheric mass (and density) during delamination and/or erosional exhumation, so it is straightforward to estimate the magnitude of surface uplift predicted for a total delamination and how its surface elevation effects would trigger an erosional response. We will demonstrate such a case in the supplemental new results. However, we note that there are many cases in nature in which surface uplift due to isostasy is not balanced by erosion in a straightforward manner, or is dependent on the geological setting (e.g., within vs. at the margins of a plateau region such as the Tibetan and Altiplano plateaus). These are challenging situations to address in a 1D model.

1. B) I also have other points that I believe need to be addressed:
   
   
   
   The timing of the various phases of each scenario appears quite arbitrary (see for example the first paragraph of section 2.2.1). Is it the case? If not the rationale should be better explained.
   
   
   
   

Each set of scenarios is designed to explore a range of conditions. For example, in 2.2.1, we present 12 erosional exhumation models to explore: high vs. average internal heat production; fast followed by slow constant erosion vs slow followed by fast; smaller (10 km) vs larger (20 km) magnitude of erosion; and short, medium and long exponential decay constants. This is intended to explore the characteristics of resulting thermal histories to better understand their responsiveness to different styles of erosion, burial and exhumation, etc. We agree that a study design section could be added to explain our rationale for these experiments, and propose to add this to the revised manuscript, also addressing Reviewer 1’s recommendation above.

1. In the FM case illustrated in Figure 2c, there is no velocity in the mantle, despite the downwards flux of material related to the subduction of the mantle part of the lithosphere as illustrated in Figure 2a; as shown in Batt and Braun (Geophysics Journal Int, 1996), this downward flux is important as it compensates almost exactly the upward flux in the overlying crust to maintain the Moho at a quasi constant temperature. Shouldn’t it be included?

This is a good point! We had designed the FM models to allow crustal uplift above a flat base, similar to the early stages of continental collision. However, the geometry presented in Figure 2a does indeed include a downward velocity that is not included in our present results. One challenge is that the rate of downward motion depends on both the dip angle of the lower plate and the partitioning of convergence between uplift in the hanging wall and subduction of the footwall (e.g., the Alpine Fault with a horizontal mantle velocity, vs. the western Andes where subduction angle is steeper). Hence, it may be difficult to provide some generic results demonstrating this effect in a similar way to the existing results. That said, we plan to run a few additional exploratory fixed-Moho models with a downward mantle velocity and include the results in the supplementary material. As this scenario is likely to produce an even larger difference between the fixed and moving Moho variants, we feel this will also highlight the importance of considering how the Moho depth evolves in thermal models.

1. 
   
   It is not clear to me whether there is any surface erosion/exhumation during the thrusting event in the TF models? Can it be clarified?

Yes, there is exhumation during the thrusting event in the TF models, as visually depicted at the surface in the schematic models shown in Fig. 2e and Fig. 5d and f, but we agree that this should be more clearly explained. We will clarify this in the revised version.

1. 
   
   In section 3.1 (and following parts of the model results section), the explanation of these results would gain a lot, in my opinion, if it used the concept of the dimensionless Peclet number. Indeed the ‘curvature’ of the thermal field due to the advective component of heat is known to be controlled by the value of the Pellet number that is proportional to the exhumation velocity and the thickness of the layer bing exhumed. The effect of heat production can be included as shown in Batt and Braun (1996) - Appendix A. 
   
   
   
   

We thank the reviewer for this insightful comment. We had avoided mentioning the Péclet number in the first draft of the manuscript because it may be less familiar to some of our intended readership. That said, we agree this would be a nice way to illustrate why the thermal histories differ for different model variants, and why the fixed versus moving Moho models have differing predicted ages. It may also simplify some of the model explanations that can be wordy in places currently. We will introduce the Péclet number and use it where relevant in the revised text.

1. 
   
   Comparing the results of the DL and EE scenarios implies that it would be very difficult to derive the existence of a delamination event from an age dataset, especially if the delamination event is “ancient”, i.e., it is followed by a long period of passive exhumation; at most, the ages would be interpreted as yielding an effective cooling rate that is slightly faster than that caused by the slow post-delamination exhumation phase; this should be stated in the discussion or in the results section.

Yes, we agree that the effect on cooling rates is an important outcome of our study, and will make sure to emphasize this result more clearly in the revised version. We are fairly convinced it is challenging to detect delamination from thermochronometer ages alone, unless they are related to a magmatic response to delamination, or erosion history is tightly constrained.

1. 
   
   In the first paragraph of the discussion, it is stated that most ages are unrelated to the timing of the tectonic event; this is, in most part, due to the imposed long-lasting exhumation phase. In the scenarios that are envisaged most ages are set during the exhumation phase implying that they cannot record the earlier tectonic event; so this conclusion is, in my opinion, a direct consequence of the amount of exhumation that is imposed in the models.
   
   
   
   

This is indeed true, and something that be challenging when interpreting the ages measured in a low-temperature thermochronology dataset, particularly when working in tectonically inactive regions. As noted in the earlier responses, we intend to provide some supplementary results where the timing of tectonic activity is better recorded in the predicted ages in order to demonstrate additional scenarios where this occurs. For example, we intend to include a model variant with a late delamination and erosional response, which may record the timing of enhanced erosion following delamination in the lowest-temperature systems.

1. 
   
   In section 4.2 of the discussion, the reader is warned that different age patterns/values can be predicted from identical scenarios but different thermal state (amount of assumed heat production for example). This, in fact, is equivalent to say that to interpret ages in terms of denudation rates, one needs to know the geothermal gradient. It has been shown that this can be avoided by measuring ages along elevation transects (Gleadow and Fitzgerald, 1987). So, in my opinion, this part of the discussion is more a reflection of the limitation of the 1D model.

We agree that vertical transects can be used to estimate rates of erosional or tectonic exhumation, so this is a good point we had not presented in our original text. However, estimation of rates of exhumation from vertical transects can still be challenging due to topographic effects, transience in rates of exhumation, and other factors (e.g., Ehlers, 2005, MSA). While vertical transect data would still present challenges for data interpretation using a 1D thermal model, there are ways to account for various topographic effects (e.g., incorporating the approach of Ketcham, 2025, GChron into Tc1D) to provide a means for data interpretation in settings where the rates of exhumation may have varied during sample exhumation. We will provide additional text in the revised version of Section 4.2 to clarify this point as it is an important consideration.

1. 
   
   The effect of a cooling and exhumed granite on the distribution of thermochronological ages in its vicinity (as discussed in section 4.3) have also been investigated by Murray et al (2018); this work could be cited here.

Good point, we will add this reference as suggested.