1 Introduction

Carbonate rocks are only a minor component of the continental crust. However, their scientific importance far outweighs their volumetric abundance. Biogenic carbonates and speleothems document the history of life, and of Earth’s climate and environment. To generate detailed time series of past changes from these carbonate archives, an accurate and precise chronological framework is essential. This framework is anchored in absolute time using radiometric dating methods, with two techniques commonly employed for this purpose. ²³⁰Th/U dating is the method of choice for young samples whose ²³⁰Th, ²³⁴U and ²³⁸U activities are out of secular equilibrium (Kaufman and Broecker, 1965; Ludwig, 2003). ²⁰⁶Pb/²³⁸U dating is the default method for older rocks, in which the secular equilibrium between ²³⁰Th, ²³⁴U and ²³⁸U has been restored (Smith and Farquhar, 1989; Roberts et al., 2020).

Ironically, the absence of detectable ²³⁴U/²³⁸U disequilibrium compromises the accuracy of the ²⁰⁶Pb/²³⁸U method. Any initial excess or deficit of ²³⁴U and ²³⁰Th affects the ²⁰⁶Pb/²³⁸U ratio and, hence, the age estimate derived therefrom. In clean, detritus-free carbonates, it is often safe to assume the absence of initial ²³⁰Th. This assumption is not valid for ²³⁴U, which can be enriched (or occasionally depleted) relative to ²³⁸U by physiochemical processes such as (1) α-recoil ejection and preferential leaching of ²³⁴U from α-damaged mineral sites, and (2) chemical fractionation between preferentially oxidised ²³⁴U⁶⁺ and ²³⁸U⁴⁺ (Fleischer, 1982; Porcelli and Swarzenski, 2003).

Corrections for initial ²³⁴U/²³⁸U disequilibrium can be done by either assuming a specific initial ²³⁴U/²³⁸U activity ratio, or by inferring the initial ratio from any measured residual ²³⁴U/²³⁸U-disequilibrium (Richards et al., 1998; Woodhead et al., 2006; Wendt and Carl, 1985; Engel et al., 2019). In Section 2 of this paper, we review the second approach using matrix exponentials. We show that initial ²³⁴U/²³⁸U disequilibrium can bias ²⁰⁶Pb/²³⁸U dates by up to 4 Myr.

Current disequilibrium correction algorithms use a ‘Monte Carlo’ approach to propagate the errors. In Section 3 we will show that this approach can underestimate the analytical uncertainties of ²⁰⁶Pb/²³⁸U dates by an order of magnitude when samples are within a few permil of secular equilibrium (which typically happens before c. 2 Ma). This observation undermines the results of several published studies (Vermeesch et al., 2025).

In Section 4 we introduce a deterministic Bayesian approach to estimate the uncertainties of disequilibrium-corrected ²⁰⁶Pb/²³⁸U dates. We use this alternative algorithm to show that beyond c. 2 Ma, disequilibrium-corrected ²⁰⁶Pb/²³⁸U dates are impractically imprecise, unless highly enriched initial ²³⁴U/²³⁸U activity ratios can be ruled out a priori. The large uncertainty associated with the ²⁰⁶Pb/²³⁸U method degrades its ability to constrain reliable chronologies for carbonates whose initial disequibrium has expired. However, a more accurate approach is available. Following the example of Richards et al. (1998), Neymark and Amelin (2008), Vaks et al. (2020) and others, Section 6 makes a case for the little-used ²⁰⁷Pb/²³⁵U clock as a replacement for the ²⁰⁶Pb/²³⁸U method.

²⁰⁷Pb/²³⁵U isochrons are, essentially, immune to the effects of initial disequilibrium (apart from a minor correction for ²³¹Pa, which becomes smaller with increasing age). In Section 7 we present examples from Siberia and Israel to show that the ²⁰⁷Pb/²³⁵U method is more accurate than the ²⁰⁶Pb/²³⁸U method, whilst being less precise for young samples. Both the Bayesian uncertainty estimation method and ²⁰⁷Pb/²³⁵U isochrons have been implemented in the IsoplotR toolbox for radiometric geochronology (Section 8).

This paper will use the following symbols and notations:

• n₃₈, n₃₄, n₃₀, n₂₆ and n₀₆: the number of atoms of ²³⁸U, ²³⁴U, ²³⁰Th, ²²⁶Ra and ²⁰⁶Pb, respectively;

• λ₃₈, λ₃₅, λ₃₄, λ₃₂, λ₃₁ λ₃₀ and λ₂₆: the radioactive decay constants of ²³⁸U, ²³⁵U, ²³⁴U, ²³²Th, ²³¹Pa, ²³⁰Th and ²²⁶Ra, respectively;

• [4∕8]_t: the present day ²³⁴U/²³⁸U activity ratio (i.e. [λ₃₄ n₃₄]∕[λ₃₈ n₃₈]);

• [4∕8]_i, [0∕8]_i, [1∕5]_i: the initial ²³⁴U/²³⁸U, ²³⁰Th/²³⁸U and ²³¹Pa/²³⁵U activity ratios, respectively;

