We thank the reviewer for their time and comments on the manuscript.

The reviewer has raised they would like to see more discussion regarding the potential limitations and influential factors in down-hole fractionation (DHF). While we agree the with the overall concept they are presenting in their commentary, we feel that this is largely beyond the current the scope of this contribution, which is to provide an improved quantitative method for the assessment of DHF fraction during LA-ICP-MS analysis, rather than conduct the assessment of these factors. The assessment of these factors using the proposed method is a worthwhile endeavour, but would be best addressed in a one or more separate studies. The factors influencing DHF are well documented in literature, although their underlying mechanisms are somewhat ambiguous still as the reviewer states.

We agree it would be worthwhile adding a brief section in the discussion or introduction that outlines the potential factors more clearly and how they are controlled for during this study and what impact they would have if neglected, but the mathematical foundation of the modelling is independent of these factors. As with all tools, user beware.

Regarding commentary about single laboratory/single ablation system, again, given the method being proposed here is the mathematical foundation for an improved model fitting algorithm, we believe addition of extra datasets from other laboratories and equipment is not necessary for this particular manuscript. It is exactly what would be needed for a large-scale assessment of DHF by the community and would be best addressed in a separate study. From the data provided in this study it is already apparent that the usual factors influencing DHF (spot geometry, material, etc) will change the resulting model fit. Laser wavelength is known to impact DHF as shown in numerous studies and so would be expected to change the resulting fit, but again not the mathematical foundation of the fitting. We propose to add older published data from the lead author that used a 213 nm laser paired with an Agilent 7500 as a supplementary figure to demonstrate this, and would add the brief section on DHF influences as outlined above.

With regard to the commentary asking how these factors were controlled for during this study. We decided to demonstrate the fitting algorithm using common reference materials for well-established methods, i.e. U-Pb geochronology of zircon, monazite, apatite etc using standard conditions for each material as accepted by the community for a given analytical setup. The references for the standard analytical methods are provided in the methods section of the paper, but for example all zircon data presented in this study was ablated using a 193 nm wavelength laser using a fluence of nominally 2.0 J cm^-2 and a repetition rate of 5 Hz. When computing orthogonal decompositions of aggregated data, the data is first grouped by fluence, spot diameter, repetition rate, and material. This ensures that only like-for-like data are being fit, and not artificially skewing results due to external factors that are known the impact DHF (e.g. fluence) not being accounted for. Again, our paper presents a methodology for quantitatively fitting the signal of LA-ICP-MS analyses to model DHF: it is independent of the specific instrument factors used during analysis.

In regard to cross-calibration due to detection modes, I am interpreting this as the change from analogue to pulse detection model. Our ICP-MS setup is set to use a pulse cut-off of ~2.5 million counts per second, and during initial tuning before an analytical session the count rates are tuned to optimise sensitivity while ensuring that counts remain in one detection mode. If an analysis did have a portion in analogue mode that switched to pulse mode due to the signal intensity decrease over time, and that calibration factor was in error, it would make the fit erroneous and require portions of the signal to be subdivided for fitting, however no such change is observed in this data set (as seen in the supplementary signal intensity plots). It would be seen as a step change in the intensity and ratio of a signal if present.

Regarding the centring of time-dependant data. The reviewer has understood correctly that only lambda 0 is sensitive to ICP tuning dependencies. Centring the data does not change the shape of the data, just the central tendency and this allows for comparison of data that have not yet been corrected against a known ratio (e.g. NIST glass, GJ). Total signal time is likely to have negligible impact on lambda 1, but may have a greater impact on lambdas 2 and higher due to the more complicated profiles these represent (sinuosity for example). This potentially would have some impact on where they plot on UMAP, but they should still plot close to analyses from the same material. Some of the variation seen for a single material (e.g. GJ1) may be explained by this effect. This can be elaborated on in the manuscript and looked at in more detail. The time component of a modelled signal is corrected to laser-fire and signal-stabilisation (after initial intensity ramping) to mitigate some of the differences that would occur from using a delayed subset of a signal (e.g. 10 seconds after the signal stabilised).

For the minor comments:

We agree that geometry is a better term as it accounts for shape, diameter, and depth.

We will correct the missing E in the scientific notation, and will change to "value*10^X" notation if that is the Geochronology formatting recommendation.

We thank Noah for his time and constructive comments on the manuscript.

Starting with the final overview comment about how was the algorithm tested. We did not create a synthetic dataset and recover the input parameters, in hindsight this would have made testing simple, and we are happy to implement this to further validate the method as it would also be useful for debugging purposes. This method was developed upon the established methodology of O'Neill (2016) and Anenburg & Williams (2022) who have used it for modelling REE patterns quantitatively. Our algorithm was tested against theirs using identical datasets to ensure consistency in results and uncertainties. Our implementation of this algorithm is much more generalised allowing it to be used for the purposes in this manuscript, but can readily be applied to any data where a (up to five order) polynomial fit is sensible.

