Climate Audit

In recent discussions of Steig’s Antarctic reconstruction, one of the interesting statistical issues is how many principal components to retain. As so often with Team studies, Steig provided no principled reasoning for his selection of 3 PCs, statements about their supposed physical interpretation were untrue and, from some perspectives, the choice of 3 seems opportunistic.

As we’ve discussed previously, the issue of PC retention was a battleground MBH issue and was never resolved in a principled way. After the fact, Mann adduced a retention strategy that was not described in the text and could not be evidenced in other networks. Wahl and Ammann invented an ad hoc argument that the MBH result “converged” (as long as the bristlecones were in); this “convergence” argument is not one that appears in any literature by known statisticians. Although Wegman described this method as having “no statistical integrity”, that did not prevent its sanctification by IPCC AR4 through their benediction of Wahl and Ammann 2007. I won’t review this long and sordid history as we’ve discussed it on many occasions in the past.

Reading some of the interesting posts by the two Jeffs and Ryan O at Jeff Id’s blog on truncated total least squares, I recalled that borehole inversion used truncated SVD. PhilB had done some work on emulating borehole inversion calculations a couple of years ago. I corresponded with him about this, but hadn’t figured out how the calculations actually worked, but recalled that results were very sensitive to retained PCs.

I thought that it would be interesting to re-visit this topic, in part to get another perspective on retention issues in a completely different context and in part to try to get a better handle on borehole reconstructions which we’ve barely touched on over the past few years.

Unfortunately it took me a while to figure out how these calculations. While I think that I can replicate some examples, I’m not 100% sure that I’ve got it right. As usual, I’ve placed materials online, in case readers can improve on what I’ve set out here.

There is an excellent data set of borehole readings at WDCP-Paleo indexed here. Data listings can be accessed from links at the bottom of the page, which in turn links to the data provided in a somewhat annoying html format, though the data can be scraped from the html pages. I’ve placed online a read-function keyed by borehole id, which will be illustrated below. This is set up for North American listings, but can be easily modified for other locations (for other locations, you’ll need to check that calls are lowercase).

Many of the borehole locations are “holes of opportunity” from hard rock mineral exploration in mining camps in northern Ontario and Quebec; many sites are familiar names to me. I might even have visited some of the sites at some point in my life.

Tamino, who stoutly defended NASA’s inability to locate post-1991 data from stations in Peru and Bolivia, had a post on boreholes a while ago, in which his type example was borehole CA-001-0 in Ottawa. This proved an odd choice as a type case as boreholes theoretically are supposed to be away from cities. For the benefit of Tamino and readers unfamiliar with Canadian geography, Ottawa is the capital of Canada. Being from Toronto, I concede urban status to Ottawa somewhat reluctantly, but must do so. Like other urban centers in North America, it’s not an obvious location for mineral exploration and Tamino surprised me a little here. However, Tamino’s readers need not fear that the Canadian federal government is going to be displaced by a mine, as I presume that either this borehole was put down for some purpose other than mineral exploration, unless the location in the database is wrong.

Borehole inversion methods have a long history. The researcher hopes to reconstruct the surface temperature history from a set of depth-temperature measurements using the heat equation. The problem turns out to be very “ill-conditioned” (the literature derives from numerical analysis, rather than statistics though the issues are similar). Therefore, various regularization methods have been adopted over the years, of which the most common appears to be truncated SVD.

Hugo Beltami’s website has many papers on the topic, most of which seem to describe the same methodology in more or less similar terms. Beltrami provides equations and examples, but it still took a while to replicate his methodology and I’m not 100% sure that I’ve got it right, though I’ve got much of it. As a benchmark case, I eventually settled on borehole ca-8813 at Belleterre, Quebec illustrated in Figure 6 of Beltrami et al 1992

The idea behind borehole inversion (expressed here in matrix terms rather than the continuous expressions of the papers):

where
is the downhole temperature,
E is a matrix derived from the heat equation (see Beltrami’s articles for the assertion) and
is the temperature history, expressed typically as a vector (step function).

Where erfc is the standard function,
is heat diffusivity here converted to years,
t is a sequence of time intervals.

The idea is to estimate the history from by matrix inversions, the mechanics of which are perhaps best seen in the script below.

First download information from the WDCO site for ca-8813 using online tools.

