## 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[,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

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

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 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 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. 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