Technical note: “U–Th Analysis” – open-source software dedicated to MC-ICP-MS U-series data treatment and evaluation

We present our stand-alone data analysis application for ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ dating with multi-collector inductively coupled plasma mass spectrometry (MC-ICP-MS). The Python-based algorithm is equipped with a graphical user interface (GUI) and comprises raw data treatment, corrections, age calculation, and error estimation. Our underlying measurement protocol employs a combination of Faraday cups (FCs) and secondary electron multipliers (SEMs), and the software allows for different detector layouts for the measurement of the least abundant isotopes ²³⁴U, ²³⁰Th, and ²²⁹Th. We especially focus on features that ensure reproducibility and enable user-friendly reanalysis of measurements such as customized calculation constants with templates. Result files are saved automatically and contain all relevant settings used. We demonstrate the relevance of adequate data outlier treatment and generally recommend using the median instead of the mean of calculated ratios. The performance of our evaluation software is demonstrated by a case study from a Puerto Rican stalagmite with growth phases from the modern era to 40 kyr old. The majority of the obtained ages reach uncertainties in the range of 0.3 %–0.6 %, underlining the capability of our measurement protocol.

1 Introduction

The U-series disequilibrium method, ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ dating, is a precise chronometer covering approximately the last 650 kiloyears and has proven to be indispensable for the age determination of marine and continental carbonate archives and their applications (Bourdon et al., 2003). The method is based on a complete disequilibrium of ²³⁴U, with its daughter nuclide ²³⁰Th, during the formation of secondary carbonates. It presumes a subsequent closed-system evolution of the activity ratio of (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{234}}\mathrm{U}$) and (${}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$) since the time of formation. Ideally, the initial ²³⁰Th activity of the material is presumed to be zero or can be estimated from the total Th concentration via an initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio. The dating applications for secondary carbonates and other appropriate materials are manifold in geochemistry, archeology, and climate science. Further development of this dating method includes improvements in both instrumentation and measurement protocols, as well as reproducible data analysis and age calculation schemes (Pourmand et al., 2014; Andersen et al., 2004; Cheng et al., 2013; Breton et al., 2015; Chiang et al., 2019; Hellstrom, 2003; Hoffmann et al., 2007; Shen et al., 2002; Shen et al., 2012; Kerber et al., 2023; Shao et al., 2019). The presently most sensitive and precise technology for high-precision U and Th isotope measurements is multi-collector inductively coupled plasma mass spectrometry (MC-ICP-MS). Recent technological advances of MC-ICP-MS include the implementation of high ohmic amplifiers, allowing the enhancement of the dynamic range of multiple Faraday cups (FCs) to 6 orders of magnitude for the simultaneous detection of very large and low isotope abundances, instead of the conventionally used combination of secondary electron multipliers (SEMs) and FCs (Breton et al., 2015). Measurement protocol updates aim at increasing measurement precision and/or decreasing input sample masses by combining new detector layouts, improving the understanding of correction factors, and ensuring a stable measurement environment (Cheng et al., 2013; Chiang et al., 2019; Shen et al., 2002; Shen et al., 2012; Hellstrom, 2003; Andersen et al., 2004; Hoffmann et al., 2007; Kerber et al., 2023; Shao et al., 2019).

We focus here on the third route for the enhancement of ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ dating, which is clear and reproducible data analysis and age calculation schemes. Up to now, only two ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ dating data analysis routines have been published (Shao et al., 2019; Pourmand et al., 2014). However, regarding the rising amount of data being produced in MC-ICP-MS laboratories, data management is becoming more and more important. For example, some samples might require later adaptation of the individual corrections of isotope ratios due to residual contamination with non-carbonate material or detection of initial ²³⁰Th from the carbonate-forming environment.

Dating young materials of only a few years to centuries in age is challenging due to the small number of counts on ²³⁰Th in particular, which implies that all correction factors including “ghost signal” corrections need to be determined very precisely (Zhao et al., 2009; Chiang et al., 2019; Kerber et al., 2023). Regarding the removal of scatter ions on the specific low-abundance masses 230 and 229 amu, Kerber et al. (2023) demonstrated an effective correction based on a linear dependence of the scattered ions on the ²³⁸U signal. Other authors separate U and Th chemically to reduce or remove the ²³⁸U beam from the low-abundance thorium isotope measurements (Chiang et al., 2019), which implies flexibility in the detector arrangement and data treatment protocol. As such scatter peaks may depend on the specific instrument or vary through time, these corrections need to be adaptable constants in the data evaluation routine. The influence on the final atomic ratio and accuracy of the ghost signals as well as typical variation in other individual corrections, such as peak tailing, mass fractionation, and isobaric interferences, are evaluated in detail in Kerber et al. (2023).

