Calculation of Uncertainty in the (U-Th)/He System

. Currently, there is no standardized approach for reporting date uncertainty in the (U-Th)/He system, partly due to 5 the fact that the methods and formulae for calculating single-grain uncertainty have never been fully described and published. This creates challenges for interpreting the expected distribution of dates within individual samples and for comparing dates generated by different labs. Here we publish two procedures to derive (U-Th)/He single-grain date uncertainty (linear and Monte Carlo uncertainty propagation), based on input 4 He, radionuclide, and isotope-specific F T (alpha-ejection correction) values and uncertainties. We also describe a newly released software package, HeCalc, that 10 performs date calculation and uncertainty propagation for (U-Th)/He data. Propagating uncertainties in 4 He and radionuclides using a compilation of real (U-Th)/He data (N = 1978 apatites and 1753 zircons) reveals that the uncertainty budget in this dataset is dominated by uncertainty stemming from the radionuclides, yielding median relative uncertainty values of 2.9% for apatite dates and 1.7% for zircon dates (1s equivalent). When uncertainties in F T of 2% or 5% are assumed and additionally propagated, the median relative uncertainty values increase to 3.5% and 5.8% for apatite dates, and 15 2.6% and 5.2% for zircon dates. The potentially strong influence of F T on the uncertainty budget underscores the importance of ongoing efforts to better quantify and routinely propagate F T uncertainty into (U-Th)/He dates. Skew is generally positive and can be significant, with ~1417% of apatite dates and ~56% of zircon dates in the data compilation characterized by skewness of 10%0.25 or greater. This outcome indicates the value of applying Monte Carlo uncertainty propagation to identify samples with substantially asymmetric uncertainties that should be considered during data interpretation. The 20 formulae published here and the associated HeCalc software can aid in more consistent and rigorous (U-Th)/He uncertainty reporting, which also is a key first step in quantifying whether multiple aliquots from a sample are overdispersed, with dates that differ beyond what is expected from analytical and F T uncertainties.


Introduction
Geochronology and thermochronology by the (U-Th)/He method was initially developed as a reliable technique 25 approximately three decades ago Wernicke and Lippolt, 1994;Wolf et al., 1996;Zeitler et al., 1987). Since that time, numerous advances such as the ability to measure the (U-Th)/He date of individual grains (e.g., House et al., 2000), improvements in kinetic models to account for the effects of radiation damage accumulation and annealing on He diffusion kinetics (e.g., Flowers et al., 2009;Gautheron et al., 2009;Guenthner et al., 2013), and the development of thermal history modeling tools that improve interpretation of these data (Gallagher, 2012;Ketcham, 2005) have led to the 30 widespread application of this technique and large amounts of data generation. However, with this progress has come recognition of the need to more rigorously and consistently report uncertainties on individual (U-Th)/He dates Ketcham et al., 2022). For example, the intra-sample variability of (U-Th)/He dates often exceeds that predicted by analytical uncertainty, both due to interpretable variation from differences in He diffusion kinetics among grains of the same sample and due to scatter from other factors (e.g., Brown et al., 2013;Fitzgerald et al., 2006;Flowers et al., 2022b). Better 35 accounting for the uncertainties of individual analyses is a key step in determining whether multiple individual analyses from a sample are actually "over-dispersed", and would help develop a more complete understanding of the causes of data dispersion by allowing the scatter attributable to analytical uncertainty to be subtracted from the overall dispersion pattern. In addition, more rigorous uncertainty reporting would improve confidence in large-N datasets, facilitate inter-laboratory data comparisons, and ultimately increase the precision and accuracy of thermal history reconstructions. 40 One current challenge to comprehensive uncertainty propagation is that, although individual laboratories have derived the methods for propagating uncertainty components into single-grain (U-Th)/He dates, these methods and the resulting formal analytical uncertainty in (U-Th)/He dates have never been described or thoroughly assessed in the literature. It also is unclear if different labs propagate uncertainties in the same manner. Uncertainty propagation in the (U-Th)/He system is complicated by the fact that the age equation has no analytical solution, precluding the direct application of typical 45 uncertainty propagation formulae that combine individual uncertainty components in quadrature through a given function. This problem may be circumvented by approximations of the He age equation that solve directly for time (e.g., Meesters and Dunai, 2005), or by the use of the general "error propagation equation" using the first derivatives of the uncertainty components with respect to time (Bevington and Robinson, 2003). However, linear uncertainty propagation methods rely on an assumption that the derivative of the first term of the Taylor series is a linear function at the scale of the uncertainties 50 being combined (Bevington and Robinson, 2003;McLean et al., 2011). As this assumption is often violated in the (U-Th)/He system, uncertainties have the potential to be skewed (i.e., asymmetric), and uncertainties propagated using standard linear uncertainty propagation may be inaccurate. Comprehensive uncertainty accounting on individual (U-Th)/He dates involves propagating not only the analytical uncertainties associated with measurements of parent and daughter amounts, but also propagating uncertainties associated 55 with alpha-ejection corrections (FT corrections, which account for He ejected from the crystal via alpha decay). The analytical uncertainty on parent and daughter amounts is generally well-characterized, and the geometric uncertainty on FT values is increasingly well-constrained (e.g., Cooperdock et al., 2019;Glotzbach et al., 2019;Zeigler et al., 2022). As FT uncertainties are better quantified, propagating both analytical and FT uncertainties into the reported uncertainty of (U-Th)/He dates is desirable (e.g., Flowers et al., 2022a). 60 Here we explain how analytical and FT uncertainties in (U-Th)/He dates may be combined to derive a single-grain (U-Th)/He date uncertainty. To address the shortcomings of linear uncertainty propagation, we adopt primarily a Monte Carlo approach to quantitatively constrain (U-Th)/He uncertainty. This procedure is both accurate and mathematically simple, and enables evaluation of asymmetric uncertainties (which linear uncertainty propagation does not provide). For completeness and to ease retrospective data comparisons, we also include a method to propagate compute date uncertainty 65 through the calculation of a (U-Th)/He date that relies on uses more traditional linear uncertainty propagation. In addition, this manuscript presents a new program written in Python 3.8 termed HeCalc (Helium date and uncertainty Calculator; Martin, 2022) that is capable of both Monte Carlo and linear methods of uncertainty propagation. We conclude by using HeCalc to reduce a compilation of real data to determine the typical contributions of each uncertainty component to date uncertainty in actual practice. 70

