Paleothermometry redux

Sir,

In “Paleothermometry by control methods” (MacAyeal and others, 1991), we presented a mathematical method for estimating past surface-temperature history from ice-sheet borehole-temperature profiles (the paleothermometry problem). Dahl-Jensen and others (1993) have suggested that our solution of the paleothermometry problem fell short of what is needed to achieve meaningful paleoclimatic inference. Naturally, we focused in our paper primarily on the virtues of the control method, perhaps to the detriment of a sufficient discussion of the vices. Therefore, we appreciate Dahl-Jensen and others’ comments for drawing attention to questions which our paper left unanswered. Answers to many of these questions, and further analysis of the Dye 3 demonstration test we presented in our paper, are provided in Firestone (1992).

The two main points we shall address in this letter are: (1) that our paper fails to demonstrate that a quantitative estimate of surface-temperature uncertainty is possible (let alone satisfactory in the examples presented in our paper); and (2) that the improved fit between calculated and observed temperature profiles that our method facilitates is illusory (i.e. deceptive). We agree with both criticisms. We disagree, however, with what we believe Dahl-Jensen and others imply: first, that failure to quantify uncertainty results from the inadequacy of our particular method alone; and secondly, that our method cannot avoid “overfitting” the data.

With respect to the first implication, we believe that all inverse methods are inadequate in quantifying uncertainty in the paleothermometry problem. A notorious difficulty of this problem is that there exists a class of possible surface-temperature histories which has no measurable effect on the borehole temperatures. The diurnal temperature cycle prior to 5000 years ago is an extreme example of one such history. It is therefore impossible for any mathematical analysis of a borehole-temperature profile, and in particular a least-squares analysis such as ours, to distinguish between these histories. It is in this sense that we agree with Dahl-Jensen and others; our method, indeed all methods, fail to quantify uncertainty.

In retrospect, we recognize that we misstated one of the conclusions in the abstract of our paper; namely, that the uncertainty of our method “can be established quantitatively”. We correct this misstatement by replacing the term uncertainty with the term resolution which, as shown below, is a more specific measure of uncertainty. Figure 9 of our paper, and the discussion surrounding equations (44) and (45), present the quantitative assessment of resolution we intended to highlight in our abstract.

With respect to the second implication, we reiterate what we stated in our paper; the control method does allow for a trade-off between fitting noisy borehole data and satisfying independent performance constraints such as estimated climate histories. The transformation of equation (7) into equation (8) demonstrates this trade-off. The point Dahl-Jensen and others made, and with which we agree, is that this trade-off should be carefully engineered to restrain the method from interpreting unmeaningful measurement noise in the borehole profile. We shall outline how this can be done.

Insofar as other points have not been discussed in sufficient detail to satisfy Dahl-Jensen and others, and conceivably other readers as well, we shall re-visit the paleothermometry problem in sufficient detail to diagnose the unsatisfactory results in the demonstration tests of our paper. The lengthy analysis that follows reflects our continuing interest in applying inverse methods to glaciological problems. In particular, we believe that the paleothermometry problem serves as a metaphor for a large class of glaciological inverse problems which are burdened by imprecise methods. An example is the problem of deducing basal traction from measurements of velocity at the surface of a glacier. The surface velocity is analogous to the borehole-temperature profile, and the basal-traction field is analagous to the surface temperature history. As demonstrated by Bahr and others (1992), this problem is ill-posed in the same sense as is the paleothermometry problem. Techniques developed here may therefore have applications that extend beyond the narrow subject of borehole-temperature analysis.

We set forth several goals to accomplish in this letter. First, we wish to show that unsatisfactory aspects of our demonstration tests do not stem from the control method, but rather from the way in which we defined our particular performance index (i.e. the way in which we defined the paleothermometry problem). Second, we wish to develop an integral-equation approach using continuous variables as a means of separating the fundamental properties of the problem from the details associated with finite-difference discretization. Third, we wish to derive a formal correspondence between the control method and other least-squares methods in common use (e.g. Anderssen and Saull, 1973; Wang, 1992). Fourth, we wish to cut through the exoskeleton of mathematical formalism that may have left some readers of our paper mystified as to what the paleothermometry problem is and how it can be solved; we re-develop our method using a familiar eigenfunction (or eigenvector) approach.

1. What is the Paleothermometry Problem?

