Speleothem Dating

Today I’m doing a short note on nonlinear equation solving of speleothem equations.

As a preface to this, I’d like to comment that I’ve not made an overall post placing speleothems in context. And they have some very interesting properties as proxies that should not be lost sight of either in the context of posts like this one or in the context of such proxies being innocent victims of data mining methods in articles such as Mann et al 2008. On a Pleistocene context, dO18 data is extremely valuable and speleothems provide such important information that is independent of other contexts. They provide information on the tropics and the monsoons – while there is still some hair on the interpretation of the proxies, it would be grossly unfair to say that current interpretation issues render the enterprise hopeless. The ability of speleothems to record interesting information on the Eemian which is independent of ice cores and foraminifera looks like a success story in my opinion, though I’ve only browsed the topic.

In another post, I queried the seeming disconnect between speleo records from China and the Dasuopu ice core record. Although some readers seemed to take this as a shot on speleos, it wasn’t really. The question is pointed much more at the dating of the Dasuopu ice core, about which I’ve had some lingering questions for some time- which have been totally frustrated by Lonnie Thompson’s failure to provide a proper archive of this important data – a failure that is now over 20 years for some cores.

Some day, I’ll try to discuss these things, but, unfortunately, I can’t cover everything. Since I have a longer term interest in the dating of the Dasuopu ice core, I want to ensure that I’m fully cognizant of the details of speleo dating methods and have therefore tried to replicate reported results using my own nonlinear solver and, at the same time, get a better appreciation of whether the sorts of potential errors in this method could be highly autocorrelated in comparisons of one core to another (a la Socotra-Qunf).

I’m going to do a couple of instalments on these methods (which will all be in this post), but today I’m only going to report on the methods of Edwards 1987, a seminal article for U-Th dating, kindly recommended to me and forwarded to me by Jud Partin.

Edwards 1987
U-Th uses a number of constants – and the constants used in a calculation are not always stated in the article. Edwards uses the following decay constants (the decay constants corresponding to half-lives): Th230: 9.195E-6; U234: = 2.835E-6; U238: 1.551E-10 and a “natural” U234/U238 ratio given by the ratio of the latter two decay constants (reported as 5.472E-5).

lambda=NULL
lambda[“Th230″]= 9.195E-6### ##text 108750 # 2 sigma error 850 years #.8%
lambda[“U234″]= 2.835E-6 ##352740 #2 sig error 710 #.2%
lambda[“U238″]= 1.551E-10 ##.4464E9 #error .0069E9
lambda[“U234:U238″] = lambda[3]/lambda[2] ## 5.472E-5 #Table 4

I manually transcribed the data in Edwards’ tables of sample calculations and uploaded it to CA, where it can be read and converted to consistent units as follows:

A=read.csv(“d:/climate/data/speleothem/edwards.1987.csv”,skip=2,header=TRUE)
A[,3]=A[,3]/10^3;A[,4]=A[,4]/10^5;A[,5]=A[,5]/10^5;A[,9]=A[,9]/10^5
A=A[1:(nrow(A)-1),] #ocean row removed

First, Edwards calculated \delta^{234}U[0] by comparing the ratio of observed U234/U238 to the “equilibrium” ratio, subtracting 1 and multiplying by 1000:

del_U234_0= (A$U234_U238/lambda[“U234:U238″] -1)*10^3
cbind(del_U234_0,A$del_U234.0,1-del_U234_0/A$del_U234.0)
#this matches up to rounding (.004%)

My replication results matched up to rounding (within 0.005%).

Next, Edwards calculated the ratio of Th230 to U238.

Th230_U238 = A$Th230_Th232 * A$Th232/A$U238
cbind(Th230_U238,A$Th230_U238,1- Th230_U238/A$Th230_U238)

For 7 of 9 samples, my calculation matched within rounding (0.005%), but for samples E-T-2 and E-L-3, the differences were about 3%. I have no idea why.

Next, Edwards calculated the activity of Th230 to U238 – the activity being the product of the concentration by the decay constant. Thus:

act.Th230_U238 = A$Th230_Th232*A$Th232/A$U238 * (lambda[“Th230″]/ lambda[“U238″])
cbind(act.Th230_U238,A$act.Th230_U238,1-act.Th230_U238/A$act.Th230_U238 )
#less than 0.005% up to 7-8 where 2-3% off

The results were driven by the previous calculation: they were within rounding for 7 of 9 samples, but, for samples E-T-2 and E-L-3, the differences were about 3%

Edwards now presents the following two equations involving [Th230/U238]_act and del_U234[0]:

Edwards observes that equation (1) yields a value for T, the age of the coral in his case, and equation (2) yields information about the starting system. Equation (1) solves easily using nonlinear methods and I accordingly wrote a short routine to solve equation (1) directly – my routine in this case minimizing the difference between the LHS and RHS of equation 1 – which will minimize at 0 for the samples at hand. (I set interval limits of 0 and 200000 for this experiment, but the optimization can be a little touchy on interval setting and this needs to be controlled in a less overt way on another occasion.) The function implementing equation (1) is shown below, with several lines differentiated to make the accounting easier:

f=function(t) {
A=(1- act)
B= exp(-lambda[“Th230″]*t)
C= del*10^-3* (lambda[“Th230″]/(lambda[“Th230″]-lambda[“U234″]))
D= (1- exp( (lambda[“U234″]-lambda[“Th230″])*t) ) #correcting per reader comment
f= (A- (B- C*D) )^2
f}

I then tried to replicate the Edwards results:

K=nrow(A); age= rep(NA,K)
for(i in 1:K) { act=A$act.Th230_U238[i];del=A$del_U234.0[i];
age[i]= optimize(f,interval=c(0,200000),tol=1E-9)$minimum}
cbind(round(age,0),A$Age,1-age/A$Age )

This yielded the following results (after I fixed the sign as described in a comment below). The calculation is done two ways – first with the values in Edwards Table and second with non-rounded values carried through.

Sample

Replication-Table

Replication (w/o rounding)

Edwards

Diff 1

Diff 2

TAN-E-1g

180

181

180

0

-1

CWS-F-1

846

849

845

-1

-4

CH-8

8293

8288

8294

1

6

AFS-12-A

121984

122429

122100

116

-329

AFS-12-B

122632

121612

122700

68

1088

AFS-12-C

124529

123516

124500

-29

984

E-T-2

129977

136896

129900

-77

-6996

E-L-3

125605

118223

125500

-105

7277

VA-1

NA

NA

NA

NA

NA