Background: uncertainty components in (U-Th)/He dates
The currently quantifiable uncertainties on single-grain (U-Th)/He dates include analytical uncertainties associated with parent and daughter isotope measurements and geometric uncertainties associated with FT alpha-ejection corrections. These are discussed in detail in Flowers et al. (2022a) and summarized more briefly here. We use the word "uncertainty" as a probabilistic statement of the distribution of repeated measurements (e.g., for a 238 U measurement of 10 ± 1 µg/g at 1σ, 75 68.27% of repeated measurements will fall between 9 and 11 µg/g), while "error" refers to the deviation of a measured value from the true value. The uncertainty in decay constants is negligible relative to other sources of uncertainty and results incauses systematic error across all (U-Th)/He measurements, and therefore is not incorporated in our uncertainty calculation methods.
In the (U-Th)/He technique, the parent nuclides (Uranium, Thorium, and Samarium; 238 U, 235 U, 232 Th, 147 Sm) are 80 typically measured using inductively coupled plasma mass spectrometry (ICP-MS), while the daughter product (Helium; 4 He) is usually measured on a quadrupole or magnetic sector noble gas mass spectrometer. Most commonly, quantification of 4 He and its parent nuclides is performed via isotope dilution to permit conversion from isotopic ratio measurements to molar amounts. Given the measurements of parent and daughter products, a (U-Th)/He date may be calculated using the equation for 4 He ingrowth 85 He= 8  (1) where each nuclide is given as an amount, t is time, and λ is the decay constant for each parent nuclide given in the subscript.
Because of the kinetic energy associated with alpha decay, individual alpha particles (i.e., 4 He nuclei) travel between 4 and 34 µm in solid matter before coming to rest, depending on the mineral density and parent or intermediate 90 daughter nuclide Ketcham et al., 2011). This redistribution of the daughter product can cause He ejection from a crystal. By assuming a homogenous parent nuclide distribution, measuring the physical dimensions of a single grain, and applying a geometric model to those physical dimensions, the proportion of alpha particles retained in a grain (the fraction trapped; FT) can be calculated for each nuclide's mean stopping distance . Determination of grain dimensions to calculate FT is usually accomplished via size measurement of individual grains using photomicrographs 95 with a calibrated digital camera (Cooperdock et al., 2019;Glotzbach et al., 2019 (2) 100 We refer to dates calculated with this correction applied as "alpha-ejection corrected" or simply "corrected" dates, while dates calculated with no correction applied using Eq.
(1) we refer to as "uncorrected" or "raw" dates. For this work, it is assumed that the amount and uncertainty of each nuclide has been constrained. The natural U isotopic ratio (137.818 ± 0.023 1s; Hiess et al., 2012) is usually used to calculate 235 U in a sample based on the 238 U amount measured. In these cases, the uncertainty in 235 U is perfectly correlated with 238 U; treatment of these uncertainties as though 105 they were independent could lead to inaccurate uncertainty calculations. Whether or not the uncertainty in the other radionuclides is correlated depends on the details of isotope spiking procedures and must be evaluated on an individual lab basis. The option to incorporate correlated radionuclide uncertainty is included in the methods for propagating uncertainty. Uncertainty and systematic error in FT values likely stems from a combination of inaccurate grain measurements, assumptions regarding the specific geometry of a given grain (i.e., deviations from the "idealized" shapes in Ketcham et al., 110 2011), and undetected parent nuclide zonation (e.g., Farley et al., 1996). When the magnitudes of these effects are constrained then the corresponding uncertainties can be propagated into the FT value uncertainty and FT values can be corrected for systematic error. Several approaches have been developed to approximate the 3D shape of individual grains to assess uncertainty associated with grain geometry estimates, generally finding that the geometric uncertainty in FT ranges between 2-8% (1s) (Cooperdock et al., 2019;Evans et al., 2008;Glotzbach et al., 2019;Herman et al., 2007;Zeigler et al., 115 2022). The magnitude of systematic error depends on grain shape and the details of the method used for FT value determination (e.g., Cooperdock et al., 2019;Glotzbach et al., 2019;Zeigler et al., 2022). Measurement of parent nuclide zonation is not currently done in typical workflows, so this source of uncertainty in FT is generally unquantified, but it could be included in FT uncertainty in the future if labs characterize zonation prior to dating, for example via grain mapping or depth-profiling of parent-isotope zonation using laser-ablation ICP-MS methods (Farley et al., 2011;Johnstone et al., 2013). 120 The isotope-specific FT values are highly correlated (Zeigler et al., 2022), so we include correlated FT uncertainty in the methods below. Several additional sources of dispersion in (U-Th)/He dates exist, including alpha implantation (e.g., Murray et al., 2014) and the influence of defects on He diffusion (e.g., Zeitler et al., 2017). These factors potentially contribute to intrasample variability, but would not cause dispersion in repeated measurements of the same grain, and thus are best considered 125 as part of multi-aliquot data compilations.

Date and uncertainty calculation methods
Here, (U-Th)/He dates are calculated by first estimating a date using an approximation of the helium age equation that solves directly for time. Using this estimate as an initial value, the exact date is then calculated iteratively using the Newton-Raphson method. We describe two independent methods (linear uncertainty propagation and Monte Carlo 130 uncertainty modeling) of calculating the uncertainty in this date given the uncertainty components described in Sect. 2 above. We exclusively use the term "linear uncertainty propagation" rather than "analytical" or "standard" propagation to avoid confusion with analytical error arising from instrument noise and standards used in analytical measurements, respectively. The linear method allows precise and repeatable calculations, while the Monte Carlo method is slightly more accurate and allows for calculation of skewed probability distributions, as discussed further in section 5. 135