source(“http://www.climateaudit.info/scripts/borehole/functions.borehole.txt”)
belle=read_borehole(id=”ca-8813″)
Data=belle$bore
z=Data$z ;range(z) # 25.98 389.17
Q=T_z=Data$T_z

In this paper, Beltrami used 20-year steps and the period 1000-2000. There has been discussion in the literature over resolution and, in later papers, Beltrami typically uses only 100-years and a 500 year interval. However, for today’s purposes, which are about truncation principles as much as anything, the denser net is advantageous. Thus:

delta=20
t=seq(0,1000,delta);K=length(t); #51 #time per page 175 Beltrami

Next we make the matrix from the heat conduction equation. This could be done in one line using the outer function, but I’ve sacrificed a little elegance for transparency:

M=length(z);K=length(t)
E=array(1,dim=c(M,K-1) )
for( k in 1:(K-1)){
for(j in 1:M) {
E[j,k]= erfc(z[j]/(2*sqrt(Kdiff*t[k]))) – erfc(z[j]/(2*sqrt(Kdiff*t[k+1])))
}}
dim(E) #[1] 44 51
dimnames(E)[[2]]=t[2:K]-delta/2;dimnames(E)[[1]]=z
E[1:3,1:5]
# 10 30 50 70 90
#25.98 -0.4646354 -0.1404778 -0.06778118 -0.04175596 -0.02900359
#34.64 -0.3295700 -0.1609920 -0.08293746 -0.05241155 -0.03690798
#43.3 -0.2229484 -0.1658662 -0.09284768 -0.06062173 -0.04344844

The second term in Beltrami’s equation includes a term for surface heat flow density. I was unable to locate any information on what value was used and the value appears to get washed out in the coefficient. However, I located an interesting map of surface heat flow density http://www.heatflow.und.edu/agus05.ppt and used a value of wm-2 from the map coloring as an approximation.

q_0=0.04 #deem inspection of map of surface heat flow density
R=q_0*z/Kcond

I did a couple of quick experiments with regression and ridge regression, both of which yielded warnings due to singularity:

A=data.frame(Q,R,E); dim(A) #[1] 44 52
fm=lm(Q~.,data=A);summary(fm) # #OLS fails
# prediction from a rank-deficient fit may be misleading
library(MASS)
fm.ridge=lm.ridge(Q~.,data=A);summary(fm)
#prediction from a rank-deficient fit may be misleading

The regression fails because of problems in inverting the matrix. Beltrami regularizes this by retaining only a few PCs; he reported the retention of all PCs with eigenvalues above 0.025 (page 173), which, for this example, proves to be, 6

MinEigen=.025 #page 173 right column
one=rep(1,M)
A=cbind(one,R,E) #this is matrix desribed on page 173 left column
svd.bore=svd(A)
U=svd.bore$u;d=svd.bore$d;V=svd.bore$v
round(d,5)
#[1] 30.18039 3.04874 0.54588 0.21885 0.07824 0.02507 0.00692
r0=sum(d>MinEigen);r0 #6

According to my original objective, I did some experiments with different numbers of retained PCs.

beta=list();Qhat=list()
for(r in 1:8) {
if(r==1) Sr=as.matrix(1/d[1:r]) else Sr=diag(1/d[1:r])
beta[[r]]=as.matrix(V[,1:r])%*%Sr %*%t(U[,1:r]) %*% Q;
Qhat[[r]]=c(A%*%beta[[r]]); # 5.779109 5.665154 5.569599 5.492003 5.429110
}

Here is the target reconstruction excerpted form Beltrami et al 1992 Figure 6.

Here is my emulation of Beltrami’s methodology On the left, I’ve plotted the original data (black points) with the fits from the first few recons; they all fit more or less perfectly with relatively few PCs. On the right, I’ve plotted the various reconstructions. The recon with 6 PCs looks pretty close to the Beltrami example, but I don’t guarantee it.

A few comments. As more PCs are added, the results start diverging wildly. I’ve shown up to 8 here, but they get more and more extreme.

Keep in mind that the underlying temperature measurements are smooth. It appears to me that the notches on the right hand side of the reconstruction are artifacts of a 1-dimensional Chladni-type situation, sort of like what we saw in a totally different context with Stahle SWM tree ring network eigenvectors. There’s a reason for this: the coefficients are obtained from the eigenvectors.

Keep in mind here the important and interesting sign change property of eigenvectors in spatially autocorrelated regions that we’ve discussed in the past. The first eigenvector has no sign change; the 2nd has one sign change and so on. I can’t see any reason for using 6 eigenvectors in this case.

At this time, I draw no morale on what any of this means. I simply was interested in seeing how the borehole recons worked and in seeing what principles, if any, were used in SVD truncation. It’s a fairly large topic, but not unconnected with issues that we’ve discussing and may perhaps give an alternative perspective on some issues.

I’m very far from a clear understanding of how much weight can be placed on these borehole reconstructions, but we’re in a much better position to experiment with these matters. If nothing else, the tools may prove helpful to some readers.

References:
Beltrami, H., A. M. Jessop, and J. C. Mareschal. 1992. Ground temperature histories in eastern and central Canada from geothermal measurements: Evidence of climate change. Global Planet. Change 98, no. 2-4: 167–184.

Click to access Beltrami_et_al_1992_GPC.pdf

Beltrami website: http://esrc.stfx.ca/pdf/