gchron-2022-23 review: Fox et al.

The manuscript “Origin of Great Unconformity obscured by thermochronometric uncertainty” by Fox et al. presents numerical models of uncertainty of He diffusion kinetic parameters from an undamaged zircon and zircon N17 (of Guenthner et al., 2013). Specifically, they consider uncertainties in the frequency factor (D0) and activation energy (Ea) for these two end-member grains, which co-vary, and they use ZRDAAM (Guenthner) with their updated kinetics for these two grains to predict zircon (U-Th)/He (zircon He) thermochronometry dates for an example deep time thermal (tT) history.

This contribution drives at a timely and important topic – the utility, and potential limitations, of zircon He thermochronometry to resolve deep time thermal histories. As the authors point out, zircon He data has been used to evaluate the timing and magnitude of erosion during the creation of the Great Unconformity, with data patterns being used to discriminate between glacial, tectonic, and geodynamic mechanisms for unroofing. However, uncertainty in He diffusion kinetic parameters, which are a function of radiation damage accumulation and annealing, have not previously been considered before and I appreciate the authors tackling this here. As the authors show with their suite of numerical models: the effect of these uncertainties is magnified in long time scales (i.e., deep time) and dispersion in a zircon He date for a given eU value and deep time tT history can be >>> 100 Ma. This suggests that, in the absence of additional diffusion experiments on zircon grains over a spectrum of effective damage, it may not be possible to resolve cooling during say the Sturtian and Marinoan glaciations, or doing a specific Neoproterozoic time period.

I think that this manuscript is suitable for Geochronology pending moderate revisions. There are elements of the manuscript that I think should and can be expanded and clarified prior to acceptance for publication:

First, the authors center their analysis on only two end-member zircon grains. One grain, “z”, is sometimes called a minimally damaged grain and elsewhere called an undamaged grain or a “low” grain. As Guenthner et al. (2013) shows, these are different. What is the source of the data for the low/undamaged grain? This is never stated in the text and should be added. (see comment below about Figure 1).

Second, and more importantly, it is not clear how the authors are taking accounting for uncertainties in intermediate or even high damage zircon grains (N17 is effectively an amorphous grain; Guenthner et al., 2013). The authors state (152-153) that they are focused on end member grains and imply their approach is an improvement on the extrapolations between high and low damage levels in Guenthner et al. (2013) (152-157). Regardless of whether the low damage end member is zero damage or Mud Tank kinetics, what happens in between the low damage and N17 is, I think, critical. What about moderately-damaged to high damage grains: how do the predicted versus observed diffusivities vary for RB140, BR231, M127, and G3? (Guenthner et al., 2013)? Note that Guenthner et al. (2013) presents data for slabs of the crystals cut parallel and orthogonal to the C axis for each grain. Related to this, looking at figure 3, I found myself wondering if D0 and Ea are correlated in the same way in moderate to high damage zircon grains, and do these correlations align with Guenthner et al. (2013), as they do for the low damage “z” grain?

The reason why this additional analysis I think matters is because the evolution of effective damage and resulting He diffusivity is complex as the authors point out (126-128). Diffusivity initially decreases from an undamaged/low damage grain to a moderate damage grain, but diffusivity then dramatically increases from a moderate to an amorphous grain. In figure 4, the greatest spread in dates when taking into account the variation in grain “z” and N17 uncertainties occurs in the range of moderate to high damage grains (1000-1500 ppm eU). If only endmembers are being considered, how are uncertainties in diffusion kinetics between these two endmembers being accounted for or propagated through a tT history? This was not clear to me.

I think the authors can address these questions in two ways: I encourage the authors to perform the same calculations that they perform on an undamaged/low damage grain and N17 on additional grains, mentioned in the paragraph above, and expand figure 2. I think showing how the D0 and Ea values compare with Guenthner et al. (2013) will be useful for the community. It was interesting to see how their models reproduced N17, but not for the un/low damage grain. At the very least, the authors should expand the text to explain how they are interpolating between un/low damage and the amorphous grain when propagating model uncertainties across a range of eU – I apologize, I could not follow how this was done with the information provided.

Third, Figure 1 would benefit from more information. Adding these details will help the reader follow the discussion in the text and also help if additional analyses are added:

-Please label each of the Arrhenius relationships for the different diffusion experiments – perhaps using different symbols and colors. I think these figures could be enlarged as well. I found myself flipping back to the original data of Guenthner et al. (2013) to identify N17 and I could not figure out the low or undamaged grain.

-I would make the trend lines a darker grey or perhaps thicker – it is hard to see them. Enlarging the figures (perhaps single column?) would help.

Additional edits and comments:

12-13: I would rephrase this topic sentence, as it is confusing as written. Perhaps something like “Thermochronology provides a unique perspective on the timing and magnitude of erosion during the generation of unconformities.”

17: I do not think However is needed at the beginning of this sentence.

23: I would add “… as a function of radiation damage accumulation and annealing.” At the end of this sentence.

25: “less well resolved” – during what time interval?

26 and 320-322: “calibration with natural laboratories” – Calibrating zircon He diffusion kinetics from natural samples that experienced deep time thermal histories is challenging, right? Because it is precisely during those long tT histories where issues such as U-Th zonation will be exacerbated. And so if you chose to advocate for an empirical approach, I would add a cautionary note that provided U-Th zonation can be quantified in analyzed grains (which is hard, right? because grains that are consumed for zircon He analyses are not available for detailed LA-ICPMS zonation investigation. And LA-ICPMS zircon He analyses also have uncertainties that would need to be taken into account in this modeling approach).

33: I would omit Keller et al., 2019, as this Great Unconformity study does not include zircon He thermochronometry data, which is the subject of the sentence.

36: I believe this is 1300 to 200 million years of missing time.

38: Here I think you could cite other papers besides McDannell et al., 2022 such as those that you site in 33 plus others, minus Keller et al., 2019.

38-40: Although I do not think you want to wade into this with this paper, I would rephrase this because the modeling approaches described in McDannell et al. (2022) reflect only one approach to this problem as I understand it.

95: “variability in radiation damage annealing parameters” – do you mean “diffusion kinetic parameters”? Although D0 and Ea depend on the effective radiation damage (i.e., accumulated and annealed damage), the work presented here does not evaluate damage annealing specifically as this phrase implies.

210: Useful to call-up the inset in Figure 4.

Figure 4: Please add tick marks for x and y axis values.

229: Grain size = equivalent spherical radius? Or?

244-245: These are incredibly precise numbers – is this level of precision warranted here? And how much does it matter?

238-256: I think this analysis would be useful for grains like RB140, BR231, M127, or G3!

Alexis K. Ault

Department of Geosciences

Utah State University

Please see attached comments

Please see attached PDF.

C. Brenhin Keller

The manuscript by Fox and co-authors re-examines certain aspects of a commonly used damage-diffusivity model for the zircon (U-Th)/He system (ZRDAAM). The authors highlight that the model, as originally derived and defined by Guenthner et al. (2013), currently lacks information on kinetic uncertainty, which could in turn influence the thermal history modeling-based outcomes that rely upon ZRDAAM. In particular, the authors argue that the lack of kinetic uncertainty could create over-interpretation of and/or overconfidence in thermal history model output as applied to a recent discussion surrounding the origin of the Great Unconformity erosion surface. I found this contribution to be of importance for three reasons: 1) it stresses the need to incorporate uncertainties into the kinetic models used by many authors and thermal history modeling software packages, 2) it shows how large zircon He datasets (~40 grains) are needed to refine and extract meaningful time-temperature information, and 3) it recognizes that ZRDAAM is a work in progress that could be strengthened by further He diffusion kinetic work. In sum, I think this manuscript is a solid contribution to ongoing discussions surrounding thermal history modeling and thermochronometric data interpretation. However, given the authors stated goals of exploring calibration issues with ZRDAAM, I think the article does not go far enough in investigating multiple aspects and nuances of the current ZRDAAM calibration. My recommendation is major revisions, which should focus on expanding the scope of the current work. I have several general comments that I would like to see the authors address before this article is published.

First, the authors focus heavily on the theoretically pristine endmember zircon, and the fact that this value is an extrapolation that, correctly, could have a large degree of uncertainty to it. It is difficult in the current preprint to determine how important this endmember is in the context of real datasets. Let me expand. When the ZRDAAM was derived, a proposed “zero-damage” grain was used to make the model setup more straight-forward, and to account for potential scenarios where users may want to model thermal histories at very short time scales. But nearly all grains, especially those used for deep-time problems, will very quickly (within ~5 myr for a grain with 500 ppm eU) accumulate damage levels in excess of the Mud Tank sample, for which empirical results do exist. So perhaps the challenge with the current ZRDAAM is less, what should we do for a zero-damage endmember, but rather, should we instead be using a sample with known kinetics as our endmember? Would the Mud Tank kinetics therefore be more appropriate here as a realistic endmember, and if so, does that reduce some of the variation in thermal history model output? Secondarily, there are both molecular dynamics (Reich et al., 2007; Saadoune et al., 2009) and empirical results from zircon-like, zero-damage orthophosphates (Farley, 2007) that attempt to place an approximation on He diffusion in a pristine zircon. How do the author’s new extrapolations match up to these results as a comparison? Finally, as Ginster et al. (2019) observed, it is very unlikely that a natural zircon would return to a pristine damage state following damage accumulation, so the likelihood of encountering a natural, zero-damage zircon is remote. I would like to see the authors acknowledge some of these points and address these questions. I think the authors have illuminated a valid and important point: that the ZRDAAM as defined by Guenthner et al. 2013 needs to be recalibrated at the low damage end and uncertainty needs to be accounted for. However, as I’ve hopefully conveyed, I would place the focus more upon using an endmember for which we can obtain kinetic information, rather than (as was done in Guenthner et al., 2013 and is done here) an idealized, extrapolated endmember that does not likely exist in nature.

Second, the authors highlight a data set from Minnesota as an example with which to test their newly constrained ZRDAAM. My concern with this modeling setup, as articulated in a line specific comment below, is that an apples-to-apples comparison to the McDannell et al. 2022 study would include many more iterations. In that study, 500,000 burn-in, and 500,000 post were used, whereas here it appears that only 100,000 burn-in and 100,000 post are used. That is, the higher amount of uncertainty as to when Neoproterozoic cooling initiates could be that the model was not run for long enough.

Third, the authors call for additional diffusion experiments and new data to be added to the ZRDAAM, which I certainly agree with. I would point out though that some of these data have already been collected by Ginster and are published in her PhD dissertation, available through multiple repositories, including the University of Arizona library. This work is yet to be published in a peer-reviewed manuscript, but the data essentially agree with the Guenthner et al. (2013) diffusion data, which suggests that the variation (and therefore the uncertainty) in measurable diffusion kinetics is lower than for extrapolated endmembers. The point should be emphasized that, from an experimental perspective, the ZRDAAM remains on solid footing, although perhaps the endmember kinetics need to be recast, as my previous comment suggests.

Line specific comments:

Line 73: Should acknowledge here that alternatives using an arguably more direct measure of alpha dose (Raman spectroscopy) have been proposed (Ginster et al. 2019) and incorporated into ZRDAAM (Guenthner, 2021). I would also elaborate the discussion here to comment on the different styles of damage that may influence He diffusion kinetics (i.e. alpha ejection, alpha recoil, fission track). The diffusion kinetics are calibrated to alpha dose (essentially alpha recoil), but the mode of annealing is debated and as yet not fully resolved. The Ginster et al. (2019) data and the model demonstration by Guenthner (2021) of these data are particularly salient given my comments above concerning an idealized zero-damage endmember.

Line 85: An additional explanation here is that, as McDannell et al. 2022 showed, these models needed to be run for many more path iterations (at least 500,000 pre and post burn-in). In the Thurston et al. (2022) study, we ran models for only 100,000 pre and post burn-in, which was admittedly likely not enough. To this point, we have re-run some of these models with the greater number of paths (currently unpublished) and indeed the earlier portions of the time-temperature history remain under-constrained by zircon He (as discussed in Thurston et al., 2022). However, the late Miocene cooling remains robust.

Line 99: The phrase “varies by hundreds of millions of years” should be qualified here. As you suggest later in the manuscript, the variation is dependent on the number of grain inputs you have and the spread in date-eU space of those inputs. Moreover, my understanding (this could be clearer, see next comment below at line 238) is that this is for scenarios that incorporate the 2 sigma from the full kinetic distribution.

Line 238: The focus in this paragraph on fixed endmembers seems out of place with one theme of this manuscript: kinetic uncertainties should be sampled in the rjMCMC approach. Some of this could be my confusion: am I correct that the models were run with Ea and D0 values that represent the 2 sigma of the kinetic distribution? If I’m not correct, then please more thoroughly explain how or why these specific kinetics were selected. If I am correct, then I understand that the authors are perhaps trying to show the worst-case, 2 sigma extremes from their distribution, but why not incorporate the full probability distribution as shown in figure 2 and sample that? The authors mention further below that computation limitations prohibit this exercise, but much of the discussion and implications seem to be cast in light of the highest possible amount of variation. If the modeling incorporates the full distribution (and samples it) is the situation as dire? I am really curious to see the outcome of a modified MCMC proposal algorithm with a selection statement that samples the kinetic distribution.

Line 249: For a better apples-to-apples comparison here with the McDannell et al. (2022) study, 500k pre and 500k post burn-in is needed. As we have seen (and learned, see the comment about Thurston up above) these deep-time problems require at a minimum 1,000,000 path iterations.

Line 294: Are the authors suggesting here that the binning and averaging approach has limitations? The point I’m most struck by, as the authors allude to, is that binning and averaging can eliminate the sensitivity of the whole data set by removing portions of date-eU space.

Willy Guenthner

UIUC