Date calculation
The initial value for iterative age calculation is obtained by calculating an approximated noniterative solution of the (U-Th)/He age equation as described by Meesters and Dunai (2005). We slightly modify the production term in this method to permit calculation of parent-specific alpha ejection-corrected effective helium production rates: = × × × 140 (3) Where pj is the 4 He production rate, N is the number of alpha particles produced by a given decay chain, and j FT, λj, and j M are the alpha ejection-correction factor, decay constant, and concentration of radionuclide j (i.e. 238 U, 235 U, 232 Th, and 147 Sm), respectively.
Using the resulting date approximation as an initial guess (t0), the (U-Th)/He date is then found using the relatively 145 simple but highly efficient Newton-Raphson method where ti and ti+1 are successive approximations of the date, and f(ti) and f'(ti) are the implicit age equation (the helium age equation set at zero; Eq. (5)) and its first derivative with respect to t, respectively. This calculation is repeated until the difference between successive iterations is less than one year. This method benefits from an accurate initial guess and a 155 quadratic rate of convergence such that generally only three to five iterations are required, though for dates >500 Ma (where the noniterative approximation produces relative errors of >0.1% ; Meesters and Dunai, 2005), as many as ten iterations may be needed.

Linear uncertainty propagation
Here we provide a method of calculating date uncertainty using linear propagation of uncertainty. We apply the 160 general formula for uncertainty propagation through a function f(a, b…z), including cross terms for correlated error where such correlations exist (Bevington and Robinson, 2003): The following equations presume that 235 U has not been measured directly, but equations that include directly 165 quantified 235 U are provided in Appendix A, and the HeCalc software released with this paper includes an option to account for either means of constraining 235 U. As an alternative to the use of ing HeCalc, these equations could be replicated in spreadsheet programs with a one-time expenditure of effort.
Applying the uncertainty propagation equation to the (U-Th)/He age equation, including potential covariance in the radionuclide and FT uncertainties (i.e., the potential that the uncertainties are not fully independent), indicates that the 170 uncertainty in a (U-Th)/He date is: (8) where, for example, σHe is the uncertainty in the 4 He measurement, and σ238-232 2 is the covariance between 238 U and 232 Th. Note that the covariance terms collapse to 0 if no correlation exists between uncertainties, while positive covariance will 175 increase the overall uncertainty. While solving the (U-Th)/He age equation for t explicitly is not possible, finding the first derivative of t with respect to each variable is possible through implicit differentiation. Specifically, where X is each variable in the (U-Th)/He age equation with an uncertainty. Using this relationship, the relevant derivatives are: Where each summation term involves addition of the four radionuclides with the same variable convention described in Sect. These equations are printed in their expanded forms, along with versions that allow for direct quantification of 235 U, in Appendix A.

Monte Carlo uncertainty 6modelling
Monte Carlo uncertainty propagation is based on the approach of combining the uncertainty in measured parameters 200 with any given probability distribution (including non-gaussian distributions as may be caused by compositional zoning; Hourigan et al., 2005) by randomly sampling each distribution a large number of times and propagating those randomly generated parameters through some function of interest (Eqs. (7) and (10); Fig. 1). This method yields a probability density histogram that describes the true uncertainty to arbitrary precision depending on the number of simulations run (Anderson, 1976;Possolo and Iyer, 2017). As such, the application of Monte Carlo techniques is mathematically straightforward, in this 205 case requiring no knowledge beyond that required to calculate a (U-Th)/He date. In addition to this benefit, Monte Carlo uncertainty analysis does not require that a function of interest have a linear first term of the Taylor series to accurately calculate uncertainty; when this assumption is violated (as in the (U-Th)/He age equation), uncertainties propagated using linear uncertainty propagation (Eq. (7)) can be inaccurate. Future work may also permit analysis of the effects of non-linear input uncertainties to date uncertainty. While the Monte Carlo method has historically been hindered by computational 210 expense, the increases in computational power in recent decades make this more accurate approach an attractive method for routine uncertainty propagation in (U-Th)/He chronology.
Here, Monte Carlo uncertainty modeling of (U-Th)/He data is performed by generating arrays of a pre-determined size N, which contain randomly generated values for each input according to the gaussian distribution described by each value's 1σ uncertainty (Fig. 1, input probability distributions). Correlated uncertainties (correlations between 238 U, 232 Th, and 215 described above using these randomly generated variables. From these arrays, 68% and 95% confidence intervals are calculated using the 15.865 and 84.135, and 2.275 and 97.725 percentiles of the samples of dates, respectively. We use 220 confidence intervals as opposed to standard deviation because some output uncertainty distributions are skewed. Although the average of the 68% and 95% confidence intervals yields the 1-and 2-standard deviation levels for reasonably gaussian (normal) distributions, this does not necessarily hold for non-gaussian (asymmetric or skewed) distributions ( Fig. 1, example output probability distributions).

Figure 1: A conceptual diagram of Monte Carlo uncertainty modeling for the (U-Th)/He system. Each independent gaussian input probability distribution (with 1σ standard deviation marked as vertical lines) is sampled at random a large number of times.
Because we assume fully correlated FT uncertainties, the isotope-specific FT distributions are represented by a single distribution where an individual percent deviation from the mean value is sampled four times. Using these randomly sampled inputs, a single date is calculated. This process is repeated until the probability distribution of interest (in this case, a skewed non-gaussian distribution 230 with the 68% confidence interval shown with vertical lines) has been sufficiently sampled, as set by the analyst.
The number of Monte Carlo simulations dictates the precision of the results because Monte Carlo analysis is a numerical approximation of uncertainty, (e.g., the lower panels in Fig. 1 become progressively smoother with an increasing number of simulations). Therefore, separate from the probability distribution describing date uncertainty, there is a predictable level of variation in uncertainty estimates and other parameters describing the probability distribution (e.g., its 235 mean) given a certain number of total Monte Carlo simulations (Wübbeler et al., 2010). Specifically, the standard error of the standard deviation of a Monte Carlo model is dependent on the uncertainty in the value itself and the number of simulations: where σµ is the standard deviation of the population mean, σt the date uncertainty, and N the number of simulations. To avoid running arbitrary numbers of simulations, we invert this equation to determine the number of iterations required to achieve a user-requested relative precision on the mean: Where N is the number of simulations to run, σt is the date uncertainty estimated by linear uncertainty propagation, x̄ is the sample mean estimated by calculation of the date using the nominal input values, and p is the user-requested precision in percent uncertainty. By using percent relative uncertainty, the value of the date itself need not be known a priori, as an 250 estimate of the standard deviation of the population mean date can be calculated using the percent relative precision and the date calculation from the input values (Eq. (18)).

