Radiometric dating of Middle Pleistocene carbonates: assessing consistency and performance of the U–Th and U–Pb dating methods

The U–Th and U–Pb dating methods are widely employed for radiometric dating of Pleistocene carbonates, such as speleothems and corals. The U–Th dating method has been progressively refined over recent decades, largely through advances in mass spectrometry, and is now capable of providing accurate and precise ages for carbonates as old as ∼ 650 000 years under optimal circumstances. Similarly, the U–Pb method, traditionally restricted to dating pre-Quaternary materials, has been adapted in recent years for dating young carbonates. As a result, there is now substantial overlap in the applicable age range of these two dating methods but, thus far, only limited assessment of their consistency and relative performance when dating samples over this shared age interval.

In this study, we conduct a systematic comparison of the U–Th and U–Pb dating methods, focusing on a significant part of their overlapping age range (approximately 630–430 ka). We achieve this by dating speleothem (secondary cave mineral deposit) samples from Corchia Cave, central Italy, using both methods and evaluate their consistency and performance in terms of age precision and other factors. We adopt analytical approaches that employ state-of-the-art multi-collector inductively coupled plasma mass spectrometry (MC-ICPMS) measurement protocols, including a U–Th measurement protocol that utilises a Faraday cup equipped with a 10¹³ Ω resistor to collect the low-abundance ²³⁴U⁺ and ²³⁰Th⁺ ion beams. This approach is particularly well-suited to dating samples approaching the limits of the U–Th method but also enables accurate determination of ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ ratios. Thus, as a secondary component of this study, we compare our ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ measurements with previously published speleothem values.

Our results demonstrate excellent agreement between the U–Th and U–Pb dating methods and suggest that both are capable of providing accurate and precise ages over this interval. We find that U–Pb age uncertainties are generally less predictable than U–Th age uncertainties but, on average, do not increase significantly over the interval considered. U–Th age uncertainties, on the other hand, tend to increase in a more predictable and approximately exponential manner. Additionally, U–Pb age uncertainties are highly dependent on the availability of sub-samples with a substantial spread in parent/daughter ratios and/or highly “radiogenic” (i.e. very low inherited Pb) material. In our dataset, U–Pb isochron age precision surpasses that of U–Th precision at ∼ 520 ka, although the exact crossover point is expected to vary for different sample types and depositional settings. Overall, these findings support the prospect of obtaining accurate and internally consistent U-series chronologies spanning the Middle Pleistocene. They also suggest that, for some carbonate samples, the U–Pb dating method may provide superior age precision to the U–Th method prior to the latter reaching its upper age limit.

Finally, our results show that most speleothems exhibit ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ ratios consistent with global carbonate values and that these ratios typically deviate from the conventional value of 137.88, widely adopted in geochronology, in agreement with previous studies.

1 Introduction

U-series dating of carbonates has provided some of the most accurate and precise age constraints on Earth processes unfolding over the Pleistocene. The U–Th disequilibrium dating technique (hereafter, U–Th method) is by far the most widely used approach for radiometric dating of Pleistocene carbonates. Its application to corals, beginning in the 1960s, and later speleothems (secondary cave mineral deposits) has led to significant advances in our understanding of Pleistocene climate. For example, U–Th dating of coral reef terraces has provided some of the most compelling evidence linking orbital insolation variations to glacial–interglacial cycles (e.g. Edwards et al., 2003), while U–Th dating of speleothems has yielded some of the most accurate and precise constraints on the timing of major climate events over the Pleistocene, including glacial terminations (Cheng et al., 2016) and millennial-scale climate oscillations associated with Dansgaard–Oeschger events (Corrick et al., 2020). U–Th dating of carbonates has also been used to address important questions in fields such as archaeology and human evolution (Hellstrom and Pickering, 2015; Pike et al., 2017) and provides data that underpin the calibration of the radiocarbon age scale (Reimer et al., 2020).

U–Th analytical methods have been continuously refined over the past few decades. The development of thermal ionisation mass spectrometry (TIMS) protocols in the late 1980s (Chen et al., 1986; Edwards et al., 1987), and later MC-ICPMS protocols, led to major improvements in U–Th age precision and sample throughput (Goldstein and Stirling, 2003). Following these major breakthroughs, U–Th analytical protocols have undergone more incremental refinement over the past two decades or so. An example of one such refinement is the advent of all-Faraday-cup analytical protocols that collect the low-abundance ²³⁴U⁺ and ²³⁰Th⁺ ion beams in a Faraday cup instead of using an ion counter (Andersen et al., 2004; Potter et al., 2005; Cheng et al., 2013). This approach allows much larger ion beams to be measured and circumvents issues associated with ion counter dead time correction and gain calibration, leading to significant improvement in age precision and an extension of the useful range of the dating method to ∼ 650 ka in ideal circumstances (Cheng et al., 2016).

The U–Pb dating method, traditionally used mostly for radiometric dating of igneous materials, such as zircon, has recently been adapted for dating Pleistocene carbonates, and it offers a valuable means of extending speleothem and coral chronologies beyond the current ∼ 650 ka limit of the U–Th technique (Richards et al., 1998; Woodhead et al., 2006; Denniston et al., 2008; Klaus et al., 2017). To date, the U–Pb method has been less widely applied to Pleistocene carbonates than U–Th dating. This is partly because it is more labour-intensive (especially when adopting high-precision isotope dilution TIMS or multi-collector inductively coupled plasma mass spectrometry – MC-ICPMS – approaches) but also because its application to Pleistocene carbonates is limited to samples with low inherited Pb content (Woodhead et al., 2012). Nevertheless, U–Pb dating is of similar utility, having, for example, provided radiometric age constraints on the timing of Pleistocene glacial terminations (Bajo et al., 2020), the age of key hominin fossils (Walker et al., 2006; Pickering et al., 2011), and chronologies for sub-arctic permafrost thawing over the past ∼ 1.5 Ma (Vaks et al., 2020). Although the U–Pb method is nominally better suited to dating older material, for which sufficient time has passed for significant accumulation of radiogenic Pb, it is also well-suited to dating Middle Pleistocene carbonates that typically have ${}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$ activity ratios (hereafter, ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$, where square brackets denote an activity ratio) that are analytically resolvable from radioactive equilibrium, thus facilitating a precise and accurate disequilibrium correction (Woodhead et al., 2006). In ideal circumstances the carbonate U–Pb dating method can also provide precise ages for materials as young as ∼ 270 ka (Cliff et al., 2010).

With recent advances in the carbonate U–Th and U–Pb dating methods, there is now considerable overlap in the applicable range of the two chronometers, especially over the interval from ∼ 650 to ∼ 400 ka. This overlap provides a valuable opportunity to assess consistency between the two methods and thus validate, to a certain extent, the assumptions underlying their use. However, since both dating approaches are based on the ²³⁸U decay series, they are not entirely independent, and for this reason, a comparison between them may not be considered to be as definitive as one involving completely independent chronometers (e.g. Min et al., 2000), all else being equal. Nevertheless, the carbonate U–Th and U–Pb dating methods do rely on different assumptions and require analysis of distinct isotopic ratios. For instance, the U–Th method relies on the precise analysis of ${[}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}]$ and is highly sensitive to any inaccuracies in this measured activity ratio. In contrast, the U–Pb method requires the analysis of Pb-based isotope ratios and often incorporates information from the ²³⁵U decay series in addition to the ²³⁸U series (although, for young samples, most of the age resolving power lies in the ²³⁸U–²⁰⁶Pb chronometer). The U–Pb dating method also relies on assumptions about the initial state of intermediate products beyond ²³⁰Th in the ²³⁸U decay series, such as ²²⁶Ra and ²³¹Pa, and closed-system behaviour with respect to these nuclides (Richards et al., 1998), which do not apply to the U–Th method.

Other independent methods of dating carbonates within this age range are available, such as electron spin resonance (ESR) dating (Radtke et al., 1988) and (U,Th)–He dating (Bender et al., 1973; Makhubela and Kramers, 2022). However, these are less well-established and tend to provide lower age precision than the U–Th and U–Pb methods, limiting their utility in assessing the accuracy of these more precise U-series-based dating approaches.

The availability of two different U-series approaches to dating carbonates over the ∼ 650–400 ka interval also raises the question of which is most suitable to a given situation. While U–Th dating is often assumed to be the most accurate and precise approach throughout most, if not all, of its applicable range, the U–Pb method may perform equally well or even better for samples approaching the upper age limit of U–Th dating. Consequently, selecting an optimal dating approach for carbonates in this age range requires careful consideration. In practical terms, U–Th dating may be considered advantageous because it typically involves analysing single samples extracted from discrete positions within a given speleothem or coral sample, making it less labour-intensive on a per-age basis. U–Pb dating, on the other hand, typically involves isochron approaches, whereby multiple sub-samples from within an individual growth domain are analysed to compute a single age. Although more labour-intensive, isochron methods also provide a first-order check on the assumption of closed-system behaviour and more accurately account for any initial (i.e. non-radiogenic) quantity of the daughter isotopes.

At present, the relative precision of each dating method is not yet well characterised for samples approaching their respective age limits (i.e. upper limit for U–Th and the lower limit for U–Pb). For “well-behaved” samples – those maintaining closed-system conditions with reasonably high U and low initial Th/Pb content – the precision of U–Pb dating is likely to surpass that of U–Th somewhere before the latter reaches its nominal upper age limit. However, the exact crossover point is difficult to predict a priori because U–Pb isochron age precision is not a simple function of isotope ratio measurement precision but also depends on various sample characteristics (e.g. Woodhead et al., 2012; Engel and Pickering, 2022). The direct comparison of samples dated using both methods therefore offers the most appropriate means of evaluating their relative performance in terms of age precision and other practical considerations.

Previous studies have compared carbonate samples dated using both the U–Th and U–Pb methods, although these have tended to include a relatively small number of samples across a restricted age range (e.g. Richards et al., 1998; Pickering et al., 2010). The most comprehensive comparison to date is that of Cliff et al. (2010), who obtained multiple U–Th and U–Pb ages from a flowstone that grew between ∼ 350–265 ka. They found good overall agreement for the youngest growth segment, although age agreement for the middle and older growth segments was more equivocal. While this study was more comprehensive than previous ones, it also had some limitations. For example, only one full isochron age was included in the comparison and there was some uncertainty in co-registering the U–Th and U–Pb ages to a common depth scale, which precluded a more statistically rigorous assessment.

Here we present a systematic comparison of U–Th and U–Pb dating methods focusing on the older part of their overlapping age range (∼ 630–430 ka) where routine application of the carbonate U–Pb dating method is a more realistic prospect. We achieve this by dating carbonate speleothem samples extracted from identical growth layers using both the U–Th and U–Pb methods and then assess consistency between them and evaluate their relative merits in terms of age precision and other practical considerations. In contrast to previous comparison studies, we compute full Tera–Wasserburg (Tera and Wasserburg, 1972) isochron ages for all samples and adopt state-of-the-art MC-ICPMS-based analytical protocols. This includes the use of an all-Faraday-cup U–Th measurement protocol that employs a ²³³U–²³⁶U double spike for mass bias correction and collects the low-abundance ²³⁴U and ²³⁰Th isotopes in a Faraday cup fitted with a 10¹³ Ω resistor, which provides superior signal-to-noise ratios for small ion beams relative to more conventional 10¹¹ Ω resistors (Pythoud, 2022). The use of this analytical protocol allows age consistency to be assessed more rigorously than in previous studies and facilitates an assessment of age precision using the most advanced U–Th measurement protocol currently available for samples in this age range.

An added benefit of adopting this high-precision U–Th analytical protocol is that it yields ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ data of sufficient accuracy and precision for assessing natural variability in this ratio. The ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ ratio is often used in U–Th and U–Pb data reduction, as well as U–Pb and Pb $/$ Pb age calculation, and has traditionally been treated as a constant in geochronology (Hiess et al., 2012; Tissot and Dauphas, 2015). However, the development of high-precision ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ measurement protocols over the past two decades has revealed variations in natural ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ ratios, typically expressed as δ²³⁸U values¹, of up to a few per mil (‰) (Weyer et al., 2008; Tissot and Dauphas, 2015). Despite these findings, only a small number of studies have examined speleothem δ²³⁸U values. As a secondary component of this study we therefore combine our high-precision δ²³⁸U measurements with pre-existing speleothem data in order to assess some of the controls on speleothem δ²³⁸U variability more broadly.

2 Study site and samples

The speleothem samples analysed in this study were collected from Corchia Cave, located in the Alpi Apuane massif of central Italy (43°59^′ N, 10°13^′ E). This cave system, which is developed in Mesozoic marbles, dolomitic marbles, and dolomites, is one of the longest (∼ 60 km) and deepest (∼ 1250 m) in Europe and has been the subject of extensive past cave exploration (Piccini et al., 2008) and scientific study (Bajo et al., 2017; Drysdale et al., 2019; Isola et al., 2019). The cave site receives the majority of its high annual precipitation (∼ 2500 mm yr⁻¹) from eastward-moving air masses originating over the North Atlantic Ocean, resulting in a coupling between large-scale surface ocean conditions in the North Atlantic and geochemical signals preserved in speleothems at the site (e.g. Drysdale et al., 2020). Owing to this linkage, Corchia Cave has yielded an abundance of palaeoclimate data of regional and global significance, including precise constraints on the timing of glacial terminations (Drysdale et al., 2020; Bajo et al., 2020) and interglacial climate variability (Regattieri et al., 2014; Tzedakis et al., 2018).

Palaeoclimate reconstruction studies undertaken at Corchia Cave have largely focussed on stalagmites collected from the large and well-decorated Galleria delle Stalattiti chamber. This chamber sits approximately 400 m below the surface at an elevation of 835 m above sea level and about 800 m from the nearest natural cave opening (Piccini et al., 2008). The microclimate of Galleria delle Stalattiti is relatively stable, typical of deep-cave environments, with a humidity consistently above 98 % and an average present-day temperature of ∼ 7.9 °C (Drysdale et al., 2019). Stalagmites from this chamber are typically well-suited to U-series geochronology with consistently high U content (typically 3–20 ppm) and very low detrital Th, as well as relatively low inherited Pb (Bajo et al., 2012). Furthermore, most stalagmites are composed of clean (i.e. free of visible detritus) compact columnar calcite, which is typically associated with the preservation of primary geochemical signals and closed-system behaviour with respect to U-series nuclides (Scholz et al., 2014). These characteristics make stalagmites from the Galleria delle Stalattiti very well-suited to U-series geochronology and, when combined with the site's suitability for palaeoclimate reconstruction, have led to a large number of U–Th and U–Pb determinations being undertaken on speleothems from the site. These data suggest that initial ${}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$ values of speleothems from this cave chamber are always less than 1 and have undergone a long-term evolution over the past several hundred thousand years towards more ²³⁴U-depleted values (Woodhead et al., 2006), with activity ratios as low as ∼ 0.65 obtained for Late Holocene materials (Zanchetta et al., 2007).

This study focusses on three stalagmites (sample IDs: CCB, CC2/CC15, and CC17-1) that were found broken near their original in situ position during field campaigns between 2001 and 2017. These stalagmites are large, with lengths of 67–130 cm and average diameters ranging from approximately 9–22 cm. Following collection, the stalagmites were halved along their vertical growth axes and polished to aid in the identification of visible growth laminations. In section, the stalagmites are white to grey in colour and predominately composed of compact translucent-to-opaque primary calcite. Exceptions are the very top of CCB, which contains some slightly porous calcite as the diameter tapers toward the tip, and the middle of CC17-1, which contains some localised, slightly porous sections associated with sharp changes in growth axis direction. These features are atypical of Galleria delle Stalattiti speleothems and were avoided during sampling.

3 Analytical methods

3.1 Sampling

Samples for both U–Th and U–Pb dating were extracted using a dental drill, targeting clean and compact calcite by carefully tracing along visible growth contours approximately 5 mm apart at the central axis but tapering along the flanks. Several solid sub-samples of ∼ 4 mm (width) × ∼ 4 mm (length) × ∼ 3 mm (depth), or 20–120 mg, were then extracted from each growth domain. For each sampling position, one sub-sample, taken from the centre of the sampling area, was used for U–Th analysis and the remaining 4–11 sub-samples were used for U–Pb analysis. The only exception to this is sample CC17-1-3, where the U–Th sample was taken from a growth layer marginally below that of the corresponding U–Pb samples. This sampling offset was negligible relative to the age and age uncertainty of these growth layers, and so for the purpose of this study the samples are assumed to have been extracted from the same growth layer.

3.2 U–Th geochronology

U and Th isotopic analysis was performed at the University of Minnesota using a modified version of the all-Faraday-cup mass spectrometry protocol of Cheng et al. (2013) but with the ²³⁴U⁺ and ²³⁰Th⁺ ion beams collected in a Faraday cup fitted with a 10¹³ Ω resistor (Pythoud, 2022). Chemistry procedures were similar to those described in Edwards et al. (1987). Briefly, solid samples of 20–110 mg (equivalent to ∼ 500 ng U) were dissolved in ∼ 1 % HNO₃ before the addition of a mixed ²³³U–²³⁶U–²²⁹Th tracer calibrated against IRMM-074/10 (Cheng et al., 2013). Approximately six drops of concentrated HClO₄ were added to each sample beaker before being capped and refluxed on a hotplate for several hours to ensure complete sample-spike equilibration. Following complete dry down, the separation and purification of U and Th proceeded in two stages. Firstly, U and Th were separated from the matrix using Fe co-precipitation in 2 N HCl by drop-wise addition of NH₄OH followed by centrifuging. Secondly, U and Th were separated and further purified via elution chromatography using BioRad AG1-X8 anion exchange resin, with the Th fraction collected in 6 N HCl and the U fraction collected in ultra-pure water. The Th fraction was then passed through a second identical column to ensure complete separation from U. Following this, the purified U and Th sample fractions were dried down and sequentially treated with HClO₄ and HNO₃ to oxidise any remaining organics before being taken up in 1 % HNO₃ for analysis.

The U and Th isotope ratios were measured separately in static mode on a Thermo-Scientific Neptune Plus MC-ICPMS (see Table S1 in the Supplement for details on the detector configuration). Samples were introduced into the mass spectrometer via a Cetac Aridus II desolvating nebulising system fitted with a C50 PFA concentric nebuliser, using an uptake rate of 30–50 µL min⁻¹. Zeros were measured for 13 min for both U and Th immediately prior to sample introduction. U measurements were then made over 5–10 min with a ²³⁴U⁺ signal size between 8–18 mV. Th measurements immediately followed U for each sample and were made over 1–5 min with a ²³⁰Th⁺ signal size between 8–20 mV.

For Faraday cups fitted with 10¹⁰ and 10¹¹ Ω resistors, amplifier gain calibration was performed via the mass spectrometer software using a standard current at the start of each run. However, the limited dynamic range of the 10¹³ Ω resistor precluded the use of this approach for the cup fitted with this resistor. Therefore, the amplifier gain on this detector was characterised by periodic (∼ fortnightly) measurement of Nd standard SRM 3135a following Pythoud (2022). Tailing of the ²³⁸U⁺ and ²³⁵U⁺ ion beams was characterised by analysis of un-spiked CRM-112A prior to each run, with the background measured at masses 233, 233.5, 234.5, 236, and 237 amu on the axial SEM with the retarding potential quadrupole (RPQ) off. A logarithmic regression fit through these data was then used to estimate tailing onto the mass 233–236 cups as a function of the signal intensity at mass 237 during measurement (Pythoud, 2022). For each U fraction, tailing was corrected by measuring the background at mass 237 on the axial SEM during sample analysis and applying the predicted $M_{\mathrm{tail}}/\mathrm{237}$ ratios (where M refers to masses 233–236). Tailing of the ²³²Th⁺ beam onto the mass 230 cup was negligible given the very high ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$ ratios of the samples analysed in this study and was not corrected for.

Table 1 U–Th dating results.

Uncertainties are given at the 2σ level. Decay constant values are ${\mathit{\lambda }}_{\mathrm{238}}=\mathrm{1.55125}\times {\mathrm{10}}^{\mathrm{-}\mathrm{10}}$ (Jaffey et al., 1971), ${\mathit{\lambda }}_{\mathrm{234}}=\mathrm{2.82206}\times {\mathrm{10}}^{\mathrm{-}\mathrm{6}}$ (Cheng et al., 2013), and ${\mathit{\lambda }}_{\mathrm{230}}=\mathrm{9.1705}\times {\mathrm{10}}^{\mathrm{-}\mathrm{6}}$ (Cheng et al., 2013). ^a ${\mathit{\delta }}^{\mathrm{234}}\mathrm{U}={[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]-\mathrm{1}$. ^b Uncertainties calculated by first-order algebraic uncertainty propagation (Appendix B). ^c Calculated by Monte Carlo simulation. More complete U–Th results are available in Tables S2 and S3 of the Supplement.

Data reduction was performed using an offline spreadsheet that implements a tau correction following Pythoud (2022) to account for the slow response time of 10¹³ Ω resistors before applying the ²³⁸U tail correction as described above and mass bias correction based on the exponential law (Russell et al., 1978) using the ${}^{\mathrm{233}}\mathrm{U}{/}^{\mathrm{236}}\mathrm{U}$ ratio measured for each sample. The mass bias factor calculated for the U fraction was subsequently used to correct the Th fraction run immediately after, based on the assumption of near-equivalent mass bias behaviour between U and Th isotopes (Potter et al., 2005). This is justified by the consistent mass bias behaviour observed on this instrument, especially over multiple successive analyses (Pythoud, 2022). Finally, corrected U ratios were normalised to CRM-112A, which was periodically analysed throughout each run, using the CRM-112A value of 52.852 ± 0.015 obtained by Cheng et al. (2013). U–Th ages and uncertainties (see Table 1 and Tables S2 and S3 in the Supplement) were then computed using a Python script that calculates ages and age uncertainties using both first-order algebraic uncertainty propagation (Appendix A1 and B) and Monte Carlo approaches (Hellstrom, 1998; Ludwig, 2003) (see the Supplement for further details).

Owing to the very high ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$ ratios of the samples analysed in this study, we opted not to correct U–Th or U–Pb ages for initial ²³⁰Th in order to simplify age and age uncertainty calculations (i.e. Appendix B). We note, however, that if such a correction were applied it would result in a negligible correction of at most a few years for both U–Th and U–Pb ages.

3.3 U–Pb geochronology

U–Pb dating was performed at the University of Melbourne. Analytical methods built on those published previously (Woodhead et al., 2012; Sniderman et al., 2016) but using the streamlined chemistry procedure described in Engel et al. (2020). Solid calcite samples were repeatedly cleaned by brief (∼ 1 min) immersion in dilute (∼ 0.01 M) multiply distilled HCl, followed by washing in ultra-pure water, before drying in a HEPA-filtered laminar flow hood. The cleaned samples were then individually weighed into pre-cleaned Teflon beakers and treated with sufficient 6 N HCl to ensure complete dissolution. A mixed ${}^{\mathrm{233}}\mathrm{U}{/}^{\mathrm{205}}\mathrm{Pb}$ isotopic tracer, calibrated against EarthTime reference solutions (Condon et al., 2015), was then weighed into the beakers which were each sealed and refluxed on the hotplate for several hours to ensure complete sample-spike equilibration. Samples were then dried down and U and Pb separated from the calcite matrix using the single-column mixed-resin ion exchange technique described in Engel et al. (2020).

Isotope ratios were measured on a Nu Plasma MC-ICPMS using a DSN-100 desolvation unit and MicroMist glass nebuliser, operating with an uptake rate of ∼ 80 µL min⁻¹. Instrumental mass bias effects were monitored and corrected using NIST SRM 981 reference material in the case of Pb and the sample's internal ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ ratio in the case of U.

U–Pb isochrons were fitted using DQPB (Pollard et al., 2023), initially employing the spine robust isochron algorithm (Powell et al., 2020). This algorithm accounts for analytical uncertainties and can accommodate datasets with weighted residuals that depart from a strict Gaussian distribution but converges to the classical least-squares solution of York et al. (2004) for “well-behaved” datasets (i.e. those where data scatter is accounted for by assigned measurement uncertainties). The spine width, s, was used to assess suitability of this algorithm for each dataset in relation to the upper 95 % confidence limit on s (here denoted s_lim) obtained by simulation of Gaussian-distributed datasets (Powell et al., 2020).

In cases where s clearly exceeded s_lim, the Robust Model 2 (RM2) algorithm presented in Pollard et al. (2023) was implemented instead. This algorithm ignores analytical uncertainties, inferring a robust dispersion scale from data scatter itself rather than assigned measurement uncertainties. Therefore, it is arguably better suited to datasets containing significant overdispersion (i.e. data scatter in excess of that attributable to assigned measurement uncertainties). However, for isochrons where s only marginally exceeded s_lim, the spine algorithm was considered likely to provide a more accurate best-fit line than the RM2 algorithm, although in this case the spine uncertainty calculations become potentially unreliable. Thus, where s exceeded s_lim by less than an arbitrary but reasonable value of 50 %, the spine regression line was retained, but uncertainties were expanded in an attempt to accommodate data scatter above what is reasonably accounted for by measurement uncertainties. This was achieved by multiplying the covariance matrix of the slope and y intercept, V_θ, by $(s/s_{\lim }{)}^{\mathrm{2}}$. Implementing the RM2 algorithm instead for these datasets does not significantly alter any of the conclusions of this study.

U–Pb ages were calculated as the intercept point between the isochron line and disequilibrium concordia curve on a Tera–Wasserburg diagram, incorporating the measured ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ value from the U–Th analysis into the iterative age-solving procedure (Appendix A2) and assuming initial ${[}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}]$, [${}^{\mathrm{231}}\mathrm{Pa}{/}^{\mathrm{235}}\mathrm{U}$], and [${}^{\mathrm{226}}\mathrm{Ra}{/}^{\mathrm{238}}\mathrm{U}$] values equal to 0 (see Fig. 1). Uncertainties were calculated using a Python script that implements the Monte Carlo simulation approach described in Pollard et al. (2023) using 10⁶ iterations (see the Supplement for further details) and also via first-order algebraic uncertainty propagation (Appendix B).

Figure 1 U–Pb Tera–Wasserburg isochron diagrams. Measured data points are plotted as 95 % confidence ellipses (grey) along with the fitted isochron lines (orange) and their 95 % confidence bands (orange shading). The dark blue lines show disequilibrium concordia curves, constructed using the measured ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ value and assuming ${[}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}]={[}^{\mathrm{226}}\mathrm{Ra}{/}^{\mathrm{238}}\mathrm{U}]={[}^{\mathrm{231}}\mathrm{Pa}{/}^{\mathrm{235}}\mathrm{U}]=\mathrm{0}$. This Concordia curve is truncated at the youngest age point associated with a physically impossible negative initial activity ratio solution following the approach of McLean et al. (2016). Provided in the text boxes is the spine width, s; the value for a spine regression fit; and the simulated upper 95 % confidence limit on s, s_lim (see text for further discussion). The MSWD (mean square weighted deviation) for a classical weighted least-squares regression (e.g. York et al., 2004) is also provided for comparison with other datasets. The regression model implemented is indicated in brackets: spine indicates the use of the algorithm described by Powell et al. (2020), and RM2 indicates the use of the Robust Model 2 algorithm (Pollard et al., 2023). The data used to produce these plots are available in Tables S4 and S5 of the Supplement.

For the purpose of determining analytical accuracy, each U–Pb session included analysis of several “synthetic zircon” standards: ET-100Ma, ET-500Ma, and ET-2Ga. These materials, produced by the EarthTime initiative, have nominal U and Pb compositions approximating zircon minerals with ages of 100 Ma, 500 Ma, and 2 Ga but without any matrix (Condon et al., 2008). While they were originally conceived for quality control during TIMS single-zircon analysis and are still routinely used for that purpose (Szymanowski and Schoene, 2020), they are also well-suited to isotope dilution (ID) carbonate U–Pb studies since the quantities of U and Pb being analysed can be controlled to approximate unknowns and the highly radiogenic nature of the Pb involved serves to rapidly identify any anomalous analytical issues, such as unexpected high blanks. In our studies, we pass aliquots of these solutions (typically equating to ∼ 5 ng of radiogenic Pb) through the same column chemistry as our samples. The measured isotopic ratios of these standards analysed during the course of this study were consistent with benchmark TIMS values (Fig. 2).

Figure 2 Radiogenic ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{238}}\mathrm{U}$ and ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ age determinations for EarthTime (ET) “synthetic zircon” standards analysed during the course of this study at the University of Melbourne (UoM) using a similar Pb load to the speleothem samples and identical column chemistry and mass spectrometry procedures. The weighted average age is plotted as the black line with grey shading indicating the 95 % confidence band. Also shown for comparison are the average inter-laboratory TIMS ages obtained by Schaltegger et al. (2021). ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ ages were calculated using ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}=\mathrm{137.818}$.

4 Results and discussion

4.1 Age consistency

We first evaluated age consistency on a pairwise basis. For each sample, we tested if the U–Th and U–Pb ages agreed within their analytical uncertainties using a hypothesis test that amounted to assessing if the age difference was statistically consistent with zero (e.g. Barlow, 1989). The test statistic was calculated as

$$\begin{array}{}\text{(1)}& z={\displaystyle \frac{\left|\mathrm{\Delta }\right|}{{\mathit{\sigma }}_{\mathrm{\Delta }}}}\end{array}$$

such that

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }=t_{\mathrm{Th}}-t_{\mathrm{Pb}},\end{array}$$

where t_Th and t_Pb represent the U–Th and U–Pb ages respectively and σ_Δ is the standard error in the age difference. The latter was calculated via first-order uncertainty propagation as

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathrm{\Delta }}=\sqrt{{\mathit{\sigma }}_{t_{\mathrm{Th}}}^{\mathrm{2}}+{\mathit{\sigma }}_{t_{\mathrm{Pb}}}^{\mathrm{2}}-\mathrm{2}{\mathit{\sigma }}_{t_{\mathrm{Th}}}{\mathit{\sigma }}_{t_{\mathrm{Pb}}}\mathit{\rho }},\end{array}$$

where ${\mathit{\sigma }}_{t_{\mathrm{Th}}}$ and ${\mathit{\sigma }}_{t_{\mathrm{Pb}}}$ are the standard errors in the U–Th and U–Pb ages and ρ is their correlation coefficient arising from their shared ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ measurement uncertainty (see Appendix B for further details). In most cases this correlation term was relatively minor, and it was negligible in cases where U–Pb regression fitting uncertainties were relatively large. We then obtained a formal p value under the null hypothesis that there is no significant difference between the ages by comparing z against the standard normal distribution.

Table 2 Summary of pairwise U–Th versus U–Pb age comparison.

^a Weighted average of U–Th and U–Pb ages. ^b Age difference calculated as t_Th−t_Pb. ^c Uncertainty (±2σ) calculated by first-order algebraic error propagation. ^d Calculated by Monte Carlo simulation. ^e The p value under the null hypothesis in which there is no significant age difference. ^f Result of the hypothesis test of age consistency (α=0.05 significance level). ^∗ The p value should be interpreted as a qualitative indicator only (see text).​​​​​​​

Strictly speaking, this approach is only appropriate if uncertainties in both U–Th and U–Pb ages conform to a Gaussian distribution. While this assumption is reasonable for the majority of ages, it is slightly inaccurate for the older U–Th ages (i.e. samples CCB-6-1, CC17-1-3, and CC15-1), where the age uncertainty distributions are positively skewed owing to the non-linearity of the U–Th age equation, which is more apparent for older ages and/or those with larger analytical uncertainties (e.g. Ludwig, 2003). Therefore, for these samples, we instead implemented a Monte Carlo procedure to simulate the age difference distribution (see the Supplement for details). Analogous to the formal hypothesis test, we concluded that the age difference was consistent with zero if the estimated 95 % confidence interval on the Monte Carlo simulated age difference overlapped zero. The results of this pairwise age comparison are given in Table 2 and Fig. 3a.

Figure 3 Comparison of U–Th and U–Pb age determinations. (a) Age difference for each sample calculated as t_Th−t_Pb. Note that the “CC” prefix of each sample ID has been omitted for brevity. (b) Collective comparison of U–Th and U–Pb ages with each age pair plotted as a 95 % confidence ellipse. The confidence ellipse for sample CCB-6-1, which was found to be inconsistent on a pairwise basis and therefore excluded from the MSWD calculation, is plotted with the dashed outline. Error bars and confidence ellipses are plotted at the 95 % confidence level.

Our results show excellent agreement between the U–Th and U–Pb ages for all samples, with the exception of sample CCB-6-1, where the U–Th age is significantly older than the U–Pb age by 41 + 34 $/$ −25 ka (95 % confidence). In this particular case, we consider the U–Pb age to be more reliable for two key reasons. Firstly, the U–Pb isochron is well formed, with good spread in U $/$ Pb ratios and data scatter that is consistent with analytical uncertainties (s<s_lim and p=0.14 for a classical least-squares fit). Secondly, the U–Pb age better conforms with the sample's stratigraphic position when ages from all three stalagmites are aligned to a common depth scale by synchronising their carbon (δ¹³C) and oxygen (δ¹⁸O) stable isotope profiles – noting that there is generally very good agreement between isotopic variations amongst coeval speleothems from the Galleria delle Stalattiti chamber of Corchia Cave (Tzedakis et al., 2018; Bajo et al., 2020). The cause of the anomalously old U–Th age is unclear, but one possibility is that this particular sub-sample was affected by localised open-system behaviour which resulted in post-depositional U loss and thus an older-than-true U–Th age determination. This is known to affect some speleothem samples even if they visually appear pristine (Bajo et al., 2016).

To evaluate agreement between the U–Th and U–Pb ages on a collective basis, we also employed a linear regression approach. If two dating methods are consistent and unbiased, we expect paired age determinations to plot along a 1 : 1 line passing through the origin within their analytical uncertainties. After excluding sample CCB-6-1, which was already found to be inconsistent on a pairwise basis, we fitted a regression line through the age pairs using the weighted least-squares algorithm of York et al. (2004), which accounts for assigned uncertainties in both variables and uncertainty correlation (Fig. 3b). The resulting best-fit regression line had a slope of 0.98 ± 0.07 (2σ), a y-intercept value of 6.0 ± 31 (2σ), and an MSWD (mean square weighted deviation) of 0.50 (n=14, p=0.74). We also assessed goodness of fit to a 1 : 1 line passing through the origin directly, obtaining an MSWD of 0.89 (n=14, p=0.55). These results suggest that the two dating methods are consistent and unbiased over the age interval considered in this study.

One limitation of this approach, however, is that – as noted above – the U–Th age analytical uncertainties for older samples CC15-1 and CC17-1-3 are slightly non-Gaussian, leading to some inaccuracy in regression results obtained by the classical least-squares approach, which assumes strictly Gaussian-distributed uncertainties. Nevertheless, these samples have relatively large age uncertainties that comfortably overlap the best-fit line. We therefore consider the results of this collective comparison to be relatively insensitive to this assumption.

4.2 Age uncertainties

The nature of age uncertainties differs considerably between the two dating methods for the age interval considered in this study. Overall, U–Th age uncertainties are more predictable than the U–Pb age uncertainties and tend to increase in an approximately exponential manner from ∼ 4.5 ka (2σ) at ∼ 450 ka to > 25 ka (2σ) at ∼ 620 ka (see Fig. 4). When U–Th age uncertainties are plotted against U–Th ages, the data tend to follow a simulation line that is based on propagating average ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ and ${[}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}]$ analytical uncertainties obtained in this study through to U–Th age uncertainties. Departures from this trend primarily reflect deviations in analytical uncertainties of one or both of the measured ratios from the average. While the exact trajectory of this line is expected to vary for different analytical protocols, the exponential form is not, provided that initial ²³⁰Th is either negligible or reasonably consistent. The trajectory of the line may also depend somewhat on initial ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values, with lower age precision generally expected to be obtained for samples with lower ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values and especially values less than unity (e.g. Meckler et al., 2012​​​​​​​), all else being equal.

Figure 4 Age precision versus age for U–Th and U–Pb age determinations The orange line in panel (a) shows the simulated U–Th age uncertainty calculated based on the average ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ and ${[}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}]$ measurement uncertainties for samples analysed in this study and an assumed ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}{]}_{\mathrm{i}}$ value of 0.72, which is typical of Corchia Cave stalagmites within this age range. The orange line in panel (b) shows the median U–Pb age uncertainty.

U–Pb age uncertainties, on the other hand, do not show an obvious trend with increasing age and do not correlate directly with average analytical uncertainties in ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{206}}\mathrm{Pb}$ or ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}\mathrm{Pb}$ isotope ratio measurements. Moreover, there is only a relatively small correlation between U–Pb age uncertainties in this dataset and ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ measurement precision or scatter of data about the isochron, quantified by either the MSWD statistic or its robust variant, s (see the Supplement). Instead, U–Pb age uncertainties appear to depend predominantly on the distribution of data along the isochron, consistent with the findings of previous studies (Woodhead et al., 2012; Engel and Pickering, 2022).

Engel and Pickering (2022) examined factors controlling U–Pb isochron age precision in Middle Pleistocene speleothems from multiple cave sites, including Corchia Cave. They found that uncertainty in Tera–Wasserburg U–Pb isochron ages generally correlates strongly with uncertainty in the isochron slope, which is in turn controlled largely by the spread in Pb $/$ U ratios along the isochron line due to variability in inherited Pb content. To quantify the spread of data points along the isochron, Engel and Pickering (2022) defined a metric termed “average isochron distance”, which quantifies the spread of data as the average Euclidean distance between each data point and the centroid of the dataset ($\stackrel{\mathrm{‾}}{x}$, $\stackrel{\mathrm{‾}}{y}$). Here we adopt an alternative metric based on the least-squares-fitted points² along the regression line rather than the nominal measured values of the sub-samples. In this context, the use of the least-squares-fitted points helps to better differentiate between the spread of data along the isochron and the scatter about the isochron (e.g. due to “geological” scatter). It also circumvents scaling issues associated with isochron diagram axes spanning vastly different orders of magnitude, as is the case for the Tera–Wasserburg diagram. The metric we adopt is termed “combined isochron spread” and is calculated as

$$\begin{array}{}\text{(4)}& d_x=\sqrt{\sum _k{\left(\widehat{x_k}-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}},\end{array}$$

where $\widehat{x_k}$ is the least-squares fitted ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{206}}\mathrm{Pb}$ value of the kth data point and $\stackrel{\mathrm{‾}}{x}$ is the centroid of the data points.

In addition to the spread of data points along the isochron, another factor that influences Tera–Wasserburg U–Pb isochron age precision is the location of the data relative to the concordia intercept point. To assess this effect on U–Pb isochron age uncertainties, we define a second metric termed the “average inherited Pb index”, which is calculated as

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{‾}}{{\mathrm{Pb}}_i}={\displaystyle \frac{\stackrel{\mathrm{‾}}{y}-y^{\ast }}{y_{\mathrm{0}}-y^{\ast }}},\end{array}$$

where y^∗ is the ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ value of the concordia intercept point, $\stackrel{\mathrm{‾}}{y}$ is the centroid of the data, and y₀ is the y-axis intercept point (i.e. the estimated inherited ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ value from the fitted isochron line). Essentially, this is the distance in y between the concordia intercept point and the centroid of the data, normalised to the maximum possible distance that could be obtained for a given age, initial ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ isotopic composition, and initial ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ value.

Consistent with Engel and Pickering (2022), we find that the spread of data along the isochron line is controlled mostly by variation in inherited Pb content rather than U content, with the latter tending to remain relatively constant across speleothem growth layers (see Fig. S4 in the Supplement). We find that U–Pb isochron age uncertainties correlate strongly with uncertainty in the isochron slope, excluding the highly radiogenic (i.e. very low inherited Pb) samples CC15-1 and CCB-C-3, which encompass multiple data points plotting close to, or overlapping, the disequilibrium concordia curve (Fig. 5a). In turn, isochron slope precision correlates with combined isochron spread, although sample CCB-C-20, which has an exceptionally low dispersion of data about the regression line (MSWD = 0.15), clearly does not conform to this general trend (Fig. 5b).

Figure 5 Factors affecting U–Pb isochron age precision. (a) U–Pb age uncertainty plotted against isochron slope uncertainty. Markers are coloured according to ${\stackrel{\mathrm{‾}}{\mathrm{Pb}}}_i$ (i.e. the “inherited Pb index”; see Eq. 5). Age uncertainty generally correlates strongly with uncertainty in the isochron slope, except where isochrons are highly radiogenic (i.e. have very low average inherited Pb). (b) Isochron slope uncertainty plotted against the reciprocal of the “combined isochron spread” (d_x) metric (see Eq. 4). Overall, isochron slope uncertainty is correlated with the inverse of combined Pb $/$ U spread. An exception to this is sample CCB-C-20, which exhibits very little data scatter about the isochron. Markers are coloured according to $\sqrt{\mathrm{MSWD}}$, employed here as an indicator of data scatter about the isochron line.

Overall, we find that U–Pb isochron ages of similar average precision can be obtained throughout the ∼ 630–430 ka interval considered in this study provided that the isochron either encompasses data points with reasonable spread in Pb $/$ U ratios or includes highly radiogenic material. In practical terms, there are a number of ways in which samples can be pre-screened for suitability in this regard. A number of options were discussed in Woodhead et al. (2012), but advances in both chemical separation procedures (allowing rapid sample throughput) and the development of in situ analytical techniques (e.g. Roberts et al., 2020​​​​​​​) have now rendered some of these approaches redundant. Currently, the most accurate and time-efficient methods for assessing suitability for U–Pb dating would appear to be either pre-screening via laser ablation ICPMS (Woodhead and Petrus, 2019), U–Th isotopic analysis using a high-throughput protocol (e.g. Hellstrom, 2003), whereby ²³²Th is employed as a proxy for inherited Pb (Woodhead et al., 2006) or a simple reconnaissance ID analysis involving a small number (∼ 3) of sub-samples taken from different parts of the growth domain.

4.3 Comparative of age precision

Given that U–Pb age uncertainties do not show an obvious trend with increasing age over the interval considered in this study (Fig. 4), it is reasonable to characterise average expected U–Pb age uncertainties using either a mean (±10.6 ka, 2σ) or median (±8.5 ka, 2σ) value. In this context, however, the median arguably provides a more useful indicator of expected U–Pb age uncertainty because it is less influenced by outliers, e.g. samples CC2-1 and CC17-1-3, which exhibit anomalously large age uncertainties relative to the other isochrons. For this dataset, median U–Pb age precision exceeds expected U–Th precision at ∼ 520 ka, suggesting that U–Pb dating may generally yield higher-precision ages beyond this point. However, caution is needed in generalising this finding to other study sites or sample types, which may be either less suitable to U–Pb geochronology (e.g. higher average inherited Pb) or more suitable to U–Th geochronology (e.g. similar U and Th content but higher initial ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values). In many cases, it is expected that this crossover point would occur at a somewhat older age.

Furthermore, even in cases where U–Pb dating delivers higher precision for a single age, there may be practical reasons to favour the use of the U–Th method. For example, U–Pb isochron ages are significantly more labour-intensive on a per-age basis (when using precise ID MC-ICPMS and TIMS approaches) compared to U–Th ages because of the need to analyse multiple coeval sub-samples (typically >4). Thus, in applications combining multiple age determination, such as in compiling a depth–age model along the axial length of a speleothem (e,g. Scholz et al., 2012; Bajo et al., 2012), the U–Th method may achieve lower final uncertainties for a given amount of analytical effort, even beyond the U–Th versus U–Pb age precision crossover point.

4.4 Inherited ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ composition

In addition to yielding an age solution, Tera–Wasserburg isochrons yield an estimate of the ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ composition of inherited Pb, assuming this to be invariant across the sub-samples used to construct the isochron (Tera and Wasserburg, 1972; Woodhead et al., 2012). The availability of multiple Tera–Wasserburg isochrons in this study, therefore, provides an opportunity to further assess the consistency or otherwise of the isotopic composition of inherited Pb at Corchia Cave. We find that inherited ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ values are reasonably consistent (weighted mean = 0.8148 ± 0.0012, MSWD = 3.3, n=14, p=0.00) with no apparent trend with age. The elevated MSWD value, however, indicates the data are overdispersed with respect to their analytical uncertainties. If we assume that this overdispersion is due to real variability in the inherited Pb isotopic composition and, furthermore, that this additional component of variability is independent and follows a strict Gaussian distribution (e.g. Ludwig, 2000​​​​​​​), we obtain a weighted average of 0.8149 ± 0.0017 (σ_excess=0.023). However, the assumption of Gaussian-distributed excess scatter is not necessarily justified by the available data here. These results are reasonably close to the average ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ value of Bajo et al. (2020), who obtained a weighted mean value of 0.8134 ± 0.0048 for Corchia Cave speleothems that grew over the interval ∼ 970–810 ka.

The relative consistency of inherited ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ values observed in this study provides an opportunity to assess the effect of fixing Tera–Wasserburg isochron slopes to a common y-intercept value. If we repeat the U–Pb age calculations with the Tera–Wasserburg isochron slopes fixed at 0.8148, all ages remain within their original unanchored isochron age uncertainties, except for sample CCB-F-1, which returns a marginally older U–Pb age, and CC2-1, which returns a significantly older age. It is worth noting that isochron CC2-1 exhibits poor spread in U $/$ Pb ratios and is therefore particularly sensitive to the effect of anchoring the y intercept. Overall, these results suggest that for isochron datasets with limited U $/$ Pb variation, anchoring y intercepts to a well-constrained common ${}^{\mathrm{207}}\mathrm{Pb}{/}^{\mathrm{206}}\mathrm{Pb}$ value may still produce reliable ages at Corchia Cave.

4.5 Speleothem ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ values

In recent years, high-precision U-isotope analyses have revealed significant natural δ²³⁸U variations across a range of different materials. These variations are attributed primarily to U fractionation associated with nuclear volume effects rather than classical mass-dependent fractionation (Fujii et al., 2009) and appear to be most pronounced during redox reactions. Notably, during the reduction of soluble U⁶⁺ to insoluble U⁴⁺, ²³⁸U is concentrated in the U⁴⁺ state, leaving the insoluble phase with a higher δ²³⁸U than the remaining dissolved U. This forms the basis for efforts to employ δ²³⁸U of marine sediments as a palaeoredox proxy (e.g. Andersen et al., 2014​​​​​​​). Processes such as adsorption (e.g. Weyer et al., 2008) and mineral leaching by strong acids (Stirling et al., 2007; Hiess et al., 2012) have also been observed to induce ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ fractionation, while the oxidation of U is not associated with significant fractionation.

Data acquired over the past ∼ 15 years suggest that the average δ²³⁸U of bulk crustal rocks (∼ −0.29 ‰) is practically indistinguishable from that of the deep mantle (∼ −0.31 ‰) (Andersen et al., 2017; Tissot and Dauphas, 2015). The upper mantle is somewhat more enriched in ²³⁸U on average due to recycling of altered oceanic crust (Andersen et al., 2015), whereas modern ocean water and marine carbonates are more depleted in ²³⁸U (∼ −0.39 ‰; Tissot and Dauphas, 2015). Accessory U minerals show more δ²³⁸U variability than bulk crustal rocks, as do materials formed in low-temperature environments (Stirling et al., 2007; Weyer et al., 2008; Hiess et al., 2012).

To date, only a small number of studies have reported δ²³⁸U values for speleothems. Stirling et al. (2007) and Cheng et al. (2013) presented δ²³⁸U values for speleothems from numerous cave systems and observed highly variable δ²³⁸U values ranging from ∼ −0.7 ‰ to +0.4 ‰, although variability within a single speleothem or cave system was significantly lower (Shen et al., 2012). Stirling et al. (2007) also identified a tentative correlation between speleothem δ²³⁸U and measured ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values, suggesting that speleothem ²³⁴U disequilibrium and ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ fractionation may be linked by a common weathering control, albeit involving distinct fractionation mechanisms.

Figure 6 Speleothem δ²³⁸U data from this study and others. Speleothem samples are grouped according to cave site where more than two analyses are available. Speleothem data analysed by Cheng et al. (2013) include SB32 from Sanbao Cave, central China (e.g. Wang et al., 2008), WM1 from Wilder Mann Cave, Northern Calcareous Alps, Austria (Meyer et al., 2009), and MO-1/2/3 from Leaena's Breath Cave, Nullarbor Plain, Australia (Woodhead et al., 2006). “Other speleothems” include the GBW04412 and GBW04413 stalagmite standards from Wang et al. (2023), the FAB-LIG sample from Stirling et al. (2007), and the Kr3 and CA-1 samples described in Cheng et al. (2013). Speleothem samples analysed by Stirling et al. (2007) include SPA-43 from Spannegal Cave, Austria (e.g. Spötl and Mangini, 2010), Argentaorola from Argentarola Island, central Italy (e.g. Bard et al., 2002), and KOZ. Also shown are the “consensus” ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ value of Steiger and Jäger (1977), the average terrestrial zircon value of Hiess et al. (2012), and the average modern ocean water (−0.39 ‰) and crustal rock (−0.29 ‰) values (Tissot and Dauphas, 2015). The global range of δ²³⁸U values for carbonates, Phanerozoic shales, and Fe–Mn crusts/Fe oxides are shown as box plots with “whiskers” indicating the 2.5th and 97.5th percentiles. These data were taken from the Uranium Isotope Database (UID) of Li and Tissot (2023). Absolute ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ values are shown on the right-hand y axis based on the CRM-112A ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ value of 137.837 from Richter et al. (2010).

The use of a high-precision all-Faraday-cup protocol for U–Th measurements in this study, employing ²³³U–²³⁶U double spike for mass bias correction, provides accurate and precise δ²³⁸U determinations. Average analytical uncertainty in measured ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ ratios was 0.0006 (2σ), which is adequate for assessing δ²³⁸U variability at the sub-per mil level (see Table 1). A comparison of our speleothem δ²³⁸U with previous studies (Fig. 6) reveals that most values fall within a range from approximately −0.6 ‰ to −0.2 ‰, consistent with the range observed for carbonate rocks globally (Li and Tissot, 2023). This supports the simple hypothesis that speleothems primarily inherit their ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ composition from carbonate host rocks without significant fractionation. However, two Nullarbor speleothem samples (MO-1/3; Woodhead et al., 2006) and the KOZ samples (Stirling et al., 2007) exhibit elevated δ²³⁸U values, suggesting either a more ²³⁸U-enriched non-carbonate source of U or the presence of an enriching fractionation process during source water U uptake and/or transport to the site of speleothem deposition. Clearly further studies, ideally involving paired speleothem–bedrock analyses, are required to investigate the cause of these high δ²³⁸U values, as well as the primary controls on speleothem δ²³⁸U values more generally.

We also compared the δ²³⁸U and measured ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values of the speleothems analysed in this study with data obtained in previous studies to further assess the relationship between speleothem ²³⁴U disequilibrium and ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ fractionation (see Fig. S5 in the Supplement). With the inclusion of additional speleothem data, we observed a weaker correlation between δ²³⁸U and measured ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values than Stirling et al. (2007), i.e. R²=0.42 versus R²=0.65, and this correlation becomes negligible if the “high-leverage” KOZ data are excluded. However, for speleothem samples that are presumed to have acted as a closed system with respect to U-series isotopes post deposition, it is more appropriate to compare δ²³⁸U values with initial ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values instead of measured activity ratios, since measured ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values also depend on sample age. When repeating this comparison using initial ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ values instead of measured values (necessarily excluding a number of samples for which no age data were available), no correlation is observed (R²=0.01). This result is inconsistent with the hypothesis that speleothem ²³⁴U disequilibrium and ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ fractionation are linked by a common weathering control.

Regardless of the dominant controls on speleothem δ²³⁸U values, the data available thus far show that the U-isotopic composition of most speleothems departs significantly from the conventional ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ value of Steiger and Jäger (1977) and also to a lesser extent the average terrestrial zircon value of Hiess et al. (2012). Using a present-day ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ value of ∼ 137.79 (the weighted average obtained for speleothem samples in this study) in place of the traditional value of 137.88 (Steiger and Jäger, 1977), for either mass bias correction of U–Pb data or in U–Pb age calculation, has a negligible effect on U–Pb ages for young (e.g. Cenozoic) samples. However, when applying the Pb $/$ Pb dating approach to older materials, the ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ value adopted can have a significant effect on age accuracy (Stirling et al., 2007; Weyer et al., 2008; Hiess et al., 2012). This is also the case in U–Th dating using an assumed ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ value for either mass bias correction or calculating ²³⁸U-based ratios from measured ²³⁵U (Stirling et al., 2007; Cheng et al., 2013), although the magnitude of the inaccuracy depends greatly on the specific analytical protocol that is adopted (Shen et al., 2012). The significant departure of speleothem ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ values from 137.88 that is observed here re-iterates the need to pay careful attention to the ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ value that is adopted when dating carbonates under these circumstances.

5 Conclusions

This comparison of the carbonate U–Th and U–Pb chronometers demonstrates excellent agreement between these dating methods over a significant part of the Middle Pleistocene (∼ 630–430 ka). In addition to establishing consistency, we identify key differences in the performance characteristics of the two dating methods. We find the following:

1. U–Th age uncertainties are reasonably predictable (for a given level of analytical precision and initial Th content) and tend to increase in an approximately exponential manner with age. U–Pb age uncertainties do not increase systematically over this interval but are highly dependent on individual sample characteristics, such as Pb $/$ U spread of data points and the availability of highly radiogenic material.

2. The average precision of U–Pb isochron ages surpasses that of U–Th ages prior to the latter reaching its upper age limit. This occurs at ∼ 520 ka for our dataset, although the exact crossover point may vary considerably across different sample types and depositional settings.

3. Even in cases where the precision of U–Pb isochron ages clearly surpasses that of corresponding U–Th age determinations, U–Th dating may still offer some practical advantages, for example, when combining multiple age determinations (e.g. in compiling speleothem depth–age models) because U–Th dating is generally less labour-intensive on a per-age basis.

Overall, these findings highlight the complementary nature of the U–Th and U–Pb dating methods and reaffirm the potential for their combined use in producing accurate and consistent Middle Pleistocene chronologies.

Finally, our assessment of speleothem ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ reveals the following:

1. The majority of speleothem ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ values are consistent with the range observed for carbonate rocks globally, although there are some exceptions.

2. Speleothem ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ values typically depart significantly from the nominal value of 137.88 (Steiger and Jäger, 1977), traditionally adopted as a constant in geochronology. This re-emphasises the need to pay careful attention to the value adopted in data processing and age calculation when it has a significant influence on calculated ages.

3. Further studies are required, ideally involving paired speleothem–bedrock analyses, to better understand the controls on speleothem ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$.

Appendix A: Age equations

A1 U–Th age equation

The standard U–Th age equation, which assumes negligible initial ²³⁰Th and negligible decay of ²³⁸U, may be written as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{a}_{\mathrm{08}}=& \mathrm{1}-e^{-{\mathit{\lambda }}_{\mathrm{230}}t}\\ & +{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{230}}}{{\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{234}}}}(a_{\mathrm{48}}-\mathrm{1})\left(\mathrm{1}-e^{({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{230}})t}\right),\end{array}\end{array}$$

where a₀₈ and a₄₈ denote measured ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ and ${}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$ activity ratios respectively (Broecker, 1963). This equation can be solved for t via iterative methods, such as Newton's method.

A2 U–Pb age equation

The Tera–Wasserburg U–Pb age equation may be written as

$$\begin{array}{}\text{(A2)}& \mathit{\beta }={\displaystyle \frac{n_{\mathrm{75}}}{U}}-\mathit{\alpha }{n}_{\mathrm{68}},\end{array}$$

where α and β are the y intercept and slope of the fitted isochron line respectively, $n_{\mathrm{75}}=\frac{{}^{\mathrm{207}}{\mathrm{Pb}}^{\ast }}{{}^{\mathrm{235}}\mathrm{U}}$ (such that the superscript ∗ denotes radiogenic Pb), $n_{\mathrm{68}}=\frac{{}^{\mathrm{206}}{\mathrm{Pb}}^{\ast }}{{}^{\mathrm{238}}\mathrm{U}}$, and U is the present ${}^{\mathrm{238}}\mathrm{U}{/}^{\mathrm{235}}\mathrm{U}$ ratio. If initial abundances of intermediate nuclides other than ²³⁴U are negligible, then Eqs. (1) and (7) in Pollard et al. (2023) simplify to

$$\begin{array}{}\text{(A3)}& n_{\mathrm{68}}=\mathit{\gamma }+{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{238}}}{{\mathit{\lambda }}_{\mathrm{234}}}}\left(\left(a_{\mathrm{48}}-\mathrm{1}\right)e^{{\mathit{\lambda }}_{\mathrm{234}}t}+\mathrm{1}\right)\mathit{\eta },\end{array}$$

where

$$\begin{array}{c}{\displaystyle }\begin{array}{rl}\mathit{\gamma }=& c_{\mathrm{1}}+c_{\mathrm{2}}{e}^{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{234}})t}+c_{\mathrm{3}}{e}^{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{230}})t}\\ & +c_{\mathrm{4}}{e}^{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{226}})t}+e^{{\mathit{\lambda }}_{\mathrm{238}}t},\end{array}\\ {\displaystyle }\begin{array}{rl}\mathit{\eta }=& h_{\mathrm{1}}{e}^{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{234}})t}+h_{\mathrm{2}}{e}^{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{230}})t}+h_{\mathrm{3}}{e}^{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{226}})t}\\ & +e^{{\mathit{\lambda }}_{\mathrm{238}}t},\end{array}\end{array}$$

and

$$\begin{array}{}\text{(A4)}& n_{\mathrm{75}}=d_{\mathrm{1}}+d_{\mathrm{2}}{e}^{({\mathit{\lambda }}_{\mathrm{235}}-{\mathit{\lambda }}_{\mathrm{231}})t}+e^{{\mathit{\lambda }}_{\mathrm{235}}}.\end{array}$$

Equation (A3) incorporates a re-arranged version of the ²³⁴U age equation that assumes negligible decay of ²³⁸U (e.g. Woodhead et al., 2006), but it is also possible to use a version of this equation that does not make this assumption (e.g. Ivanovich and Harom, 1992, p. 780) instead. In this study it makes little difference which form of the equation is adopted, so the simpler version is used. The coefficients c_i, h_i, and d_i are “Bateman coefficients” given by the following:

$$\begin{array}{c}{\displaystyle }{c}_{\mathrm{1}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{234}}{\mathit{\lambda }}_{\mathrm{230}}{\mathit{\lambda }}_{\mathrm{226}}}{({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{238}})({\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{238}})({\mathit{\lambda }}_{\mathrm{226}}-{\mathit{\lambda }}_{\mathrm{238}})}}\\ {\displaystyle }{c}_{\mathrm{2}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{238}}{\mathit{\lambda }}_{\mathrm{230}}{\mathit{\lambda }}_{\mathrm{226}}}{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{234}})({\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{234}})({\mathit{\lambda }}_{\mathrm{226}}-{\mathit{\lambda }}_{\mathrm{234}})}}\\ {\displaystyle }{c}_{\mathrm{3}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{238}}{\mathit{\lambda }}_{\mathrm{234}}{\mathit{\lambda }}_{\mathrm{226}}}{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{230}})({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{230}})({\mathit{\lambda }}_{\mathrm{226}}-{\mathit{\lambda }}_{\mathrm{230}})}}\\ {\displaystyle }{c}_{\mathrm{4}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{238}}{\mathit{\lambda }}_{\mathrm{234}}{\mathit{\lambda }}_{\mathrm{220}}}{({\mathit{\lambda }}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{226}})({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{226}})({\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{226}})}}\\ {\displaystyle }{h}_{\mathrm{1}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{230}}{\mathit{\lambda }}_{\mathrm{226}}}{({\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{234}})({\mathit{\lambda }}_{\mathrm{226}}-{\mathit{\lambda }}_{\mathrm{234}})}}\\ {\displaystyle }{h}_{\mathrm{2}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{234}}{\mathit{\lambda }}_{\mathrm{226}}}{({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{230}})({\mathit{\lambda }}_{\mathrm{226}}-{\mathit{\lambda }}_{\mathrm{230}})}}\\ {\displaystyle }{h}_{\mathrm{3}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{234}}{\mathit{\lambda }}_{\mathrm{230}}}{({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{226}})({\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{226}})}}\\ {\displaystyle }{d}_{\mathrm{1}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{231}}}{({\mathit{\lambda }}_{\mathrm{231}}-{\mathit{\lambda }}_{\mathrm{235}})}}\\ {\displaystyle }{d}_{\mathrm{2}}={\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{235}}}{({\mathit{\lambda }}_{\mathrm{235}}-{\mathit{\lambda }}_{\mathrm{231}})}}\end{array}$$

Having fitted a suitable isochron line and obtained a ${[}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}]$ measurement that is representative of the material sampled for Pb $/$ U analysis, Eq. (A2) may also be solved for t via iterative methods, such as Newton's method.

Appendix B: Age uncertainties and uncertainty correlation

U–Th and U–Pb age uncertainties and their uncertainty correlation may be computed via first-order uncertainty propagation. A convenient way to set up these calculations is via the matrix approach outlined in McLean et al. (2011). Starting with Eq. (63) in McLean et al. (2011) we have

$$\begin{array}{}\text{(B1)}& \left[\begin{array}{cc}{\mathit{\sigma }}_{t_{\mathrm{Th}}}^{\mathrm{2}}& \mathrm{cov}(t_{\mathrm{Th}},t_{\mathrm{Pb}})\\ \mathrm{cov}(t_{\mathrm{Th}},t_{\mathrm{Pb}})& {\mathit{\sigma }}_{t_{\mathrm{Pb}}}^{\mathrm{2}}\end{array}\right]={\mathbf{J}}^{\top }\mathbf{VJ},\end{array}$$

where J is the Jacobian matrix of derivatives and V is the covariance matrix of the age equation input variables. These matrices may be written in full as

$$\begin{array}{}\text{(B2)}& {\displaystyle }\mathbf{J}=\left[\begin{array}{cc}\frac{\mathrm{d}{t}_{\mathrm{Th}}}{\mathrm{d}{a}_{\mathrm{48}}}& \frac{\mathrm{d}{t}_{\mathrm{Pb}}}{\mathrm{d}{a}_{\mathrm{48}}}\\ \frac{\mathrm{d}{t}_{\mathrm{Th}}}{\mathrm{d}{a}_{\mathrm{08}}}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{d}{t}_{\mathrm{Pb}}}{\mathrm{d}a}\\ \mathrm{0}& \frac{\mathrm{d}{t}_{\mathrm{Pb}}}{\mathrm{d}b}\end{array}\right],\text{(B3)}& {\displaystyle }\mathbf{V}=\left[\begin{array}{cccc}{\mathit{\sigma }}_{a_{\mathrm{48}}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{\sigma }}_{a_{\mathrm{08}}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}& \mathrm{cov}(\mathit{\alpha },\mathit{\beta })\\ \mathrm{0}& \mathrm{0}& \mathrm{cov}(\mathit{\alpha },\mathit{\beta })& {\mathit{\sigma }}_{\mathit{\beta }}^{\mathrm{2}}\end{array}\right].\end{array}$$

The derivatives in the first column of J are given by

$${\displaystyle \frac{\mathrm{d}{t}_{\mathrm{Th}}}{\mathrm{d}{a}_{\mathrm{48}}}}=-{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{1}}}}\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{230}}}{{\mathit{\lambda }}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{234}}}}\right)\left(\mathrm{1}-e^{({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{230}})t}\right)$$

and

$${\displaystyle \frac{\mathrm{d}{t}_{\mathrm{Th}}}{\mathrm{d}{a}_{\mathrm{08}}}}={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{1}}}}$$

such that

$$\begin{array}{}\text{(B4)}& k_{\mathrm{1}}={\mathit{\lambda }}_{\mathrm{230}}\left(e^{-{\mathit{\lambda }}_{\mathrm{230}}t}+(a_{\mathrm{48}}-\mathrm{1})e^{({\mathit{\lambda }}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{230}})t}\right).\end{array}$$

These expressions are equivalent to those given in Ludwig and Titterington (1994). The derivatives in the second column of J are given by

$$\begin{array}{c}{\displaystyle \frac{\mathrm{d}{t}_{\mathrm{Pb}}}{\mathrm{d}{a}_{\mathrm{48}}}}=-{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{238}}}{{\mathit{\lambda }}_{\mathrm{234}}}}\right)\mathit{\alpha }\mathit{\eta }{e}^{{\mathit{\lambda }}_{\mathrm{234}}t},\\ {\displaystyle \frac{\mathrm{d}{t}_{\mathrm{Pb}}}{\mathrm{d}a}}=-{\displaystyle \frac{n_{\mathrm{68}}}{k_{\mathrm{2}}}},\\ {\displaystyle \frac{\mathrm{d}{t}_{\mathrm{Pb}}}{\mathrm{d}b}}=-{\displaystyle \frac{\mathrm{1}}{k_{\mathrm{2}}}},\end{array}$$

where

$$\begin{array}{}\text{(B5)}& k_{\mathrm{2}}=\mathit{\alpha }{\displaystyle \frac{\partial n_{\mathrm{68}}}{\partial t}}-{\displaystyle \frac{\mathrm{1}}{U}}{\displaystyle \frac{\partial n_{\mathrm{75}}}{\partial t}}.\end{array}$$

The partial derivatives in k₂ are

$$\begin{array}{rl}{\displaystyle \frac{\partial n_{\mathrm{68}}}{\partial t}}=& {\displaystyle \frac{\partial \mathit{\gamma }}{\partial t}}+{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{238}}}{{\mathit{\lambda }}_{\mathrm{234}}}}({\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}\left((a_{\mathrm{48}}-\mathrm{1})e^{{\mathit{\lambda }}_{\mathrm{234}}t}+\mathrm{1}\right)\\ & +\mathit{\eta }(a_{\mathrm{48}}-\mathrm{1}){\mathit{\lambda }}_{\mathrm{234}}{e}^{{\mathit{\lambda }}_{\mathrm{234}}t}),\end{array}$$

where $\frac{\partial \mathit{\gamma }}{\partial t}$ and $\frac{\partial \mathit{\eta }}{\partial t}$ are trivial to compute since γ and η are simply composed of sums of exponential terms, and

$${\displaystyle \frac{\partial n_{\mathrm{75}}}{\partial t}}=({\mathit{\lambda }}_{\mathrm{235}}-{\mathit{\lambda }}_{\mathrm{231}})d_{\mathrm{2}}{e}^{({\mathit{\lambda }}_{\mathrm{235}}-{\mathit{\lambda }}_{\mathrm{231}})t}+{\mathit{\lambda }}_{\mathrm{235}}{e}^{{\mathit{\lambda }}_{\mathrm{235}}t}.$$

Having calculated the U–Th and U–Pb age uncertainties, ${\mathit{\sigma }}_{t_{\mathrm{Th}}}$ and ${\mathit{\sigma }}_{t_{\mathrm{Pb}}}$, and the covariance between them, cov(t_Th, t_Pb), via Eq. (B1), the correlation coefficient is simply given by

$$\begin{array}{}\text{(B6)}& \mathit{\rho }={\displaystyle \frac{\mathrm{cov}(t_{\mathrm{Th}},t_{\mathrm{Pb}})}{{\mathit{\sigma }}_{t_{\mathrm{Th}}}{\mathit{\sigma }}_{t_{\mathrm{Pb}}}}}.\end{array}$$