the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 28 Aug 2025
Modelling diffusion, decay and ingrowth of U–Pb isotopes in zircon
Abstract. Understanding the thermal evolution of geological terranes provides essential insights into tectonic processes, crustal evolution, and mineral resource formation. Zircon U–Pb geochronology is widely used to date geological events, yet these dates are altered by a wide-range of processes, including diffusion of radiogenic isotopes at high (>800 °C) temperatures. This study utilises the Underworld3 numerical code to couple diffusion processes with radioactive decay and ingrowth in two-dimensions. We assess the numerical solutions against a series of benchmarks to test the implemetation, and apply the models to examine lead-loss due to thermal events and complexities that arise from multiple zircon growth episodes. Our approach bridges analytical U-Pb isotope measurements with a diffusion-decay-ingrowth numerical model, providing insights into how the thermal evolution of a region alters zircon U–Pb isotope ratios. We apply the methodology to the Trivandrum block in southern India—a region characterised by a prolonged high-temperature event—comparing multiple temperature–time paths with analytical U–Pb isotope data to provide constraints on the thermal evolution of the region. The modelling framework can be easily modified to investigate diffusion-decay-ingrowth across various minerals and isotopic systems, providing a tool to decipher the thermal history of a region recorded in isotopic data.
- Preprint
(6459 KB) - Metadata XML
-
Supplement
(16412 KB) - BibTeX
- EndNote
Status: final response (author comments only)
-
CEC1: 'Comment on egusphere-2025-2278', Juan Antonio Añel, 22 Sep 2025
-
CC1: 'Reply on CEC1', Ben Steven Knight, 22 Sep 2025
Hi Juan A. Añel,
Sorry about that. When I submitted the manuscript underworld3 hadn't had a proper release yet. It since has, and the details are as follows:
The paper:Moresi, L., Mansour, J., Giordani, J., Knepley, M., Knight, B., Graciosa, J.C., Gollapalli, T., Lu, N., Beucher, R., 2025. Underworld3: Mathematically Self-Describing Modelling in Python for Desktop, HPC and Cloud. JOSS 10, 7831. https://doi.org/10.21105/joss.07831
The software version on zenodo:
Moresi, L., Mansour, J., Giordani, J., Knepley, M., Knight, B., Graciosa, J.C., Gollapalli, T., Lu, N., Beucher, R., 2025. Underworld3: Mathematically Self-Describing Modelling in Python for Desktop, HPC and Cloud. https://doi.org/10.5281/zenodo.16838572
Thanks,Ben
Citation: https://doi.org/10.5194/egusphere-2025-2278-CC1 -
CC2: 'Reply on CEC1', Ben Steven Knight, 25 Sep 2025
Hi Juan Antonio Añel,
Apologies, please find an updated code availability statement. I have updated the ms accordingly, which will be included in the future re-submission.
Thanks,
Ben
Code availability:
underworld3 v0.99.1 is available from Moresi et al., 2025, licensed under LGPL Version 3. Scripts to replicate the research are available from Knight (2025), licensed under Creative Commons Attribution 4.0 International. Figures were made with Matplotlib v3.10.1 (Matplotlib, 2025), available under the Matplotlib license at https://matplotlib.org/.
References:
Ben Knight: bknight1/diffusionDecayIngrowth: GMD submission (GMD_submission). Zenodo, 2025. https://doi.org/10.5281/zenodo.16524625
Moresi, L., Mansour, J., Giordani, J., Knepley, M., Knight, B., Graciosa, J.C., Gollapalli, T., Lu, N., Beucher, R. Underworld3: Mathematically Self-Describing Modelling in Python for Desktop, HPC and Cloud. Zenodo, 2025. https://doi.org/10.5281/zenodo.16838572
The Matplotlib Development Team: Matplotlib: Visualization with Python (v3.10.1). Zenodo, 2025. https://doi.org/10.5281/zenodo.14940554
Citation: https://doi.org/10.5194/egusphere-2025-2278-CC2
-
CC1: 'Reply on CEC1', Ben Steven Knight, 22 Sep 2025
-
RC1: 'Comment on egusphere-2025-2278', Anonymous Referee #1, 03 Oct 2025
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.
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?
- D is defined as chemical diffusivity and as thermal diffusivity, which is also defined as k.
- 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
- 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
- Fig4 – the colors chosen for the different models are almost identical and hard to distinguish
- 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.
- 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.
- alphaSiO2 is not defined
Citation: https://doi.org/10.5194/egusphere-2025-2278-RC1 -
RC2: 'Comment on egusphere-2025-2278', Anonymous Referee #2, 12 Oct 2025
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.
l:19: “ingrowth of lead” -> “ingrowth of radiogenic lead”
l:21: please add “e.g.” in these very general references.
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.
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)
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)
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.
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.
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.
l:86:87: you could combine these in 1 sentence.
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.
Eq. (4) You missed a minus sign before lambda
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.
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.
Equations 5 to 7: Please cite, or derive these equations.
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).
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.
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).
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.
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.
L 151: Equation (9) is not an error function. Please simplify, i.e. “The analytical solution of Eq. (9)…”
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.
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.
L 162: “Isotropic”: what is this supposed to mean. This is not a title.
The choice of colors in fig. 4 is very bad. I can hardly see the differences.
L 198: diffusion of what is slow
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.
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.
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)
Citation: https://doi.org/10.5194/egusphere-2025-2278-RC2
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 1,620 | 63 | 19 | 1,702 | 29 | 13 | 10 |
- HTML: 1,620
- PDF: 63
- XML: 19
- Total: 1,702
- Supplement: 29
- BibTeX: 13
- EndNote: 10
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Dear authors,
Unfortunately, after checking your manuscript, it has come to our attention that it does not comply with our "Code and Data Policy".
https://www.geoscientific-model-development.net/policies/code_and_data_policy.html
You have archived the underworld3 code on GitHub. However, GitHub is not a suitable repository for scientific publication. GitHub itself instructs authors to use other long-term archival and publishing alternatives, such as Zenodo. Therefore, the current situation with your manuscript is irregular. Please, publish your code in one of the appropriate repositories and reply to this comment with the relevant information (link and a permanent identifier for it (e.g. DOI)) as soon as possible, as we can not accept manuscripts in Discussions that do not comply with our policy.
Additionally, you mention that you have used Matplotlib for part of your work. However, you do not provide the information on the version of the library that you have used, which is important to ensure the replicability of your work. Please, do it.
Therefore, you must reply to this comment with a modified 'Code and Data Availability' section that addresses and solves the mentioned issues, and include it in a potentially reviewed manuscript, containing the information of the new repositories.
I must note that if you do not fix this problem, we cannot accept your manuscript for publication in our journal.
Juan A. Añel
Geosci. Model Dev. Executive Editor