Helium date and uncertainty Calculator (HeCalc) code
In this section we describe the implementation of the above methods of date and uncertainty calculation in the new HeCalc software (Martin, 2022). For ease of access and to best provide this software as a resource to the (U-Th)/Hhe lium 255 community, HeCalc is available both as a standalone program with a graphical user interface (GUI) and as a package in Python 3 available for download from the Python Package Index (PyPI) via pip commands. The descriptions below apply specifically to the GUI version of the software and the main_hecalc() function in the Python package; those interested in writing their own code and incorporating the component functions provided in HeCalc may consult the associated documentation for more detailed programming considerations. 260

Input
The input for HeCalc is designed to be straightforward and flexible (Table 1). Input files may be in Excel (.xls/.xlsx), comma separated value (.csv), or tab-delimited text (.txt) format. In addition to data input through a file, HeCalc users may manually input values to calculate a date and uncertainty for a single set of data by clicking on the "Manual" tab. If importing data through a file, the file must contain columns for sample name, U, Th, Sm, He, and all FTs with the headers 265 Sample, mol 238U, mol 232Th, mol 147Sm, mol 4He, 238Ft, 235Ft, 232Ft, and 147Ft (Table 1). Although "mol X" is required as the input column header for U, Th, Sm, and He, the actual units of the input data may be any unit of quantity (e.g., atoms, mol/g, etc.) as long as they are identical. The 1σ uncertainty for each value, in the same units, must be included in the column following each respective value, even if the applied uncertainty is 0 (e.g., for FT values with unknown uncertainty); there is no naming requirement for these headers. If 235 U was measured directly, columns for this measurement 270 and its uncertainty should also be present. Correlated uncertainty between the radionuclides and between the isotope-specific FT values can be input using their Pearson correlation coefficient, which is related to the covariance as: where σab 2 is the covariance between variables a and b. The correlation coefficient is preferable to inputting covariance 275 directly as it has the intuitive meaning of being in the range of [-1, 1] where 1 is perfectly correlated and 0 is fully uncorrelated, while numerical covariance is generally unintuitive. These values may be included in the input file using headers with the naming convention "r 238U-235U" and "r 238Ft-235Ft" (Table 1); either ordering of the correlated uncertainties in the header (i.e., "r 238U-235U" vs. "r 235U-238U") is permitted. Uncertainties are assumed to be uncorrelated unless these columns are explicitly included. Example input files both with and without correlated uncertainty 280 are provided as templates in the code's repository (see the Code availability Section for a direct link). The order in which these columns appear is unimportant as long as the uncertainty associated with each value follows that value. Extraneous columns with differing headers also will not interfere with the code's execution. Additionally, if an input Excel file has multiple sheets, the first sheet will be read in by default. If this sheet does not contain the required column headers, the program will ask for the name of the sheet to use instead. In this way, HeCalc ideally allows for input of 285 any given lab's standard data reduction spreadsheet or other typical data product with no or minimal alteration, allowing it to be integrated seamlessly into a lab's existing workflow.
In addition to data input, several further options are provided. The number of decimals included in the output is determined by the user (this option affects only output and does not impact the statistical aspects of the code). The user can also select whether to perform linear uncertainty propagation, Monte Carlo uncertainty propagation, both, or neither. If 290 Monte Carlo uncertainty propagation is selected, the desired precision of the mean is specified in percent as described above.
In practice, the precision of the mean date need be no better than the number of significant figures present the in data; for common (U-Th)/He analyses, this equates to a precision of ~0.01%, which generally requires on the order of 10 4 -10 5 simulations. The program also contains the ability to generate histograms using the Monte Carlo results. If this option is chosen, this histogram may be parameterized as a skew-normal distribution (Azzalini and Capitanio, 1999;O'Hagan and 295 Leonard, 1976 1.00.9 r 238Ft-147Ft 1.00.9 r 235Ft-232Ft 1.00.9 r 235Ft-147Ft 1.00.9 r 232Ft-147Ft 1.00.9 a All values in this table are purely for illustration and do not reflect actual data. b r values are uncertainty convariance as given by their Pearson coefficient (Eq. 19). If no values are provided, uncertainties are assumed to be uncorrelated with r = 0.

Output
There are two main outputs from HeCalc: the results of the date calculation and uncertainty propagation, and the 300 histograms of the Monte Carlo results for each sample (Table 2). At a minimum, the sample name, raw date, and corrected date are saved to an Excel sheet titled "Uncertainty Output" that includes a header with the input file's name and directory. The raw and corrected dates in these columns are calculated using each exact input value (e.g., mol 238U = 1 in Table 1); we refer to these dates as "nominal dates" below. The selection of linear uncertainty propagation causes columns to be added titled "Linear raw 1σ uncertainty", "Linear raw 2σ uncertainty", "Linear corrected 1σ uncertainty", and "Linear corrected 2σ 305 uncertainty", with "raw" indicating no alpha ejection correction. If Monte Carlo error propagation is selected, a header line specifying the user-requested precision is added, and the columns "MC average 68% CI, raw", "MC +68% CI, raw", "MC -68% CI, raw", "MC average 95% CI, raw", "MC +95% CI, raw", "MC -95% CI, raw", and the corresponding values for FTcorrected dates (titled with "corrected" instead of "raw") are included along with a column giving the number of Monte Carlo simulations run. The confidence intervals are reported as the 15.865 and 84.135 percentiles (the 68% confidence 310 interval) and the 2.275 5 and 97.725 5 percentiles (the 95% confidence interval) of the Monte Carlo results, converted to uncertainty values by reference to the nominal date. Throughout this manuscript, the asymmetry of the confidence intervals will be calculated with respect to the nominal date calculation. It is worth noting that the nominal date does not strictly correspond to the mode of the histogram, and instead falls toward the skewed side, meaning that the skew calculations presented here are a slight underestimate of the actual asymmetry in the distribution. 315 If the user chooses to include histograms in the output, an Excel sheet titled "Histogram Output" is added to the workbook, with columns for the center of each histogram bin (i.e., the individual intervals in the histogram) and number of simulations in that bin as x-and y-values for the both the raw and FT-corrected dates. Four total columns are therefore present for each sample. The number of bins is equal to 1/1000 th the number of simulations run or ten bins, whichever is greater. If parameterization is selected, the histogram is fit to a skew-normal distribution. Although this distribution does not 320 perfectly replicate the histograms generated by HeCalc, it allows for first-order interpretations using continuous probability distributions. Columns are appended to the end of the "Uncertainty Output" sheet titled "Hist raw fit a", "Hist raw fit u", "Hist raw fit s", and the corresponding values for FT-corrected calculations. These parameters correspond to the shape ("a", the skewness), location ("u", a measure of central tendency), and scale ("s", the width of the distribution) parameters for a skew-normal probability distribution function (Azzalini, 1985;O'Hagan and Leonard, 1976). 325 Hist corrected fit s 10.531 Parameterization selected (requires Monte Carlo) a The header for the file will contain a line for the file path of the input file and (if Monte Carlo propagation is selected) the user-requested precision.
b "Hist fit a" refers to the shape or skewness of the histogram. "Hist fit u" is a measure of central tendency of the histogram. "Hist fit s" is a measure of the width of the distribution. (Azzalini, 1985;O'Hagan and Leonard, 1976).

Discussion
Here we use the methods described above to calculate the dates and uncertainties for a compilation of real apatite 330 and zircon (U-Th)/He data, and then examine the overall uncertainty budget and skew in this dataset. In the appendices we additionally explore the influence of theoretical input uncertainties on date uncertainty, on skew in Monte Carlo date distributions, and on the differences in dates derived from the linear and Monte Carlo methods (Appendices C-E).

Uncertainty budget in real data
To assess the uncertainty budget in real (U-Th)/He data, we compiled 1,978 apatite and 1,753 zircon analyses that 335 were acquired with typical and consistent (U-Th)/He methods and instrumentation (quadrupole noble gas mass spectrometer and quadrupole ICP-MS) in the University of Colorado Thermochronology Research and Instrumentation Laboratory (CU TRaIL). These data were measured from October 2017 to March 2020 following the methods described in Sturrock et al. (2021) for apatite and Peak et al. (2021) for zircon. Figure 2 shows the histograms of percent relative uncertainty in the absolute amounts of 4 He, 238 U, 232 Th, and (for apatite) 147 Sm for this dataset, while Table 3 lists the median value and 68% 340 confidence interval calculated using the percentile approach (Sect. 3.3) for these distributions. Analytical uncertainties are typically higher for apatites than for zircons due to the lower 4 He and radionuclide amounts for apatites relative to zircons, which causes apatite measurements to have lower count rates that are more impacted by blank and background uncertainties. For both apatite and zircon analyses, radionuclide uncertainties are higher than 4 He uncertainties. For example, for apatites, the percent uncertainties in 238 U, 232 Th, and 147 Sm amounts are 3.2%, 2.8%, and 2.8%, respectively, which are 3-4x less 345 precise than the uncertainty in the 4 He amount of 0.86% ( Fig. 2a; Table 3). For zircons, the uncertainties of 238 U and 232 Th are 1.8% and 2.2%, respectively, which are 5-6.5x less precise than the 4 He uncertainty of 0.34% (Fig. 2b, Table 3).  N = 1,753). Note that the y-axis scales differ for these plots. We analyze the uncertainty in these data with and without propagating FT uncertainty using HeCalc (Fig. 3). For illustrative purposes, we explored two different scenarios assuming 2% and 5% uncertainties in the isotope-specific FT values and fully correlated FT uncertainties. These uncertainties are based on those reported by Zeigler et al. (2022), which 355