The ice-sheet heat-transfer problem is non-linear because the thermal properties of ice are temperature-dependent and because the equations of heat transfer couple to the equations of momentum conservation through thermal advection and the variation of ice flow with temperature. For the purpose of answering the questions posed by Dahl-Jensen and others, however, we shall consider the simpler linear problem. In addition, since they are not germaine to the questions, we shall adopt a homogeneous initial condition and a homogeneous boundary condition at the base of the ice sheet (i.e. zero geothermal-heat flux). In this circumstance, we may immediately write the relationship between the known borehole-temperature profile, θ(z), and the unknown surface-temperature history, T_a(t), as

(1.1)

where the variables are defined as in our paper, and G(z,t) is a Green’s function (Carslaw and Jaeger, 1988, chapter XIV). In our paper, we used a discrete, finite-difference treatment of the heat-transfer problem. Our method was thus described in terms of matrix equations. For clarity, here we shall adopt a continuous view of the heat-transfer problem; this leads to integral equations such as Equation (1.1).

As defined in Equation (1.1), θ(z) is a linear functional of the surface-temperature history; in particular, it is an integral transform of T_s(t). The Green’s function, G(z,t), describes the physics of the problem and determines the mathematical properties of the function θ(z), such as that it be analytic. Comments by Dahl-Jensen and others concerning our treatment of vertical ice velocity and accumulation rate are simply details absorbed into G(z,t); they need not concern us here.

Before outlining our approach to the paleothermometry problem, it is important to distinguish between three related inverse problems. The mathematical apparatus used to solve the third problem is the subject of our paper. We introduced the other two, and in particular the first, to emphasize the fact that the way in which the problem is defined strongly determines the properties of the resulting solution (e.g. oscillations).

PROBLEM 1. Given Borehole Data, θ_B,(z), a Definition of the Heat-Transfer Physics, Boundary and Initial Conditions, and G(z, t),Find T_S(t) Using Equation (1.1)

This is not the problem we solved in our paper. In fact, it may not be possible to solve this problem in general because borehole data, θ_b(z), may not satisfy the mathematical properties required of solutions of the heat-transfer equations (i.e. analyticity). In addition, Equation (1.1) is an integral equation of the first kind, that is, T_s(t) appears only under the integral sign (Courant and Hilbert, 1953, p. 159). This kind of integral equation is notoriously difficult to solve and we did not do so in our paper.

One of the properties of integral equations of the first kind that is germaine to the question of uncertainty is that they can have a null space, i.e. there might be a function (or functions) such that

(1.2)

(By setting the righthand side of Equation (1.2) to zero, we mean 0 K; in other words, the history has no effect on the borehole-temperature profile.) If a solution to Equation (1.2) exists, then the surface-temperature history which solves problem 1 can be accurate only to within a multiplicative constant of this solution. On intuitive grounds, we believe that a non-trivial solution to Equation (1.2) is physically unlikely; i.e. it is unrealistic to think that there are climate histories which have literally no effect on borehole temperatures at time t_f .

What is intuitively more reasonable to us, as we stated in our paper, is that finite-precision arithmetic and limited borehole-measurement sensitivity make the possibility of non-trivial solutions to Equation (1.2) a certainty in practical situations. Finite-precision arithmetic (such as that which is performed on a computer) suggests that some temperature histories cannot be operated on by the integral transform expressed by Equation (1.1) without serious arithmetic error. In fact, this arithmetic error may make the associated borehole-temperature profile that results from to be zero. More importantly, limited measurement sensitivity implies that there are an infinite number of borehole-temperature profiles θ(z) that cannot be distinguished from zero. The associated with these profiles, which cannot be distinguished in the solution, constitute the null space of the problem in Equation (1.1).

In this circumstance, no inverse method can guide us in the choice between the infinite number of equally satisfactory solutions to Equation (1.1) which differ only by solutions of Equation (1.2); or, in other words, by their projection into the null space. Typically, we choose to accept a solution which has a zero projection. In the Dye 3 demonstration of our paper, the projection is determined by the climatology η(t) we imposed in equation (8) of our paper. A description of the temporal structure which has been excluded from our solution serves as the best statement we can provide about the uncertainty of the solution. Since any function that belongs to the null space can be added to our solution without affecting the fit to borehole data, it is not possible to represent this uncertainty quantitatively.

PROBLEM 2. Approximate the Solution of Equation (1.1) in a Least-Squares Sense

Here we recall the definition of the performance index, equation (7) in our paper, but make use of the expression in Equation (1.1)

(1.3)

Henceforth, we shall drop the subscript from θ_b(z) to show that the measured profiles no longer need to satisfy the stringent mathematical properties of exact solutions to the heat-transfer problem. This is indeed the main benefit provided by adopting a least-squares approach.