• [4∕8]_m and s[4∕8]_m: the measured ²³⁴U/²³⁸U activity ratio and its standard error, respectively.

• 123.456 ± 0.012 represents a(n approximate) 95% confidence interval; whereas 123.456(78) is a more succinct notation that indicates a value of 123.456 with a standard error of 0.078.

2 Disequilibrium corrections in a nutshell

Even though the U–Pb decay systems consist of numerous steps (14 for the ²³⁸U→²⁰⁶Pb chain and 11 for the ²³⁵U→²⁰⁷Pb chain), conventional U–Pb geochronology ignores this complexity and the method is mathematically treated as a set of simple parent-daughter pairs. This simplification is justified once a state of secular equilibrium is established between all the intermediate daughter products in the decay chains. Such secular equilibrium is practically reached after 1 to 2 million years. As mentioned in Section 1, any disequilibrium that might exist prior to this secular equilibrium can be used as a chronometer in its own right.

Initial disequilibrium of the U-decay series affects the accuracy of the U–Pb method. For example, ignoring any initial excess ²³⁴U results in an overestimated ²⁰⁶Pb/²³⁸U age, and ignoring any initial ²³¹Pa deficit results in an underestimated ²⁰⁷Pb/²³⁵U age. Therefore, initial disequilibrium is one mechanism to produce discordant U–Pb results. Extreme ²³⁴U-enrichments have been observed in places such as South Africa ([4∕8]_i < 12; Kronfeld et al., 1994), Siberia ([4∕8]_i < 6; Vaks et al., 2020) and Japan ([4∕8]_i < 12; Kuribayashi et al., 2025). Using the [4∕8]_i = 12 value as an upper bound, the maximum effect of initial ²³⁴U/²³⁸U-disequilibrium can be approximated as follows:

(1)

where λ₃₈ = 0.155125(83),Gyr^-1 (Jaffey et al., 1971), λ₃₀ = 9.1705(16)Myr^-1 and λ₃₄ = 2.82206(80)Myr^-1 (Cheng et al., 2013). For old carbonates (> 100 Ma, say), a 4 Myr bias may be inconsequential. However, for young carbonates, the relative effect of initial disequilibrium can result in order-of-magnitude levels of bias. A disequilibrium correction is needed to remove this bias.

If the intermediate daughter is sufficiently long lived and the sample is sufficiently young (t < 5∕λ, say) to retain some of its disequilibrium, then the activity ratios can be back-calculated to the time of isotopic closure (assuming subsequent closed-system behaviour). This strategy applies to ²³⁴U/²³⁸U-disequilibrium and, for very young samples, to ²³⁰Th/²³⁸U-disequilibrium. The complex evolution of the U-decay series was first described by Bateman (1908) and subsequently applied to U–Pb geochronology by Ludwig (1977), Wendt and Carl (1985) and Engel et al. (2019). Here we opt for an alternative formulation, using matrix exponentials (Albarède, 1995). For example, the ²³⁸U→²⁰⁶Pb decay chain can be expressed in matrix form as follows:

(2)

where λ₂₆ = 0.4332(19)kyr^-1 (Audi et al., 2003), and the shortest lived intermediate daughters (< 1 kyr half lives) have been omitted. The solution to Equation 2 is a so-called matrix exponential:

(3)

which expresses the present day amounts of the different isotopes as a function of the initial amounts. An interesting result is obtained by setting t = ∞ in Equation 3 to estimate the activity ratio under secular equilibrium:

(4)

Note that this activity ratio is not exactly equal to unity. This is because 0.0055% of ²³⁸U is lost during ²³⁴U’s mean lifetime of 1∕λ₃₄ = 354kyr. Equation 3 can also be inverted to express the initial amounts as a function of the present day amounts:

(5)

Equations 3 and 5 can be used to construct a concordia diagram in the presence of disequilibrium. If measured activity ratios are used to infer the initial conditions, then the concordia line terminates where those inferred activity ratios reach unrealistic values (e.g., [4∕8]_i = 500 and [0∕8]_i = 0; Figure 1a). Beyond ten or so ²³⁴U half lives, it becomes very difficult to estimate [4∕8]_i from [4∕8]_m, and it is even more difficult to quantify the analytical uncertainty of the disequilibrium correction. In Section 3 we review the current ‘Monte Carlo’ approach to uncertainty estimation for initial disequilibrium correction and in Section 4 we propose an alternative ‘Bayesian’ approach, which offers significant advantages for samples that are close to secular equilibrium.

3 ‘Monte Carlo’ uncertainty estimation

Existing data processing software for disequilibrium-corrected ²⁰⁶Pb/²³⁸U geochronology, such as DQPB (Pollard et al., 2023), estimate the uncertainty of the disequilibrium correction by Monte Carlo simulation. Given a linear array of isotopic data in Tera-Wasserburg space (i.e. n₀₇∕n₀₆ vs. n₃₈∕n₀₆), paired with an activity ratio measurement [4∕8]_m with standard error s[4∕8]_m, this approach works as follows:

1. Draw a random value [4∕8]_t from a normal distribution with mean [4∕8]_m and standard deviation s[4∕8]_m.

2. Fit a straight line to the U–Pb data and find the [4∕8]_i-value and isochron age (t) that are consistent with both the U–Pb measurements and [4∕8]_t. In other words, use Equation 5 to estimate [4∕8]_i from [4∕8]_t and repeat this for different values of t until the linear fit to the U–Pb data is optimised.

3. Repeat steps 1 and 2 until the entire distribution of [4∕8]_t-values has been sampled.

4. If step 2 fails or produces physically impossible results (e.g., t < 0), then ignore the corresponding [4∕8]_t-value. Otherwise add the [4∕8]_i and t-values to a list of acceptable results.

5. Use the spread of the acceptable [4∕8]_i and t-values to quantify their respective uncertainties.

For the purpose of the present study, we have implemented our own version of this algorithm, using R and IsoplotR (Vermeesch, 2018). The only major difference between our code and DQPB is that it does not sample the [4∕8]_m-distribution randomly, but uses a targeted approach to sample [4∕8]_m as a sequence of regularly spaced normal quantiles. This is faster and produces deterministic results that do not depend on the seed of a random number generator. Figure 1 summarises the application of this approach using the ‘Corchia’ dataset of Pollard et al. (2023), producing identical results to DQPB. To reflect this equivalence of outcomes, we will refer to our version of the algorithm as a ‘Monte Carlo’ method, despite the fact that it does not actually use a random number generator.

Next, let us apply the same approach to older materials such as sample AV03 (Bolt’s Farm, South Africa) of Pickering et al. (2019). The uncorrected U–Pb isochron age for this sample is 5.6 ± 0.9 Ma, which is 22 half-lives of ²³⁴U. Consequently, the measured present-day ²³⁴U/²³⁸U activity ratio is statistically indistinguishable from secular equilibrium, at [4∕8]_m = 1.0046 ± 0.0063. Despite the lack of measurable disequilibrium, the ‘Monte Carlo’ approach appears to have successfully applied a disequilibrium correction, resulting in a corrected age that is less than half the uncorrected age, with a precision of better than 12% (Figure 2). How is this possible? The answer lies in the rejected solutions (step 4 of the algorithm), which are marked in black in Figure 2b. Ignoring these ‘physically impossible’ initial ratios suppresses the equilibrium solutions and skews the distribution of Monte Carlo solutions towards high [4∕8]_i-values (Figure 2c) and young ages (Figure 2d).

To demonstrate that the result of Figure 2 is wrong, let us replace the measured ²³⁴U/²³⁸U activity ratio with the equilibrium ratio (Equation 4):

Plugging this value into the ‘Monte Carlo’ algorithm yields an impossible result (Figure 3). It has applied a disequilibrium correction without any actual disequilibrium, by ignoring exactly half of the [4∕8]_t distribution (Figure 3b). This was necessary because, for this old sample, essentially any [4∕8]_t value that is less than the equilibrium ratio would require a negative [4∕8]_i ratio, or a negative isochron age t.

4 A Bayesian approach

The previous section showed that the ‘Monte Carlo’ algorithm produces incorrect results for samples whose [4∕8]_t-distributions fall within, say, three standard errors (i.e., 3 × s[4∕8]_m) from secular equilibrium. One way to address this issue is for the ‘Monte Carlo’ algorithm to issue a warning when the [4∕8]_m value is close to [4∕8]_∞. This is the approach taken by DQPB (Pollard et al., 2023). In this section we introduce an alternative approach that automatically handles problematic cases, without the need to define a nominal ‘applicability cutoff’. Given any values of [4∕8]_m and s[4∕8]_m, our new ‘Bayesian’ algorithm proceeds as follows:

1. Define a prior distribution for [4∕8]_i. In a first instance, we will use a uniform distribution that stretches between a minimum [4∕8]_i-value m = 0 and a maximum [4∕8]_i-value M = 20. However, this uniform distribution can easily be replaced by a more informative prior. One flexible way to capture a diversity of prior information is the logistic normal distribution:

   (6)

   where μ and σ are the location and dispersion parameters of the distribution, respectively. An application of this informative prior will be given in Section 7.

2. Draw a random sample from the prior distribution, carry out a constrained isochron regression and register the resulting age (t) and corresponding [4∕8]_t-value. In other words, use Equation 3 to estimate [4∕8]_t from [4∕8]_i, and repeat this for different values of t until the linear fit to the U–Pb data is optimised. Register the likelihood of this linear fit using the same algorithm as used for regular U–Pb isochron regression (Ludwig, 1998; Vermeesch, 2020).

3. Calculate the likelihood of the inferred [4∕8]_t-values under a normal distribution with mean [4∕8]_m and standard deviation s[4∕8]_m. Combine with the likelihood of the linear fit (obtained in step 2) to produce the ‘posterior’ probability of initial ratios.