350
found that the FT uncertainty depends partly on grain geometry. Figure 3 shows the distributions of percent relative uncertainty calculated by the Monte Carlo method for the apatite and zircon corrected (U-Th)/He dates. Table 4 lists the median value and 68% confidence interval for these distributions. Propagating only analytical uncertainties (on radionuclides and 4 He) yields median date uncertainties of 2.9% and 1.7% for apatite and zircon, respectively. With uncertainty in FT included, the date uncertainty increases substantially. For apatites, the uncertainty value increases from 2.9% for analytical 360 uncertainties alone to 3.5% and 5.8%, respectively, when FT uncertainties of 2% or 5% are also propagated. For zircons, the uncertainty increases from 1.7% to 2.6% and 5.2%, respectively. The addition of FT uncertainty most heavily impacts the analyses with more precise radionuclide and 4 He measurements because FT uncertainty comprises a correspondingly larger proportion of the uncertainty budget. Given that FT geometric uncertainty estimates are on the same order of magnitude asand in some cases larger than-typical analytical uncertainties in (U-Th)/He dating, additional efforts to constrain FT 365 uncertainty (i.e., on minerals other than apatite, such as zircon) are important to rigorously calculate uncertainties in individual (U-Th)/He dates.

1,753) analyses. The top panels show input uncertainties for only analytical uncertainties ( 4 He and radionuclides), while the lower panels additionally include 2% or 5%geometric uncertainties in FT. Calculations assume fully correlated isotope-specific FT uncertainties (r = 1) based on Zeigler et al. (2022)
, and assume fully uncorrelated radionuclide uncertainties for these radionuclide data. Note that the y-axis scales differ for these plots.