Problem 3. Approximate the solution of Equation (1.1) in a least-squares sense; but add subsidiary performance conditions, such as a cost function that depends on deviations between T_s(t) and a climatology η(t)

In this circumstance, the performance index J of problem 2 is modified as follows

(1.4)

To address better the questions Dahl-Jensen and others have about the choice of the parameter e in our paper, we have replaced e with a new mixture parameter α ∈ [0, 1] which describes the trade-off between the fit to borehole data and the fit to climatology. The above performance index is the same as that defined in equation (8) of our paper when α/(1 – α) = ϵ. We deviate from the original notation to explain better the role of preconceived climatology and to indicate better what considerations must be used to choose α (formerly ϵ). Observe that when α = 0, borehole misfit is the only consideration used to determine the surface-temperature history (problem 3 becomes problem 2); and, when α = 1, borehole misfit becomes irrelevant; only deviation from climatology is important. (We note that α could be defined as a function of time. For simplicity, we have refrained from doing so here.)

To minimize J defined in Equation (1.4), we follow the course of action described in our paper and consider the conditions to be satisfied when the variation of J is zero, δJ = 0:

(1.5)

If we define

(1.6)

and

(1.7)

then Equation (1.5) reduces to finding the solution, T_s(t), of

(1.8)

Equation (1.8) is an integral equation of the second kind (Courant and Hilbert, 1953, p. 112), and the kernel K(t,t’) is symmetric. The solution of Equation (1.8) constitutes what we define here to be the least-squares solution to the paleothermometry problem.

2. A Family of Paleothermometry Problems

In shifting our attention from problem 1 to problem 3, we have changed the mathematical and physical nature of the paleothermometry problem. This shift is necessary because it allows us to avoid serious mathematical difficulties. Problem 1 leads to an integral equation of the first kind, which may not have a solution due to the incompatibility between the data, θ_b(z), and solutions of heat-transfer equations (i.e. analyticity). Problem 3, however, leads to an integral equation of the second kind which always has a solution, often an infinite number of solutions. Problems 1 and 3 also differ in that K(t,t’) is a symmetric kernel whereas G(z, t) is not. In addition, the data in problem 1, θ_b(z), can be rough; whereas, the data in problem 3, Θ(t), are smoothed by the integral transform implied by Equation (1.7). From what we know about singular-value decomposition (SVD), we expect the kernel K(t,t’) to have eigenvalues that are near zero. As we shall show later, this motivates the need to introduce the climatology and a non-zero a. Even with all the advantages gained by abandoning problem 1 in favour of problem 3, we must still deal with the possibility that the kernel K{t,t’) may have singularities (i.e. may not be square-integrable), and so may be very difficult to work with in a specific application.

Efforts to define a workable problem, both here and in our paper, call attention to the fact that there is no universally accepted definition of the paleothermometry problem. If problem 1 were workable, methods for its solution could be compared in direct and absolute terms. In reality, problem 1 is unworkable (ill-posed); thus, in comparing methods, we must remember that there can be subtle differences in the problems to which the methods are applied. In other words, differences between our method and that used, for example, by Dahl-Jensen and Johnsen (1986) are obscured by the fact that the two methods may not have been applied to precisely the same problem.

3. Is the Control Method Relevant?

Thus far, we have said nothing about control methods. We shall next argue that the control method we used in our paper to solve the paleothermometry problem is irrelevant from the standpoint of Dahl-Jensen and others’ commentary. In other words, this method does nothing to determine the properties of the solution to the paleothermometry problem. (Heat-transfer physics determines these properties.) The control method is simply one of many strategies for obtaining the conditions which minimize J defined by Equation (1.4). In other words, the control method provides an alternative avenue to the derivation of Equation (1.8), and nothing more.

To show the correspondence between the control method and the least-squares problem represented by Equation (1.8), we develop the Green’s function which solves the forward problem represented by equations (1)–(4) of our paper. This forward problem for T(z, t) is restated here with a homogeneous boundary condition at z = H (the bed of the ice sheet) and a homogeneous initial condition to simplify the analysis

(3.1)

for 0 < z < Η and 0 < t < t_f, and

(3.2)

(3.3)

(3.4)

Here, and elsewhere, subscripts (t and z) denote partial differentiation by the subscripted variable. The solution, T(z, t) to Equations (3.1)– (3.4) may be written in terms of a Green’s function, which satisfies the following equations:

(3.5)

for 0 < z < Η and 0 < t < if, and

(3.6)

(3.7)

(3.8)