4. Repeat steps 2 and 3 to constrain the posterior distributions of [4∕8]_i and t. This can either be done using a Markov chain, or with a targeted approach of appropriately spaced [4∕8]_i values.

Applying this method to the 580 ka Corchia example (Figure 4) yields essentially identical results to the ‘Monte Carlo’ algorithm (Figure 1).

However, when the Bayesian approach is applied to the older Bolt’s Farm data (Figure 5), it produces a very different result than the ‘Monte Carlo’ approach (Figure 2). The posterior distributions for Bolt’s Farm sample AV03 (shown in Figure 5c and d) still have maxima at [4∕8]_i = 9 and t = 2.6 Ma, just like the ‘Monte Carlo’ distributions (Figures 2c and d). But unlike the ‘Monte Carlo’ solution, the result of the Bayesian approach also assigns a significant probability to older ages, including the uncorrected U–Pb date of 5.6 Ma. The similarity of this posterior distribution to the prior distribution reflects the fact that the measured ²³⁴U/²³⁸U activity ratio contains relatively little information. The resulting uncertainties are large, but correctly reflect our ignorance about the true extent of the disequilibrium in this case.

Finally, changing the [4∕8]_m-ratio to the equilibrium value (Figure 5) produces a posterior distribution that is nearly identical to the prior distribution. This means that the likelihood function contains almost no information. In other words: the measured ²³⁴U/²³⁸U activity ratio does not tell us anything about the initial disequilibrium (except that [4∕8]_i < 10). If it cannot be ruled out that the sample may have experienced extreme ²³⁴U/²³⁸U disequilibrium, then it is not possible to undo the effects of this disequilibrium using the modern (measured) ²³⁴U/²³⁸U activity ratio.

5 The case against ²⁰⁶Pb/²³⁸U-dating of old carbonates

The applicability range of the ²⁰⁶Pb/²³⁸U method depends on the [4∕8]_i-ratio and on the precision of the [4∕8]_m-measurements. Although sub-permil level analytical uncertainties can be routinely achieved for individual [4∕8]_m-measurements, the external reproducibility is likely worse than this in most samples. This is due to a combination of two competing factors. First, [4∕8]_i is negatively correlated with U concentration (Osmond et al., 1976; Zhou et al., 2005; Kuribayashi et al., 2025). Second, the U concentration must exhibit large variations to form a statistically robust ²⁰⁶Pb/²³⁸U isochron.

The combination of these two effects has the potential to cause intra-sample variations in [4∕8]_m exceeding the analytical uncertainties. The exact magnitude of the dispersion is unknown because most speleothem U–Pb dating studies report only one or a few [4∕8]_m-values per isochron. Here we will assume a conservative value of s[4∕8]_m = 2o, based on a collection of twelve speleothems analysed by Walker et al. (2006).

Disequilibrium corrections using measured ²³⁴U/²³⁸U activity ratios are only feasible if those activity ratios are statistically distinguishable from secular equilibrium. To turn these conclusions into quantitative guidelines, let us define “statistically distinguishable” as “at least 3 × s[4∕8]_m removed from secular equilibrium”. Using this definition, a sample with [4∕8]_i < 2.7 would become indistinguishable from secular equilibrium after c. 2 Ma. In other words, it is impossible to correct a 2 Ma sample whose [4∕8]_i = 2.7, say. The uncorrected ²⁰⁶Pb/²³⁸U isochron age of such a sample would be 2.6 Ma, corresponding to a bias of 30%. Table 1 shows the outcomes of the same exercise for a range of other ages.

Figure 6 presents a more extensive exploration of the magnitude (panel a) and precision (panel b) of disequilibrium-corrected ²⁰⁶Pb/²³⁸U geochronology assuming the aforementioned 2o reproducibility. Alternative versions of this diagram can be generated by modifying the reproducibility value in the R-code provided in the supplementary information (Vermeesch, 2025).

Based on these considerations, we judge carbonate ²⁰⁶Pb/²³⁸U geochronology to be unreliable beyond c. 1.5 Ma, and impossible beyond c. 2 Ma unless initial ²³⁴U/²³⁸U disequilibrium can be confidently ruled out (Figure 6). However, there is a solution to the conundrum of ²³⁴U/²³⁸U-disequilibrium. This solution is the ²⁰⁷Pb/²³⁵U isochron method (Richards et al., 1998; Engel et al., 2019; Vaks et al., 2020; Vermeesch et al., 2025).

6 A potential ²⁰⁷Pb/²³⁵U fix to ²⁰⁶Pb/²³⁸U’s problems

In the previous section, we showed that the ²⁰⁶Pb/²³⁸U method’s accuracy is hampered by the extreme enrichment (up to double the equilibrium value or more) of ²³⁴U observed in certain groundwaters (e.g., Osmond et al., 1976; Kronfeld et al., 1994; Kuribayashi et al., 2025). This problem can be solved by avoiding ²³⁴U altogether and sidestepping the ²³⁸U–²⁰⁶Pb decay chain in favour of the ²³⁵U–²⁰⁷Pb decay chain (Neymark and Amelin, 2008).

