Authors have developed a numerical code to model the effects of thermal histories on U-Pb zircon ages in 2D. Their code is done within the framework of Underworld3 modelling package that is very versatile and used for a wide range of problems on geodynamics. The authors basically show several benchmarks of simple cases using analytical solutions to diffusion equation or finite differences in 1D. All the used equations are well established and the methodologies also. The authors apply the code to a case study of rocks from Trivandum block.

I have found that the paper is in general clear and well organized. It is more a technical contribution that a research paper. Unfortunately, the code used for the modelling has to be run within the environment of Underworld3 which is very specialized and difficult to run for non-experts. It is also confusing to give the same name to the specific code the authors have developed as the Underworld3 modeling environment which is much of wider applications. Thus I was not able to test the code and if it is properly running. Moreover, I wonder really how many researchers will be able to use the code given the complexity of running it. This is a pity because it significantly decreases the impact of the contribution.

• Thanks for pointing this out. We have developed a binder image for others to run the models on the web, but may be slow due to limitations of running models in the cloud. We have also added instructions, based on the underworld3 installation guide, on how to install and start running models quickly on your own laptop. We have also developed a package, UWDiffusion, to make setting up and running the models much easier.

Aside from this, there are quite a few details that need to be fixed before publication. Many figures should be improved for clarity. The same symbol is used for different variables and there are some symbols that are apparently not well define. Please find below some examples:

• Line 125- “one-dimensional heat pipe” what is the reference for such solution?

We have updated the reference for the solution in the manuscript, thanks for pointing this out.

• D is defined as chemical diffusivity and as thermal diffusivity, which is also defined as k.

Thanks for pointing this out, we have switched k to D to be consistent throughout (L162 & L168).

• In some modeling results no units are given for the variables (e.g space) so that the results seem normalized but in others (and in the plots) units are explicitly given

Thanks for pointing this out. We have updated the figures to make sure they are labelled correctly.

• Fig3- because the differences are small, they are hard to see-it would be better to also have a plot of the difference between the models rather than only the models

We have elected to modify the plot to show the error across the domain for two selected models. The location of the error is the same, but errors are higher at the interface when the cell size is higher/resolution is lower.

• Fig4 – the colors chosen for the different models are almost identical and hard to distinguish

The figure is very poor, thanks for point this out. We have changed the colours to be more clear.

• I understand that the code can also do anisotropic diffusion but this is not applicable to the case being discussed and can lead to confusion of the U-Pb systematics for the reader. I would remove this part and just mention that the code can also deal with diffusion anisotropy.

Thanks for this comment, we have removed the anisotropic model as it is confusing.

• Diagrams of Fig 5 are too small to the difference between the different cases
• Fig 6.Calling profiles to paths is confusing with diffusion profiles.

Thanks for pointing this out, it is confusing. We have changed from profiles to paths to distinguish between T-t paths and profiles across the model space.

• alphaSiO2 is not defined

Thanks for pointing this out. We have updated the text to make it clearer alphaSiO2 and alphaTiO2 are the activity of each element in the Ti-in-zircon calculation.

Dear editor, I have read the manuscript entitled “Modelling diffusion, decay and ingrowth of U-Pb isotopes in zircon” by Knight and Clark. The authors show a newly developed tool for the decay and growth of isotopic elements in zircon. Such tools are welcome and useful to our community. However, I found that the documentation of the methodology was not optimal. For example, in many places, initial and boundary conditions are neglected from the problem configuration. Thus, this manuscript reads well for potentiual users but is not detailed enough for model developers that need to reproduce the results. I strongly suggest that the authors revise their manuscript, and I have provided detailed comments below.

l:17: I suppose that researchers used zircon for the same property much before Rubatto, (2017). Please add some basic reference or add (e.g.) in the citation.

• Thanks for this comment, we have updated the text accordingly.

l:19: “ingrowth of lead” -> “ingrowth of radiogenic lead”

• we have changed the text, thanks for pointing this out

l:21: please add “e.g.” in these very general references.

• we have changed the text, thanks for pointing this out

l:25: “Temperatures that exceed a mineral’s closure..” You should mention before: i) what is the closure temperature, ii) that this interpretation assumes the diffusion within a pristine (undamaged) zircon.