Here’s a version of the same calculation using more recent lambda values (from Cheng 2000 as reported in Cobb et al 2003. The constants in these calculations do matter – so one needs to know what were used. (They are not always reported.)

Sample

New Constants

Edwards

Difference

TAN-E-1g

174

180

6

CWS-F-1

816

845

29

CH-8

7956

8294

338

AFS-12-A

112844

122100

9256

AFS-12-B

112138

122700

10562

AFS-12-C

113786

124500

10714

E-T-2

125257

129900

4643

E-L-3

109192

125500

16308

VA-1

NA

NA

NA

The Edwards calculations are of historical interest only. I’ve got some work in progress on more recent studies, where the procedures tend to be a little more poorly documented and harder to figure out, but I’m working through them.

PS. I looked quickly at Isoplot documentation. Here’s what puzzles me right now – Edwards equation (1) solves without getting involved with isochrons, so I don’t understand why they would enter the picture for this particular system.

Reference:
Edwards, R. L., J. Chen, and G.J. Wasserburg. 1987. 238U-234U-
230Th-232Th systematics and the precise measurement of
time over the past 500,000 years. Earth and Planetary Science Letters 81: 175-192. url


23 Comments

  1. David L Hagen
    Posted Dec 12, 2008 at 11:06 AM | Permalink

    There appears to be a trend in the ratios. This suggests that identifying that trend may improve the equation, or help identify the difference between Edward’s method and his equations.

  2. mugwump
    Posted Dec 12, 2008 at 11:24 AM | Permalink

    There’s an exact solution to equation (1) if you make an approximation (that I don’t know is valid):

    Let A be the LHS of (1).
    B = lambda_230.
    C = lambda_234.
    D = delta_234 U(0) / 1000.
    K = D * B / (B-C)

    With these definitions, after some manipulation of equation (1) I get:

    log(A + K) = -BT + log(1 + K exp(CT))

    If K exp(CT) >> 1 then log(1 + K exp(CT)) ~ logK + CT and we can solve for T:

    T ~ log(1 + A/K) / (C-B)

    Is K exp(CT) >> 1 ?

    [If 1 >> K eCT we can use log (1 + x) ~ x. But I have no idea about the relative magnitudes of these numbers.]

  3. mugwump
    Posted Dec 12, 2008 at 11:25 AM | Permalink

    Hmm: sup tags got lost. “eCT” in the previous post is exp(CT).

  4. Steve McIntyre
    Posted Dec 12, 2008 at 11:28 AM | Permalink

    Here’s one way of getting a value around 184 for the first data point. In my f -function, I changed the tolerance from 1E-9 to 500, and then ran

    optimize(f,interval=c(0,200000),tol=500)

    this gave 184.785. But the other results were still way different than Edwards.

  5. Nick S
    Posted Dec 12, 2008 at 12:20 PM | Permalink

    Steve,

    You state that:
    “Equation (1) solves easily using nonlinear methods and I accordingly wrote a short routine to solve equation (1) directly – my routine in this case minimizing the difference between the LHS and RHS of equation 1 – which will minimize at 0 for the samples at hand”.

    I assume by this that you mean that you know that the global minimum of the function you wrote has an objective value of 0. How do you arrive at this conclusion?

    Furthermore, assuming that you know the above to be true, check that the solution that you’re getting is indeed satisfying the original problem (i.e. that the LHS = RHS). The reason I say this is because the optimization problem that you’re solving is a non-convex optimization problem, meaning that many local optima exist. I don’t know what nonlinear solver you’re using to optimize your function, but I strongly suspect that it’s gradient based (as most nonlinear solvers are). If so, then there’s no theoretical guarantee that the point you end up with is the global optimum (i.e that your objective function optimal value is equal to 0, and thus that the LHS = RHS in your original problem).

    Finally, it may be the case that even if your solution does indeed satify the condition that the LHS = RHS (i.e. that your optimal objective value is equal to 0), you may still get a different solution to your problem (that that reported in the literature for instance) because your function may take multiple optimal solutions (i.e. many solutions that have an optimal objective function value of 0). Just by looking at your function, I don’t know if its global minimum is unique.

    Steve: I realize that the function is non-convex but it is convex in the interval that I’ve specified. I noted this issue in the thread itself. So it’s an issue for a more general formulation of the problem but not for the testing done here.

  6. The walrus
    Posted Dec 12, 2008 at 12:21 PM | Permalink

    Hmm..
    I might be missing something but the line
    D= (1- exp( -(lambda[“U234″]-lambda[“Th230″])*t) )
    seems to have an extra minus
    Shouldn’t it be
    D= (1- exp( (lambda[“U234″]-lambda[“Th230″])*t) ?

  7. SB
    Posted Dec 12, 2008 at 12:47 PM | Permalink

    There’s an easier approximation for the first couple of table entries. You can use Exp(k*T) = 1 + k*T when k*T is small. Since the magnitude of k is about 10-5 in both exponentials, for T much less than this you get a nice straight line. Try plotting RHS-LHS against T to confirm.

    It might be handy if you gave the values for del and act that gave the table entries, rather than to have to download/calculate them.

    I did wonder if he was using more precision than he had put in his tables, but that doesn’t seem likely.

  8. The walrus
    Posted Dec 12, 2008 at 3:07 PM | Permalink

    For what it’s worth. After removing the minus from the exponential, the code in this post produces the following output:

    [,1] [,2] [,3]
    [1,] 180 180 0.0001052673
    [2,] 846 845 -0.0014062701
    [3,] 8293 8294 0.0001219000
    [4,] 121984 122100 0.0009533586
    [5,] 122632 122700 0.0005561474
    [6,] 124529 124500 -0.0002318603
    [7,] 129977 129900 -0.0005933254
    [8,] 125605 125500 -0.0008329258
    [9,] 200000 NA NA

  9. Steve McIntyre
    Posted Dec 12, 2008 at 4:17 PM | Permalink

    #8. Aha, nice to have such excellent readers. That’s exactly right. On the plus side, this emulation now works pretty well without needing anything convoluted from isoplot. I’ve amended my tables to reflect this debugging.

  10. Steve McIntyre
    Posted Dec 12, 2008 at 4:51 PM | Permalink

    Now that this calculation is properly benchmarked (thanks to walrus), I’ve re-run the Edwards table with the constants reported in Cobb et al 2003, which turn out to make a noticeable difference. One needs to know the constants in order to compare things.

  11. Steve McIntyre
    Posted Dec 12, 2008 at 4:58 PM | Permalink

    I collated Cobb et al 2003 Table 2 and was able to replicate the results pre-Thorium correction with the following script (it looks like they used 2000 as benchmark year to convert from age to yea.

    ##input constants
    lambda=NULL
    lambda[“Th230″]= -log(.5)/75690 # 9.157711e-06 # 0.5=exp(-k*245240);
    lambda[“U234″]= -log(.5)/235240 # 2.946553e-06 #
    lambda[“U238″]= 1.55125E-10
    lambda[“Th232″]= -log(.5)/ 1.41E10 # 4.915937e-11 ;
    natural=NULL
    natural[“Th230/Th232″]= 4.4E-6#Cobb table 2 note i Th232:U238 = 3.8

    ##read data
    A=read.table(“http://data.climateaudit.org/data/speleothem/cobb03.csv”,skip=1,header=TRUE)
    A[,c(3,4)]=A[,c(3,4)]/10^3;A[,6]=A[,6]/10^6
    A$age=2000-A$Date_uncorr

    ##replicate age calculation
    A$del_U234_0= (A$del_U234_T)* exp(- lambda[“U234″]*A$age)
    ##this works backwards using Equation (2) to get del_U234_0
    range(A$del_U234_0) # 142.2236 148.2506

    f=function(t) {
    A=(1- act)
    B= exp(-lambda[“Th230″]*t)
    C= del*10^-3* (lambda[“Th230″]/(lambda[“Th230″]-lambda[“U234″]))
    D= (1- exp( (lambda[“U234″]-lambda[“Th230″])*t) )
    f= (A- (B- C*D) )^2
    f}

    K=nrow(A); age=rep(NA,K)
    for(i in 1:K) { act=A$act.Th230_U238[i];del=A$del_U234_0[i];
    age[i]= optimize(f,interval=c(0,200000),tol=1E-6)$minimum
    }
    names(age)=A$Sample

    ##report the results
    test=data.frame(names(age),round(age,0),A$age,round(age- A$age,0) )
    dimnames(test)[[2]]=c(“Sample”,”Replication”,”Cobb”,”Diff”)
    test

  12. Geoff Sherrington
    Posted Dec 12, 2008 at 6:39 PM | Permalink

    and a “natural” U234/U238 ratio given by the ratio of the latter two decay constants (reported as 5.472E-5).

    It is a fair assumption, because we know that a large work input is required for enrichment, i.e. to change the natural ratio of uraium isotopes (in the enrichment case, U238:U235). Plus, both are uranium and geochemically, they behave quite similarly.

    It is less certain that the U234:Th230 ratio can be assumed invariant, but I lack examples where it has differed. The reason is that uranium and thorium act differently in the geochemical environment. As a specialist case, the gross ratio of U:Th in the natural mineral monazite varies all over the place, but then most of its thorium typically comes from a separate decay chain to most of the uranium. There are other acounts of chemical U:Th being quite variable in nature, but without further digging I cannot find an example U234:Th230 disequilibrium – and they both come from the same decay chain.

    The event causing a different ratio would be the distance separation over time of U234 from Th230, a partitioning effect. If an amount of U234 was in the sample aeons ago, but then it gradually was removed, while the Th230 stayed in situ, the relation breaks down. (Or the reverse case). My past work did not need to pick up fine distictions of this type so we did not go lookng for them. I hope that there has been some experimental confirmation that U234 and Th230 do not partition in the natural environment over time, to an extent that would affect age dating. Unfortunately, usually you need to know the age of the specimen to determine this, so circular arguments are possible.

    This is hair splitting at fundamental level, but then audits have to go fundamental sometimes.

  13. Steve McIntyre
    Posted Dec 12, 2008 at 9:21 PM | Permalink

    I’ve seen some comments about preferential leaching of U^{234} in diagenesis, which I’ll discuss. The formulas shown so far do not allow for non-U238 originated Th230. The adjustment formula used in Cobb et al 2003 is:
    (\frac{Th^{230}_{rad}}{U^{238}})_{act}=  (\frac{Th^{230}_{meas}}{U^{238}})_{act}- (\frac{Th^{232}_{rad}}{U^{238}})_{act}(\frac{Th^{230}_{nr}}{Th^{232}})_{act} exp(-\lambda_{230}T)

    I was able to replicate the “minimum” corrected results in Cobb et al 2003 within a year or two as follows:

    act.Th230nr_Th232=natural[“Th230/Th232″]* lambda[“Th230″]/lambda[“Th232″];
    # 0.8196591 #this is a constant
    A$act.Th232_U238= (A$Th232/A$U238)*lambda[“Th232″]/lambda[“U238″];
    # [1] 1.162356e-04 6.804403e-05 1.871566e-05 9.592091e-05 7.704772e-05 …

    f=function(t) { #modified from prior function to implement Cobb equation 3
    act= act.Th230_U238-act.Th232_U238*act.Th230nr_Th232*exp(-lambda[“Th230″]*t )
    A=(1- act)
    B= exp(-lambda[“Th230″]*t)
    C= del*10^-3* (lambda[“Th230″]/(lambda[“Th230″]-lambda[“U234″]))
    D= (1- exp( (lambda[“U234″]-lambda[“Th230″])*t) )
    f= (A- (B- C*D) )^2
    f}

    ## calculate for Cobb data
    K=nrow(A); age= rep(NA,K)
    for(i in 1:K) {
    del=A$del_U234_0[i];
    act.Th230_U238=A$act.Th230_U238[i]
    act.Th232_U238=A$act.Th232_U238[i]
    age[i]= optimize(f,interval=c(0,200000),tol=1E-6)$minimum}
    names(age)=A$Sample

    ##report
    test=data.frame(names(age),round(age,0),2000-A$Date_mincorr)
    test=cbind(test,test[,2]-test[,3])
    dimnames(test)[[2]]=c(“Sample”,”Replication”,”Cobb”,”Diff”)
    test

    So far, so good. With only a few lines of code, we’re getting accurate replication of the mechanics and hopefully can test some permutations and combinations.

    • Geoff Sherrington
      Posted Dec 13, 2008 at 3:32 AM | Permalink

      Re: Steve McIntyre (#13),

      Bother, seems I can’t write short blogs.

      Physics and mathematics of age dating.

      All uranium isotopes are radiogenic. All thorium isotopes are radiogenic. That is, they are unstable and they break down to lighter particles, releasing radiation of some type. There are three main natural decay series, now headed by U238 (half life 4.5 Gy), U235 (half life 705 My) and Th232 (half life 14 Gy).

      One can easily tell which series is mentioned, because as a rule of thumb, all daughters in the U238 series have atomic weights which, divided by 4, leave 2 as a remainder. The members of the U225 series are always odd numbers and the Th232 series are all neatly divisible by 4.

      So, U238 will not decay into (say) Th232 because there is a difference of 6 in the atomic weight, when 4,8,12 etc are the only allowed increments within a series (because of alpha decay, atomic weight 4).

      Thus in #13:

      The formulas shown so far do not allow for non-U238 originated Th230.

      There is no natural Th230 that does not originate from U238.

      It is conventional that age dating is performed within one of the 3 decay chains, not between them.

      The U238 chain.

      For accuracy, the chosen isotopes need to have half lives comparable to the ages seeking determination. Thus, U238 itself is too long for most uses. Its daughter U234 (half life 245,000 years) and in turn its daughter Th230 (half life 75,000 years) form a pair which cover geological ages of interest here.

      Another common method in the U238 chain uses U/Pb isotopes, one variation using stable Pb206, the end point. A disadvantage is that between U and Pb there is radon gas Rn222 (half life 3.8 days) and it can escape a closed system and thus affect all isotopes below it. This is an example of radioisotopic disequilibrium. It is best to avoid gas intermediates.

      The choice of isotopes is also affected by the ease of measurement. In general, with radiometric emission methods, gamma is easier than beta is easier than alpha because of a fall-off in penetrating power. Sometimes one has no choice because of the mode of decay.

      With mass spectrometer methods, the preparation of the specimen has to present 100% of each isotope, that is, there should be no losses in dissolution or ionization or whatever preparation is used. Mass spectrometry also works best when the abundance of the sought isotopes is not wildly different, that is, the half lives are similar. (Long half lives mean abundance, short mean very small concentrations).

      The Th232 chain.

      Apart from Th232 itself and the final stable Pb208, the half life of the longest daughter is 6 years for Ra228. Thus,Th232/Pb208 dates are possible. The Th232 stays about constant but the Pb208 builds up quickly with time. Th232 methods are not preferred – radon is another intermediate.

      Geochemistry.

      “Dirty” samples.

      Specifically for speleothems, production is assumed to be from precipitation from dissolved solids in dripwater. Now, if the water contained suspended matter as well as dissolved, the composition of the suspended matter could affect the results.

      See, for example, Mebus Geyh:

      hhttp://www.geochronometria.pl/pdf/geo_20/geo_20_02.pdfhttp://www.geochronometria.pl/pdf/geo_20/geo_20_02.pdf

      Two prerequisites for using the U234/Th230 method discussed here are that –

      1. Terrestrial material (e.g. speleothem) contained solely uranium immediately after its formation and was free of (.detrital.) thorium.
      2. The dated material had behaved as a closed system for uranium through the time of aging.

      The geochemistry of thorium, as a generalization, favours detrital species more than the geochemistry of uranium. So, if a bulk chemical analysis was performed on a dirty speleothem sample, it would probably indicate more total thorium (rather than uranium) than was useful for the radiometric dating. Therefore, empirical/observational/guesswork equations such as those of Cobb et al 2003 in #13 are used to correct for the dirt. But these equations fail theoretically because no geochemical law says that the Th232/Th230 ratio should be constant and be able to be used as a correction.

      There is thus an uncertainty in U234/Th230 dating introduced through the use of Th232, which comes from another decay chain and whose initial abundance might have been rather different from the U238 heading the other decay chain. This in turn urges prudence before detailed analysis is taken too far.

      • Jud Partin
        Posted Dec 13, 2008 at 8:27 AM | Permalink

        Re: Geoff Sherrington (#18),

        But these equations fail theoretically because no geochemical law says that the Th232/Th230 ratio should be constant and be able to be used as a correction.

        then measure more 232Th/230Th ratios in a sample. And between samples. See Partin et al. 2007 (SOM) for some examples of this “failure”.

        • Geoff Sherrington
          Posted Dec 13, 2008 at 7:44 PM | Permalink

          Re: Jud Partin (#19),

          Analogy. You have a Pb/Zn mine with 1% Pb and 10% Zn on average. At another hypothetical mine, you have 5% Pb. Can you calculate the Zn content? No. No known law or generalisation says what the ratio should be. The problem cannot be overcome by running more analyses.

          This isotope correction equation problem is that Th232 and Th230 come from different decay chains, almost certainly have different (unknown) histories from place to place and cannot be related except by empiricism. There is no absolute argument in the correction method.

          The paper by Frank et al 2000 had much discussion about errors. This is understandable. The spread of ratios of U234/U238 in table 1, subscript m, is from 1.903 to 2.076. Intuitively, one would expect no measurable variation in a closed system over a short term wrt the half lives. Alternatively, one would expect the authors to give an explanation for this wide divergence. Therefore, one has to question the estimates of error for each determination. In the same table, the range of values of Th230/U234 is from 0.630 to 0.717. These are some of the most critical values for calculating age. Their range is about the same as U234/U238 (which should be near-constant). The ranges of both are so wide that doubt is cast on the measurement accuracy, especially of Th, which is 750 times less abundant than U in the spring waters. The authors note that in these waters “the 230Th/234U activity ratio is rather small, amounting to 0.0004 =/- 0.00012.” This is another of the crucial measurements for age dating.

          In laboratories, it is not uncommon to confuse error estimation on synthetic specimens with measurement on complex actual samples. (e.g. this paper has a note on background correction of mass spectrometry for scattered atoms, which in some spectroscopy like XRF requires non-linear multiple regression analysis, but let’s not go down that road again).

          To their credit, the authors of Frank et al have taken many precautions but when one deals in nanograms-picograms per gram, errors are a major complication. And that’s before approaching the problems of closed systems and contamination of the speleothem.

          Re: Scott Gibson (#17),

          Well said.

  14. Steve McIntyre
    Posted Dec 12, 2008 at 10:15 PM | Permalink

    OK, just in case Dr Fleitmann is ever in need of a quick algorithm to calculate ages from the data provided in his table, here’s how (to a very close approximation):

    Set the constants – here I’ve used the constants from Cobb et al 2003 as Fleitmann 2007 did not mention their constants (that I could notice anyway):

    lambda=NULL
    lambda[“Th230″]= -log(.5)/75690 # 9.157711e-06 #
    lambda[“U234″]= -log(.5)/235240 # 2.946553e-06 #
    lambda[“U238″]= 1.55125E-10
    lambda[“Th232″]= -log(.5)/ 1.41E10 # 4.915937e-11 ;
    lambda[“U234:U238″]=lambda[3]/lambda[2]

    Input the information from Table 1 (collated from the pdf by pasting to Notepad) and tidying manually:

    A=read.csv(“http://data.climateaudit.org/data/speleothem/fleitmann.2007.csv”,sep=”\t”)
    #Th, U in ppb

    Calculate the information needed for Edwards Equation (1). These can be obtained from the information in Fleitmann Table 1 format by simple operations:

    K=nrow(A); age=rep(NA,K)
    A$act.Th230_U238= A$act.Th230_U234* A$act.U234_U238;A$act.Th230_U238
    # [1] 0.08534431 0.08867938 0.10596096 0.11616480 0.11570169 …
    A$U234_U238= A$act.U234_U238 * lambda[“U238″]/lambda[“U234″];A$U234_U238
    # [1] 7.896412e-05 7.966958e-05 8.490788e-05 8.493947e-05 ….
    A$del_U234_0= (A$U234_U238/lambda[“U234:U238″] -1)*10^3

    Use the same nonlinear function as for Edwards and Cobb (w/o thorium adjustment):

    f=function(t) {
    A=(1- act)
    B= exp(-lambda[“Th230″]*t)
    C= del*10^-3* (lambda[“Th230″]/(lambda[“Th230″]-lambda[“U234″]))
    D= (1- exp( (lambda[“U234″]-lambda[“Th230″])*t) )
    f= (A- (B- C*D) )^2
    f}

    K=nrow(A); age= rep(NA,K)
    for(i in 1:K) { act=A$act.Th230_U238[i];del=A$del_U234_0[i];
    age[i]= optimize(f,interval=c(0,200000),tol=1E-6)$minimum
    }
    names(age)=paste(A$id,A$depth,sep=”-“)

    Report:

    test=data.frame(names(age),round(age,0),A$age*1000,round(age- A$age*1000,0) )
    dimnames(test)[[2]]=c(“Sample”,”Replication”,”Fleitmann”,”Diff”)
    test

    This yields the following, which is pretty close for most years:

    peleo-depth
    Replication
    Fleitmann
    Diff

    H5-5
    6377
    6414
    -37

    H5-25
    6572
    6582
    -10

    H5-65
    7389
    7391
    -2

    H5-98
    8121
    8119
    2

    H5-135
    8157
    8153
    4

    H5-160
    8087
    8083
    4

    H5-175
    8789
    8788
    1

    H5-205
    8227
    8224
    3

    H5-235
    8379
    8381
    -2

    H5-255
    9165
    9164
    1

    H5-295
    9398
    9405
    -7

    H5-325
    9603
    9605
    -2

    H5-350
    9869
    9867
    2

    H5-400
    9922
    9920
    2

    H5-435
    10144
    10147
    -3

    H5-457
    10086
    10091
    -5

    H5-480
    10003
    10011
    -8

    H5-498
    10253
    10250
    3

    H12-85
    230
    232
    -2

    H12-151
    449
    447
    2

    H12-221
    635
    633
    2

    H12-338
    1296
    1298
    -2

    H12-372
    1505
    1505
    0

    H12-391
    1649
    1647
    2

    H12-400
    1638
    1633
    5

    H1-450
    1793
    1791
    2

    H12-571
    2594
    2597
    -3

    H1-585
    2605
    2598
    7

    H12-595
    5337
    5333
    4

    H12-638
    5887
    5887
    0

    H12-660
    5864
    5862
    2

    H1-665
    5876
    5850
    26

    H12-744
    5887
    5887
    0

    H12-760
    6152
    6157
    -5

    H1-780
    6370
    6341
    29

    Q5-2
    405
    403
    2

    Q5-40
    1496
    1400
    96

    Q5-62
    3046
    3049
    -3

    Q5-140
    3791
    3789
    2

    Q5-197
    3847
    3845
    2

    Q5-259
    4437
    4441
    -4

    Q5-340
    4438
    4435
    3

    Q5-428
    4974
    4969
    5

    Q5-503
    5836
    5838
    -2

    Q5-565
    6473
    6470
    3

    Q5-574
    6635
    6634
    1

    Q5-627
    6948
    6960
    -12

    Q5-702
    8028
    8030
    -2

    Q5-762
    8879
    8880
    -1

    Q5-780
    8995
    8996
    -1

    Q5-860
    9904
    9903
    1

    Q5-903
    10207
    10210
    -3

    Q5-961
    10629
    10623
    6

    S3-149
    416
    418
    -2

    S3-245
    745
    746
    -1

    S4-10
    9166
    9163
    3

    S4-75
    9755
    9674
    81

    S4-89
    9277
    9283
    -6

    S4-120
    9657
    9661
    -4

    S4-121
    9659
    9660
    -1

    S4-135
    9706
    9670
    36

    S4-168
    9802
    9760
    42

    S4-205
    9943
    9908
    35

    S4-220
    10089
    9964
    125

    S4-230
    10075
    10077
    -2

    S4-395
    11183
    10864
    319

    S4-415
    10959
    10742
    217

    S4-420
    10823
    10681
    142

    S4-421
    10762
    10761
    1

    D1-312
    1483
    1480
    3

    D1-432
    1882
    1890
    -8

    D1-565
    2182
    2170
    12

    D1-665
    2661
    2670
    -9

    D1-810
    3142
    3140
    2

    D1-875
    3659
    3490
    169

    D1-960
    3716
    3720
    -4

    D1-1040
    4292
    4290
    2

    D1-1115
    4529
    4530
    -1

  15. Steve McIntyre
    Posted Dec 12, 2008 at 11:02 PM | Permalink

    I can’t tell what Fleitmann et al 2007 did on their Th adjustment. Here’s what they say (and, in my opinion, this does not count as an adequate operational description):

    In some cases, detritus corrections were needed when 230Th/232Th ratios were high (see Table 1). If necessary, ages were corrected assuming a 232Th/234U ratio of 3.8 (see Frank et al., 2000 for details). Detritus corrections were needed for stalagmites S3 and S4 from Kahf Defore (Southern Oman) and stalagmites H5 and H12 from Hoti Cave (Northern Oman). We must, however, emphasize that their chronologies are either confirmed by annual layer counts (S3 and S4, Figs. 3d, e) or by Th–U ages (H5 and H12) for which no or minor detritus correction was needed (see Table 1 for details).

    As an experiment, I used the adjustment for non-radiogenic Thorium of Cobb et al 2003 and ended up with what looks like a 75 year displacement to the present of Q5 – pretty much what I thought might be the case from the correlation analysis. Here’s what I did (And I’ll probably post this up separately since it’s sort of interesting).

    Calculate various activities to go from Fleitmann 2007 Table 1 format to the format of Cobb et al 2003 Equation (3):

    act.Th230nr_Th232=natural[“Th230/Th232″]* lambda[“Th230″]/lambda[“Th232″]; act.Th230nr_Th232
    # 0.8196591 #this is a constant
    A$act.Th232_Th230=1/A$act.Th230_Th232 ;A$act.Th232_Th230
    ##act.Th232_Th230= (Th232/Th230)*lambda[“Th232″]/lambda[“U238″]=1/act.Th230_Th232
    # [1] 0.1941747573 0.0224567707 0.0049637645 0.0073168947 ….
    A$act.Th232_U238= A$act.Th232_Th230 * A$act.Th230_U238; A$act.Th232_U238
    # 1.657171e-02 1.991453e-03 5.259653e-04 ….

    Set up the nonlinear solver as with Cobb 2003 case:

    f=function(t) {
    act= act.Th230_U238-act.Th232_U238*act.Th230nr_Th232*exp(-lambda[“Th230″]*t )
    A=(1- act)
    B= exp(-lambda[“Th230″]*t)
    C= del*10^-3* (lambda[“Th230″]/(lambda[“Th230″]-lambda[“U234″]))
    D= (1- exp( (lambda[“U234″]-lambda[“Th230″])*t) )
    f= (A- (B- C*D) )^2
    f}

    Calculate for the entries in Fleitmann Table 1 and report:

    K=nrow(A); age= rep(NA,K)
    for(i in 1:K) {
    del=A$del_U234_0[i];
    act.Th230_U238=A$act.Th230_U238[i]
    act.Th232_U238=A$act.Th232_U238[i]
    age[i]= optimize(f,interval=c(0,200000),tol=1E-6)$minimum}
    names(age)=A$Sample
    test=data.frame(paste(A$id,A$depth,sep=”-“),round(age,0),1000*A$age)
    test=cbind(test,test[,2]-test[,3])
    dimnames(test)[[2]]=c(“Sample”,”Replication”,”Fleitmann”,”Diff”)
    test

    As an indication of the average non-radiogenic Th adjustment for Q5 and D1 in the overlap period 2700-4700 BP (and this is an indication only), calculate the mean adjustment for each core:

    temp1=(A$id==”Q5″); range(test$Replication[temp1])
    temp2=(A$id==”D1″); range(test$Replication[temp2])
    temp=(test$Replication>=2700)&(test$Replication< =4700)
    mean(test$Diff[temp&temp1]) # -89.8
    mean(test$Diff[temp&temp2]) # -16.4

    Q5 is brought forward about 75 years more than D1. This is the first time that I’ve done this sort of calculation so these results have to be taken with caution.

    But the result is definitely intriguing as this adjustment for non-radiogenic Th ends up reversing the previously reported “anti-correlation” between D1 and Q5 in the overlap period – and is convincing evidence to me that relative dating errors may be highly autocorrelated (as I surmised before).

    This also indicates a refutation for those critics who argued against even experimenting with displacing one core against another.

    • Jud Partin
      Posted Dec 13, 2008 at 1:06 AM | Permalink

      Re: Steve McIntyre (#15), Fleitmann’s description sounds fine to me. The 230Th/232Th correction yields younger ages as some of the 230Th is attributed to detrital material and not to U decay. The Isoplot software documentation (and references therein specifically Ludwig and Titterington 1994) should help you out. While Fleitmann assumes a ratio based on literature, the 232Th/230Th ratio for a sample can be measured using isochrons. But this is an expensive and difficult measurement.

      Also, displacing one core relative to another is completely justifiable b/c of age errors. What I didn’t agree with was a wholesale moving of the whole record by a fixed amount. The section between each date should be viewed as independent.

  16. Scott Gibson
    Posted Dec 13, 2008 at 2:48 AM | Permalink

    Many years ago, I had some experience with radiometric dating procedures. As I remember, the assumptions made for uranium-series dating of speleothems were as follows:

    1. The speleothem is usually calcite crystallized from precipitation from groundwater entering an air-filled cave.

    2. Uranium is soluble in the near-surface groundwater, whereas thorium is not, so initial thorium concentration is assumed to be negligible. As a result, all thorium in the sample is assumed to be a decay product of uranium. It is assumed that once precipitated, uranium and thorium are not removed or added to through chemical processes; only radioactive decay changes the composition.

    3. There are several relevant decay series going on, all at the same time, complicating the issue. Simplified, the most important decay series can be stated to be: 238U -> 234Th -> 234Pa -> 234U -> 230Th -> 226Ra etc. To date the sample, the intitial 230Th/234U ratio must be known. Because the thorium concentration during calcite deposition is assumed to be negligible, this ratio is assumed to be 0. Also, the ratio of 238U/234U must be known, because some 234U (and hence 230Th) comes from the decay of 238U. In my opinion, this assumption is the most difficult part of uranium-series dating, because the ratio of natural 238U/234U changes with time.

    Because of all of the assumptions, I have always been somewhat skeptical of uranium-series dates. A glance at the literature suggests that some consistent dates have been obtained, however we know how data can appear consistent without actually being correct. I would like to see other opinions on this.

  17. Steve McIntyre
    Posted Dec 13, 2008 at 11:08 AM | Permalink

    In some cases, detritus corrections were needed when 230Th/232Th ratios were high (see Table 1). If necessary, ages were corrected assuming a 232Th/234U ratio of 3.8 (see Frank et al., 2000 for details).

    I’ve posted up a copy of Frank et al 2000 here. I don’t see anything in Frank et al 2000 involving a 232Th/234U ratio of 3.8. It describes adjustments using the assumption that modern spring water was relevantly similar (wrt Th) to ancient spring water.

    Also what is the criterion involved in “in some cases”? This is not an operational criterion. In what cases? What’s the hurdle?

  18. Posted Jan 6, 2009 at 2:56 AM | Permalink

    I recently came accross your blog and have been reading along. I thought I would leave my first comment. I dont know what to say except that I have enjoyed reading. Nice blog. I will keep visiting this blog very often.

    Kate

    http://educationonline-101.com

  19. Posted Apr 26, 2009 at 1:05 PM | Permalink

    Can you recommend some further reading for the beginner as this is interesting stuff

One Trackback

  1. […] Life is easy, math is hard.  I’m not even sure what a speliothem is, nevermind that they’re dating already. […]

Follow

Get every new post delivered to your Inbox.

Join 3,382 other followers

%d bloggers like this: