General comments

The authors present a new, open-source Python package designed to forward model large numbers of low-temperature thermochronometric ages from time-temperature paths generated by geodynamic numerical models. The manuscript, and Python package, comprise timely and potentially substantial contributions to the fields of geochronology and geodynamics – as the authors mention in the paper, advances in geodynamic numerical modeling have enabled increasingly high-resolution simulations of geological and thermal events across various tectonic settings and thermochronology data provide quantitative constraints against which such models can be validated. The GDTchron Python package presented here expands on thermokinematic modeling approaches (which predict time-temperature histories and low-temperature thermochronology ages based on deformation kinematics and topographic evolution but do not directly incorporate dynamics; Braun, 2003; Braun et al., 2012; Almendral et al., 2015) and previous strategies for integrating time-temperature information from geodynamic models with thermochronology data (which has involved exporting t-T paths from geodynamic models and importing them into thermal history modeling software like HeFTy or QTQt; Ketcham, 2005; Gallagher, 2012) in that it allows for extraction of time-temperature histories, and thermochronology age predictions, for a large number of particles from multiple steps in a geodynamic model.

The paper is well written, the accompanying Python notebooks are well documented/commented and readily accessible via Github, and the figures are generally illustrative. However, this research has several areas that I consider in need of major improvements before manuscript publication and widespread implementation of the Python package by users in the geodynamics and thermochronology communities:

Specific comments

1. The evolution of topography can profoundly impact thermal histories and predicted thermochronology cooling ages, particularly for the low-temperature systems that this contribution focuses on (e.g., Braun, 2005; Gilmore et al., 2018). The authors briefly mention this and caution users to implement geodynamic models that account for erosion and sedimentation – and even make changes to the parameters used to approximate topographic evolution to “provide more realistic surface processes” in their rift-inversion orogenesis example model. However, particularly considering potential readers/GDTchron users who may have specific expertise in either geodynamic modeling or thermochronology, the authors need to include more detailed information on i) how surface processes and other parameters that affect thermal evolution, like heat production, are incorporated into geodynamic numerical models, ii) the extent to which those choices may impact t-T histories and predicted cooling ages, and iii) related recommendations or considerations when integrating geodynamic models and thermochronology datasets. A more in-depth comparison of the strategies for incorporating surface processes in the Vasey et al. (2024) rift-inversion model, versus the new rift-inversion model showcased in this study, would provide a valuable, and cautionary, case study. Both of these models incorporate the same geodynamic evolution, but produce notably different low-T cooling ages – a change solely related to the mechanism of topographic evolution. These results raise questions regarding the sensitivity of predicted cooling ages, and measured thermochronology data, to geodynamic mechanisms versus surface processes that need to be addressed in this manuscript to help educate potential GDTchron users and establish best practices.
2. Though the authors highlight advances to understanding of thermochronometer systems and diffusion/annealing kinetics in their motivation for this work, they have not incorporated the most widely used and up-to-date kinetics for predicting (U-Th)/He or fission track cooling ages or track length distributions (e.g., Enkelmann and Garver, 2016; Flowers et al., 2009; 2023a; 2023b; Guenthner, 2021; Guenthner et al., 2013, 2014; Ketcham, 2019; Ketcham et al., 2007; Tamer & Ketcham, 2020) in their Python package. This is flagged as a potential future direction, but until these updates are made, users will either be limited in comparisons of geodynamic model predictions to thermochronology data or existing thermal history models, or will need to follow the same workflows of previous studies that exported t-T paths predicted by the geodynamic models for interpretation using other software like QTQt or HeFTy (as performed in Ternois et al., 2021). I encourage the authors to incorporate apatite fission track annealing kinetics after Ketcham et al. (2007), and options to use radiation damage and annealing models for the apatite and zircon (U-Th)/He systems (e.g., Flowers et al., 2009; Guenthner et al., 2013; or more recent updates), in this publication. Baseline parameters for U and Th concentration, grain size, Dpar, and other parameters could draw from standard values– though users would benefit from options to include specific values at specific sample locations to compare with thermochronology data. Routines for the implementation of these kinetic models, and grain parameters, in the PecubeGUI (Bernard et al., 2025) may provide a useful framework for such updates.
3. While related references are briefly mentioned and cited (e.g., Ketcham, 2005; 2024; 2025; Gallagher, 2012; Braun, 2003; Braun et al., 2013; etc.), I’d like to encourage the authors to include a more in-depth discussion of how the GDTchron package fits into the framework of related work, including 3-D and 4-D thermokinematic and 2-D (time-depth and time-temperature) thermal history modeling. This may be particularly relevant to discussions and recommendations related to recommendations/best practices for incorporating crustal thermal parameters, surface processes, and disentangling the effects of each system on predicted cooling ages.
4. Finally, the tendency of the implemented interpolation scheme to assign new syndeformational deposits the same thermal history as adjacent, pre-existing particle paths may be a real limitation in areas that involve significant accumulation of such deposits during polyphase tectonic and geodynamic histories (like the Pyrenees studies cited as examples in this manuscript, Capaldi et al., 2022; Curry et al., 2021). Are there any changes to the interpolation program, as implemented, that would help set up for future changes in ASPECT or other geodynamic models that might include better tracking and ID of particles related to syndeformational sedimentation?