• Thanks for pointing this out, we have added in the definition of what the closure temperature is and that this can be decreased due to radiation damage.

l:53-58 (and equations 1-2): From what I read before, there are two isotopes for Uranium (and Lead) that you want to consider. These isotopes (for U) have different decay constants. Thus, equations (1-2) is strictly not accurate since we do not know if this is for one isotope, or the sum of them. One way around is to add an index if you refer to any of the two isotopes, or otherwise, mention that C is the total concentration (I do not think that you need the total, but just to be clear)

• That is correct, we do consider the different isotope systems with differing decay constants. We have presented eq. 1 and eq. 2 in the generic form (for the element) and had hoped this would be enough to distinguish the different decay constants for the different isotope systems outlined in Table 2. We understand this is confusing and have changed the text and formulas to reference parent and daughter isotopes rather than the elements to make it clearer. Thanks for pointing this out.

l 58-61: since these symbols are introduced for the first time. Please mention their units (or if they are somehow normalized in the code)

• Thanks for pointing this out, we have added in the units for each variable/symbol

l:73: please change this notation. Since you used already partial derivatives before, please use a partial derivative or the dot product of an outward vector with the gradient of C. Also please define n (outward vector) in any case.

 Thanks for pointing this out, we have re-worded the sentence

Eq. (3): The units in equation 3 are wrong. The activation energy and the gas constant must have both J or kJ otherwise the result inside the exponential term is not consistent. I suppose that it is programmed correctly, and this is just an oversight.

• Thanks for pointing this out, we have corrected the units in equation 3. It is implemented correctly in the code but incorrectly stated in eq. 3.

Update formula in l. 84: The update formula is a bit of a mixture of notation. On one hand you are using delta t (for the discretized time increment). On the other hand, you use the laplacian of concentration. This does not show how you do the discretization of the concentration in space. Since you are using the FEM approach you should mention:

• What kind of shape functions you use?
• Is your grid rectangular or it uses triangles?

then, the discretized form could be shown using operators.

• Thanks for pointing this out, the mixture of notation was confusing. We have updated the text to use consistent notation. We have also updated the text to highlight we use a simplex mesh with a polynomial degree of 1 (linear shape functions) (L95-97).

l:86:87: you could combine these in 1 sentence.

• We have updated the equation but left each variable on separate lines. 

l:95-97: mentioning half-life here has nothing to do with the equations below. Please remove it. The decay constant is enough to mention that you have a decay. However, if you would like to discuss half-lives anyway, please mention it in a separate sentence after the ODE.

• We have modified the sentence to make it clearer, removing a lot of the text. Thanks for pointing this out.

Eq. (4) You missed a minus sign before lambda

• Thanks for pointing this out, we have added it in

Eq. (4) Are you assuming that the system is closed with respect to diffusion, I suppose so in this example. In this case you should use total derivatives.

• Yes, eq. 4 and 5 are just showing the change in isotope over time due to decay and ingrowth.

l:105: “both required” -> “both are commonly used”. In fact, only one is required from a mathematical sense. It is just that we use two for better constraints when we measure it.

 Thanks for pointing this out, we have updated the text accordingly.

Equations 5 to 7: Please cite, or derive these equations.

 The calculations are written out explicitly in the text, and are derived from the half-life of each isotope.

Eq. (8) Please mention that this comes from dividing eq. 6 to eq. 5.

l:113-120: please add an appropriate reference (or more) for the concordia plot interpretations (Wetherill etc).

 We have updated the text to reflect we use the TW plot and added in references. Thanks for pointing this out.

Figure 2 and lines 120-123: It is not clear how you plotted this diagram. I.e. are you solving just the ODEs (Eqs. 5-7), or is this the result of the FEM 2d calculation? If figure 2 is the result of solving just the ODEs then this figure is very trivial and has no purpose. If, on the other hand, you want to demonstrate that at the limit of negligible diffusion (i.e. when D is set to 0) your code reproduces this figure; then this should be mentioned together with initial, boundary conditions etc.