In addition, the correction for initial Th may cause large age corrections and propagated uncertainties, in particular since adequate initial Th values based on the ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$ ratio may be variable and difficult to detect (Hellstrom, 2006; Wenz et al., 2016; Wortham et al., 2022). There are different methods to estimate the initial Th isotope ratio: first, isochrons can be used to determine the isotopic composition of the detrital component in the carbonate (Ludwig and Titterington, 1994; Wenz et al., 2016; Stinnesbeck et al., 2020; Töchterle et al., 2022). Secondly, analyses of modern drip waters or recent carbonate deposits allow estimation of the value and sources of initial Th (Wortham et al., 2022; Li et al., 2022). In some cases the “true” age of a stalagmite can also be inferred from other dating methods, such as radiocarbon (Akers et al., 2019; Huang et al., 2024) or the stratigraphic order (Hellstrom, 2006). Also, several approaches can be combined (Warken et al., 2020; Akers et al., 2016; Roy-Barman and Pons-Branchu, 2016).

Another aspect involves updating half-lives (such as from Cheng et al., 2000b, to Cheng et al., 2013), which makes re-evaluation of previously measured data necessary. These tasks are error-prone, in particular when they require copying and pasting data in, e.g., spreadsheets. Also, a clear and unified documentation of the applied constants and the way of saving data is desirable. Additionally, the statistical methods, for example for outlier correction, should undergo clear documentation. Altogether, this helps to report $\mathrm{Th}/\mathrm{U}$ ages in a standardized way (Dutton et al., 2017).

In this study, we present our user-friendly graphical user interface (GUI) and the underlying algorithm for data treatment and age calculation. The software is currently optimized for Thermo Fisher Neptune MC-ICP-MS instruments, but the open-source code in principle allows adaptations to other setups and instruments. Methods to treat outliers in measurement data are particularly highlighted. As a case study, we present newly obtained ages from a stalagmite from Larga Cave, Puerto Rico, which shows a modern growth phase, as well as continuous deposition during the last glacial into the deglaciation, thus demonstrating the performance of our method for both very young and older sample materials. Our protocol enables a precise determination of speleothem growth rates, which allows a comparison to a coevally deposited stalagmite from the same cave, highlighting the influence of in-cave processes on speleothem growth rates. In particular, this dataset showcases how initial ²³⁰Th correction models can be easily tested with our software and GUI presented here and how those influence speleothem chronologies.

2 Methods

2.1 Standards and reference materials

We use our in-house triple-spike solution (TriSpike) with a ²³³U concentration of (0.038556 ± 0.0000009) ng g⁻¹, a ²³⁶U concentration of (3.86778 ± 0.00009) ng g⁻¹, and a ²²⁹Th concentration of (0.018055 ± 0.000008) ng g⁻¹ (2 standard error of the mean) (Kerber et al., 2023). For standard bracketing, we employ Harwell uraninite 1 (HU-1) as a reference material. Its activity ratios (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$) and (${}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$) are presumed to be 1, as it is a secular equilibrium material. Abundance sensitivity and hydride correction are determined by measuring CRM-112A U reference solution and an in-house ²³²Th standard. The CRM-112A gravimetric standard solution has a ²³⁸U concentration of (4.3021 ± 0.0015) µg g⁻¹, while the in-house ²³²Th standard calibrated with TriSpike has a ²³²Th concentration of (505.8 ± 1.02) ng g⁻¹ (2σ uncertainties). CRM-112A solution is also used to track the values of the two ghost signal constants, k₂₂₉ and k₂₃₀ (Kerber et al., 2023). For k₂₂₉, it is measured without the addition of TriSpike, while in the case of k₂₃₀, the spiked CRM-112A solution is employed. For age determination, the ²³⁰Th and ²³⁴U decay constants determined by Cheng et al. (2013) are used. Ages are reported with 2σ statistical standard mean error but do not include half-life uncertainties.

2.2 Chemical preparation and instrumentation

The chemical preparation of carbonate samples includes sample dissolution in ultra-clean nitric acid, spiking with TriSpike, and two steps of wet-column chromatographic ion exchange separation of U and Th from matrix elements using Eichrom UTEVA resin (Douville et al., 2010; Wefing et al., 2017; Matos et al., 2015). Chemical blanks are commonly below 0.4 fg for ²³⁴U and 0.04 fg for ²³⁰Th, and Ca matrix concentrations are required to be below 10 ppm. For the mass-spectrometric measurement, samples are dissolved in 1 % HNO₃ and 0.05 % HF. All samples were measured by MC-ICP-MS (Thermo Fisher Neptune Plus) at the Institute for Environmental Physics, Heidelberg University (Germany). The desolvating system CETAC Aridus II is used as an inlet.

The mass spectrometer is equipped with Faraday cups (FCs) and a central secondary electron multiplier (SEM). The central detector can be selected between the SEM and an FC connected to a 10¹³ Ω amplifier. ²³⁸U is measured on a 10¹⁰ Ω amplified resistor. All other FCs are connected to 10¹¹ Ω amplifiers. The cup setting to collect isotope signals on masses 238 to 229 is shown in Table 2 in Kerber et al. (2023). The first cycle collects all U isotopes for 2 s, with ²³⁴U on the central detector (FCs and SEM). The second and third cycles collect the Th isotopes for 2 s integrations time, with ²³⁰Th and ²²⁹Th on the central SEM. These cycles are repeated for an optimal number for each measurement.

A measurement sequence starts with the determination of abundance sensitivity and tailing on two different solutions for both uranium and thorium. Each sample and standard measurement is preceded by a procedural blank measurement to ensure that the background signal has gone back to a clean state. CRM-112A measurements are carried out to track ghost signal values at the beginning and end of a measurement sequence. Samples are bracketed with HU-1 as a reference material. Samples, standard, and procedural blanks are measured with the same configuration.

Mass fractionation (or mass bias) is corrected via the natural ratio of ${}^{\mathrm{235}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$ due to the lack of natural Th isotopes. In our setup, the ratio of the artificial isotopes ²³³Uand²³⁶U in the spike is monitored to double-check the mass bias correction. The calibration of FC gain and SEM yield is described in detail in Kerber et al. (2023): while there is an internal electronic calibration function for the calibration of 10¹⁰ Ω and 10¹¹ Ω amplifiers, the 10¹³ Ω amplifier in our setup is calibrated manually. For this, the gain factor is determined regularly by measuring ²³⁵U alternately on the 10¹³ Ω and 10¹¹ Ω amplified cup. In an analogous manner, the SEM yield is routinely determined by measuring ²³⁵U alternating on the SEM and on a 10¹¹ Ω FC at a signal intensity of ∼ 5 mV. Since HU-1 standards are measured with the same detector configuration in standard bracketing mode, the observation of the (${}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$) values of HU-1 measurements allows monitoring and manual optimization of gain and yield values for each measurement sequence in the data analysis scheme.

In its current version, the GUI is written for this type of measurement protocol. Many adaptions, such as fewer procedural blank measurements or other isotope ratios for mass fractionation correction, require small changes in the code but are easily feasible for users sufficiently competent in Python. Numerous MC-ICP-MS instruments, measurement protocols, and cup configuration settings can be used for such isotope measurements.

Here, we have developed our data treatment and GUI for a quasi-simultaneous measurement of U and Th isotopes as detailed previously (Kerber et al., 2023). We do not need to specify the number of cycles or the exact scheme of the sequence as the software will extract this information from the data. However, the cup configuration matters. Our method implies that U and Th isotopes are present in the final solutions, which requires a U isotope scattering correction. This is not necessarily needed when measuring U and Th isotope solutions independently. Nevertheless, the open-source software presented here can easily be adopted to other cup configurations or raw data outputs from different instruments.

2.3 Speleothem sample description

Stalagmite B1 was collected in 2019 in Larga Cave, Puerto Rico (18°19^′ N 66°48^′ W; 350 m.s.l.; Fig. S1a in the Supplement), from a passage in the deep part of the cave connected to the “collapse room”. The host rock overburden at the location of the sample is about 40–60 m. It is in total 60 cm long and has an average diameter of 15 cm (Fig. S1b). The drip site was still active and was monitored with spot measurements over several years, revealing varying drip intervals between 2 and > 120 s. A total of 37.7 mL of water from the drip site of stalagmite B1 was analyzed for its U and Th activity ratios. Samples for ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ dating of the speleothem with typical input masses of 100–150 mg were cut using a diamond-wire saw along the growth axis. Chemical preparation, mass-spectrometric measurements, data treatment, and evaluation of drip water and the speleothem samples followed the methods described in Kerber et al. (2023) and in this study.

Larga Cave is located in the north central karst region of Puerto Rico (Fig. S1a). Previous work including extensive cave air and drip monitoring has demonstrated that the cave is a valuable location to study of the influence of changing climate on past rainfall patterns in the western tropical Atlantic (Vieten et al., 2018a; Warken et al., 2022b; Vieten et al., 2018b). In particular, the main passage of Larga Cave is subject to seasonally varying ventilation, which results in pCO₂ values of 600 ppm close to atmospheric values during winter and higher values up to 1800 ppm in summer (Vieten et al., 2016). In contrast, in the deep part of the cave, where stalagmite B1 was also collected, ventilation is strongly muted, and cave-air pCO₂ values are higher, with values up to 2300–3600 ppm (Vieten et al., 2016). As a result of this ventilation regime, growth rates are expected to vary both seasonally and between different locations inside the cave (Vieten and Hernandez, 2021). So far, two speleothem records from Larga Cave have been published; the most recent covers the past 500 years (Vieten et al., 2024), and the second stalagmite grew during the period of 46.2–15.3 ka with a hiatus from 41.1 to 35.5 ka (Warken et al., 2020). For ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ dating of Larga speleothems, high initial Th contents have to be considered – a phenomenon that regularly occurs in speleothem records from the Caribbean and Central American region (Fensterer et al., 2010; Steidle et al., 2021; Moseley et al., 2015; Schorndorf et al., 2023; Stinnesbeck et al., 2020; Beck et al., 2001; Akers et al., 2016; Rivera-Collazo et al., 2015).

3 Data treatment and analysis procedures

The whole analysis procedure from raw data treatment to age calculation is conducted in one GUI featuring three tabs: Input for isotopic ratio calculations, Inspect for outlier correction of the signal, and Analysis for age calculation. The source code is accessible at https://github.com/puahd/UTh_Analysis (last access: 11 November 2024). It is based on the open-source PyQt5 Python library (https://pypi.org/project/PyQt5/, last access: 11 October 2024). To execute the GUI, the user has to run the file “main.py”. The folder https://github.com/puahd/UTh_Analysis/dist (last access: 11 November 2024) also contains a compiled .exe file for the GUI (“U–Th Data Analysis.exe”) as well as default configuration files (“constants – coral.cfg” and “constants – stalag.cfg”. The input and output format of files is .csv or .xlsx. The GUI consists of three consecutive tabs, the functionalities and the underlying calculations and processes of which will be described in the following.

3.1 Input tab

In Input, as presented in Fig. 1, the user can navigate to the folder containing the raw mass spectrometer data and start the calculation of corrected isotopic ratios (Box 1 in Fig. 1). All tab screenshots present data from stalagmite B1. Prior to the calculations, a configuration file containing all necessary constants used in the calculations needs to be loaded (also Box 1). This file contains constants and correction factors used for evaluation of the activity ratios and ages, such as mass fractionation coefficients, decay constants, the exact masses of the isotopes, and the values applied for initial ²³⁰Th correction model. All constants can be edited manually either in the configuration file directly or within the GUI using the button “edit”. An exemplary configuration table is also provided in the Supplement (Fig. S2). To apply a ²³⁰Th correction model a value can be set for the activity ratio and uncertainty of the contaminating material (“A²³⁰Th232Th Init.”). The conventional approach to account for initial Th would be an activity ratio of 0.75 ± 0.38, assuming an upper continental crust ${}^{\mathrm{232}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ weight ratio of 3.8 (Taylor and McLennan, 1985) with an uncertainty of 50 % (Ludwig and Paces, 2002) and ²³⁰Th, ²³⁴U, and ²³⁸U in secular equilibrium for the detrital material. Exemplary templates for corals and speleothems with conventionally used correction models are provided. For speleothems, a typical activity ratio of (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$)_ini/detr of detritus is estimated to 0.75 ± 0.38 (see above). Nevertheless, this ratio may require adjustment according to local conditions. The coral template assumes as a default value an activity ratio of 8 ± 4, which is estimated for corals dwelling in waters of the northeast Atlantic upper thermocline (Wefing et al., 2017). For one data series, only one correction constant, the (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio of the contamination, can be added to the calculation. Hence, in the case that several factors need to be explored, the data series requires repeated treatment.

Figure 1 shows the layout of the GUI Input tab. Once the constants are implemented and the input data are selected, it is optional to choose an output path to store the analysis output (Box 1). If no path is specified, the results will be stored in the raw data folder. When clicking the “settings” button next to the output path (highlighted in red in Box 1), a menu opens in which the following parameters of the sample can be noted: denomination, type of archive, lab numbers, geographic origin, and a general description. The first and last laboratory numbers are automatically read out from the raw data. The final output result files will then be saved in a newly created folder under the name _[labnumber₁-labnumber_n] in the directory chosen before. The metadata information transferred through the GUI dialogue window is stored in a .json file in the respective folder. In the “custom constants” panel (Box 2), some settings can be selected: for example, if the blank has already been subtracted in the mass-spectrometric software or not. Next to this panel, an overview of the files read in from the folder is shown. After running the evaluation script with the loaded data and adjusted settings (Button “run”, highlighted in red in Box 1), the results of tailing and hydride correction, respectively, as well as the calculated ratios are displayed in the tab (Box 3). In addition, four excel .xlsx output files are created by the software at this stage and stored in the directory path folder: Ratios.xlsx, Tailing.xlsx, PrBlank.xlsx, and Intensities.xlsx. Ratios.xlsx contains all calculated ratios and their errors as also presented in the GUI (Box 3). Tailing.xlsx summarizes the U and Th tailing values (in counts per second (cps) per volt (V) ²³⁸U) for each mass. In PrBlank.xlsx, the average values for each mass of the procedural blank measurements before each standard and sample are presented. Intensities.xlsx contains the full data tables, with the signals in cps or V for each mass over all cycles. Every standard or sample has its own sheet.

Figure 1 Screenshot of the Input tab of “U–Th Analysis”: (1) in the top part the data folder is selected (“directory path”), and the “constants” file (“constants”) can be loaded (“load”), edited (“edit”), or created (“new”). In addition, it is possible to set an “output path”. Red boxes show the “settings” button to enter metadata for saving, as well as the “run” button to start the analysis. Box (2) shows the custom “constants” box as well as the file overview for the selected folder. In Box (3), the plots on the top left show the interpolated tailing. On the top right, numerical values of U and Th tailing and hydride correction are presented. The calculated ratios are shown in the bottom panel.

The algorithm of the Input tab starts by reading in the .exp measurement files for sample and standard measurements, process blank (instrumental background), and uranium and thorium abundance sensitivity measurements. The lines for all cycles for all isotopes are imported into a pandas data frame. Firstly, matrices for tailing, hydride, and process blank correction are produced that are later subtracted from the isotopic masses used for ratio building. The individual steps are carried out as follows.

• Tailing. Uranium tailing is determined by measuring the off masses 228.5, 233.5, 236.5, 236.7, 237.05, and 237.5 before a measurement sequence starts. The first half-mass can be changed between 228.5 and 227.5 as we observed a scatter peak around this mass that switched its exact position every few months. Thorium tailing off masses are 227.5, 228.5, 229.5, 230.5, and 231.5. For interpolation to full masses, we use piecewise cubic Hermite-interpolating polynomial fits (Kerber et al., 2023). The masses that undergo ²³⁸U tailing correction are ²³³U, ²³⁴U, ²³⁵U, ²³⁶U, ²²⁹Th, ²³⁰Th, and ²³²Th, while ²³²Th correction is applied to ²²⁹Th and ²³⁰Th.

• Hydride isobaric interference. Hydride correction is determined by measuring 239 amu for UH⁺ and 233 amu for ThH⁺ during the abundance sensitivity measurements. The instrumental background (or memory) is referred to here as a process blank. It is measured between all sample and standard measurements for 70 s. Typical blank levels afterwards are 0.5 cps for ²³⁰Th and 6 cps for ²³⁴U. The matrices from these three corrections are then used for data reduction of each isotope.

• Detector setting. Three main different detector layouts are possible and are detected automatically by the software: (1) all isotopes on the cup; (2) ²³⁴U, ²³⁰Th, and ²²⁹Th on the SEM; and (3) ²³⁴U on the FC as well as ²³⁰Th and ²²⁹Th on the SEM. In normal operation, options (2) and (3) are used, depending on the ²³⁴U concentration of the respective samples. ²³⁴U signals above 2 mV are measured on the center FC, which is the case for the absolute majority of samples.

The sample and standard data, corrected for yield, gain, and mass fractionation, are now used for the calculation of all relevant isotopic ratios followed by subsequent outlier tests, as described in Sect. 3.2.

3.2 Inspect tab

Following the initial raw data treatment in the Input tab, the Inspect tab (presented in Fig. 2) allows visualizing and retreating the data prior to final age calculation. In particular, the settings for the outlier test can be adapted.

Figure 2 Screenshot of the Inspect tab of “U–Th Analysis”: (1) table of ratio results (from Input tab), (2) overview of measurement files in the folder, (3) metadata of a selected file and signal over the measurement cycle number for one isotope (which can be selected in step 4), and (4) option selection panel for the signal plotting.

The Inspect tab allows the user to plot the signal data points over the measurement cycle number for all isotopes in the individual measurement files of the sequence. At the top of the tab (1), the table of ratio results from the Input tab is presented. On the left (2), the list of measurement files (.exp) is shown. By clicking on a specific file, the metadata and the signal plotted over the measurement cycle number are presented (3). On the bottom left (4), four dropdown menus are available: the first one, “Isotope”, allows selecting one isotope from all of the isotope species measured. “Mean” offers switching between the mean and median of the signal. The “Deviation” menu provides three options for the assessment of data dispersion: standard deviation, median absolute deviation, and interquartile range. By setting “Scaling” to absolute or relative, the y axis of the plot on the right can be changed between signal intensities in V or cps and relative values. Any selection in the dropdown menus leads to an automatic update of the plot on the right. The mean and median, as well as the dispersion ranges, are presented as dashed blue lines. Data points outside of the dispersion range are marked in red as outliers.

The “Calculate ratios” button provides the option to recalculate the ratios using the updated mean and deviation selection for all isotopes. The default settings are median and standard deviation. However, these updated options are then used to exclude outliers from the ratio arrays, not the signal intensity arrays themselves. This means that exactly the same data points are not necessarily marked as outliers in the signal intensity plots and will be excluded but rather the ones for which signal ratios of two isotopes are outside of the accepted deviation range. The option selected in the Mean menu will then also be used to calculate the average of the isotope ratios. The method of calculating the uncertainty of outlier-corrected isotopic ratios via the standard error, however, is fixed. In total, the software provides three different options for dispersion, including (i) the standard deviation (s), (ii) the interquartile range (IQR) (Tukey, 1977), and (iii) the median absolute deviation (MAD) (Leys et al., 2013; Huber, 2004; Rousseeuw and Croux, 1993). For the calculation of the MAD we assume normally distributed data. The treatment of the means of ratios may have undesirable statistical properties for low or fast-changing signals (Ogliore et al., 2011; McLean et al., 2016), which could be taken into consideration when updating the software.

3.3 Analysis tab

In a last step, age calculation is carried out in the Analysis tab presented in Fig. 3. Here, additional input data are necessary from the sample weight tables (1). There are several ways to import these tables: either by clicking “load” and navigating to the respective folder or by manually creating the table directly in the GUI (“create”). An exemplary weight table is provided in the Supplement (Fig. S2). In the panel “Metadata history”, the previously loaded sample weight tables in the directory path folder are shown and can be directly imported. “Start analysis” starts the data analysis and calculates the ages. Outputs are both presented in the GUI as a result table and stored in a Results.xlsx file. In the case that an output path was specified, Results.xlsx is created in both the output and the directory path folder. If the output path is missing, the file is only saved in the directory path folder. If an output directory has been created for specific lab numbers, all following analysis of these same files will be written to the same output directory but not overwrite earlier Results.xlsx. The Results.xlsx has five sheets: Inputs, Calc, Results, Constants, and Options. Inputs presents sample weight and metadata as well as the calculated ratios. In Calc, all steps of the age calculation such as concentrations and activity ratios are shown. Results is a summary of the most important calculation steps and final age values, and the same table as is presented in the GUI as the table of results in Fig. 3. Constants contains the whole list of values from the (potentially edited) .cfg file. In Options the average and dispersion measure options are stored.

Figure 3 Screenshot of the Analysis tab of “U–Th Analysis”: (1) load sample weight tables (metadata files). The bottom panel lists the history of previously loaded tables. The button highlighted in red starts the analysis to calculate ages (“Start analysis” button). The panel in Box (2) displays the table of results.

The equations for activity ratios to calculate ages are implemented according to Ivanovich and Harmon (1992), with the following.

$$\begin{array}{}\text{(1)}& {\displaystyle }\left({\displaystyle \frac{{}^{\mathrm{234}}\mathrm{U}}{{}^{\mathrm{238}}\mathrm{U}}}\right)\left(t\right)=\left({\left({\displaystyle \frac{{}^{\mathrm{234}}\mathrm{U}}{{}^{\mathrm{238}}\mathrm{U}}}\right)}_{\mathrm{init}}-\mathrm{1}\right)\cdot e^{-{\mathit{\lambda }}_{\mathrm{234}}\cdot t}+\mathrm{1}\text{(2)}& {\displaystyle }\begin{array}{rl}\left({\displaystyle \frac{{}^{\mathrm{230}}\mathrm{Th}}{{}^{\mathrm{238}}\mathrm{U}}}\right)& =\mathrm{1}-e^{-{\mathit{\lambda }}_{\mathrm{230}}t}+{\displaystyle \frac{\mathit{\delta }{}^{\mathrm{234}}\mathrm{U}}{\mathrm{1000}}}\\ & \cdot \left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{230}}}{{\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{234}}}}\right)\cdot \left(\mathrm{1}-e^{-({\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{234}})t}\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(3)}& \mathit{\delta }{}^{\mathrm{234}}\mathrm{U}=\left({\left({\displaystyle \frac{{}^{\mathrm{234}}\mathrm{U}}{{}^{\mathrm{238}}\mathrm{U}}}\right)}_{\mathrm{meas}}-\mathrm{1}\right)\cdot \mathrm{1000}\left(\mathrm{‰}\right)\end{array}$$

To obtain ages corrected for initial and detrital ²³⁰Th, the ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ activity ratio used in Eq. (5) is corrected using the initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$)_ini/detr ratio and

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\left({\displaystyle \frac{{}^{\mathrm{230}}\mathrm{Th}}{{}^{\mathrm{238}}\mathrm{U}}}\right)}_{\mathrm{corr}}& ={\left({\displaystyle \frac{{}^{\mathrm{230}}\mathrm{Th}}{{}^{\mathrm{238}}\mathrm{U}}}\right)}_{\mathrm{meas}}-{\left({\displaystyle \frac{{}^{\mathrm{232}}\mathrm{Th}}{{}^{\mathrm{238}}\mathrm{U}}}\right)}_{\mathrm{meas}}\\ & \cdot {\left({\displaystyle \frac{{}^{\mathrm{230}}\mathrm{Th}}{{}^{\mathrm{232}}\mathrm{Th}}}\right)}_{\text{ini/detr}}\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{230}}}{{\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{234}}}}\right)\cdot e^{-{\mathit{\lambda }}_{\mathrm{230}}\cdot t}.\end{array}\end{array}$$

These equations need to be solved numerically. For the determination of age uncertainty, the usual approach is to repeat the numerical determination of the age for several thousand runs in a Monte Carlo simulation while randomly sampling the input ratios from a normal distribution, with μ corresponding to the ratio's value and σ corresponding to the uncertainty in this parameter.

4 Example dataset: stalagmite B1

To demonstrate our data evaluation tool, we present newly obtained ages of stalagmite B1 from Larga Cave, Puerto Rico. The results of activity ratios and calculated ages can be accessed in Table S1 in the Supplement. Analysis of the speleothem samples reveals moderate U concentrations in the range between 300 and 600 ng g⁻¹ and minor detrital ²³²Th contamination, with (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios of typically > 300. However, in both the top 20 mm and around 450 mm distance from the top (dft) lower (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios of ca. 40–125 are measured. U isotopic composition varies between 450 ‰ and 640 ‰ of δ²³⁴U values. Uncertainties of the uncorrected ages are typically in the range of 0.2 % to 0.6 % (Table S1). Drip water shows a high U concentration of 0.825 ng g⁻¹ and elevated initial Th concentrations, with an activity ratio of K = (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) = 11.1 ± 0.1. We have used the software to test how the chronology changes to assess the influence of a varying initial Th activity ratio. For this, we used three different correction models, including the measured initial Th ratio of the drip water (K = 11.1 ± 0.1), the detrital correction value of K = 0.75 ± 0.38 derived from the bulk Earth crust chemical composition, and a value of K = 23.7 ± 7.5 as previously determined using isochrons on speleothem PR-LA-1 from the same cave (Warken et al., 2020). Figure S3 in the Supplement shows the ages corrected for initial ²³⁰Th using the three different models. Only the initial ²³⁰Th value measured in the drip water yields a stratigraphic order of the corrected ages supporting the use of this value. Residual variability around the mean chronology increases and age inversions appear in the record when using a different value of K. Figure S4 shows different growth models for stalagmite PR-LA-B1. Growth rates vary between ca. 10 and 150 µm a⁻¹, with the highest values during the warm Bølling–Allerød period at ca. 13.97 ± 0.051 and 13.114 ± 0.073 ka BP as well as the late Holocene growth phase after 0.277 ± 0.008 ka BP. The lowest growth rates occur during the final stage of Heinrich Stadial (HS) 1 (16.23 ± 0.082 to 13.97 ± 0.051 ka BP), HS3 (31.02 ± 0.10 to 29.38 ± 0.12 ka BP), and HS4 (40.81 ± 0.16 to 39.12 ± 0.12 ka BP).

5 Discussion

5.1 Outlier correction

Outlier correction is carried out automatically by the software by adapting the dispersion measure of the raw data, and in the following we argue that generally means should be replaced by medians. Shao et al. (2019) addressed this problem by implementing manual outlier removal by comparison to box plots based on interquartile ranges. We opted for the automatic version as this is more time-efficient for large datasets. The different dispersion measure options described in Sect. 3.2 are relevant because measurements are not always ideal cases with normally distributed data and thus outliers. During measurements, short-term system instabilities occur for a variety of reasons, such as varying gas flow in the inlet system, plasma instabilities, and varying size of sample aerosols causing outliers in the signal intensities. Even though only the ratios between the different isotopes are of interest, strong changes in signal intensity may lead to varying isotope ratios as a result of changing variance. Such a difference may be amplified by the use of different detectors or with respect to different magnetic field settings, which are not necessarily responding at exactly the same amplitude. Moreover, signal decreases (detuning events and temporal clocking) cause the statistical variance to increase locally.

Figure 4 Measurement example. (a) ²³⁸U signal intensities in volts over measurement cycles for a carbonate sample during routine lab measurements. (b) Corresponding uncorrected ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ ratio. The mean and median as well as the three different dispersion measures are plotted.

Figure 4 shows an example: the upper panel displays periodic dips in the ²³⁸U signal intensity during a measurement. In the lower panel of Fig. 4, the uncorrected (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$) activity ratio for the same measurement is plotted. For both curves, the different measures to calculate dispersion are shown. It is clearly visible that the median agrees much better with the majority of signal intensity values than the mean, which is strongly influenced by the periodic dips due to the asymmetry in the statistical distribution. Such an obvious difference is not visible in the isotope ratio but within the resulting uncertainty. Consequently, we propose generally using the median instead of the mean by default. This is more accurate in the case of asymmetric small-scale oscillations inside the non-outlier interval and has no disadvantages.

Applying the standard deviation as a dispersion measure in Fig. 4 does not cover most of these outliers due to their large number and relatively small deviation. Thus, applying another dispersion measure for outlier removal is necessary here and in addition more robust and easier to accomplish than manual deletion of all of the outliers. It is important, however, to stress that the outlier correction using the selected dispersion option is run on the calculated ratios after correction, not on the signal intensities themselves. This implies that when all isotopes are affected in the same way, they pass the outlier test. This is, however, unlikely, at least for ratios of isotopes measured in different magnetic field settings. The dispersion measure of the outlier-corrected ratio array is the same in every case, as described above.

5.2 Detrital thorium correction

Thorium correction is often crucial for studying carbonates where the correction is significant but the initial ²³⁰Th value is unknown, potentially variable, or when studying “dirty” carbonates such as tufa and travertine (Mallick and Frank, 2002; Hellstrom, 2006; Wenz et al., 2016). Several studies have shown that this correction is particularly important for speleothem records from the Caribbean and Central American region, where values were found including 2 ± 1 (Schorndorf et al., 2023) or 14 ± 4 (Moseley et al., 2015). In Larga Cave, initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) ratios are presumably even higher where Warken et al. (2020) obtained a value of 23.7 ± 7.5 using isochrons on a stalagmite from the cave. Besides the terrestrial regime, this aspect is also relevant for marine archive such as corals, where studies propose a large range of seawater (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratios. While Cheng et al. (2000a) set the range to 80 ± 80 for deep-sea solitary corals and Frank et al. (2004) calculated 10 ± 4 from seawater in the eastern North Atlantic deep sea, values between 0.4–3.1 were determined for tropical corals (Shen et al., 2008). The range of both absolute values and uncertainties for these widely studied archives is hence enormous, and the choice of the appropriate correction model becomes particularly important when (i) samples are very young and have generated only small amounts of ²³⁰Th from U decay or (ii) when ultrahigh precision is at play since any possible correction of the data contributes to the final age uncertainty. In our case study, we have run the correction of the ages of stalagmite B1 using three different correction models (Table S1). The resulting differences are visualized in Fig. S3 and demonstrate the significant impact not only on absolute corrected ages, but also their uncertainties. For the young age at 7 mm dft (0.0466 ± 0.0045 ka BP), the difference in the absolute corrected age when using a correction factor other than the drip water value of K = 11.1 ± 0.1 is ca. ±50 years, which corresponds to a relative difference of the order of 100 % (compare Table S1). Another example is the sample at 554 mm dft (37.81 ± 0.14 ka BP for K=11.1 ± 0.1), for which the other correction models also lead to substantially different ages of 41.37 ± 0.19 ka BP (K=0.75 ± 0.38) and 33.3 ± 2.7 ka BP (K=23.7 ± 7.5); hence, the differences are still in the range of ca. 10 %. Notably, the low relative error of the initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio of the drip water results in equally low uncertainties of the corrected age in the range of 0.4 %. In contrast, the relative uncertainty of the age corrected with K=23.7 ± 7.5 increases to 8 %. Our GUI permits an easy adjustment of the initial (${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$) activity ratio for Th correction, which allows a direct assessment of the resulting corrected ages and uncertainties and thus provides a convenient basis for further comparisons of the data. The use of standardized software instead of handmade tuning reduces the susceptibility to potential errors, e.g., from copying and pasting, and ensures reproducibility in the case that a re-evaluation of the data is required at a later stage.

5.3 In-cave comparison of speleothem growth rates

The high number and high precision of ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ ages of speleothem B1 allow investigation of growth rate (GR) changes. Comparison with northern hemispheric climatic changes suggests that speleothem B1 growth is sensitive to prominent millennial-scale temperature variability, with a higher growth rate during warmer phases and vice versa. In particular, during the cooler and drier Heinrich Stadials (Warken et al., 2022a), growth rates are reduced. In addition, the results allow a comparison of the two coeval stalagmites from Larga Cave as shown in Fig. S4. Overall, GRs of PR-LA-B1 are about 5 times lower than observed for PR-LA-1, where average annual growth rates are up to several millimeters per year. The difference in mean GR is also reflected in the shape of both speleothems, with PR-LA-1 exhibiting a large and variable diameter between ca. 15 and 35 cm (Warken et al., 2020), while B1 is thinner with a diameter of 10–15 cm (Fig. S1b). Differences in speleothem growth rates and the shape of a stalagmite may result from temperature, carbonate saturation, drip rate, and carbon dioxide contrast between cave air and saturation concentration of drip water (Merz et al., 2022; Skiba and Fohlmeister, 2023; Kaufmann, 2003; Dreybrodt, 1999). Ca concentrations in Larga Cave show no significant differences between drip sites (Vieten et al., 2018a, 2018b, Warken et al., 2022b). Therefore, the amplified GR and larger diameter of PR-LA-1 could be the result of the considerably lower pCO₂ values in the main passage than compared to the back part of the cave (Vieten et al., 2016; Vieten and Hernandez, 2021), which facilitates enhanced oversaturation of the drip water with respect to calcite and, hence, stronger degassing of CO₂ and speleothem growth (Merz et al., 2022). Hence, the two stalagmites reveal growth differences potentially related to ventilation conditions.

6 Conclusions

We provide an algorithm combined with a user-friendly GUI application for the treatment of ²³⁰$\mathrm{Th}/\mathrm{U}$ MC-ICP-MS data obtained by Thermo Fisher Neptune instruments and subsequent age calculation and correction. The two programs published so far, which are explicitly aimed at ${}^{\mathrm{230}}\mathrm{Th}/\mathrm{U}$ dating data reduction and age calculation, are both written for Thermo Fisher Neptune instruments as well. Pourmand et al. (2014) described a Mathematica routine, distributed as a Computable Document Format (.cdf) file, while Shao et al. (2019) published a MATLAB algorithm with GUI. We have chosen to use Python for our algorithm and GUI to keep it open-source. The advanced user might want to change settings, which makes an open-source language and libraries a major advantage. However, the stand-alone executable .exe format of the GUI also allows user-friendly handling for non-programming experts. Our program supports multiple types of detector configurations: the FC–FC-based approach as well as FC–SEM combining protocols. It is, however, adapted for combined Th and U measurements in three magnetic field lines (compare Kerber et al., 2023), but other methods (such as separate solutions for Th and U) can be implemented with small changes in the code. Furthermore, we offer the first-order Taylor derivation as a time-saving option for uncertainty calculation of final ages. Our application is especially designed to take reproducible and clear data management into account by using a collection of methods: this enables automatic creation of folders containing the “results” files and information on the sample metadata, as well as .xlsx output files that automatically contain all constants used for calculation, in addition to the settings for outlier correction. Manually changing input constants, e.g., correction of initial and detrital Th, does not require going to the code directly. So, the whole analysis scheme does not require any copying and pasting from one excel table to the other, and the constants used for calculation are easy to update.

Lastly, we demonstrated our protocols and data analysis scheme by accurately measuring and evaluating 30 speleothem ages from Larga Cave, Puerto Rico. Analyses of the growth rates and comparison with a coevally growing stalagmite from the same cave highlight the importance of in-cave processes for speleothem deposition rates.