Technical corrections

Line 54: Specify that the simplified model involves block uplift here?

Line 55: Because routines incorporating erosion and sedimentation are so important in predicting cooling ages for low-T thermochronometers, I strongly recommend providing more information about what these routines were, & what geological or other evidence supports their use, here in the intro in addition to later in the paper.

Line 109 / Section 3.1: What kind of computational resources are required to run GDTchron? Can you provide more information here regarding the types of scenarios that might run on a personal computer/workstation, versus HPC? (I recognize that related information is documented within the Github & other package resources, but would be valuable to include for readers in this section of the paper too)

Line 141: As a reviewer with expertise in thermal history and structural/thermokinematics modeling but no hands-on experience with geodynamic software like that listed above, I am curious what common inputs and outputs to these software are, particularly as they relate to the temperature field, partical motion, and the development/decay of topography. Even though your target audience may already be familiar with these programs, I recommend including a more comprehensive background/context section that helps frame geodynamic model capabilities, and limitations, for a broader audience. E.g., how is topography typically handled in these systems, and what are recommendations/considerations for evaluating incorporation of topographic change and sedimentation in geodynamic models? In addition, there may be readers with expertise in thermochronology, but not ASPECT etc. modeling, who would like to use GDTchron to extract t-T paths and predict cooling ages based on published geodynamic models to help inform their sampling strategies. Providing a little more context for the geodynamic modeling component would help clarify the intended workflow and capabilities of the GDTchron package.

Line 169: Choose a different word than “flow” when discussing how material exits the top of the box? “Flow” has implications for rheology/deformation style that do not necessarily apply here.

Line 170: By “temperature structure of the box remains constant,” do you mean that there is no advection or change in isotherms due to particle uplift & exhumation? It appears so, from the model diagrams? If not, please consider revising this simple example to allow for isotherm advection, which will produce a more geologically accurate case study?

Line 170: Specify "particle / rock uplift" or "exhumation" here, as there is no surface/topographic uplift in this scenario that maintains flat topography (e.g., has complete erosional efficiency).

Lines 176-177: Language related to Fig. 5 here is a little confusing – could be interpreted to mean that the timesteps were generated with GDTchron, but they were defined in ASPECT, right? Then the outputs of that model were processed in GDTchron?

Lines 179-180, sentence starting with “Ages young with depth…”: Phrasing makes it seem like 0 Ma is reached spatially below the PRZ/PAZs, but the unit is time, not space. Change to "Ages young with depth through the partial retention/annealing zones of each thermochronometric system. At the base of these zones, cooling ages are 0 Ma because diffusion and annealing outpace He and fission track production" or similar?