There are two kinds of ²⁰⁷Pb/²³⁵U-isochrons. The simplest kind plots ²⁰⁴Pb/²⁰⁷Pb-ratios against ²⁰⁴Pb/²³⁵U-ratios, defining the following linear relationship:

(7)

where ²⁰⁴Pb is used as a proxy for common Pb. Alternatively, one can also use ²⁰⁸Pb to fulfil this role. This gives rise to a ²⁰⁸Pb_i∕²⁰⁷Pb vs. ²³⁵U/²⁰⁷Pb isochron, where ²⁰⁸Pb_i is the non-radiogenic ²⁰⁸Pb-component (with the decay products of ²³²Th removed). Equation 7 is then replaced with:

(8)

where

(9)

in which λ₃₂ = 0.0495(25)Gyr^-1 (Le Roux and Glendenin, 1963) and the ²³²Th/²⁰⁷Pb-ratio can be obtained from the product of the ²³²Th/²³⁸U, ²³⁸U/²³⁵U and ²³²Th/²⁰⁷Pb ratios. Because Th is insoluble in water, radiogenic ²⁰⁸Pb is often absent from carbonates. Therefore, it is generally safe to assume that ≈.

The ²³⁵U→²⁰⁷Pb method has just one long-lived intermediate daughter, ²³¹Pa, which requires a correction. Due to ²³¹Pa’s short half-life of 32.65 kyr (λ₃₁ = 21.158(71)Myr^-1; Audi et al., 2003), it is generally not possible to measure any remaining disequilibrium in the Myr time range where the ²³⁵U→²⁰⁷Pb method offers a tangible advantage over the ²³⁸U→²⁰⁶Pb method. Therefore, the [1∕5]_i-value must be assumed.

Pa is chemically similar to Th, and insoluble in water. Therefore, ²³¹Pa is always depleted relative to ²³⁵U in carbonates. So whereas [4∕8]_i can vary anywhere between 0 and 12 (Osmond et al., 1976), [1∕5]_i is always less than 1 and can be safely assumed to be zero. In a worst case scenario, in which one assumes [1∕5]_i = 0 but the true activity ratio is [1∕5]_i = 1, this would only bias the ²³⁵Pb/²³⁵U age by a relatively small amount (Table 2).

The degree of potential bias of the ²⁰⁷Pb/²³⁵U method decreases with increasing age, unlike the ²⁰⁶Pb/²³⁸U method, whose bias increases with age (Table 1). In this sense, the ²⁰⁷Pb/²³⁵U and ²⁰⁶Pb/²³⁸U methods are complementary to each other. The ²⁰⁷Pb/²³⁵U is most accurate for samples older than 1 Ma, whereas the ²⁰⁶Pb/²³⁸U is more accurate for samples younger than 1 Ma. Note that the latter is similar to the applicability range of the ²³⁰Th/U method. So one could argue that disequilibrium-corrected ²⁰⁶Pb/²³⁸U dating is of limited use to carbonate U–Pb geochronology (except to infer [4/8]_i; Engel et al., 2019). Although the ²⁰⁷Pb/²³⁵U method outperforms the ²⁰⁶Pb/²³⁸U method at c. 1 Ma in terms of accuracy, its poorer precision means that its potential benefits do not materialise until c. 2 Ma. In the next Section, we will demonstrate this by applying both methods to three different case studies.

7 Case studies

Having made a largely theoretical case against ²⁰⁶Pb/²³⁸U dating and for ²⁰⁷Pb/²³⁵U dating of old carbonates that are suspected to have experienced initial ²³⁴U/²³⁸U disequilibrium, we will now compare and contrast the two chronometers using three practical case studies. The first example will demonstrate the accuracy of the ²⁰⁷Pb/²³⁵U method by showing its consistency with disequilibrium-corrected ²⁰⁶Pb/²³⁸U dates of young (< 2 Ma) samples.

The second example uses ID-TIMS data to serve two purposes. First, it will show that the ²⁰⁷Pb/²³⁵U method produces more accurate, more consistent, and more precise results than the ²⁰⁶Pb/²³⁸U method for older (> 2 Ma) carbonates. Second, it will demonstrate how the Bayesian framework can use prior information to overcome the the inaccuracy of the ²⁰⁶Pb/²³⁸U method.

In the third case study we apply the ²⁰⁷Pb/²³⁵U method to LA-ICP-MS data using ²⁰⁸Pb as a proxy for common Pb. In addition to highlighting a successful application where the ²⁰⁷Pb/²³⁵U method produces demonstrably superior results to the ²⁰⁶Pb/²³⁸U method, this dataset also illustrates limitations of the ²⁰⁷Pb/²³⁵U approach with a sample that yields a precise ²⁰⁶Pb/²³⁸U isochron and an unusably imprecise ²⁰⁷Pb/²³⁵U isochron.

7.1 ID-ICP-MS data from Siberia