While the reviewer has raised a similar concern to reviewer 1 about the results being applicable to only our analytical setup, we again believe this is not a limitation of the current scope of the study but would be happy to add older published data from the lead author that used a 213 nm laser paired with an Agilent 7500 as a supplementary figure. The scope of the study is to provide a tool that can reliably and quantitatively model DHF patterns to enable studies that assess the differences between the factors that influence DHF and differences between analytical setups, laboratories, etc., rather than make a full assessment of these. The assessment of these factors, and results from different laboratories would best be done in a separate (or several) community-led studies involving many institutes.

Fluences in appendix A were chosen based on the analytical methods cited in the methods section. The chosen fluences for a given material have been used consistently for several years in the host laboratory with numerous studies published using those conditions. The initial studies aimed to optimise signal sensitivity, count rate, signal complexity, and pit depth while minimising DHF for a given material. As the reviewer has rightly stated, this will vary between laboratories but this is where the proposed method can be of significant use in enabling quantitative comparison between setups so long as the appropriate metadata (laser and tuning conditions) are known for a given fit. We can add a sentence or two to the manuscript highlighting this point.

For the reviewer's overall commentary on the clarity, we will use their input to refine the manuscript to improve this aspect. One of the co-authors who is more removed from the mathematical and technical knowledge on this manuscript will look over the manuscript again with the reviewer's feedback in mind to assist further with this.

Regarding the discussion of applications of the method. We will aim to clarify this further in a revision. Our intent is to provide the method so it could be used as a drop-in replacement for standard polynomial fitting used in current data reduction software (Iolite, LADR) for DHF correction whilst then providing the quantitative parameters and enable DHF correction uncertainty propagation. For our purposes, it is used to demonstrate the methods validity for down hole modelling, and we will use it to quantitatively assess the appropriateness of using one material to perform down-hole correction on another.

The proposed method, could however be used to quantitatively model any geological process which has a (mathematically) linear relationship between x and y (for example, REE abundance ionic radii, as has been done previously).

With specific regard to line 293, while we have not investigated how such an implementation would work in practice, we believe it would be possible to remove the dependence of the "down-hole method" (Paton el al.) on a homogeneous material for DHF correction by fitting each unknown directly. There are numerous parameters of fit that could numerically assess the sensibility of a fit to an individual signal, and if the signal were particularly noisy (due to things like inclusions) and the parameter of fit suggested such it could fall back to the standard "down-hole model" of a reference material. We may remove this sentence at the reviewer's recommendation, given the unexplored nature of this implementation. Regardless, we think it worthwhile elaborating on how the proposed method differs from the existing down-hole modelling methods as they are only briefly discussed in the manuscript's introduction.

Figure 2: We would happily make modifications to the caption and figure to address the reviewer's comments here. For figure 2b we could remove λ₄ as it is effectively zero, and add a supplementary figure that displays the shape of each component (λ 0 - 4) more clearly. The dot between λ_j and variable x is intended to signify that is it λ_j multiplied by the orthogonal function φ(j) of x. More concretely, λ₁* φ₁(x) is the linear orthogonal polynomial component of x. This needs to be revised on the figure as they are in the wrong spot, it should be λ_j⋅φ(x).

Line 81: It currently represents an arithmetic mean due to fitting on centred (i.e. geometric mean subtracted) ratios. These by definition end up centred around 0 with both positive and negative ratios causing log transformations to not be possible. This is more a limitation of how the data is initially ingested and the order in which operations are done. This could be overcome with a code rework. We don't see this as an issue currently, as we are only attempting to show this method is possible, and can make improvements to the code in future (contributions are welcome) to further improve the accuracy of the method from a compositional data perspective.

Line 115: We'll work to revise this paragraph for clarity, as it is rather complicated as currently written.

An example of what this paragraph was attempting to highlight: for session aggregated fits lambda 1-4 values for a non-centred and centred fit are slightly off from each other.

For example, lambdas 1-4 for 91500 on 2020-02-24 are 0.00065959, -3.5747e-5, 4.6642e-7, 4.1489e-8 when not centred, and 0.00066039, -3.5474e-5, 5.2599e-7, 5.198e-8 when centred. After investigating this further, it appears to be an artefact of the data not being calibrated against a known reference material and their being some drift in the session (which has not been corrected for). This would cause underlying data to have slightly greater spread along the y-axis and would then impact λ₁ (the slope), which then impacts λ₂, and so on with each higher term being partially dependant on the term before it. Using a synthetic dataset and introducing "drift" across it will help confirm this.

Importantly, individual analysis fits have identical lambda values for centred and non-centred datasets, implying that the algorithm operating correctly and this is an input data issue unique to session aggregated data.

This further highlights the need to centre data when comparing different analyses or aggregating data to remove the influences of machine drift and tuning.

Figure 4: We'll take this feedback on-board. The figure is there to demonstrate that we can compare data from different sessions without needing to first calibrate them to a known-ratio.