• Thanks for pointing this out, we have updated the figure and included how they are calculated, preferring to not include the decay-ingrowth in the equation and instead calculate the amount numerically over a timestep. We test the decay and ingrowth equations as source terms but the accuracy was very dependent on the timestep.

l:126-128 (and eq. 9). Firstly, equation 9 is one of the most classic analytical forms of diffusion solutions. This equation is definitely not rooted in Zhang et al. (2022). Any classic textbook (e.g. Crank 1975 – The mathematics of diffusion – page 15) has this solution. Secondly, this is actually the solution of an infinite slab, therefore I would not call it a pipe (pipes can also be cylindrical etc).

• Thanks for pointing this out, we have updated the text to remove the ‘heat pipe benchmark’ and referred to it as ‘rectangular pulse benchmark’ now instead.

L: 136: There is no k in this solution. Only diffusivity. Also, you mention U_a is a concentration but has not defined it before. Furthermore, you can start from the distribution of a square pulse at t=0. The initial form of this analytical solution, although discontinuous, is well known and will in any case be smoothed by your spatial discretization.

• Sorry about that, we have updated the text to change k to D to be consistent throughout that D represents diffusivity.

l:141-142: You have not mentioned the boundary conditions of the numerical calculation. In reality the analytical solution has boundary conditions at infinity, therefore, I guess it should be ok to assume either BC only for very limited timescales. These details need to be mentioned (also the total size of your model) for the reproducibility of your results.

• Thanks for pointing this out, we have updated the text to be explicit with the BCs used.

L 151: Equation (9) is not an error function. Please simplify, i.e. “The analytical solution of Eq. (9)…”

• Thanks for pointing that out, we have updated the text accordingly.

L: 53: It does contain boundary conditions (BC); they are just at infinity. You either need to expand your domain or reduce your diffusion timescale and use the results as an approximation. Alternatively, you need to derive another analytical solution that have the same BC as you have.

• Thanks for pointing this out, we have updated the text to be explicit with the type of BCs used.

l 150-160: This benchmark is not needed if you limit the timescale as I mentioned before. Also, this has now conductivity where you need diffusivity. These variables have different units; it’s not just a naming issue. Please fix the code, the text and mention explicitly the initial and the boundary conditions of this benchmark. Furthermore, please mention the CFL condition for this case. To my opinion, this benchmark serves no purpose but if you want to keep it, please document the things I mentioned above to make your model reproducible.

• Thanks for pointing this out, we have updated the code accordingly.

L 162: “Isotropic”: what is this supposed to mean. This is not a title.

• Apologies, it should have been ‘isotropic diffusion’ to represent no changes in diffusivity due to the planes.

The choice of colors in fig. 4 is very bad. I can hardly see the differences.

 Thanks for pointing this out, we have updated the figure accordingly.

L 198: diffusion of what is slow

• We have updated the text to show highlight both elements (U and Pb) diffuse slowly at 750C.

l:231: These are not really pulses since they can be 400Myr long (many orogenies fit into that time). I would call them thermal perturbations.

• Thanks for pointing this out, we have renamed them to 'events'

l:249: What is “temperature profile four”. I am lost please refer to a figure and explain in detail. This whole paragraph was difficult to follow. Also, how you perform the different “growths” is not clear. Are you adding layers with given isotopic ratios? If yes, please mention it explicitly in the model description. Now it is too late.

• Thanks for pointing this out, we have updated the text. Our second growth is ‘instantaneous’, where the U isotopic ratios are calculated at that time, and updated in the outer domain.

l: 330-331: This approach was suggested in the past for diffusion in garnet. However, without accurate estimates of time duration, this method has large errors. In other words, if you do not know in advance the duration of your diffusional reset, you cannot propose that this is a viable method. If, on the other hand, you have independent estimates for the duration of the “pulse”, then yes, I see your point. (thus this comment is relevant for lines 333-335)

• Thanks for pointing this out, timescales derived by diffusion will have errors associated with them due to the nature of obtaining diffusion coefficients. However, they provide an additional constraint that can be used to constrain the thermal evolution of a region. If used in conjunction with diffusion/discordance of other minerals and elements that have different closure temperatures, then the errors can be reduced by combining independent observations. They can also be further reduced if used with other, independent constraints.