where δ(z – ξ) is the delta function. Notice that Equation (3.5) represents the adjoint form of Equation (3.1). The expression which gives T(z, t) in terms of G(z, t; ξ) is found by manipulating the following integral

(3.9)

Performing integration by parts, we get

(3.10)

Consequently, using boundary conditions,

(3.11)

Thus, ft,

(3.12)

which is a restatement of Equation (1.1) with kG_z (0, t; ξ) identified with G(ξ,t).

To minimize the performance index expressed in Equation (1.4) using the control method we solve the adjoint trajectory problem given by equations (17)–(20) in our paper for λ(z,t) (the motivation for deriving the adjoint trajectory problem is outlined in our paper). The solution is expressed in terms of the same Green’s function G(z, t; ξ) used to solve the forward problem in Equations (3.1)– (3.4),

(3.13)

To determine the surface-temperature history, T_s(t), we simply solve equation (21) of our paper, which is restated below (recall that ϵ = α/(1 – α)

(3.14)

Making use of Equation (3.13), and recalling that G(0, t; ξ) = 0, this expression becomes

(3.15)

Replacing T(ξ, t _f) with the expression given by Equation (3.12), we obtain

(3.16)

which is a restatement of Equation (1.8).

We come now to a crucial point in our response to Dahl-Jensen and others’ commentary. The correspondence between Equations (1.8) and (3.16) shows that the control method developed in our paper serves only to assist in the solution of problem 3. The surface-temperature history produced by the control method is identical to the least-squares solution (i.e. the solution to Equation (1.8)). The fact that a control method was used in our paper to obtain the least-squares solution is thus a minor detail, and is therefore irrelevant to our discussion of Dahl-Jensen and others’ comments concerning the mathematical properties of the solution to our Dye 3 demonstration problem. In what follows, we shall describe these properties without making further reference to the control method.

4. The Least-Squares Solution: Continuous Case

We are now ready to develop a general solution T_s(t) to problem 3 by solving Equation (1.8). If the symmetric kernel K(t,t’’) were square-integrable (finite when t = t’ = t_f , if, for example), we could exploit the fact that K(t, t’) has eigenvalues and eigenfunctions (Courant and Hilbert, 1953, p. 122) to describe T_s(t). By square integrability we mean

(3.1)

where M is a finite bound. We state without proof, however, that the K(t,t’) associated with an inverse linear heat-diffusion problem is assured not to be square integrable (this will be demonstrated explicitly for the heat-conduction problem in a semi-infinite solid below). Strictly speaking, we are therefore not able to use an eigenfunction analysis comparable to the eigenvector analysis made possible by adopting discrete, finite-difference versions of the variables.

The lack of square-integrability is yet another indication of the great difficulty of the paleothermometry problem. Even with all the benefits gained by adopting a least-squares approach, the problem is still intractable owing to the singularity in K(t, t’) which, stated without proof, occurs at t = t’=t_f. Once a numerical method is adopted, such as in our paper, the kernel in Equation (1.8) is converted to a square, symmetric matrix, and the question of square-integrability is converted to a question of how to cope with an ill-conditioned matrix (a matrix with a near-zero determinant). Eigenvalues and eigenvectors can be found for the matrix representing the discrete, numerical problem even when the matrix is ill-conditioned (using, for example, singular value decomposition). We shall thus proceed with an eigenvalue and eigenfunction approach to solving the continuous problem by replacing the kernel K(t,t’) with an (unspecified) approximation K(t, t’) that is square integrable. We must not forget, however, that error is introduced into the solution of Equation (1.8) when we make this replacement. We shall not attempt to characterize this error since it is not pertinent to the analysis of the Dye 3 demonstration test we discuss below.

Keeping the above remarks in mind, we carry on with the formal solution of Equation (1.8) subject to the approximation of K{t,t’) by K(t,t’) Let and denote the eigenvalues and eigenfunctions associated with (Courant and Hilbert, 1953, p. 122):

Furthermore, let form a complete set, i.e. arbitrary functions can be expressed as linear combinations of these functions. Then, if we expand T_s , Θ and η as follows

(3.3)

(3.4)

(3.5)

Equation (1.8) reads

(3.6)

and the solution T_s(t) becomes

(3.7)

We have now produced a description of the formal solution to the least-squares inverse problem with climatology (problem 3). This formal solution is abstract in the sense that we do not actually specify the eigenfunctions or eigenvalues. The advantage gained by expressing T_s(t) as an infinite sum of eigenfunctions is that we can discuss the general properties of uncertainty, and the role of α, without the distractions of a specific application.

5. Least-Squares Solution: Discrete Case

Equation (1.8) and its solution in Equation (4.7) are the continuous versions of the matrix equation and the finite-difference version of T_s(t) derived in our paper. To display the parallel between the continuous and the discrete solutions to the paleothermometry problem, we re-express the solution derived in our paper as follows. First, we recall that in the above derivation, we have assumed the initial and basal boundary conditions to be homogeneous; thus, we take T(l) = 0 and G = 0, where T(l) and G are the discrete analogs of the initial and basal boundary conditions as defined in our paper. Secondly, we use the description of the adjoint trajectory, equations (36) and (37) of our paper, to derive the expression for T _s = (T _s(Δt), T _s(2Δt), …, T _s(nΔt ), …, T _s(N – 1)Δt)) ^T ∈ R ^{N – 1} in terms of T(N) – Θ, where T(N) = (T(0), T(Δz), T(2Δz), …, T(nΔz ), …, T((M – 1)Δz)) ^T ∈ R^M is the discrete analog of the temperature profile at time t_f = (N - l)∆t, and Θ is the discrete analog of the borehole data. Thirdly, we use equation (27) of our paper to express T(N) in terms of the T_s. Fourthly, we invoke ϵ = α/(1 – α). The result of these four steps is the discrete analog to Equation (1.8):

(5.1)

where the vector H ∈ R ^{N – 1} is the discrete analog of the climatology η(t), and the vector Γ ∈ R ^{N – 1} is the discrete analog of the continuous function θ(t), and is given by

(5.2)

where the rectangular Μ×(N-1) matrix g is the discrete analog of the Green’s function G(z,t), and is given by

(5.3)

The (N - 1) × (N - 1) matrix Κ is the discrete analog of the kernel K(t,t’) appearing in Equation (1.8), and is given by

(5.4)

and where A, Β and C are matrices or vectors defined in our paper, and the superscript Τ denotes the transpose. We remind the reader that A, Β and C are determined by the physics of the heat-transfer process only; thus, Κ and the properties of its eigenvectors are independent of the climatology or borehole data. (Note two typographical errors in the definition of A in equation (28) of our paper. The terms which read (2κ)/(2Δz ²) should not have the 2 in the denominator; and the term in the last column of the second-to-last row should have a minus sign where a plus sign appears.)

The solution to Equation (5.1) can be expressed formally in terms of the eigenvalues and eigenvectors associated with the symmetric matrix Κ : R ^{N – 1} → R ^{N – 1}. Let and respectively, denote the eigenvectors and eigenvalues of K. In other words, let

(5.5)

The symmetry of Κ assures us that the eigenvectors can be made orthonormal, and will span the vector space R^N-1. In parallel with the continuous version of the problem, we expect the matrix k to be ill-conditioned in the sense that some of its rows may have vanishingly small (near-zero) elements, and this may lead to a near-zero determinant. This is the discrete analog of the difficulty associated with the non-square-integrability of the continuous kernel K(t,t’). As a result of this ill-conditioning, it may be very difficult to compute all the eigenvalues and eigenvectors of K. In this circumstance, singular-value decomposition (Press and others, 1989, p. 52) can be used to determine s of the eigenvalues and eigenvectors. The remaining Ν - 1 - s eigenvalues can be taken as zero, and the remaining N-1- s eigenvectors can be chosen via an orthonormalization process to ensure that the set spans the vector space R^N-1.

We may expand T_s , Γ and H as follows

(5.6)

(5.7)

(5.8)

where the coefficients and are given by

(5.9)

(5.10)

(5.11)

Substituting Equations (5.6)– (5.8) into Equation (5.1), and making use of Equation (5.5), jnves a relation between the unknown coefficients and the known coefficients and

(5.12)

for i = 1, ..., {Ν - 1). The solution of the discrete form of the paleothermometry problem is thus

(5.13)

The advantage gained by expressing T_s in terms of the eigenvectors of Κas opposed to the equally valid expression

(5.14)

is that the mathematical properties of T_s can be immediately appreciated. The fact that are independent of the borehole data and climatology tells us how the solution depends on the physics of the heat-transfer process. The eigenvectors depend solely on the matrix Κwhich depends, in turn, on a statement of the discrete form of the forward problem.

We note in passing that Equation (5.1) is readily solved by standard linear-algebra techniques and the result is the same as that which comes from solving equation (40) in our paper for T(N) and then using the adjoint equations, equations (36)–(38) of our paper, to get T_s. We chose to work with equation (40) of our paper instead of Equation (5.1) because the matrix in equation (40) is dimensioned by the number of vertical grid points (101 for the Dye 3 example), whereas the matrix in Equation (5.1) is dimensioned by the number of time steps (1199 for the Dye 3 example).

The close parallel between the continuous and discrete versions of the solution to problem 3 reminds us that the basic properties of the solutions are independent of the discretization. We thus expect the physics of the problem to have the greatest effect on the intrinsic properties of the solution. Details of the numerical method used to find the solution are important, as Dahl-Jensen and others have pointed out, but are not the sole cause of what Dahl-Jensen and others cited as unsatisfactory in our Dye 3 surface-temperature history.

6. Regularization and the Role of Climatology

We now answer the question concerning trade-off between the fit to borehole data and the fit to climatology. For both the continuous and discrete versions of the least-squares solution, the mixture parameter, α, determines the trade-off between the coefficients of the climatology, c_i , and the coefficients of the borehole data b_i. This is why we emphasize that our method does indeed allow for imprecision in the borehole data.

Whether we failed to select a value of the trade-off parameter that would take adequate account of borehole-data error in our Dye 3 demonstration is an issue we are not able to resolve without further work. We anticipate that an objective choice of α will depend on an evaluation of the relative confidence ascribed to the borehole-temperature measurements (which determine the 6_is) and the climatology (which determine the c_is). For example, one might choose , where is a statistical measure of temperature-measurement uncertainty and σ_c ² is a measure of the uncertainty in the climatology.

There is another important role played by the mixture parameter α which must be considered when selecting its value. It serves to regularize the problem of determining T_s(t). Even without specifying the approximation K(t,t’) to the symmetric kernel K(t,t’) in Equation (1.6), or the symmetric matrix K in Equation (5.4), we can anticipate that the eigenvalues μ_i will tend to converge to 0 as i → ∞ in the continuous case (Courant and Hilbert, 1953, p. 130) or as i → (N – 1) in the discrete case. Thus, when α = 0 (i.e. when we force our retrodiction to depend only on achieving the best fit between predicted borehole temperatures and data), the expressions in Equations (4.7) and (5.13) could become divergent because of zeros, or extremely small numbers, in the denominator of the expansion coefficients. By specifying α > 0, the denominator remains finite when μ_i . The series expressions thus become well-behaved from an arithmetic standpoint. The benefit of regularizing the problem is offset by the fact that, as the problem becomes better conditioned from a mathematical standpoint, the retrodicted surface-tern perature history depends more on the climatology η(f) and less on the borehole-temperature data.

We were well aware of the important role of α and the preconception η(t) in regularizing the Dye 3 demonstration problem (see figure 11 of our paper). The value of ϵ = α/(1 – α) we used in our demonstration was chosen solely on this basis (regularization). We urge future researchers to take note of the fact that an objective choice of α (or the original ϵ) depends on two needs: (1) the mathematical need to regularize an otherwise ill-conditioned problem (as represented by the series-expansion representations of the solutions to the continuous and discrete versions of the problem), and (2) the strategic need to optimize the balance between fitting borehole data and matching preconceived climatology. The need to regularize the problem can be met by consideration of the physics only (physics determines the eigenvalues and eigenfunctions). The need to optimize the balance requires an understanding of the relative accuracy of the borehole data with respect to the climatology.

7. Why does the Solution Oscillate?

The form of the surface-temperature solution given by Equations (4.7) and (5.13) allows us to explain why the solution produced in the Dye 3 demonstration oscillates. Although we have not specified the exact form of the or we can anticipate that they are functions or vectors which exhibit oscillatory behavior. The surface-temperature history, T_s(t) or T_s , may thus be expected to oscillate too.

To demonstrate this point, we plot in Figure 1 several of the eigenvectors of the matrix K defined in Equation (5.4) using the A, Β and C developed for the “demonstration using synthetic data” described in our paper. (In this circumstance Κ is a 49 × 49 matrix, and can be decomposed into eigenvalues and eigenvectors on a Macintosh computer. The matrix Κ from our Dye 3 demonstration is 1199 =× 1199, too big to fit into our computer.) As seen in the figure, the eigenvectors are oscillatory, and this imposes oscillations of the kind seen in our demonstration analysis of Dye 3. Noteworthy properties of the eigenvectors are: (1) that they decay as t = 0 is represented by 5000 years BP in Figure 1), (2) they oscillate with increasing frequency as t → t _f, and (3) the most oscillatory eigenfunctions are associated with the smallest eigenvalue. This is why the history produced in our Dye 3 demonstration (see the caption for figure 6 of our paper) displayed what we referred to as “insignificant” (meaningless) oscillations as t → t _f if.

Fig. 1. Five eigenvectors associated with the largest eigenvalues of the 49 × 49 symmetric matrix Κ in Equation (5.4). Most of these vectors are oscillatory, and the amplitude of their oscillation increases as t → t _f.

The eigenvector analysis presented in Figure 1 confirms Dahl-Jensen and others’ concern about the formulation of our problem as a source of oscillation. Oscillations are indeed a natural consequence of solutions of problem 3 (the least-squares problem). The physics of the heat-transfer process (the Green’s function G(z, t) or the matrix g) and the least-squares formulation determine the properties of the kernel and matrix K. These properties, in turn, determine the oscillatory nature of the eigenfunctions and eigenvectors from which the solutions, T_sub(t) and T_s respectively, are composed. The distinction should be made, however, that the oscillations DahlJensen and others found unacceptable are inherent in all least-squares solutions of our particular paleothermometry problem. Indeed, oscillations appear to be inherent in least-squares solutions of problems similar to ours such as the rock-borehole paleothermometry problem presented by Wang (1992). Had we defined our problem differently, perhaps by penalizing the time derivative of T_s(t) in the performance index, then we would have obtained a less oscillatory solution.

8. The Rock-Borehole Paleothermometry Problem

Dahl-Jensen and others have pointed to the literature concerning paleothermometry problems involving temperature profiles in rock, and suggested that a comparison between our method and those described in this literature would be a “useful first step”. We have taken this suggestion and have made two comparisons. First, we asked a student at the University of Chicago to repeat the analysis of an observed temperature profile in permafrost from northern Alaska (site AWU) reported by Lachenbruch and Marshall (1986). The student used the techniques described in section 5 of this letter. (We are indebted to Ms H. Wang for performing this analysis.) Secondly, to examine correspondence with the theoretical development presented in this letter, we found the exact expressions for G(z,t) and Κ(t,t’) associated with the thermal conduction problem in a homogeneous half-space Earth.

The results of these two comparisons are interesting, but shed no new light on the comments made by Dahl-Jensen and others. In short, the analysis of the northern Alaskan temperature profile produced a surface-temperature history that oscillates with increasing frequency as t → t _f. As with our Dye 3 demonstration, the fit to borehole data could be improved with seemingly unlimited precision (depending on the choice of a). The cost of this improvement is an increase in the complexity of the surface-temperature history. Thus, while our method can fit the data presented by Lachenbruch and Marshall (1986), indeed with greater precision than is justified by the data, the resulting surface-temperature history can be so complex (depending on the choice of a) that it fails to confirm the conclusion reached by Lachenbruch and Marshall (1986); namely, that the local climate has warmed by 2–4 deg in the last few decades to a century.

The fact that two so widely different surface-temperature histories can explain the same borehole temperature profile (as is the case for our Dye 3 demonstration and the analysis of Dahl-Jensen and Johnsen (1986)) troubles us. Can we conclude that the simplest of the two histories is a better representation of the true climate history than the other? Perhaps this question may have an affirmative answer; however, the analysis presented here and in our paper cannot provide the answer to this question.

For the second comparison, involving the heat-conduction problem in a semi-infinite Earth, we display the Green’s function, G(z, t), and symmetric kernel, K(t,t’), introduced in section 1. The temperature profile in the semi-infinite region z < 0 is assumed to satisfy the following conditions

(8.1)

(8.2)

(8.3)

(8.4)

where subscripts denote pardal differentiation by the subscripted variable, and where k is assumed constant. A non-homogeneous initial condition and a steady, nonzero geothermal flux at z → –∞ need not concern us here.

The solution T(z,t) to Equations (8.1)– (8.4) is expressed as the convolution of the Green’s function, G(z, t), with the surface-temperature history, T_s(t), as in Equation (1.1). The Green’s function is readily determined (Carslaw and Jaeger, 1988, p. 62);

(8.5)

This Green’s function is used to obtain the symmetric kernel using Equation (1.6),

(8.6)

We remark that K(t, t–) for this problem is not square integrable. It blows up as t, t′ → t _f. This reflects our physical intuition that surface-temperature history in the most recent part of the past has a stronger influence on the borehole-temperature profile than the history in the more distant past.

9. Uncertainty

We now address Dahl-Jensen and others’ main question: is our method able to quantify uncertainty in the surface-temperature history? There are many factors which contribute to uncertainty. The emphasis in our paper was on the uncertainty which stems from the physics of heat diffusion; in other words, from the dissipation of thermal memory. This source of uncertainty is independent of borehole-temperature data, and its quantification yields the lower bound on uncertainty which would be valid if the data were error-free. As demonstrated in our paper, the model-resolution matrix provides a convenient quantitative description of this lower bound. (A careful consideration of the model-resolution matrix, incidentally, will answer the questions Dahl-Jensen and others raised about our demonstration test with synthetic data.)

A second contribution to error discussed in our paper concerns the preconceived climatology. Again, this contribution is independent of measurement error in the data. We refrained from analyzing this source of error for the Dye 3 demonstration because our choice of climatology was intentionally oversimplified.

We are unsure of how to describe quantitatively the uncertainty generated by borehole-data error. Dahl-Jensen and others’ reminder that this error could have coherent spatial structure associated with borehole-fluid convection leads us to believe that this error will translate into a non-uniform temporal distribution of surface-temperature uncertainty. We lack sufficient information about the data to proceed with a formal analysis of this error. We therefore offer merely a suggestion as to how one might approach such an analysis in future work.

We suggest that the observed borehole-temperature profile be expressed as the sum of the actual temperature profile and an error profile. In this situation, each expansion coefficient of Equation (4.4) becomes the sum of two terms

(9.1)

where , is the expansion coefficient for the actual borehole-temperature profile, and e_i is the expansion coefficient for the error profile Ɛ(z). Recalling Equation (1.7), [G(z,t)S(z)dzdt. (9.2)

(9.2)

Taking α = 0 to simplify the discussion, we can write Equation (5.7) as

(9.3)

where is the uncertain surface-temperature history, and where the summation has been split into two parts to highlight the effect of truncation. The second term on the righthand side of Equation (9.3) represents the part of the measurement error which can be quantified, provided the e¿s, for i = 1, ..., Ν, can be estimated from information about the borehole-data collection process and anticipated convection cells. The last term on the righthand side of Equation (9.3) represents the part of the measurement error which cannot be quantified because of our inability to sum the series to infinite terms. We can anticipate that as i → ⇮ the eigenvalues μ_i → 0, and that the eigenfunctions φ(t) contain smaller time-scale structure. The last series on the righthand side of Equation (9.3) may therefore diverge when the e_is for i > Ν are incorrectly estimated. Contributions to the surface-temperature history due to this error would have very short time-scales but very large, or infinite, amplitude. Physically, this amounts to saying that the uncertainty associated with a temperature cycle of short duration, a diurnal cycle in the ancient past for example, cannot be quantitatively assessed.

To quantify the second term on the righthand side of Equation (9.3) in a manner such as would yield error bars, or confidence intervals, on the retrodicted surface-temperature history, we suggest that a family of coefficient sets be estimated. This could be done using information about the accuracy of the measurement devices and the anticipated spatial structures of the convection cells. The family of coefficient sets can then be transcribed into a family of surface-temperature histories using the eigenfunction expansion in Equation (9.3). The envelope which bounds this family of histories would then constitute an error-bar representation of measurement error.

10. Conclusion

We believe that the analysis presented above clarifies our position concerning the quantification of uncertainty and the trade-off between fitting noisy borehole data and satisfying secondary performance constraints, such as simplicity. As to the first position, no absolute measure of uncertainty is possible in the paleothermometry problem due to the inherent memory loss of heat diffusion. What is possible to quantify, however, is the resolution of an inverse method such as ours. As to the second position, a tendency to overfit measurement noise in borehole-temperature data can be restrained through the imposition of precisely described secondary constraints on the surface-temperature history. These secondary constraints determine the functional properties of the surface-temperature history, such as the degree to which it oscillates. There are no strict rules which define these secondary constraints; thus, it is perfectly natural to expect two or more widely differing surface-temperature histories to result from the analysis of the same borehole-temperature data.

Concerning the question of whether solutions of borehole paleothermometry problems could serve as a check on the interpretation of ice-core isotope stratigraphy, we direct the reader to Firestone (1992). To our understanding, the interpretation of isotopic stratigraphy depends as much on an empirical relationship between present-day surface temperature and present-day isotopic compositions of precipitation as it does on complex theories of the hydrologie cycle. We are aware of no proof that this empirical relationship cannot change with time. We thus re-affirm our interest in developing better methods to solve the paleothermometry problem. Furthermore, we suggest that this development is no less justified than the continued development of ice-core geochemical methods for inferring the paleoclimate.

The accuracy of references in the text and in this list is the responsibility of the authors, to whom queries should be addressed.