Skew in real data
We also analyze the compilation of real data in Figure 2 for skew in Monte Carlo-generated date probability distributions (i.e., asymmetric uncertainty). Skew refers to the extent of asymmetry in the "tails" of a distribution 380 (Fig. 1, D1). Although "Sskewness" is a statistical concept that strictly referring refers to the third standardized moment of a population, which is a this metric is unitless and generally unintuitive metric. H, so here we additionally report a percent skew to more intuitively convey the distribution asymmetry, calculated by taking the percent difference between the positive and negative 68% confidence intervals with respect to the date. This asymmetry would most accurately be reported as separate positive and negative uncertainty values referring to the 68% confidence interval rather than the more typical 1σ 385 uncertainty reporting of symmetric uncertainty. As an example, for 10% skew, a 100±6.4 Ma date would be represented as 100 [+6.7, -6.1] Ma at the 1s level. Analysis of theoretical data reveals that skew increases with relative input uncertainty and varies with age (Appendix D; Fig. D1-D3).
In our data compilation, positive skew is common and can be significant (Table 3). This is consistent with 390 aAnalysis of theoretical data that reveals that skew increases with relative input uncertainty and varies with age (Appendix D; Fig. D1-D3). , as would be predicted based on the theoretical data analysis and the real data uncertainties in Table 3. For the real dataset, wWith only analytical uncertainty included, the median skew in apatites and zircons is 4.4% 0.05 and 3.2%0.02 (or 4.4% and 3.2% skew), respectively. The inclusion of FT uncertainty increases the skew (Fig. 4a-b, Table 4). For apatites, 395 skew rises to 6.0%0.13 and 11.4%0.33 (or 6 and 11.4% skew) for 2% and 5% uncertainty in FT, respectively (Fig. 4a, Table  4). For zircons, the same combinations of uncertainty yield skews of 5.2%0.11 and 11.0%0.32 (or 5.2 and 11% skew; Fig.  4b, Table 4). For a 2% FT uncertainty, ~1417% of apatite data and ~65% of zircon data have a a percent skewness of 10%0.25 or greater. As an example, for 10% skew, a 100±6.4 Ma date would be represented as 100 [+6.7, -6.1] Ma at the 1s level. 400 General practice in (U-Th)/He dating has been to report symmetric uncertainties. Our analysis reveals that for most cases this is appropriate, and averaging asymmetric uncertainties in data reporting is unlikely to substantially impact interpretations. However, for highly asymmetric uncertainties, it may be appropriate to report positive and negative uncertainties separately, and only combine the reported uncertainties if they are indistinguishable within the appropriate number of significant figures. Our results suggest that skew may be an important consideration when interpreting (U-Th)/He 405 data with less precise 4 He and radionuclide measurements, because these data generally have date uncertainties with greater skew. In these cases, asymmetric uncertainties may be important for determining whether a dataset is consistent with a given hypothesis within uncertainty. However, a challenge to interpreting dates with asymmetric uncertainties is that no widely used inverse thermal history modeling software for (U-Th)/He data permits asymmetric uncertainty input. Future work implementing skewed probability distributions in such software may enhance interpretation of the subset of (U-Th)/He data 410 characterized by highly skewed uncertainties.

Comparison of linear and Monte Carlo uncertainty propagation results for real data
Finally, we compare the uncertainties derived from linear uncertainty propagation with the averaged 68% confidence intervals from Monte Carlo propagation for the compiled dataset. For nearly all analyses, the uncertainties yielded by the two methods are within 1% of each other, regardless of the amount of FT uncertainty included (Fig. 5a, Table 4). Thus, error due to linear uncertainty propagation (i.e., inaccurately calculated uncertainty resulting from an assumption of linearity) is 425 largely insignificant in this dataset, and the uncertainties yielded by the two methods are interchangeable in most circumstances. In contrast, skew, discussed in the previous section, is only revealed by the Monte Carlo method using the equations presented here and included in the HeCalc software. The asymmetric uncertainties of some samples, specifically those with atypically large input uncertainties, is more important for accurate uncertainty analysis than error introduced in the uncertainty calculation due to linear uncertainty propagation. 430  (N = 1,978) and (b) zircon (N = 1,753) analyses. Input uncertainties include analytical uncertainty only ( 4 He and radionuclides), and analytical uncertainties propagated with 2% or 5% geometric uncertainty in FT, assuming fully correlated individual FT uncertainty values (r = 1). The outlines of 435 covered histograms are included to show detail for each. Note that the y-axis scales differ between the two plots.

Conclusions
Here we publish fully traceable end-to-end calculations of uncertainty in (U-Th)/He dates, including the 440 propagation of uncertainties in FT values. We also provide a software package, HeCalc, to do these calculations explicitly and to perform more accurate Monte Carlo propagation of these uncertainties. Using a compilation of apatite and zircon (U-Th)/He analyses, we find that for a common instrumental setup (quadrupole noble gas and ICP mass spectrometers), uncertainty in radionuclide quantification is generally 3-6.5x larger than the uncertainty in 4 He measurement. When only 4 He and radionuclide uncertainties are propagated, the typical alpha ejection corrected (U-Th)/He date uncertainty is 2.9% for 445 apatites and 1.7% for zircons. The inclusion of 2% and 5% geometric uncertainty in the FT values yields greater date uncertainty of 3.5% and 5.8% for apatites and 2.6% and 5.2% for zircons. For the compiled dataset, the asymmetry in the 68% confidence interval can be significant, especially for dates with less precise input uncertainty. With 2% uncertainty included in FT, 174% of all apatite and 65% of all zircon analyses have a skewness of greater than 10%0.25. The results of linear uncertainty propagation for these data agree with those from Monte 450 Carlo uncertainty propagation to within ~1%, indicating that this error is negligible. Given that Monte Carlo uncertainty propagation permits calculation of skewed probability distributions and does not make an assumption of linearity in the (U-Th)/He age equation, we propose that this method should be preferred for uncertainty calculation in (U-Th)/He data. However, the current lack of a means of including asymmetric uncertainty in thermal history modeling, and the roughly equivalent symmetric uncertainty values from Monte Carlo and linear uncertainty propagation methods, indicates that the 455 results are likely interchangeable for common workflows, pending advancements in the (U-Th)/He method and interpretative models. The methods presented here allow for more rigorous inter-laboratory data comparisons and retrospective data analyses by providing a consistent means of quantifying the uncertainty budget of a given (U-Th)/He analysis. Further developments of the (U-Th)/He technique are also facilitated by this study. In particular, this work suggests that continued 460 refinement of FT uncertainty is warranted, and provides a framework into which those developments may be placed. Using the Monte Carlo results, asymmetric uncertainty may also be quantified, and could potentially be included in future versions of thermal history modeling software. Finally, fully accounting for analytical and geometric uncertainties will better isolate the magnitude of overdispersion and promote more effective examination of its causes. (a411)

