Inverting Boreholes « Open Mind

Inverting Boreholes

June 7, 2007 · 8 Comments

I’ve been studying borehole temperature inversions lately, so I thought I’d post some ideas on the topic. The idea is to measure the temperature at depth below the surface of the ground; this is a borehole temperature profile. From this information, we attempt to reconstruct the temperature history at the surface; this is the ground surface temperature history (GSTH). The process of translating the temperature profile to the surface temperature history is referred to as inversion of the borehole temperature profile.

My exposition may not be clear enough, so I encourage questions, and may end up editing this post to make things clearer. This is, however, another of those mathematical posts, so if that sort of thing doesn’t frighten you, read on…

In an earlier post, I gave the 1-dimensional heat equation:

${\partial u \over \partial t} = \lambda {\partial^2 u \over \partial z^2}$ ,

where $u(t,z)$ is the temperature at time $t$ and depth $z$ , and $\lambda$ is a constant called the thermal diffusivity of the material we’re studying (in this case, the rock beneath our feet). The function $u(t,z)$ gives the temperature at all times, at all depths. This particular form of the equation assumes a few things, namely that the thermal properties of the rock (in particular, its diffusivity $\lambda$ ) is constant throughout the rock, and that there is no significant heat generation within the rock (which can happen from the decay of radioactive material). If these conditions do not hold, we modify the heat equation, but for the purposes of discussion I’ll stick to the simple form.

I’d like to show how we can solve the heat equation, to invert a borehole temperature profile. We will take great advantage of the fact that the heat equation is a linear equation. This means that if we find two different solutions $u_1(t,z)$ and $u_2(t,z)$ to the heat equation, then the function which is their linear combination

$au_1(t,z) + bu_2(t,z)$ ,

where $a$ and $b$ are constants, is also a solution to the heat equation. Thus we can build up more complex solutions from linear combinations of simpler solutions. So, we’ll seek simple solutions, then combine them linearly to make more intricate ones.

First let’s see what happens if the temperature never changes. In that case

${\partial u \over \partial t} =0$ .

Then, from the heat equation, we have

${\partial^2 u \over \partial z^2} = 0$ .

This has the unique set of solutions

$u = a + bz$ ,

where $a$ and $b$ are constants. If we were to plot the temperature in a borehole as a function of depth, for this case we would get a straight line:

At the surface, the temperature is 5.2 deg.C, 500 m below the surface, it’s 12.3 deg.C, hence it’s getting hotter the deeper we go into the earth. This is actually the case for our planet; the heat from deep inside the earth is flowing outward toward the surface, creating a nearly linear depth-temperature profile in most boreholes.

But a perfectly linear temperature-depth profile only happens when the surface temperature is constant. If $u(t,z)$ is the actual temperature function, it is a solution of the heat equation. Also, the function $a+bz$ is a solution. Hence (because the heat equation is linear) their difference

$\theta(t,z) = u(t,z) - a - bz$ ,

is also a solution of the heat equation, i.e.,

${\partial \theta \over \partial t} = \lambda {\partial^2 \theta \over \partial z^2}$ .

Now let’s seek solutions for this equation which are separable, i.e., which have the form

$\theta(t,z) = T(t) Z(z)$ .

Essentially, we’re seeking a function $\theta(t,z)$ which is the product of a function $T(t)$ of $t$ alone with another function $Z(z)$ which is a function of $z$ alone. The heat equation becomes

${dT \over dt} Z = \lambda T {d^2 Z \over dz^2}$ ,

or simply

$T' Z = \lambda T Z''$ .

Dividing both sides by $TZ$ , we get

${T' \over T} = \lambda {Z'' \over Z}$ .

The left hand side is a function of $t$ alone (being the ratio of two functions of $t$ alone), the right-hand side is a function of $z$ alone. The only way this can happen, is for both sides to equal a constant, which I’ll call $\alpha$

${T' \over T} = \lambda {Z'' \over Z} = \alpha$ .

This is actually two equations, one for $T$ and another for $Z$ . We can solve the equation for $T$ directly

$T(t) = e^{\alpha t}$ ,

and for $Z$ as well

$Z(z) = e^{\sqrt{\alpha / \lambda}~z}$ .

Now we have a problem. For real $\alpha$ , the function $e^{\alpha t}$ explodes to infinity quite rapidly as $t$ goes to infinity. If $\alpha > 0$ , then the function explodes as $t \to + \infty$ , while if it’s less than zero, then the function explodes as $t \to - \infty$ . But we want our functions not to do that; the temperature history of the earth should be bounded. We can solve this problem by choosing the constant $\alpha$ to be imaginary, i.e.

$\alpha = i \omega$ ,

where $\omega$ is real, and $i$ is the square root of -1. Then we have

$T(t) = e^{i \omega t}$ ,

with $\omega$ real.

This introduces another problem, that then the function $T(t)$ will be complex, while the actual temperature anywhere in the earth should be a real number. We’ll solve that problem later, by combining solutions linearly in such a way that their imaginary parts cancel, leaving a real-valued solution.

What is the square root of $\alpha / \lambda = i \omega / \lambda$ (for $\omega$ real)? If $\omega$ is positive, there are two solutions,

$\sqrt{i\omega / \lambda} = \pm (1 + i) \sqrt{\omega \over 2\lambda}$ ,

while if $\omega$ is negative, we have the two solutions

$\sqrt{i\omega / \lambda} = \pm (1 - i) \sqrt{-\omega \over 2\lambda}$ .

How do we choose between the plus and minus sign? We do so by noting that with the plus sign, the function $Z(z)$ explodes rapidly as $z \to \infty$ , but with the minus sign it doesn’t. So we’ll choose the minus sign in both cases. We don’t have to worry about exploding to infinity as $z \to - \infty$ , because the earth doesn’t extend to negative depth (that would be above the surface). Hence for any real constant $\omega$ , we have a solution of the form

$\theta(t,z) = T(t) Z(z) = e^{i \omega t} e^{-(1 \pm i)(\sqrt{|\omega| \over 2\lambda})z}$ .

Just for convenience, I’ll write

$k = \sqrt{|\omega| \over 2\lambda}$ ,

and I can write our solution as

$u(t,z) = e^{i\omega t - (1 \pm i) kz}$ ,

where we choose the + sign (in $\pm$ ) when $\omega$ is positive, the – sign when it’s negative.

Now we have enough solutions (by choosing different values of $\omega$ ) to construct the actual physically meaningful solution we want. We could, for instance, combine solutions for $\omega = \omega_1$ and $\omega = \omega_2$ , with constant coefficients $F_1$ and $F_2$ , to get the solution

$\theta = F_1 e^{i\omega_1 t - (1 \pm i) k_1 z} + F_2 e^{i\omega_2 t - (1 \pm i) k_2 z}$ .

But we want more than that; we want to combine solutions for all possible values of the parameter $\omega$ . We can do this with an integral

$\theta(t,z) = \int_{-\infty}^\infty F(\omega) e^{i\omega t - (1 \pm i) kz} ~d\omega$ .

where we now have a coefficient for each $\omega$ value which defines a function $F(\omega)$ . For any choice of this function, we get a valid solution to the heat equation. Note: for each $\omega$ , the parameter $k$ is calculated as $k = \sqrt{|\omega| / (2\lambda)}$ .

What is the value of this function at the surface ( $z=0$ )? Suppose the actual temperature at the surface (the ground surface temperature history) is $f(t)$ . Then, putting $z=0$ in our solution, we get

$f(t) = \theta(t,0) = \int_{-\infty}^\infty F(\omega) e^{i\omega t} ~d\omega$ .

This may be familiar to you as the Fourier series form of the function $f(t)$ . The coefficient function $F(\omega)$ which defines the function $f(t)$ is called the Fourier transform of the function. It is well-known how to compute the Fourier transform, namely,

$F(\omega) = {1 \over 2\pi} \int_{-\infty}^\infty f(\tau) e^{-i\omega \tau} ~d\tau$ .

If we substitute this expression for $F(\omega)$ is the formula for our heat-equation solution, we get

$\theta = {1 \over 2\pi} \int_{-\infty}^\infty \int_{-\infty}^\infty f(\tau) e^{i\omega(t-\tau) - (1 \pm i)kz} ~d\tau ~d\omega$ .

I have expressed this as an integral over $d\omega$ (as well as $d\tau$ ), with a $\pm$ to indicated that the imaginary part of the exponent factor of $z$ is positive when $\omega$ is positive, negative when negative, and that the factor $k=\sqrt{|\omega| / (2\lambda)}$ . I will now change variables, so as to make the expression a bit simpler. The integration variable $\omega$ will become $x$ , and I won’t provide all the details, only the final result:

$\theta = {4 \lambda \over \pi z^2} \int_{-\infty}^\infty \int_{0}^\infty f(\tau) x e^{-x} \cos [x + (2\lambda (\tau - t)/z^2) x^2] ~dx ~d\tau$ .

I can also write this as

$\theta = {4 \lambda \over \pi z^2} \int_{-\infty}^\infty \int_{0}^\infty f(t+\tau) x e^{-x} \cos [x + (2\lambda \tau/z^2) x^2] ~dx ~d\tau$ .

Finally, I will write it as

$\theta = {4 \lambda \over \pi z^2} \int_{-\infty}^\infty f(t+\tau) G(2\lambda \tau / z^2) ~d\tau$ ,

where

$G(q) = \int_{0}^\infty x e^{-x} \cos [x + q x^2] ~dx$ .

The function $G(q)$ is referred to as a Green’s function.

We see that our solution for the borehole temperature profile involves the surface temperature history $f(t)$ for all times (including the future!), because the variable $\tau$ is integrated from $-\infty$ to $+\infty$ . But the temperature profile in a borehole at time $t$ should not depend on the surface temperature in the future, only the surface temperature in the past! If we actually compute the value of the Green’s function $G(q)$ (which I did numerically), we get the (not so) surprising result that for $q>0$ , the value is $G(q) = 0$ ! Hence in fact, the borehole temperature profile does not depend on the surface temperature in the future, only the surface temperature history (in the past).

All this is fine; it enables us to translate the surface temperature history $f(t) = \theta(t,0)$ to the borehole temperature profile $\theta(t,z)$ . But that’s not what we wanted to do! We wanted to transform the temperature profile $\theta$ at some particular time $t$ into the surface temperature history $f(t)$ (for $t$ in the past). This is the inverse problem.

The way it’s done is to assume a form for the temperature history, then use our formula to compute what the borehole temperature profile would be, given that form. A common form assumed is $f(t) = T_1$ during time interval 1, $f(t) = T_2$ during time interval 2, etc., where the time intervals run up to the present day. Before time interval 1, we assume the temperature is $f(t) = T_0$ . Plugging that form of temperature history into our formula, we geta a borehole profile that depends on the parameters $T_0,~T_1,~T_2,~...,~T_n$ (where we use $n$ time intervals, the last of which ends at the present). We also must remember that the actual temperature profile in a borehole includes the linear part, the $a + bz$ , so we have to estimate the constants $a$ and $b$ , which are also parameters of our model.

We then seek the parameter values which give the best match between the computed borehole profile, and the observed one. Then, we can approximate the actual temperature history by our “synthetic” history. Of course, it’s only an approximation. The actual temperature history at the surface does not consist of a series of episodes of constant temperatures $T_0,~T_1,~...$ . That’s why most reconstructions (like the one shown for borehole CA-001-0 in the previous post) are rather coarse. But if we choose our time intervals wisely, our approximation will be at least an indication of the truth; in particular, we’ll get a good idea whether or not the average temperature during each episode was hotter or colder than during the previous, and successive, episodes, and we’ll get an estimate of how much hotter or colder.

There’s another difficulty with this procedure. Solving for the “best-fit” values of our parameters generally boils down to the problem of inverting a large matrix. The trouble is, this matrix is ill-conditioned. That means that small errors in our borehole measurements can translate into large errors in our estimated parameters. Hence the sensitivity of this reconstruction procedure is far from ideal.

All of which is why I’ve been working on alternative methods to reconstruct the surface temperature history from a given borehole temperature profile. That will be the subject of future posts.

Categories: Global Warming · climate change · mathematics

8 responses so far ↓

Anonymous (kind of) // June 8, 2007 at 6:17 am | Reply

Although I didn’t read it carefully, looks like a nice derivation.

Let me guess…. you make it probabilistic and put in a prior for T. A proper prior would be either smooth (composed of smooth basis functions, such as low-frequency wavelets), or smooth but allowing almost discontinuous changes now and then. The latter would be computationally much harder, so you are not probably trying it. And maybe you don’t go the bayesian way after all, but stick with some novel basis that (optimally?) takes into account the loss of accuracy in the far past.

Right? :)
cilibrar // June 18, 2007 at 5:58 am | Reply

sounds like a job for a kernel method to me.

[Response: Indeed it is. In fact, the function which one integrates with $f(t)$ is referred to as the heat kernel. But I wanted to show a derivation.]
George // June 19, 2007 at 5:09 pm | Reply

A statement was made in the previous post (Notes from the Underground”) that

“If one keeps smaller singular values in the psuedoinverse than the residuals get smaller but one no longer gets a hockey stick shape. Isn’t the ill-conditioning of the matrix suggesting that there isn’t a unique answer to be obtained? Aren’t the physics equations telling us something?”

I have a question about this:

How much does the ill-conditioning issue really affect the conclusion that there has been warming over the 20th century relative to the previous few centuries (ie, a hockey stick shape for the temperature time graph)?

First, it is not clear to me why one needs to do a full blown inversion if one is simply interested in picking out a basic temperature trend over the last few centuries — eg warming of the past century relative to the previous 4 centuries.

For a profile like the one you provided in the previous post which is pretty much linear from a depth of about 100m downward and sharply hooked to the right from 100m up to the surface, it would seem that you can just infer the recent warming trend by visual inspection.

Second, I understand that the inversion does not give a unique answer, but it seems to me that an essential physical reason for this is that you can add (almost arbitrarily) various high frequency components with relatively large amplitude (eg, for daily and seasonal temperature variation) to the overall solution without appreciably affecting the outcome.

In other words, there is a large number of different solutions that involve various combinations of such high frequency components that produce essentially the same profile at greater than a few tens of meters depth , since the average contribution of these high frequency components will be zero over the longer term.

What this indicates to me is that one has to be careful about reading too much into the ill-conditioning ( “non-uniqueness”) issue. It may or may not be a significant problem, depending on the question one is trying to answer.

[Response: Indeed one can add high-frequency components, because they only propogate to a very shallow depth. In fact such components are present, due to both the annual temperature cycle of the seasons (high frequency) and daily temperature cycle (*very* high frequency), but these are damped out in the first few meters of the borehole. Most reconstructions don't include data from the top several meters, so the effect of such high-frequency fluctuations is not detectable.

Getting a "rough" reconstruction (like the one given above) is a reliable and accurate procedure, but for getting a detailed reconstruction with high resolution in time, the ill-conditioned nature of the matrix is just too much. The problem with ill-conditioned matrices is that they magnify the errors in the data, and when a matrix is sufficiently ill-conditioned, the errors get magnified until they're bigger than the signal.]
George // June 20, 2007 at 3:15 pm | Reply

“The problem with ill-conditioned matrices is that they magnify the errors in the data”

Isn’t it possible — at least in some cases — to tell the difference between “errors” in the data (from noise in a temperature sensor) and the actual temperature signal?

The extreme case would be the “spike” on the profile at a certain depth. This is clearly non- physical, since heat flow alone would not produce such a concentration. Only a local source at depth would do so.

That is just the extreme example, but it seems that one can infer a lot from the way the profile (including its derivatives — eg, second derivative [curvature], which actually appears in the heat equation) change from one depth to another.

The profile at a given depth in a borehole has to be consistent with the profile at other depths (particularly the surrounding depths).

Unless “errors” (eg, instrumental noise) just happened to induce a change in the profile at a given depth that was consistent with what was going on at the surrounding depths (ie, consistent with surrounding profile and its derivatives), one should be able to distinguish it from the actual signal.

If there is a “jump” in either the profile or its derivatives, it would be an indication that something unusual was going on (and most likely not in the temperature history).

The obvious example is a break in the profile, but the existence of a large positive curvature in the profile at one depth with a large negative curvature at a nearby depth (producing a “kink” in the profile) would also most likely indicate something weird was going on — and that the weirdness was not in the temperature history itself.

I’d also have to say that by looking at the curvature of the profile alone, one would be able to glean an awful lot of information about the temperature history, since the curvature at each depth tells you how fast the temperature is changing at that depth and whether it is increasing or decreasing. This information is useful for making decisions about excluding noise), but my guess is that it also may be useful for reconstructing the surface temperature history.
george // July 6, 2007 at 6:54 pm | Reply

By the way, I am still curious about your alternate (non-inversion) method for reconstructing the surface temperature from a bore-hole profile.

I realize you have been posting on a lot of other topics, but I found this one particularly interesting.
Heiko Gerhauser // August 21, 2007 at 6:45 pm | Reply

Hello Tamino,

I wonder whether you might know how the borehole data are used to obtain world averages for temperature?

http://www.geo.lsa.umich.edu/climate/core.html

http://www.globalwarmingart.com/wiki/Image:Global_Warming_Map_jpg

What you’ll notice in the above graphic is that the coverage is quite spotty. I checked the respective IPCC AR4 chapter, and they basically give the same map (except that a few more proxies are added in, which makes it harder to read).

In particular, the oceans are not covered at all (for obvious reasons), and neither are vast tracts of land, while some small areas have huge numbers of boreholes.

The “good” side of there being vast areas with little coverage is that it’s easy to check them out, for example north of 60 degrees there are only 6 bore holes in Canada. There’s only one in North Africa and a few in South America and Eastern Siberia.

I looked at them, and saw vastly different graphs, eg for Canada north of 60 degrees from a recent 2.5C fall between 1800 and 1950 to a 3C rise between 1500 and 1800.

I don’t understand how those kinds of data can be meaningfully averaged out.

Am I failing to see something obvious (misreading the data?)?

http://www.geo.lsa.umich.edu/climate/RECONSTRUCTION/CA-066-0.html

http://www.geo.lsa.umich.edu/climate/RECONSTRUCTION/CA-290-1.html

Are there more boreholes than indicated on the map? Are these uncorrected data that need some sort of adjustment?

[Response: I'm not expert in this, but I'll give you my opinions.

As far as I know(!), borehole data are averaged using area weighting, the same way thermometer readings are. It wasn't always this way; area weighting has brought borehole results into closer agreement with other proxy reconstructions. But given the uneven area coverage, and the fact that ground surface temperature is not surface air temperature, there are definitely issues involved which make it "iffy" to rely on boreholes alone for proxy reconstruction. Nonetheless, they don't require calibration and of respond to temperature alone (tree rings, e.g., indicate temperature, water availability, and other factors), so they have advantages as well as disadvantages.

The sometimes drastic differences between results on the same continent indicates that there can be dramatic differences from one location to another, even for "nearby" locations. I don't know whether borehole results (ground surface temperature) exhibit the same high degree of spatial correlation that surface air temperature does. It certaintly drives home the need to use averages over as many samples as possible, and that single-borehole results cannot be extrapolated to large areas with as much confidence as air-temperature indicators.

The strongest indication of borehole reconstructions for past temperature is that they agree so well with other proxy reconstructions. When multiple lines of evidence converge, the probability is much greater that they are correctly characterizing what's happened in the past. It can happen that multiple lines of evidence are wrong *and* wrong in the same way (witness the problems with satellite and balloon-borne temperature measurements), but that's pretty rare; when one or both lines of evidence is wrong it's very unlikely that they'll agree on the same wrong result.

Hence boreholes are an excellent "check" for the last 400-500 years. I'd guess that that's one of the main reasons the NAS (National Academy of Sciences) report on paleoclimate reconstructions expressed higher confidence in the last 400 yrs than further back in time. Unfortunately beyond 500 yr ago, the inherent uncertainty in borehole reconstructions gets much bigger, and not all boreholes go deep enough to reach back that far in time anyway. Also, further back in time not only the temperature resolution, but the *time* resolution deteriorates rapidly.

So, for the last 400 - 500 yr I'd take borehole data, areally averaged, as good confirmation of the trend indicated by proxy reconstructions in general. Beyond that, they're "fuzzy" enough that I wouldn't rely on them in isolation.]
Heiko Gerhauser // August 21, 2007 at 7:57 pm | Reply

Thanks, that was helpful.
Phil. // October 5, 2008 at 9:50 pm | Reply

Tamino, when developing the method for inverting light scattered data into particle size for distributions of particles I found that the method of using a distribution functional form and fitting worked very well. However, I also found that the method of using a non-negative least squares worked well since the true equation includes a non-negativity constraint. I don’t know if such an approach would be helpful in this case?

	suricat on Vapor Lock
	Timothy Chase on Key Messages
	suricat on Vapor Lock
	Jon P on Open Thread #14
	dhogaza on Open Thread #14
	Lazar on Key Messages
	george on Key Messages
	Hank Roberts on Open Thread #14
	t_p_hamilton on Open Thread #14
	Chris S. on Embarrassing Questions
	Igor Samoylenko on Embarrassing Questions
	Chris S. on Vapor Lock
	Horatio Algeranon on Open Thread #14
	Timothy Chase on Key Messages
	Frank O'Dwyer on Key Messages

S	M	T	W	T	F	S
« May				Jul »
					1	2
3	4	5	6	7	8	9
10	11	12	13	14	15	16
17	18	19	20	21	22	23
24	25	26	27	28	29	30

Open Mind