Line 190 / Fig. 5: I strongly recommend devising a way to show the movement/change in particle positions in this figure - as vectors, a box with no fill showing before/after positions, positions of partial annealing/retention zones before and after uplift, etc. For Y-axis for Figs. 4 and 5, I strongly recommend flipping the axis such that 0 km is at the surface, and values increase with depth (or negative elevation). This recommendation uses the same perspective as structural cross sections and other geodynamic models, and I think would make both thermal structure and predicted cooling age plots more intuitive for readers and potential GDTchron package users.

Lines 194-196, sentence starting with “One could explore the effects…”: Predicted cooling ages for all of these scenarios also depend on the mechanism for incorporating or approximating topographic evolution – which may have just as much impact on the distribution of predicted cooling ages as the aforementioned factors, especially for low-T thermochron data in regions of relatively slow rock uplift/exhumation (Gilmore et al., 2018; Ketcham, 2025; & references therein).

Lines 200-204: Completely agreed – though I think it would be helpful to articulate how these factors are incorporated into geodynamic models to help educate readers and ensure potential GDTchron users make effective choices when coupling geodynamic models and interpretations of thermochronology data.

Lines 214-215, sentence starting with “Although thermokinematic models…”: Consider including a brief discussion of the similarities & differences of each approach, e.g., Pecube versus GDTchron? What is each good at (e.g., Pecube incorporates a sophisticated framework for estimating topographic evolution, Pecube-D tracks thermal histories of particles involving large amounts of horizontal shortening/extension and complex faulting, statement on how these relate to the conditions of most geodynamic models?)

