This paper presents a comprehensive statistical meta-analysis of CA-ID-TIMS zircon age data, identifying differences in the skewness of age distributions across different magmatic environments. The analysis suggests that volcanic and porphyry age distributions often show left-skewed patterns, whereas plutonic distributions tend to be right-skewed. While there is considerable overlap between the age distributions, the overall trend appears consistent and aligns with intuition based on cooling rates. The empirical confirmation of such expectations is valuable.

The paper’s strengths lie in its rigorous statistical approach, including quantitative assessments of potential biases and uncertainties. Although not new, the effect of subsampling of zircons is studied, which is useful in a practical sense. Most of the relevant literature is cited, the methods are clearly explained, and figures are well prepared making it straightforward to follow the manuscript.

The authors extend their analysis by proposing an explanation for the observed skewness differences, using a thermodynamic model combined with zircon saturation and bootstrapping techniques. Based on this the authors suggest that eruptive sampling of a monotonously cooling magma body is an important process in controlling natural zircon age distributions, with implications for using these distributions to infer magmatic fluxes from thermal modelling. However, I find this aspect of the study problematic for several reasons:

• It is well established based on thermal modelling that simple, monotonous cooling of a magma body will not produce zircon crystallization durations like those shown in your compilation. Most plutons are assembled incrementally, with even small, shallow laccoliths exhibiting zircon age distributions that are inconsistent with simple, homogeneous cooling at the emplacement level. Thus, attributing a significant role to eruptive subsampling in a monotonously cooling magma reservoir has limited relevance for explaining the observed zircon age distributions.
• This issue is even more problematic when considering that heat transfer modelling shows that monotonous cooling of an intrusion does not necessarily produce a right skewed zircon age distribution in the first place. Using thermodynamic modelling and bootstrapping to obtain synthetic zircon ages, you assume that the thermal evolution can be described by a single time-temperature curve. If the heat equation is used instead to model cooling of an intrusion in 3D, we’ll find that the magma spends the longest time at temperatures with smallest temperature gradient to the surroundings. This shifts the synthetic zircon distribution to the left depending on the conditions. When it comes to calculating synthetic zircon age distributions, the spatial and temporal dynamics of heat transfer are hugely important. Given that heat transfer is ignored in your modelling, I don’t think that the presented cooling, recharge and subsampling scenarios are particularly insightful.
• The study’s words of caution against using zircon age distributions to infer magma fluxes via thermal modelling in volcanic settings may mislead readers. Specialists in this method may find this feedback irrelevant as heat transfer is not considered, while non-experts might be encouraged to discredit the results of these studies. Rather than discouraging the approach, the study should acknowledge the limitations of its current modelling and emphasize that magma flux quantification requires thermal models accounting for complex magmatic processes and 3D heat transfer.

I do not suggest re-running the entire study with thermal modelling but rather recommend a more detailed discussion of the assumptions and limitations inherent in the thermodynamic modelling approach used here. Additionally, based on the comments above,  revisiting the conclusions on the relevance of the proposed truncation process and on magma flux quantification, as well as exploring other potential causes for skewness differences—such as emplacement depth—would enhance this paper. For instance, could left-skewed plutonic distributions correspond to shallower intrusions with rapid cooling rates?

I would be happy to support publication of this article if the authors address the issues raised above.

Other comments:

L12-14  At this stage I’ve got lost without reading the rest of the paper. Could you already here briefly explain how you go from bootstrapping to inferring the processes controlling the difference between volcanic/porphyry zircons on one side and plutonic zircons on the other.

L14-15  Here as well, how do you find that higher fluxes can contribute to the skew?

L34-36 I don’t think we can say this with any certainty at this stage. Yes, the early work (Caricchi et al. 2014,2016) seemed to suggest that there are higher fluxes at work in large eruptive systems compared to plutons, but it turned out that the originally proposed methodology, based on similarity of zircon age distributions (the topic of this work), is too much impacted by statistical biases to reconstruct magma fluxes. That is why other authors used a different sort of approach based the age span, temperature distributions and geological constraints to recover fluxes. A second complication in all this is that different authors use different measures of fluxes such as volume fluxes km3/yr and others area normalized fluxes km3/km2/yr. If we consider Toba as an example (Liu et al. 2021), reconstructed volumetric fluxes are one of the highest ever quantified, while the area normalized flux is less than a typical stratovolcano. This makes a huge difference and is also not considered in the early papers that claim systematically different fluxes for plutons and caldera forming systems.

L47-49 I have a somewhat strong opinion on this but to me this assumption does not make sense at all. Grounding the comparison of zircon age distributions on autocryst-antecryst differences is rather misleading. There is plenty of evidence that magmatic systems are built incrementally by repeated injection of magma batches. Thermally it is quite difficult to produce any long ranging zircon age populations without incremental assembly, given that large single reservoirs would lose heat too quickly. A single magmatic system can therefore contain multiple subsystems that may or may not show some chemical differences in zircon chemistry given variations in recharge magma compositions, differences in partitioning of trace elements in different parts of a reservoir, depth-dependent differences in crustal melting efficiency, and many other factors that would lead to distinct differences in the chemistry of zircons that then could be interpreted as antecrysts. This does not mean that we can’t make sense of the age distributions or have to exclude these crystals. None of this is a real problem if the goal is to reconstruct magmatic fluxes or to understand the general dynamics of magmatic systems.

L78 Antecrysts cannot be more easily identified with this approach but xenocrysts can.

L84-88 The term “magmatic event” is not right here given that you consider the entire evolution or at least large parts of the history of a magmatic system e.g. over 500 ka in the example you give above.

L85-89 See discussion about antecrysts above.

L92 “making the autocryst–antecryst divide ambiguous, and possibly detrimental to the understanding of the underlying system”. Okay, couldn’t have said it better.

 L106 “… ECDF older than the two youngest dates (to ignore age gaps at the young end of the distribution which are not considered here):”. That is fine but I wonder why only the two youngest dates? Don’t you have any gaps that could be beyond that threshold in the young part of the distribution?

L150-153 A minimum of 10 zircons to resolve the shape of the distribution seems rather optimistic given the bootstrapping you show later in Fig. 7.

L253-254 “overall divide” is somewhat of an overstatement given that most of the volcanic and plutonic distributions overlap. I agree that there is a tendency for the volcanics and porphyry to be left skewed and right skewed for the plutonics.

L359-374 1) Given that there are only a few, if these are the part of second older age population could be tested by simply looking at the relevant datasets. 2) Your bootstrapping plot (Fig. 7) shows quite some variability also for the n=30 case. My suspicion is that you could also explain the extreme cases if you slightly change the shape of the theoretical sampling distribution and/or just use more iterations.

L338-444 FYI. The effect of the specific zircon saturation curve in use is quite large when calculating synthetic zircon distributions based on thermodynamic modelling. However, when using heat transfer the effect of the temperature distribution through time is much more important. In this case, switching the zircon saturation model around has only a minor effect on what the synthetic zircon distribution looks like, so this is secondary at best.

References:

Liu, P. P., Caricchi, L., Chung, S. L., Li, X. H., Li, Q. L., Zhou, M. F., ... & Simpson, G. (2021). Growth and thermal maturation of the Toba magma reservoir

General Comments:

This manuscript addresses a critical issue in modern igneous petrology:  How best to understand the measured timescales of crystallization in igneous rocks.  The lens here is an important subset of this problem, zircon crystallization timescales, when measured using whole grain analyses at precisions that are much higher than the crystallization timescales.

In this work, the authors produce a large number of carefully curated U-Pb datasets, present a quantitative model for comparing their distributions (Wasserstein difference) and reducing the degrees of freedom in them (PCA), extract trends in the data that correlate broadly with the type of igneous system, and then compare these trends with predictions they make based on simple forward models.  I found it to be a very clearly written paper, with good explanations of the quantitative techniques and explanations of how they related to the underlying datasets.  The figures are high quality, the writing is excellent, and the choice of topic is well within the remit of Geochronology.