A rich dataset of 72 speleothem dates is available from the Botovskaya and Ledyanaya Lenskaya (LLC) caves in Siberia. Vaks et al. (2020) used the U–Pb method to extend an important palaeoclimatological archive from these caves that was previously dated using the ²³⁰Th–U disequilibrium method (Vaks et al., 2013a). Samples were analysed by isotope dilution ICP-MS and were found to exhibit a significant level of ²³⁴U/²³⁸U disequilibrium. U–Pb ages were estimated using a two-step procedure. First, the common-Pb contribution was removed by two-point isochron regression through an inherited composition that was inferred by inspection of apparent linear trends in Tera-Wasserburg concordia space. Second, a ²³⁴U/²³⁸U disequilibrium correction was applied to the radiogenic end-member composition, using the procedures described in Section 2. This correction combined the measured ²³⁴U/²³⁸U activity ratios with an assumed absence of initial ²³⁰Th and ²³¹Pa.

The inferred [4∕8]_i-values were c. 2 and 3–5 for LLC and Botovskaya cave, respectively. This corresponds to an age correction of 15% for LLC and 60% for Botovskaya cave (Figure 7). Uncertainties were estimated using the Bayesian procedure of Section 4, making the optimistic assumption that the analytical uncertainty of the [4∕8]_m-measurements faithfully captures all sources of dispersion. The scatter of the [4∕8]_i-values suggests that this may not be the case. This caveat notwithstanding, the disequilibrium-corrected ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U ages overlap within uncertainty in all but four of the samples.

The ²⁰⁷Pb/²³⁵U age uncertainty is invariably larger than the ²⁰⁶Pb/²³⁸U age uncertainty. In fact, below c. 1 Ma, it could be argued that the ²⁰⁷Pb/²³⁵U age uncertainties are unusably imprecise (s[t]∕t > 50%). However, above c. 1 Ma, the uncertainty reduces to acceptable levels (s[t]∕t < 5%). Extrapolating this trend further into the past confirms the earlier assertion that beyond c. 2 Ma, the ²⁰⁷Pb/²³⁵U method outperforms the ²⁰⁶Pb/²³⁸U method in both accuracy and precision.

7.2 ID-TIMS data for ASH-15

ASH-15 is a carbonate U–Pb dating reference material sourced from a flowstone in Ashalim cave of southern Israel (Nuriel et al., 2021). 37 ID-TIMS measurements were obtained from the flowstone, including 12 from horizon D and 25 from horizon K. The latter are shown in Figure 8. Nuriel et al. (2021) report an uncorrected semitotal-Pb/U isochron age of 2.965 ± 0.011 Ma for ASH-15.

No ²³⁴U/²³⁸U activity ratio measurements are available for ASH-15D and ASH-15K. However, two other horizons of the same flowstone (ASH-15A+B and ASH-15C1) are characterised by [4∕8]_m-values of 0.99939 ± 0.00108 and 0.99925 ± 0.0015, indistinguishable from secular equilibrium. Younger flowstones in Ashalim cave yield an average [4∕8]_i-value of 1.0470 with a standard deviation of 0.01492 (Vaks et al., 2010, 2013b). A wider survey of 904 speleothem samples dated in southern and central Israel by Chaldekas et al. (2022) have average [4∕8]_i-values of 1.081 ± 0.138. Despite this lack of observable disequilibrium, Mason et al. (2013) suggest a [4∕8]_i-value of 1.5–2.0 to explain the minor degree of discordance of the common-Pb corrected Tera-Wasserburg ratios.

To investigate the effect of initial U-series disequilibrium on ASH-15, Figures 8a and b apply the Bayesian inversion algorithm to the ASH-15K data, using the [4∕8]_m-value of ASH-15C1 and assuming that [0∕8]_i = 0 (i.e. no initial ²³⁰Th). In a first attempt, we will use the same uniform prior from m = 0 to M = 20 as before. This results in a disequilibrium-corrected ²⁰⁶Pb/ ²³⁸U-isochron age of 3.47 + 0.053∕ - 0.794 Ma, and an inferred [4∕8]_i ratio of 0.067 + 2.23∕ - 0.052. As expected, initial equilibrium is very plausible with P([4∕8]_i > 1.0) = 0.47. The initial activity ratio preferred by Mason et al. (2013) cannot be ruled out either, but is less likely, with P([4∕8]_i > 1.5) = 0.25 and P([4∕8]_i > 2.0) = 0.09.

In a second attempt, we used the [4∕8]_i-values of Chaldekas et al. (2022) to construct an informative prior, using the logistic normal formulation of Equation 6 with m = 1, M = 3, μ = 1.081 and σ = 0.2. This produces a posterior distribution for [4∕8]_i that is, essentially, identical to the prior, confirming that the [4∕8]_m data carry virtually no additional information. The corresponding disequilibrium-corrected ²⁰⁶Pb/ ²³⁸U-isochron age is 3.107 ± 0.065 Ma.

In contrast with the widely varying scenarios for the ²⁰⁶Pb/ ²³⁸U method, the ²⁰⁷Pb/ ²³⁵U isochron age calculation is straightforward:

1. An uncorrected ²⁰⁷Pb/ ²³⁵U isochron age of 3.039 ± 0.068 Ma, assuming [1∕5]_i ≈ 1.

2. A corrected ²⁰⁷Pb/ ²³⁵U isochron age of 3.086 ± 0.068 Ma, assuming [1∕5]_i = 0.

which are nearly identical to the ²⁰⁶Pb/²³⁸U date using the informative [4∕8]_i-prior. We would like to conclude the discussion of ASH-15 by remarking that disequilibrium issues do not affect the suitability of this sample as a reference material for in-situ U–Pb geochronology. This is because standardisation is done relative to the uncorrected isotopic composition (Horstwood et al., 2016).

7.3 LA-ICP-MS data from Siberia

For the final case study, we return from Israel to the Botovskaya cave deposits in Siberia. Section 7.1 and Figure 7 show abundant evidence that this cave is strongly enriched in initial ²³⁴U, with [4∕8]_i-values ranging from 3 to 5 according to the results of Vaks et al. (2020). The effect of this strong disequilibrium can be confidently undone for the young (< 500 ka) speleothems shown in Figure 7. However, Botovskaya cave also contains speleothems that are considerably older than this, going all the way back to the Plio-Pleistocene. By now it should be clear that the ²⁰⁶Pb/²³⁸U method is ill suited to unlock this older archive. Figure 9 summarises some preliminary U–Pb results from two of these older cave deposits (sample SB-1625-22 and sample SB-72-8) obtained by LA-ICP-MS.

As mentioned before under the discussion of the Bolt’s Cave data, the ²⁰⁴Pb measurements produced by this technique are imprecise and potentially inaccurate. In this case, ²⁰⁸Pb was measured and can be used as a substitute for ²⁰⁴Pb. Th/U ratios (also measured by LA-ICP-MS) were extremely low, allowing us to ignore the radiogenic ²⁰⁸Pb contribution.

The uncorrected ²⁰⁶Pb/²³⁸U isochron age of sample SB-1625-22 is 2.66 ± 0.10 Ma and exhibits significant overdispersion with respect to the analytical uncertainties (MSWD=36). Model-3 isochron regression (sensu Vermeesch, 2024) indicates that this excess scatter is equivalent to an age dispersion of 104 ± 28 kyr (Figure 9a). In reality, the excess dispersion around the isochron is unlikely to reflect diachronous isotopic closure. A more likely explanation for the scatter of the ²⁰⁶Pb/²³⁸U data is spatial variability of the [4∕8]_i-values, as discussed in Section 5. In summary, the actual dispersion estimate probably has no physical meaning, but the isochron age should be as accurate as mathematically possible. The reduction of the scatter around the Botovskaya isochrons (Figures 9a and c) towards the y-intercept also suggests that the common Pb-ratio is not substantially correlated with the postulated heterogeneous initial disequilibrium.

Given the antiquity of the sample and the difficulty of measuring [4∕8]_m by LA-ICP-MS, no initial disequilibrium measurement was made. Switching from ²⁰⁶U/²³⁸U to ²⁰⁷Pb/²³⁵U isochron space lowers the age to 1.60 ± 0.10 Ma whilst reducing the dispersion of the data around the isochron line (MSWD=1.5, Figure 9b). To verify the accuracy of this result, it is useful to point out that a [4∕8]_i-value of 3.9 would bring the corrected ²⁰⁶Pb/²³⁸U isochron in alignment with the ²⁰⁷Pb/²³⁵U isochron. Such a value is consistent with the initial activity ratios of the more recent Botovskaya deposits (Figure 7b). This not only supports the accuracy of the ²⁰⁷Pb/²³⁵U isochron results, but also suggests that the [4∕8]_i ratios have remained stable over hundreds of thousands of years.

We would like to conclude the results section by drawing attention to the fact that the ²⁰⁷Pb/²³⁵U method is not always successful. Figures 9c and d show that Botovskaya sample SB-72-8 produces a well defined linear array in ²⁰⁶Pb/²³⁸U isochron space, but fails to do so in ²⁰⁷Pb/²³⁵U isochron space. Such cases are not rare. The ²⁰⁷Pb/²³⁵U approach only works in samples that are sufficiently rich in U and sufficiently poor in common Pb. Speleothems from Botovskaya are rich in U (30–170 ppm for SB-72-8), so that here the problem seems to originate from common Pb (0.4–4 ppm for SB-72-8).

8 Implementation in IsoplotR

All the methods described in this paper have been implemented in the IsoplotR toolbox for geochronological data processing (Vermeesch, 2018). The matrix exponential disequilibrium correction method of Section 2 has been part of IsoplotR since version 3.0, whereas the deterministic Bayesian uncertainty estimation routine of Section 4 was introduced in version 5.2. At the time of writing, IsoplotR (version 6.7) supports twelve different U–Pb data formats. Disequilibrium corrected U–Pb isochron regression is available for all these formats, in different forms.