Line 220 / Fig. 6: Please clarify that ages are relative to the start time/duration of the model. Consider flipping Y axes for both the model grid (so that 0 km is at Earth's surface and depth/negative elevation increase downwards) and the Surface Ages plots (so that 0 Ma is at the bottom, and older ages are on top - this is the framework commonly used for Pecube-D and other thermokinematic models). What do you mean by "Maximum Age (Ma)" as Y-axis label on bottom charts? Seems like these are just predicted cooling ages? Please consider adding some vertical exaggeration to the middle plots, and perhaps intermediate ticks/values to the Plastic Strain and AHe Age gradient keys – predicted cooling ages only use a third of the space, and it is challenging to read these gradient values without zooming in 200+% into the figure.

Line 229, sentence starting with “This 2D model…”: Which model? That in Vasey et al. (2024) using hillslope diffusion to approximate surface processes, or the revised model in this study that uses Fastscape? Here, it would be valuable to use this example, and the changes made to "provide more realistic surface processes", as a framework for discussing what changes were made to the mechanisms of approximating topographic evolution and their influence on predicted cooling ages.

Line 231, sentence starting with “At 16 Myr…”: Time reference is unclear here - do you mean  after 16 Myr into the model run, or at 16 Ma? At the end of the sentence, it would be helpful to add “… with a total model run time of 36 Myr" or similar.

Lines 236-237, “resulting in slightly younger AHe (~10 Ma) and AFT (~14 Ma) ages…”: Here, can also emphasize that these predicted ages are partially reset with respect to inherited model age?

Line 240: Update “…lithospheric thickening and uplift” to “lithospheric thickening, uplift, and erosion”

Lines 259-261, “Implementing these alternative kinetic models…”: Why not incorporate the most widely used kinetic models in this first release, especially if kinetic models are styled after those used in Ketcham (2005)? Doing so would significantly expand the potential for comparisons between predictions based on geodynamic models and thermochronology data for specific study areas.

Lines 263-265, “Importantly, the ages produced when applying outputs from geodynamic models are dependent on parameters that would be highly variable among samples within real systems (e.g., U and Th concentrations, grain size, Dpar)…”: Here, viable solutions might include incorporating values based on best practices (e.g., for grain sizes for He analyses) or values based on standards, then allowing for users to adjust GDTchron sample parameters to best match thermochronology data of interest? This suggestion is inspired by similar options in thermokinematic modeling programs including PecubeGUI (Bernard et al., 2025).

References not included in the manuscript:

Almendral, A., Robles, W., Parra, M., Mora, A., Ketcham, R. A., & Raghib, M. (2015). FetKin: Coupling kinematic restorations and temperature to predict thrusting, exhumation histories, and thermochronometric ages. AAPG Bulletin, 99(8), 1557-1573.

Bernard, M., van der Beek, P., Colleps, C., & Amalberti, J. (2022, May). PecubeGUI: a new graphical user interface for Pecube, introduction and sample-specific predictions of apatite (U-Th)/He and 4He/3He data in the Rhone valley, Switzerland. In EGU General Assembly Conference Abstracts (pp. EGU22-2277)., 524–525, 1–28, https://doi.org/10.1016/j.tecto.2011.12.035, 2012.

Gilmore, M. E., McQuarrie, N., Eizenhöfer, P. R., & Ehlers, T. A. (2018). Testing the effects of topography, geometry, and kinematics on modeled thermochronometer cooling ages in the eastern Bhutan Himalaya. Solid Earth, 9(3), 599-627.

Ketcham, R. A. (2025). Incorporating topographic deflection effects into thermal history modelling. Geochronology, 7(3), 449-458.

This manuscript presents a new open-source Python package for predicting cooling ages for AHe, ZHe, and AFT chronometers from time-temperature (tT) histories extracted from the results of geodynamic numerical simulations. It particular, it addresses the challenge of the large number of tT paths that are produced by such simulations that need to be forward modeled in order to validate geodynamic simulations with thermochronologic observations. Using established approaches for numerically integrating diffusion and annealing behaviors along a particle tT path, this modeling package can take in outputs from common geodynamic modeling software. It also addresses some of the known problems with extracting specific particle tT paths from such numerical results using an interpolation scheme. This contribution explicitly presents an open-source starting point for developing tools that are capable of jointly leveraging the power of geodynamic simulations and thermal history analysis of low-temperature thermochronologic data. As such, some features (such as the implementation of thermochronometer kinetics models) are overly simplified, with the invitation for community-generated innovation and updates that can implement more complex and realistic behaviors as needed.

The central challenge that this new Python package is designed to address is interesting, important, and relevant within the scope of GChron and the Special Issue it was submitted to. The accompanying online documentation and demonstrations are a laudable companion resource. Creating accessible tools to bridge the divides between the geodynamic and thermochronologic communities has the potential to support advances in multiple geoscience disciplines. However, as a thermochronologist, there are a few things about the design and description of the proof-of-concept examples presented here that I think would benefit from some adjustment, which I describe below. I think with revision, this will be an excellent contribution.

MAJOR COMMENTS

1. Before walking through specific examples, I have one general comment for all figures and corresponding text: the axes on many graphs are presented in the opposite sense from what is typical for thermochronological data and exhumation studies, which overcomplicates the presentation of the numerical results. In most cases, this is because the depth (y position) and time (elapsed model time) information is being directly exported from the simulations without being converted to and presented in a geological reference frame. These conversions are simple, but for thermochronological data, this matters! It matters because depth below the surface corresponds to closure behavior, and cooling ages are in Ma (millions of years before present) not model time (Myr elapsed). Users (and readers) need to be able to efficiently evaluate if the numerical results “make sense” thermochronologically. I provide some details below, but overall I suggest that unless it is essential to present the modeling reference frame, all the results for each example be reframed as depth = depth below the surface and time = geological time (Ma). I realize that following these suggestions may require parallel updates to the very nice online documentation the authors have prepared; however, I also anticipate that just a few simple lines of code are required to flip axes or convert model depth to geological depth.

2. Forward modeling example (Section 2.3)

   1. The forward modeling demonstration has a strange starting condition from a thermochronologic perspective that is not explained. The modelled tT path starts at 100 C, but it is used to predict ages for thermochronometers that are completely (ZHe), partially (AFT), or not at all (AHe) closed at 100 C. The manuscript should either: be explicit about the assumptions of such a model design and describe clearly how the consequences of these assumptions produce the model results; or, have a starting condition that is hotter than the temperature sensitivity of all systems being modelled. This is not sufficiently discussed in the current text. For example, at line 95, the text mischaracterizes the reason why the ZHe age is predicted to be 34.8 Ma. “For the zircon (U-Th)/He (ZHe) system, this path results in an age of 34.8 Ma, corresponding to nearly the full duration of the 35 Myr history, given that the sample remains at a temperature largely below the ZHe partial retention zone throughout this history.” (emphasis added). The starting condition (and isothermal hold T) is not “largely” colder than the ZHe partial retention temperatures for the Reiners et al. 2004, it is much much colder (at least 50C colder) than even the coldest part of the PRZ for these kinetics (and for the current ZRDAAM kinetics too, for such as tT history). From a thermochronological perspective, this zircon grain simply “appears” at 100C at 35 million years ago, and accumulates He that whole time. Having the model start in the midst of the AFT PRZ is even more complicated and potentially problematic from my perspective. Unless this is intentional? If so, why?
   2. When I put this tT history into HeFTy (v 2.3.1), using the Ketcham 1999 annealing model, a Dpar of 1.75, and c-axis projection (and Ketcham 2003 c-axis projection) it predicts an AFT age of 19.8 Ma, not 21.3 Ma. (The predicted track-length distribution (Fig. 1c) is not reported in numbers so I cannot compare that to HeFTy’s prediction.) Section 2.2 describes that Ketcham’s 2005 approach (i.e., HeFTy) is being used for the AFT system here, but was this new code benchmarked against HeFTy? If I’m not setting up the AFT correctly in HeFTy to mirror what is being done here, then more information is needed. Such a difference in ages may not be geologically important, but from a numerical calculation perspective it is essential to know why this discrepancy arises. Starting the history from within the AFT PRZ makes this discrepancy additionally difficult to diagnose.
3. Simple “uplift”
   1. This scenario is repeatedly described as “uplift,” but thermochronologic ages never directly document uplift. They only document rock cooling, which in this scenario is driven by rock exhumation that happens when the imposed rock uplift occurs without any corresponding surface uplift (in other words, rock exhumation). Thus, this is a “simple exhumation” scenario, and should be described as such, starting in the abstract. I suggest never using the word “uplift” without specifying “surface” or “rock”, to avoid this confusion. See England and Molnar, 1990 (Surface uplift, uplift of rocks, and exhumation of rocks: Geology, v. 18, p. 1173–1177)
   2. The model design description at lines 168-170 states: “At 10 Myr, the bottom of the box is pushed upwards at a rate of 1 mm/yr, allowing material to flow out the top of the box while the temperature structure of the box remains constant,” If I understand correctly, in the geodynamic model used to generate the tT paths modeled here, a simple linear geothermal gradient was held constant despite a massive transient change in rock exhumation rate, from 0 km/Myr to 1 km/Myr and then back to 0 km/Myr. By imposing rock uplift and exhumation without any corresponding change in the geothermal gradient, this example is an extraordinarily unrealistic departure from how we know the Earth works. Especially considering the imposed rapid exhumation rate: the Peclet number for this scenario [Pe = (exhumation rate * lithosphere thickness)/thermal diffusivity, which describes the competition between heat advection and conduction] is close to 1 (assuming a standard thermal diffusivity of 25 km2/Myr), and so this is nearly an advection-dominated system that would not even be well represented by using a steady-state solution of the conductive geotherm for a 1 km/Myr exhumation rate. These considerations and their consequences for cooling age patterns have been quantified and evaluated by the thermochronologic community for decades (see for example the “Quantitative Thermochronology” book by Braun et al published in 2006). Given (my admittedly under-informed sense of) the abilities of geodynamic models, I’m surprised at this oversimplification, and surprised that the thermal limitations of this geodynamic result (and the corresponding age predictions) are not acknowledged. The main point of this manuscript is to demonstrate how the new tool addresses some of the major challenges of transferring particle-path tT predictions from geodynamical tools to thermal history analysis tools. And there will necessarily be some simplifications. However, given that the audience of this Technical Note spans both thermochronolgists and geodynamicists, I think it is essential that the examples use simplified, but still physically reasonable, thermokinematic scenarios.
   3. Is 10 Myr sufficient for these thermochronometers to achieve equilibrium at partial retention/annealing temperatures for each system? In other words, how arbitrary is this choice of duration, and does that matter for this demonstration?
   4. The text should be more clear about reporting ‘model time’ (Myr elapsed since the start of a simulation) vs. ‘geological time’ (Ma, what the cooling ages are documenting and how we think about geodynamic histories); see also Major Comment #1 above. The use of Ma vs Myr is consistent and correct, but the practical difference between these concepts is not explained. For example, Figure 5 is mixing these two opposite concepts. My suggestion: Unless it is essential for the reader to see the ‘model time’ framing of the design of these numerical simulations, I suggest the authors simply convert everything into ‘geological time’ for us. This matters because thermochronologists always start examining data (whether synthetic or real) by comparing cooling ages to intervals of geological time that are of interest. In this simple example, I’m looking to see how the predicted cooling ages at the end of the model run (“20 Myr – After Quiescence” in Fig. 5) document (or not) the imposed exhumation event in order to assess the rigor and utility of this approach. But, rather than just looking at Figure 5, I have to do the extra work of converting the geodynamic model time to geological time, or vice versa. (Geologic time = How long ago did the exhumation event start and end?) I also have to know, as a reader, that this conversion is necessary. 
4. Complex model, Figure 6
   1. For the “AHe Ages” plots, why is the result shown for 150 km thickness? It seems like He ages are only calculated for the top 50 km, and even that is an order of magnitude greater lithospheric thickness than these ages are sensitive to. The relatively near-surface details (between the surface and the effective closure of each system) is what is important for assessing the utility of this tool, but that cannot be seen here.
   2. For the “Surface Ages” plots, why is age increasing down the y-axis? This is the opposite of how such data are usually visualized, including in similar tools (see for example, Fig. 6 from this paper by McQuarrie et al: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2018TC005340)
5. Description and simplification of the thermochronometer kinetics models
   1. Sections 2.1 and 2.2 would benefit from reporting the nominal closure temperature ranges for each system, given the kinetics being used in this version of GDTchron.
   2. Simplifying the diffusion and annealing kinetics models being implemented in this first version of GDTchron makes sense, and I support the approach taken here as a first step (so long as it is not the last). But, I think this choice warrants more discussion up front, in the methods section of the paper, to help convince thermochronologists of the basic utility of the current GDTchron package despite the limited kinetics, highlight its potential for future adaptation, and describe the consequences of simplifying kinetics in this way for those who are not familiar with chronometer kinetics. For He diffusion kinetics, for example, the 2000 AHe kinetics parameters and 2004 ZHe kinetics parameters being used here are extremely out of date; from the perspective of thermal history analysis, these kinetics models have been superseded by the radiation damage accumulation and annealing models and the old “static” kinetics parameters should never be used to interpret real thermochronologic data, especially for scenarios with certain styles of thermal histories. But of course, the choice of kinetics model only ends up being important in some contexts. So, I suggest that the authors move much of the information presented at the beginning of section 5 (lines 250-260) into section 2.1 and 2.2. The authors could support this choice by discussing how other numerical tools (e.g., Pecube, age2exhume) have navigated this challenge, and in what types of thermal histories / geological scenarios not using the radiation damage kinetics models matters most. 

MINOR COMMENTS

Figure 1b: Can the user adjust the grain size and composition? Can different kinetics be chosen?