Line 132: We'll change this to "ratio" for clarity. Fitting signal intensity rather than ratio might provide us with greater insight on how specific elements are behaving during ablation (as a result of ionisation efficiencies, volatilities etc) and could be a potential application of the method. We have not explored this though.

Line 136: N is the total number of points. I can change this to "for all i from 1 to N" to improve clarity.

Eqs 7-9: The hat is removed from y in an already revised version as this is a typographic error.

The reviewer is correct in Eq. 9 being the variance-covariance of the lambda parameters, not the measurements in y, and the sigma term being a vector rather than a matrix. These are typographic errors. Eq. 9 is used to calculate the lambda values.

These equations and corresponding text will be changed to reflect this.

For Eq 8, notably the Omega term. The input uncertainties come from prior data ingestion steps, but they are not fully quantified with parameters like dead-time uncertainty, dwell time uncertainty etc. They are currently calculated from the underlying gas-blank subtracted count per second data for the ratio numerator and denominator intensities via the function load_agilent. These uncertainties are used as a weighting for the fit.

The reviewer is also correct that the exact current implementation is a weighted least squares due to the diagonal Omega matrix, but this is a specialised case of generalised least squares. The intention of having the algorithm set up like this currently is to allow for true covariance-variance input (accounting for dwell time and dead time) in future (Omega would then not be a diagonal matrix).

We welcome future, external discussion with the reviewer on this topic to refine the code accordingly but this does not impact the methodological foundations of the proposed fitting method as it gives results consistent to the REE fitting algorithm of Anenburg & Williams (2022) on their REE data.

The plotted lambda uncertainties in this paper are all calculated by the fitting algorithm using input uncertainties as weights for the fit, and not the default unit weights that are provided in the algorithm as these are only a fall back if the user does not provide uncertainties for weighting.

Figure 5, line 212: We would happily work through these comments to improve clarity here.

For reference, (a) are the fits for an aggregated dataset for a specified sample from all sessions (e.g. GJ, centred, all sessions), (b) are the fits for an aggregated dataset for a specified sample within a session (e.g. GJ, centred, single session). x is λ₁, y is λ₂, of an orthogonal polynomial fit to U/Pb (y) against time (x). The data have been lined up to a common time datum, laser on using the automatic times algorithm implemented in the code.

The reduced chi squared values are quite small, indicating either over-fitting of the data or the error variance has been overestimated. This is likely related to the set up of input and output analytical uncertainties, and as the reviewer has stated will lead to smaller uncertainties with increasing number of data points.

I've implemented a normalised root-mean-square-error (NRMSE, 0 being a "perfect fit" and 1 being a "nonsensical fit" of the data). For the worst case sample using the least appropriate order (lowest AIC weight), this value is around 0.25 and for an appropriate fit, the value is usually between 0.08-0.16 suggesting good-quality fits.

NRMSE values from the sample aggregated data (i.e. all of 91500), session fits (91500 per session) and analysis fits (91500 per analysis) very consistent across these models. I would expect to see greater variability in the NRMSE across the different aggregations if the input data were increasing the overall scatter due to them not following the same linear (x to y) relationship as each other.

These measures of fit aren't directly reported in the manuscript as they are directly available from the provided code and data, and it would create an excessively large table that is quite difficult to understand.

There is room for improvement in the uncertainty quantification and use of these measures of fit, with particular consideration of some of the pitfalls in reduced chi squared for linear regression (Andrae et al. 2010, preprint, http://arxiv.org/abs/1012.3754). However, this is will not change the overall conclusions drawn from this work, nor does it prevent the use of the proposed method, as the underlying fitting algorithm is well established in literature and as the reviewer has confirmed, is mathematically sound. The code is open source, and we welcome further discussion and contributions from the reviewer (and any other interested party) to improve the uncertainty quantification of the lambda parameters in the code.

We'll take on the feedback to highlight where a user should utilise analysis or aggregated data lambda values. For individual down-hole fitting, analysis fitting is what should be used (i.e. assessing how well modelled a specific analysis is), session based aggregation would be used to check within lab variability, and sample based aggregation would be used for between lab variability/long term assessment of reference materials or to choose a suitable reference material to use as a calibration for an unknown.

Line 215: We will use this feedback to make improvements here. From my current understanding, the lambda values are unit less. I.e. λ₀ is derived from the mean cps/cps (units cancel) and all other terms are multipliers of functions, that themselves have units, e.g. λ_j*φ_j(x) where, for example, φ₂ = (x_i - γ₁)(x_i - γ₂) and has units of x² (seconds squared).

Line 227: We will take this feedback on board as it also assists with the prior two comments. In essence λ₀ is the only value sensitive to ICP-tuning and drift, thus if data are first centred to account for these factors the values of λ₁ and higher, should only represent the shape parameters of an analytical signal and not change with respect to a centred and non-centred dataset of single analysis fits. For aggregated datasets, if the values of λ₁ and higher are not identical for centred and non-centred datasets this suggests an external influence is not yet accounted for (i.e. machine drift within a session or intersession tuning).