Appendix B: Implications of gaussian input uncertainties in HeCalc
Negative dates are permitted in the probability distributions produced by HeCalc; this is because the input distributions are presumed to be gaussian, meaning that if the input variables have high relative errors, negative molar 500 amounts of U, Th, Sm, and He are possible. This behavior is formally correct for gaussian uncertainties, albeit non-physical. For low count rates associated with high relative uncertainty, a Poisson distribution (rather than gaussian distribution) is appropriate, and would prevent negative input values. However, high relative input uncertainties are generally due to a measurement being near or below background rather than low count rates where the underlying poisson distribution of the data is not well approximated by a gaussian. As a result, there are potential instances of negative molar amounts included in 505 the Monte Carlo calculations. In some rare instances when a negative amount of a given parent nuclide is produced in the generation of random data, the (U-Th)/He date equation may have multiple or no solutions. In these cases, the result is simply removed from the sample of calculated ages. The total number of such removals is tracked, and if the proportion removed exceeds the requested precision level, all results associated with the Monte Carlo simulation is reported as NaN (i.e., "not a number") and only the linear uncertainty propagation results are returned. For typical inputs of routine analyses 510 with a few percent relative uncertainty (Sect. 5), the impact of this phenomenon is negligible.

Appendix C: Uncertainty in date as a function of input uncertainties
We examined the overall behavior of date uncertainty from 0 to 4.6 Ga as a function of relative input uncertainties of 1%, 5% and 20% on 4 He (Fig. C1a), radionuclides (Fig. C1b) and isotope-specific FT values (Fig. C1c). This range of dates was generated by fixing the 238 U and 232 Th values while varying 4 He values (no 147 Sm was included because of its 515 generally negligible influence on apatite and zircon results). Th/U ratios representative of a typical zircon (based on the Fish Canyon Tuff zircon reference standard), a typical apatite (from a compilation of apatite data; Sect. 5), and the Durango apatite reference standard (0.6, 1.25, and 16.1, respectively) were used. For all calculations, an isotope-specific FT value of 0.7 was applied to all isotopes to permit comparisons between raw and FT-corrected dates (while isotope-specific values will differ in real data, we simplify these to a single value here). We initially explored the influence of individual uncertainties on 520 the date by varying the relative uncertainty of one input parameter ( 4 He, radionuclides, or FT) while fixing all other uncertainties at 0 (Fig. C1). We then evaluated how combinations of input uncertainties can influence the date (Fig. C2), although this is more fully evaluated in practical terms using real data, as in Sect. 5. For these exercises, we use the results from Monte Carlo uncertainty propagation, as this technique is in theory fully accurate (see Appendix E for further discussion). We used a constant number of simulations set at 10 8 to provide precise estimates of skew and comparisons 525 between the Monte Carlo and linear uncertainty propagation methods. This number of simulations corresponds to a minimum precision of the mean date of ~0.0002% (2 µg/g). For individual input uncertainties, at young dates the input and output relative uncertainties are similar. If all uncertainty is in the 4 He value or correlated FT values, the relative date uncertainty is equivalent to the input uncertainty at zero age (Fig. C1a&c). For uncertainty in the radionuclides, the relative date uncertainty at zero age is approximately 80% 530 the magnitude of the relative input uncertainty (a 4:5 ratio). The exact scaling between input and output uncertainties is dependent on the Th/U ratio (Fig. C1b&c). The relative date uncertainty decreases with increasing absolute date for constant relative input uncertainties. For example, while uncertainty in 4 He has a one-to-one relationship with date uncertainty at zero age, at 4.6 Ga the date uncertainty is approximately half that of the input 4 He uncertainty (Fig. C1a). The same phenomenon is observed for uncertainty in radionuclides and FT (Fig. C1b&c). The relative extent of decreasing uncertainty as a function of increasing 545 date is dependent on Th/U ratio and is independent of the magnitude of input uncertainty (i.e., the three vertically stacked panels in Fig. C1a, 3b, and 3c are identical aside from the scale of their y-axes).
Decreasing uncertainty with increasing date is also observed for multiple input uncertainties. Figure C2 illustrates examples of combining uncertainties with differing magnitudes in quadrature. At zero age, the uncertainty on the date introduced by each input parameter combines roughly in quadrature to provide the uncertainty on the date. For example, at 550 zero age, a 5% uncertainty in all parameters ( 4 He, radionuclides, correlated FT uncertainty), each of which alone introduces a date uncertainty of 5%, 4%, and 5%, respectively, together yields a date uncertainty of 8.1% ( √ 0.05 2 + 0.04 2 + 0.05 2 ≅ 0.081; solid line, Fig. C2b). Alternatively, if input uncertainties have highly differing magnitudes, the larger uncertainty will dominate and the resulting combined uncertainty will be approximately equal to the larger uncertainty. As an example, a 10% and 1% uncertainty combined in quadrature will result a 10.05% uncertainty. This behavior suggests that reducing the 555 magnitude of the largest input uncertainty will be the most effective means of reducing overall date uncertainty. The phenomenon of decreasing uncertainty with increasing date is caused by the "roll over" of the helium ingrowth curve due to its nature as an exponential function. Because of this roll over, constant uncertainty in the independent variable (i.e., 4 He or the radionuclides) will correspond to smaller uncertainty in the dependent variable (the date) as the value of the 565 dependent variable increases. Fig. C3 is a schematic showing log plots of date vs. 4 He and date versus 238 U that illustrates this phenomenon. For young dates, the exponential term in the date equation (Eq. (2)) approaches zero, meaning that the relative uncertainties input to the He age equation will be roughly reflected in the output uncertainties (Sample 1, Fig. C3). For older dates, this exponential term becomes increasingly large, resulting in roll over of the ingrowth curve and reducing the date uncertainty relative to the inputs (Sample 2, Fig. C3). The exact form of this roll over is dictated by the relative 570 abundance of each radioisotope, resulting in the variations observed in Fig. C1 for differing Th/U values. (Sample 1, orange line) and one old (Sample 2, green line) are provided with gaussian constant percent relative uncertainty (1σ, depicted by the shaded region). The apparent asymmetry in the uncertainty along the x-axis is a result of the logarithmic plot. The non-linearity of the (U-Th)/He age equation is exaggerated for this schematic by decreasing the uranium decay constant to improve visibility of its effects. Note that in log-log space, the spread (i.e., uncertainty) in the x-axis is constant for constant input uncertainty, but the resulting uncertainty on the y-axis shrinks with increasing date. 580