Formats 1–3 contain neither ²⁰⁴Pb nor ²⁰⁸Pb. Therefore, isochron regression for these formats must be done by semitotal-Pb/U regression in Tera-Wasserburg concordia space. Formats 4–6 include ²⁰⁴Pb as a common Pb tracer. These formats permit the calculation of both ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U isochrons, either jointly (by three-dimensional total-Pb/U isochron regression; Ludwig, 1998) or separately. To take full advantage of the ²⁰⁷Pb/²³⁵U method’s superior accuracy, it is recommended to use the two-dimensional option. Formats 7 and 8 use ²⁰⁸Pb as a common Pb tracer. They are also amenable to both ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U isochron regression, either jointly (by total Pb/U–Th regression; Vermeesch, 2020) or separately. Formats 9–10 and 11–12 are simplified versions of formats 4–6 and 7–8, respectively, which only permit two-dimensional regression. Formats 9 and 11 are meant for ²⁰⁶Pb/²³⁸U isochron regression, whereas formats 10 and 12 are meant for ²⁰⁷Pb/²³⁵U isochron regression.

The disequilibrium corrections can be accessed from IsoplotR’s GUI (either online or offline) by using the ‘isochron’ function and ticking the ‘apply disequilibrium correction’ check box in the options menu. Alternatively, the same functionality can also be accessed from the command-line API. Bayesian uncertainty estimation is possible using either interface, but visualising the posterior distributions of the parameter space is currently only possible from the command line.

The online supplement provides all the R code that was used to reproduce the figures in this paper (Vermeesch, 2025).

9 Conclusions

In this paper, we presented a critical appraisal of carbonate U–Pb geochronology, and proposed three improvements to the technique. First, we introduced a matrix exponential solution to the initial disequilibrium problem, extending the work of Albarède (1995). This formulation produces identical results to the conventional solution by Engel et al. (2019), but can be written out more succinctly, and can more easily be modified to suit other problems. For example, the matrix exponential approach can be adjusted to calculate disequilibrium-corrected U–Th–He ages (Farley et al., 2002; Danišík et al., 2017). Second, we presented a deterministic Bayesian algorithm to quantify the statistical uncertainty associated with the disequilibrium correction. This algorithm was used to demonstrate that, for samples older than c. 2 Ma, disequilibrium-corrected ²⁰⁶Pb/²³⁸U geochronology is unreliable. Third, we advocated the use of the ²⁰⁷Pb/²³⁵U isochron method as a more accurate alternative to the ²⁰⁶Pb/²³⁸U method.

Although our findings are most relevant to young carbonates, the inaccuracy of the ²⁰⁶Pb/²³⁸U method equally applies to old samples. Only the relative difference between the ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U ages reduces with time. The absolute difference remains constant at up to 4 Myr (Equation 1). The corresponding systematic uncertainty cannot be removed without making unverifiable assumptions about the initial ²³⁴U/²³⁸U activity ratio. For samples older than > 100 Ma, say, the systematic error caused by initial disequilibrium is generally smaller than the random errors associated with the isotope ratio measurements. However, given a sufficiently precise set of isochrons, it is theoretically possible to reconstruct the U-disequilibrium conditions at the time of isotopic closure from the difference between the ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U clocks.

Engel et al. (2019) advocate using the same procedure in Quaternary studies. They propose a two-step procedure, whereby the difference between the ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U dates is used to estimate [4∕8]_i; and this [4∕8]_i-value is then used to calculate a corrected ²⁰⁶Pb/²³⁸U-age. IsoplotR implements a one-step algorithm that achieves the same goal using the total-Pb/U algorithm of Ludwig (1998) and the total-Pb/U–Th algorithm of Vermeesch (2020). However, we would like to add a note of caution about the usefulness of this joint regression procedure. Beyond c. 2 Ma, all the age-resolving power of the paired ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U approach resides in the ²⁰⁷Pb/²³⁵U clock, so the ²⁰⁶Pb/²³⁸U data add no value.

Code availability

IsoplotR is free software released under the GPL-3 license. The package and its source code are available from https://cran.r-project.org/package=IsoplotR (last access: July 21, 2025, Vermeesch, 2018). The test data used in Section 8 are provided in the supplementary information (Vermeesch, 2025)

Author Contributions

NM developed the matrix exponential solution to the initial disequilibrium problem (McLean et al., 2016); PV developed the deterministic Bayesian inversion algorithm; RP proposed the ²⁰⁷Pb/²³⁵U fix to the ²⁰⁶Pb/²³⁸U problems; AV, TG and SB provided the Siberian data; PV wrote the paper, with feedback from the other authors.

Competing Interests

PV and NM are Associate Editors of Geochronology.

Acknowledgements

This research has been supported by the Natural Environment Research Council (grant no. NE/T001518/1, awarded to PV) and the Leverhulme Trust (grant no. RPG-20202-334, awarded to SFMB). This paper benefited from thorough reviews by Perach Nuriel, Timothy Pollard and Robyn Pickering.

References