Specific Comments:

Like I said above, this is a very good paper.  The only real issue is that the authors don’t address that there are important differences between their forward models and real processes that drive crystallization and preservation of crystallization age distributions.  As I’m writing this I note that this is something that the first reviewer also pointed out, so I won’t spend a lot of time on it, but I think this is probably worth addressing qualitatively, for completeness.  Both in terms of what the T-t history is like, and additional complications such as those highlighted by Klein and Eddy, and biases in the measured zircon record.  I appreciate that it may not be possible to include those complications at this time in a quantitative model, but they are both real and important.  

Specific comments/technical corrections

L43-44: It’s worth noting that sampling biases may also play a role.  I can’t speak for every analyst, of course, but typically the largest and highest quality crystals are selected from a population.  The size bias alone may be significant, as is the fact that whole grains modified by CA bias U-Pb dates by U concentration, dissolve components that are accessible to HF during the partial dissolution phase of CA, and integrate the remaining grain.

L120:  This is an excellent summary of the method.

L134-137:  I think it’s fair to stick to one type of dataset but the objections to 230/238 geochronology aren’t strong.  Variable and asymmetric uncertainties aren’t difficult to deal with, and initial [230/238] are determined in many high quality datasets (and doesn’t vary infinitely).  There are good reasons that plutonic datasets are rare (or absent?) and that’s because the technique is exclusively applicable to young rocks, and that’s a reason enough alone to exclude it.

L159:  This is very good.

L201-206:  I’m a bit confused by this paragraph and unsure as to what I’m supposed to take away from it as a reader.  It describes a technique that is not used, justifies why it is suboptimal, but then in the last sentence tells me that the results are the same?  I have no specific objection here but I think the authors might want to clarify for a reader what they should take away from this.  For me personally, as a geochemist who has used MDS in the past, I’d be just as happy if this paragraph didn’t exist?  I don’t think many readers will look at this and wonder why they didn’t use MDS.

L220-221:  Minor point but mixing and matching confidence/coverage intervals is confusing to a reader (the parameter fit uncertainties are at a coverage factor of 1 and the result is at 2).

L221:  This is a really useful result, and I could see it being something used by others in future work – it might be worth adding a bit more context for those who do?  Can you specify the time range that this is applicable to, and if the fit parameters uncertainties are correlated? In the introduction you list 0.1-4567 Ma, but that’s probably not relevant here (if I do the calculation right it’s ~40% at 0.1 Ma.)  It’s also worth noting that these uncertainties differ markedly from those listed rather optimistically on line 138 (e.g., if I did the calculation correctly, it bottoms out at 0.11% in the Neoproterozoic).  I’d recommend revisiting line 138 to bring it in line with the actual uncertainties parameterized in this equation – that’s more of a boilerplate ID-TIMS boast than a realistic assessment of actual typical precisions from young zircon populations, where for even low-blank labs, the precision ends up being limited by the blank isotopic composition variability.

L235-240:  It’s worth revisiting the issues raised by Klein and Eddy at this point.  Not that they need to be incorporated but just as a caveat for the reader.

L268:  Please be cautious here in assuming that natural processes won’t result in gaussian distributions.  The central limit theorem applies to natural processes as well – the samples we recover average a range of disparate pieces of data and proceses, and averaged samples from any distribution are gaussian.

L305:  Can you compare this to other estimates if you are aware of any?  (This isn’t a leading question, I’m can’t find one myself)

L412:  This seems like an unnecessary swipe at LA-ICPMS labs.  Using data that is accurate (including accurate uncertainties) is obviously critical and there’s no reason to cast aspersions on a specific technique.  Peak hopping TIMS Pb data on Phoenix instruments have uncertainties on isotope ratios underestimated by ~1.7x due to autocorrelation (and some other factor on any other TIMS), but that fact is equally out of place here.  LA-ICPMS is an easy target because of the large number of practitioners and therefore the larger n of outliers.