595
The magnitude of skew correlates directly with the magnitude of input uncertainty (Fig. D1-D2). For low percent input uncertainties on all parameters, the magnitude of skew is low. For example, uncertainties of 1% for all inputs yield ≤2%0.06 skewness for dates from 0 to 4.6 Ga ( Fig. D2a-c, top panels). Only when the input uncertainties are larger does the effect of skew on the dates become substantial (Fig. D2a-c, middle and bottom panels). In the case of larger uncertainty in 4 He (Fig. D2a, middle and bottom panels), skew increases from zero to progressively larger negative values at older dates. 600 The inverse is true for uncertainty in the radionuclides and FT; skew is highest when uncertainty in these parameters is high for young dates and decreases with increasing age (Fig. D2b-c, middle and bottom panels). Note that although asymmetric uncertainties as high as ~40% a skewness of 0.85 can be yielded by radionuclide uncertainties of 20%, such large uncertainties are anomalous and do not typify most high-quality (U-Th)/He datasets (Sect. 5). When uncertainty is included in multiple input parameters, the overall skew is a combination of the skew resulting from individual input uncertainties (Fig. D3). Unlike date uncertainty, which combines individual inputs in quadrature, the 615 combination of skew from individual inputs does not follow an easily predictable trend. For input uncertainties of 1% and 5% for 4 He and the radionuclides only, the skew is largest at zero age (~1% and 5%, respectively), and declines with increasing age (dashed lines in Fig. D3). Because uncertainty in 4 He generates negative skew at older dates (Fig. D2a), the skew from these combined uncertainties becomes negative at dates ⪆3 Ga as the skew resulting from 4 He uncertainty overwhelms the skew from radionuclides, which has the opposite sign and is greatest at young dates (Fig. D2b). Similarly, 620 for input uncertainties of 1% and 5% for all parameters (including FT), the skew is largest at zero age (~2% and ~11% respectively) and declines with decreasing age (to ~0% at 4.6 Ga; solid lines in Fig. D3). Figure D3: The skew in the date probability distribution resulting from combining multiple input uncertainties. Skew is shown as a percent difference between the 68% confidence intervals. This figure shows each possible combination of (a) 1% 625 uncertainty and (b) 5% uncertainty in 4 He, radionuclides, and fully correlated FT (r = 1) for a Th/U ratio of 1.3 (the green curve, derived from the Fish Canyon Tuff zircon reference standard).

Appendix E: Comparison of linear and Monte Carlo uncertainty propagation results
To compare linear and Monte Carlo error propagation derived uncertainties, we average the two 68% confidence 630 intervals to determine uncertainty from both methods at the 1σ level. For data with high skew, this method provides a means of comparing the scale of these two differing output distributions directly. The magnitude of the error in uncertainty estimation from linear uncertainty propagation due to nonlinearity in the date equation is proportional to the magnitude of the input uncertainties. As shown in Fig. E1a, for uncertainty in He alone, the Monte Carlo and linear methods yield identical results at younger dates, with linear uncertainty propagation beginning to underestimate the true uncertainty values at older 635 dates as the absolute magnitude of 4 He uncertainty increases (reaching a maximum discrepancy of ~2% for input uncertainties of 20% at 4.6 Ga). Uncertainty in radionuclides and FT have the opposite effect; the discrepancy between the Monte Carlo and linear methods is greatest (~3% for input uncertainties of 20%, dependent on Th/U ratio and correlation in FT uncertainties) at zero age and decreases with increasing date. This small extent of error indicates that Monte Carlo and linear methods are in general agreement. 640 Figure E1: The percent error introduced by using linear uncertainty propagation instead of Monte Carlo uncertainty propagation for dates from 0 to 4.6 Ga. Error is shown for uncertainty in (a) 4 He, (b) radionuclides, and (c) FT using the percent difference between the Monte Carlo and linear uncertainty propagation results, with the average 68% confidence intervals used to represent a single uncertainty value for the Monte Carlo results. Relative input uncertainties of 1% (top panels), 5% (middle panels), and 20% 645 (bottom panels) were applied. The line colors correspond to the Th/U ratio for typical zircon (0.61, green dashed curve, derived from the Fish Canyon Tuff zircon reference standard), for typical apatite (1.25 orange solid curve; derived from apatite data compilation), and for the Durango apatite reference standard (16.1, blue dash-dotted curve). Note that unlike Fig. C1, the y-axis scale is different for each panel.
As linear uncertainty propagation relies on an arithmetic calculation rather than random sampling, this method 650 provides predictable and repeatable results for uncertainty calculations and is amenable to encoding in spreadsheet programs, facilitating the inclusion of the equations provided in Sect. 3.2 in existing spreadsheet-based workflows. However, the presence of skew in (U-Th)/He date uncertainties and the inaccuracies in uncertainty calculation induced by non-linearity in the (U-Th)/He age equation indicate that the more accurate Monte Carlo uncertainty propagation method is more universally applicable. Although in the past computational (in)efficiency was generally considered the weakness of Monte Carlo 655 methods, running as many as one million Monte Carlo simulations in HeCalc takes less than one second on a modern computer for a typical sample. This number of random samples provides a sufficiently large population that output histograms are relatively smooth, yielding accurate uncertainty calculations with sufficient significant figures that the modelto-model variation induced by random sampling is negligible. As the Monte Carlo method in HeCalc is not excessively computationally intensive and provides both skew and accurate uncertainties, we suggest that the Monte Carlo method is 660 preferable to linear uncertainty propagation in the (U-Th)/He system.