Introduction
Recent revisions of carbon-14 dating chronology and eustatic sea-level records suggest that oxygen-isotope profiles in Greenland ice cores may not provide reliable records of past surface temperature (Reference FairbanksFairbanks, 1989). In particular, the isotopic signal previously attributed to the Younger Dryas cold period 10 000 years BP alternatively may indicate isotopic modificationof the oceanic reservoir from which Greenland precipitation is derived. Reference FairbanksFairbanks (1989) identified glacial meltwatcr run-off as the likely cause of this modification. This interpretation has generated a controversy concerning the areal extent, timing and cause of the Younger Dryasevent (Reference Broecker and DentonBroecker and Denton, 1989). In an effort to resolve this controversy and to verify the interpretation of ice-core isotope stratigraphy, we explore the use of ice-sheet temperature profiles as paleothermometers.
The interpretation of temperature-depth profiles to infer past climate is the subject of a large body of literature (see, for example, Reference Paterson and ClarkePaterson and Clarke, 1978; Reference Budd, Young and RobinBudd and Young, 1983; Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen, 1986; Reference RitzRitz, 1989). Sophisticated numerical models of ice-sheet heat transfer are common and few improvements in their accuracy or fidelity to physical processes remain to be developed. Our paper does not concern these models but rather the way in which they are used to estimate past climate.
In this paper we introduce control methods as a means of recovering paleoclimatic information from ice-sheet temperature profiles. Control methods grew out of the calculus of variationsduring the last 30 years. The typical problem that these methods address consists of optimizing the terminal state of a system which is constrained to evolve according to a prescribed set of differential equations. Control methods might be used, for example, to determine the best time to launch a spacecraft in order to achieve a desired planetary rendezvous. As this example suggests, the control variable is often an initial or boundary condition.
We wish to exploit the fact that control methods deal with optimal boundary conditions, together with equations of evolution as constraints, to recover the surface temperature from ice-sheet boreholes. The present temperature-depth profile represents the terminal constraint on athermal evolution contraled by surface-temperature history. We thus appeal to control methods to give us the history that is most compatible with the observed terminal state. Here, we provide a tutorial on how to use control methods to determine the surface-temperature history fromborehole temperatures. As a demonstration, we repeat the analysis of the Dye-3, Greenland, temperature profile conducted by Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen (1986).
Notation
-
A Snow-accumulation rate (ms−1)
-
c Heat capacity (J kg−1 ° C−1)
-
∈ Adjustable penalty parameter (m s−1)
-
G Basal heat flux (W m−2)
-
H Ice thickness (m)
-
η Pre-conceived surface-temperature history (° C)
-
θ Observed temperature-depth profile (° C)
-
J Performance index (° C2 m)
-
k Thermal diffusivity (m2s−1)
-
k Thermal conductivity (W° C−1 m−1)
-
Elevation where vertical strain rate changes from linear function of -
z to constant value (m)
-
ρ Density (kg m−3)
-
λ Adjoint trajectory (Lagrange multiplier) (° C)
-
T Temperature (° C)
-
t Time (s)
-
tf Terminal time (present) (s)
-
T 0 Initial temperature profile (° C)
-
Ts Surface temperature (° C)
-
w Vertical velocity (m s−1)
-
z Vertical coordinate, positive downwards (m)
Concepts of Paleoclimate Inference
Reconstructing the past climate from borehole temperatures involves three conceptual problems: the forward problem, the inverse problem and the least-squares inverse problem. In the forward problem, the objective is to determine the present-day temperature profile by solving the heat equation forward in time with specified initial and boundary conditions. The forward problem is parallel to the thermal evolution that produced today’s temperature-depth profile. While forward problems are easy to solve, the particular forward problem representing Nature’s evolution through history cannot be solved. We are not privy to the boundary and initial conditions.
In the inverse problem, the objective is to determine the boundary or initial conditions by solving the heat equation backward in time starting with the observed, present-day temperature profile. To solve the inverse problem, the observed temperature profile must satisfy special mathematical conditions. Under practical circumstances, these conditions are too stringent; thus, the inverse problem may be impossible to solve.
The least-squares inverse problem has the same objective as the inverse problem, but does not impose stringent mathematical constraints on the functional representation of the temperature data. The least-squares inverse problem can be solved in all practical circumstances. This solution will be a non-unique surface-temperature history that produces a least-squares best fit to the observed temperature-depth profile.
Forward Problem
The forward problem represents the mathematical description of the heat-transfer process that takes an initial temperature field at time t = 0 and yields a terminal temperature field at time t = t f. We identify time t f with the present, and time t = 0 with a time in the past. For demonstration purposes, we adopt a one-dimensional heat-transfer equation
Notation has already been defined in a separate section, and t and z derivatives are denoted by subscripts.
If the basal boundary condition, G(t), surface boundary condition, Ts(t) and initial condition, T 0(z), are known, a solution of Equation (1) gives the terminal temperature state, T(z,tf ).
Inverse Problem
The forward problem described by Equations (1)–(4) may be written symbolically as
where R[] is a functional of a continuous function Ts(t), and T(z, t f) is the terminal temperature state. If, at time t f the temperature-depth profile is measured and is denoted by θ(z) then the unknown surface-boundarycondition T s(t) can be written in terms of the inverse of R[]:
Regardless of whether the heat-transfer operator R[] can be inverted, there are reasons to believe that Equation (6) cannot be evaluated in practice. The function of depth θ(z) is likely to be contaminated with structure that is physically incompatible with solutions of the forward problem. If, for example, θ(z) is an interpolation polynomial, it will lack the infinite differentiability associated with solutions of the forward problem. Alternatively, if θ(z) is a truncated Fourier series, it may lack terms representing wavenumbers found in the initial condition. The overriding concern is that measurement error may simply make the observed profile θ(z) too rough to be compatible with the diffusive heat-transfer process.
Least-Squares Inverse Problem
Having established that θ(z) may be incompatible with solutions of the forward problem, we resort to approximation methods. In this situation, our goal is to estimate the surface-temperature history, T s(t), which will yield a terminal temperature state, T(z,t f), that closely resembles the observed temperature-depth profile, θ(z) without necessarily being equal to θ(z) at all points. Stated differently, we wish to find a (possibly non-unique) T s(t) which will minimize the following performance index:
At this stage, the performance index J may not capture all of the possible goals of our analysis. In particular, we may choose to trade off some precision in the fit between T(z,t f) and θ(z)to impose independent constraints on T s(t). Such constraints could arise, for example, from a desire to minimize the climatic swings called on to fit borehole-temperature data. In this case, a term is added to J which accounts for deviations of T s(t) from a predefined reference temperature history η(t):
The value of ∈ in Equation (8) controls the trade-off between the fit to data and the conformance to a preconceived reference history. In practice, ∈ should be sufficiently small to avoid an inferred T s(t) which simply tracks the reference history. If ∈ is too small or zero, singular value decomposition (SVD) must be used to invert the matrix generated by the numerical implementation of the control method (see the Appendix for the circumstance when ∈ is identically zero). SVD is necessary where matrix eigenvalues are computationally indistinct from zero such as when there are limitations on the ice-sheet thermal memory (sec Reference Press, Flannery, Teukolski and WetterlingPress and others, 1989, chapter 2). A minimum-norm condition is imposed on the inferred T s(t) when SVD is used, and this is comparable to the imposition of a η(t) =0 pre-conception over ancient parts of the climate history.
Equations of Evolution as Constraints
A great difficulty in finding the surface-temperature history that minimizes J is the implicit relationship between the terminal temperature state T(z,t f) and T s(t). This crucial constraint can be conveniently enforced by using a Lagrange-multiplier function, λ(z,t). (Lagrange undetermined multipliers are commonly used in calculus of variations problems to impose constraints. One can learn about Lagrange multipliers from Reference Thacker and LongThacker and Long (1988).) We introduce Equations(1), (3) and (4) as constraints on the relation between T(z,t f) and T s(t) by adding three terms to the definition of J:
The Lagrange multiplier function, λ(z,t) is an abstract variable introduced to enforce a constraint. We shall soon see, however, that λ(z,t) has an important intuitive significance. In particular, λ(z,t) is the gradient of J ' with respect to T s(t) (Reference Thacker and LongThacker and Long, 1988; Reference WunschWunsch, 1988). This gradient proves useful when searching for the minimum of J '. If one were to guess a particular T s(t) for example, then one would use λ(z,t) to determine how best to improve this guess in order to move down the gradient towards the minimum of J '.
Another benefit provided by the introduction of the Lagrange multiplier function is that T(z,t),T s(t) and λ(z,t) can be treated as independent unknowns in the minimization of J ' The constraints represented by Equations (1)–(4) are automatically satisfied by a minimum value of J ' because
Before determining the conditions under which J ' is minimized, it is useful to perform several integrations by parts to swap partial derivatives acting on T for partial derivatives acting on λ. This step serves to establish a partial differential equation and boundary conditions for λ(z,t) The result of these integrations by parts is shown below:
Minimization and the Adjoint Trajectory
Having defined the performance index, we now identify the conditions for its minimization:
These conditions yield a system of equations which determine T s(t), T(z,t) and λ(z,t). Equations (1)–(4) belong to this system, so are not rewritten here. The additional equations are
Equations (17)–(20) are the adjoint form of Equations (1)–(4), and, as a result, the Lagrange multiplier function λ(z,t) is also called the adjoint trajectory. These adjoint equations represent a heat-transfer process which runs backward in time (the diffusive term appears with a change in sign). The starting condition for the adjoint trajectory is applied at the final time, t f; and this avoids the instability problems associated with reversal of the diffusion process. In finding the adjoint trajectory, one is able to propagate backward in time information about the mismatch between the terminal state and the borehole data. This information is then used to establish whether the gradient of J ′ with respect to T s(t) is zero or, if not, to determine how best to correct the T s(t) to move down the gradient of J ′ towards the minimum.
Demonstration Using Synthetic Data
We demonstrate the least-squares inversion procedure by generating a synthetic temperature-depth profile and then using the adjoint-trajectory technique derived above to recover the forcing. We begin by running a finite-difference model of Equations (1)– (4) forward in time to specify a synthetic θ(z) Next, we solve Equations (17)–(21) to derive a surface-temperature history using the synthetic data θ(z). Finally, we compare the estimated and the known surface-temperature histories and evaluate the adjoint-trajectory inversion algorithm.
Finite-Difference Forward Problem
We generate a synthetic temperature-depth profile by solving the following finite-difference representation of Equations (1)–(4) with a simple, arbitrary specification of boundary conditions and parameters.
In the above equations, index i refers to the spatial grid points (where i = 1 is the surface and i = 25 is the base), and index n refers to the time level (where n = 1 is the initial time and n = N is the terminal time).
In this demonstration, we disregard the thermal inertia of bedrock below the ice, and we choose simple representations of the thermal diffusivity and vertical velocity profile (constantvertical strain rate). We set k = 1.4 × 10−6m2s−1, k = 2Wm −1° C−1, ρc = 1.43 × 106J m−3° C−1, A = 0.25 m a−1 and G = 0.05 W m−2. We represent the ice column using 25 grid points, which are indexed from i = 1 (the surface) to i = 25 (the base). We make the ice 2000 m thick. For our initial condition, we use the temperature profile given by Equation (25). To generate the synthetic temperature profile, we step Equations (22)–(24) through 49 time steps of 100 a duration while applying a surface-temperature history given by
where ω = 7.96 × 10−11 s−1 represents a 2500 year climate cycle.
The final synthetic temperature-depth profile (Fig. 1) has no obvious features that indicate it was the result of an oscillatory surface-temperature history with a 2500 year period. Nevertheless, we shall see that the least-squares inversion procedure does quite well in recovering the amplitude and period of the oscillation.
Fig. 1. Initial condition and terminal state for the demonstration of adjoint trajectory methods with synthetic data. The least-squares inuerse problem is to determine the surface-temperature history that changes the temperature profile from its initial state to its terminal state in 5000 years.
Finite-Difference Formulation
Finite-difference Equations (22)–(24) may be written in compact notation to make the formulae derived below more readable.
Here A and B are 25 × 25 matrices, and C, G, T(n+1) and T(n) are column vectors of dimension 25.
where the vertical velocity at grid point i is
We now attempt to recover the surface-temperature history from the synthetic temperature profile generated in the forward problem. We accomplish this by choosing T s(n),n = 2, …,N, to minimize the following finite-difference version of the performance index given by Equation (13).
where the measured borehole temperatures are represented by
In this demonstration, we take ϵ = 10−8 and η(t) = −25° for all values of n. (When choosing ϵ, one must take into consideration the number of time steps N). Minimizing 
yields the following finite-difference equations for the adjoint trajectory λ(n):
where
Equations (36)–(38) can be used to express T s(n) in terms of T(N) – Θ and η(t) The result can be substituted into Equation (27) which is then applied recursively to express T(N) in terms of T(1), Θ, G and the η(n)′ s. This result represents a linear system of equations which can be solved for T(N).
and
Equation (40) is solved for T(N). Equation (36) is then solved by stepping backward in time using Equation (37) as the initial condition. These steps yield the adjoint trajectory, λ(n). Equation (38) is then used to determine T s(n) for all values of n between 2 and N = 50.
Result
The inferred surface-temperature history is compared with the known surface-temperature history (Equation (26)) in Figure 2. Figure 3 shows the mismatch between the observed (synthetic) temperature-depth profile and the terminal temperature state associated with the inferred history. This mismatch is no greater than 10−3° C and displays small spatial scales characteristic of the effects of spatial discretization.
Fig. 2. The surface-temperature history inferred from the adjoint trajectory method compared with the known surface-temperature history used to fabricate the observed terminal state.
Fig. 3. Mismatch between the observed terminal state synthesized from the known surface-temperature history and the terminal state derived from the inferred surface-temperature history.
The Greenland Temperature Record
To compare our method with those used previously, we repeat the analysis of the Dye-3, Greenland, borehole-temperature profile conducted by Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen (1986). In their study, a surface-temperature history was derived using sound physical intuition as the primary guide for minimizing the mismatch between the observations and a computed temperature-depth profile. Differences between the results derived here and those of Dahl-Jensen and Johnsen display the algorithmic advantages (disadvantages) gained from using control met-hods. The physical assumptions and intuition used to define the Dye-3 analysis problem are not questioned here because our comparison is intended to emphasize the algorithmic aspects of control methods.
To overcome the competing demands of high temporal resolution and long time span, Dahl-Jensen and Johnsen partitioned the surface-temperature history of Dyc-3 into long-term and short-term parts. The long-term part was estimated from geochemical properties of the Dye-3 ice core,and was intended to capture the gross features of the most recent part of the 100 000 year glacial cycle and is shown in Figure 4. (Dahl-Jensen and Johnsen also allowed the surface-accumulation rate to vary over the glacial cycle,) The temperature-depth profile that would result from the long-term part of the history alone docs not precisely fit the observations. Deviations between this reference profile and the borehole measurements, shown in Figure 5, are corrected by the short-term part of the surface-temperature history which resolves the last 25 000 years with greater detail. The short-term part thus represents a schedule of temperature corrections that must be made to the long-term history to best explain the borehole temperatures.
Fig. 4. Recent part of the long-term history of surface-temperature representing the gross features of the 100 000 year glacial cycle (Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen, 1986).
Fig. 5. Temperature deviations between the observed borehole temperature profile and a reference profile associated with a long-term climate history determined by Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen (1986). Observations of borehole temperature have been reported by Reference Gundestrup and HansenGundestrup and Hansen (1984) Also shown is the calculated temperature-depth profile associated with the inferred schedule of surface-temperature fluctuations derived from the control method.
Dahl-Jensen and Johnsen determined the short-term history of temperature corrections by a trial-and-error technique designed to correct mismatch between the reference profile and the observed profile. (They referred to a least-squares performance index similar to ours; but their method of finding its minimum differs from ours because they did not make use of the adjoint trajectory.)
We focus our analysis entirely on determining the short-term history of surface-temperature fluctuations needed to correct the long-term reference history defined by Dahl-Jensen and Johnsen. The deviations between the reference profile (the borehole temperatures computed from the long-term history) and the observed profile are defined to be θ(z) in this exercise, and are shown in Figure 5.
To interpret these deviations using control methods, we modify the finite-difference formulation discussed previously to account for several additional physical processes associated with the heat flow at Dye-3: (i) the existence of temperature variations in 3000 m of bedrock below the ice, (ii) a more complex variation of vertical velocity through the ice column, and (iii) depth-dependent thermal diffusivity (which is variable in the ice, but constant in the bedrock). Temperature variation in the bedrock is accounted for by continuing the finite-difference grid below the ice/bed interface. Since the borehole measurements, θ(z) do not extend into the bedrock, we multiply T(N) in the first term of Equation (34) by a rectangular truncation matrix so that the result may be compared directly with Θ. (This truncation matrix has N rows and M columns where N is the number of finite-difference levels within the ice and M is the number of finite-difference levels within the ice and bedrock. The elements of this matrix are all zero except for Is along the diagonal.) The ice column (2000 m) and bedrock (3000 m) are resolved by 101 finite-difference grid points with a 50 m spacing.
The vertical velocity, w(z), assumes that the ice thickness at Dye-3 has not changed, that the accumulation rate has been constant, and that the vertical strain rate is constant in the upper 1700 m of the ice column and decays linearly to zero in the lowest 300 m of the ice column. This variation of vertical strain rate is motivated by Reference ReehReeh’s (1989) analysis of the age-depth relationship for Dye-3 and by previous considerations by Reference Dansgaard and JohnsenDansgaardand Johnsen (1969). An expression for this velocity profile is written as follows:
where 
is the elevation above the ice/bedrock interface (300 m) below which the vertical strain rate is a linear function of z.
The accumulation rate is assumed constant at 75% of its present value of 0.49 m a−1. Dahl-Jensen and Johnsen suggested that A was constant at 0.245 m a−1 prior to 15 000 years BP, increased linearly with time to its present value (0.49 m a−1) between 15 000 and 10 000 years BP, and remained constant since 10 000 years BP. For simplicity (and because Dahl-Jensen and Johnsen did not consider variable A in determining the short-term history), we hold A constant. Variable accumulation rate can be accommodated through our method; but we choose not to consider such variation here because it would necessitate establishing a relation between A(t) and T s(t) which was not considered by Dahl-Jensen and Johnsen. We demonstrate the effect of our assumption in the discussion of uncertainty.
Depth-dependent thermal diffusivity is treated to allow for the weak variation of thermal conductivity and heat capacity within the ice and for the change in thermal properties across the bedrock interface. The thermal diffusivity in the ice is determined from the observed temperature profile and the sensitivity of diffusivity to temperature summarized by Reference RitzRitz (1989). The resulting variation of diffusivity with depth representing today’s conditions is assumed to hold for all time. We do not correct the diffusivity for changes in temperature at a given level over time. The thermal diffusivity of bedrock is assumed constant at 1.4 × 10−6 m2 s−1 (this is the value used by Reference Firestone, Waddington and CunninghamFirestone and others (1990)). Strictly speaking, the diffusivity of Greenland bedrock below Dye-3 is not well established. Here we use a value representative of ancient igneous rock.
We insulate the bottom of the ice and bedrock column (3000 m below the ice/bedrock interface). (The gcothermal heat flux, 38.7 m W m−2, is accounted for in the reference temperature profile associated with the long-term history, and is not needed in calculating temperature deviations about this reference profile.) We start our calculations at 30 000 years BP with the initial condition Ti (1) = 0 (no deviations from the reference profile are assumed to exist initially). We take time steps, Δt of 25 years. We take ∈ = 0.25 × 10−8 and η(n) = 0° C, based respectively on our experience and on our desire to minimize the surface-temperature fluctuations necessary to explain the data. We discuss the effects of a larger ∈ in the discussion of uncertainty.
Results for Dye-3
The inferred record of surface-temperature fluctuations needed to correct the long-term history Fig. 4 is shown in Figure 6. Mismatch between the observed deviation field (Θ) and the calculated terminal temperature state (T(N)) is shown in Figure 7. This mismatch is smaller than that reported by Dahl-Jensen and Johnsen in achieving their fit to the observed deviation field (< 0.03° C). In this sense, the control method has found a better match to the borehole observations.
Fig. 6. History of surface-temperature fluctuations inferred by control methods. The oscillations which affect the last 1500 year period are severe, and may not be entirely reliable. They can be reduced or eliminated by either increasing the spatial and temporal resolution of the finite-difference grid or by a more judicious choice of ∈ and η(t). (We recommend that ∈ be defined as a function of time, and that its value be increased over the period of the undesirable oscillations.)
Fig. 7. Mismatch between the calculated and observed temperature deviations in the borehole Fig. 5. This mismatch is lomer than that reported by Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen (1986) using trial-and-error methods to specify Ts(t).
A striking characteristic of the inferred record shown in Figure 6 is the increase in amplitude and frequency of surface-temperature fluctuations over the most recent part of the inferred history. This increase does not appear to be associated with a similar characteristic seen in the borehole observations (Fig. 5), thus the oscillations occurring during the most recent part of the history may be insignificant. To understand the significance of these oscillations and to determine how best to control them, we consider the finite-difference implementation of the control method.
Under typical circumstances, the observed and calculated borehole-temperature profiles will differ by a mismatch profile which oscillates over very short (grid point-to-grid point) vertical scales. This mismatch profile is the initial condition for the adjoint trajectory, which, in turn, determines the optimal surface-temperature history. The questionable surface-temperature oscillations seen over the recent past are thus directly related to spatial oscillations in the mismatch profile.
The backward diffusion process which governs the adjoint trajectory ensures that the short-scale (grid point-to-grid point) oscillations in the mismatch profile are damped rapidly as the adjoint trajectory evolves into the past. High-frequency oscillations in the surface temperature are thus restricted to the most recent past. This explains the appearance of the most recent part of the inferred surface-temperature record (Fig. 6). We conclude that the short-term wiggles seen in this record have little significance because they are generated by short-scale wiggles in the mismatch to borehole data which are not significant.
Three steps can be taken to control the wiggles: (1) the spatial and temporal resolution of the finite-difference grid can be increased (to restrict better the mismatch profile to smaller spatial scales which will be more strongly damped), (2) the initial condition specified for the adjoint trajectory can be spatially filtered to remove insignificant small-scale information, and (3) the trade-off parameter ∈ can be increased. In experiments where the time-step size was increased from 25 to 100 years, the short-term oscillations increased dramatically (over some parts of the recent history, the ice-sheet surface temperature would exceed the melting temperature of ice).
In experiments where e was increased (see the discussion of uncertainty following the next section), the amplitude of the oscillations decreased. This benefit was counterbalanced by: (1) diminished fit between calculated and observed borehole-temperature deviations, and (2) reduced detail in the ancient parts of the inferred climate history. (The ability to resolve the inferred Younger Dryas event, for example, disappears when e is increased above about 10−6). For this reason, it may be best to avoid using a non-zero ∈ in future implementations of the control method.
Dye-3 Climate Record
The surface-temperature history of the Holocene is determined by adding the inferred record of fluctuations (Fig. 6) to the reference climate (Fig. 4) derived by Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen (1986). This history is shown in Figure 8. From 30 000 years BP to about 10 000 years BP our inferred surface temperature shows a gradual rise from −29.6° to −22.7° C. This rise and the inferred temperature minimum at 15 000 years BP are features inherited from the long-term reference climate shown in Figure 4. Beginning at about 10 000 years BP, there is a rapid cooling event in our inferred history which suggests that the Younger Dryas affected Greenland. The timing and duration of this cold period is determined by the structure of the observed temperature deviations in the lowest 1000 m of the borehole (Fig. 5). The minimum temperature for this apparent Younger Dryas event is −24.35° C and is attained at 7900 years BP. As shown by Figure 4 and 6, the sharp onset of this event at 10 000 years BP is a result of the break in slope of the long-term reference climate. Inferred fluctuations about this reference climate (Fig. 6) are smooth and gradual during the period, so the rapid onset of cooling at 10 000 years BP in the inferred climate is not well constrained by our analysis.
Fig. 8. Inferred surface-temperature history for Dye-3, Greenland. This history is constructed by adding the infeired schedule of surface-temperature fluctuations Fig. 6 to the long-term reference climate Fig. 4. Also shown is the climate history inferred by Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen (1986) using the trial-and-error method.
After the apparent Younger Dryas event, our inferred history displays a period of rapid warming leading to a temperature maximum of −13.5° C at 4125 years BP. Following this temperature maximum, the inferred history displays stronger and more rapid oscillations. The temperature minimum (−24.1° C) at 2475 years BP is probably reliable. As discussed previously, the temperature oscillations following this minimum are insignificant.
The difference between our inferred climate history and that determined by Dahl-Jensen and Johnsen (Fig. 8) can be attributed to the more precise fit between observed and calculated borehole temperatures achieved by the control method. Much of the detail in the inferred surface-temperature history depends on relatively small improvements in the fit to data. Our results suggest that an improvement of this fit from approximately 0.03° C (the maximum difference reported by Dahl-Jensen and Johnsen) to approximately 0.01° C achieved here is necessary to bring out the detail of the paleoclimatc. The sensitive dependence of the inferred climate on the fit to data suggests that borehole-temperature measurement accuracy should be improved over the 0.01° C level commonly cited in ice-sheet research.
Uncertainty
How well do borehole data constrain the inferred surface-temperature history? To answer this question, we ad-dress three sources of uncertainty: (1) thermal diffusion causes information to dissipate with time (not all parts of the inferred history are equally well constrained by the borehole data), (2) incorrect description of the heat-transfer process, and (3) the influence of the preconception, η(t), in defining the minimum of the performance index. Uncertainty in the measured borehole temperatures constitutes another source of uncertainty in the inferred surface-temperature history, but we do not address this issue here.
One way to estimate the uncertainty imposed by thermal diffusion is to investigate the model-resolution operator (see Reference MenkeMenke, 1989, chapter 4). This operator is defined for any inverse method, and is independent of the actual data being inverted. For the discrete least-squares inverse problem in which ∈ = 0 (see Appendix), this operator is a matrix defined by
where K is the matrix defining the forward heat-transfer problem in Equation (A3) of the Appendix, and K−g is the generalized inverse of K derived from SVD.
The significance of R is readily seen when one considers the following thought experiment. If the true surface-temperature history were known and designated in its discrete form by Ts t then R would transform Ts t into the inferred surface temperature Ts inf:
The extent to which Ts t is degraded through this transformation is a measure of the best performance possible from the inverse method.
If the inversion procedure determining Ts inf from data Θ were perfect, R would be the identity matrix (is on the diagonal and zeros on all off-diagonal elements). In most circumstances, however, the non-zero matrix elements of R are spread in a band about the diagonal; and this reflects the degree to which Ts t is filtered in the determination of Ts inf.
As can be seen in Figure 9, R computed for the Dye-3 analysis presented above(with ∈ = 0) exhibits an off-diagonal spread of non-zero matrix values, and this suggests the extent to which thermal diffusion degrades the ability to resolve the past climate. Near the corner of the matrix which represents the mapping of the most recent history on to itself, R closely resembles the identity matrix. This indicates that the recent climate history will be well resolved by the inverse method. Elsewhere, the matrix is diagonally dominant, but the non-zero matrix elements are spread a substantial distance from the diagonal. This indicates that, even with the most precise borehole data, the inferred history is at best a temporally averaged representation of the true history.
Fig. 9. The model-resolution matrix R. The topography of the surface displayed on the left represents the values of the matrix elements seen as a function of position within a square dimensioned by the number of matrix rows and columns. A cross-section of this topograpliy (from A to A′) which represents the resolution of events occurring at about the time of the Younger Dry as is shown on the right. The degree to which the true surface-temperature history is degraded by the heat-diffusion process is measured by the spread of the non-zero matrix elements off the diagonal of the matrix. For the model-resolution matrix displayed here, the off-diagonal spread suggests that, if the true Younger Dry as event lasted exactly 1000 years, the inferred surface-temperature history would at best show a cool event which lasts about 5000 years.
To estimate the uncertainty caused by the second problem, inadequate physics, we ran a scries of sensitivity tests (Fig. 10) to display how the inferred history responds to possible errors in the specification of accumulation rate and vertical strain rate. To investigate accumulation rate, we performed the analysis once with A held fixed at the estimatedice-age value (50% of today’s value) and once again with A held fixed at today’s value (0.49 m a−1) (Reference Dahl-Jensen and JohnsenDahl-Jensen and Johnsen, 1986). (To save computer time, we ran all sensitivity tests with a time step of 100 years and ∈ = 10−8. The effect of this change was to increase the amplitude of short-term climate oscillations in the recent part of the inferred history.) The range of inferred climate variation associated with the changes in accumulation rate is quite striking. When A is reduced, the local temperature minimum which we associate with the Younger Dryas event virtually disappears. When A is increased, the local minimum is enhanced; but its timing shifts towards the present. We can conclude from these sensitivity tests that uncertainty in the snow-accumulation rate is an important source of error in the inferred climate history. In future studies, it may be possible to use ice-core geochemistry to pin down the variation of snow accumulation through time.
Fig. 10. Sensitivity of the inferred surface-temperature history to changes in accumulation rate and vertical strain rate. The curves labeled “ice-age accumulation rate”, “present accumulation rate”, and “constant vertical strain rate” represent, respectively, the result of reducing the snow-accumulation rate to 0.245 m a−1, the result of increasing the snow accumulation rate to 0.49 ma−1 and the result of adopting Equation (33) for the vertical velocity field. Short-term oscillations over the last 1000 years of the history are intensified because of the reduced temporal resolution (100 year time steps) adopted for the sensitivity tests.
We investigated the effect of vertical strain-rate variation with depth by changing the specification of w from Equation (43) to Equation (33). The effect of this modification is comparable to that associated with the un-certain snow-accumulation rate.
The uncertainty in snow-accumulation rate and vertical strain rate suggests that the timing and duration of the Younger Dryas cold event derived by our analysis can be questioned. Decreased snow accumulation and the variable vertical strain rate both reduce the inferred surface temperature at 10 000 years BP. These changes also cause the age of the earliest temperature minimum (identified here with the Younger Dryas event) to increase.
Uncertainty in the inferred Dye-3 temperature history may also arise from the fact that the ice column at Dye-3 is not located at the summit of the ice sheet, thus there may be importanttwo-and three-dimensional flow characteristics which affect local heat transfer. In addition, the geometry of the ice sheet may have changed significantly over the course of the last 30 000 years. Our implementation of the control method deals with one-dimensional heat transfer within a time-independent geometry. We are thus unable presently to estimate the full range of uncertainty implied by these possibilities.
To demonstrate the uncertainty associated with the choice of e, we show the inferred climate histories and fits between observed and calculated borehole temperatures associated with a variety of ∈ values (10−8,10−4,10−2) in Figure 11 and 12. Short-term climate oscillations can be made to disappear by increasing ϵ; but the costs of this benefit are: (1) less detail in the early part of the climate history, and (2) greater mismatch between calculated and observed borehole-temperature deviations.
Fig. 11. Sensitivity of the inferred surface-temperature history to changes in ∈. As ∈ is increased, short-term oscillations over the recent part of the history diminish, hut detail in the ancient part of the history is lost.
Fig. 12. Sensitivity of the calculated borehole-temperature deviations to changing ∈.
Conclusion
Control methods provide a superior alternative to traditional trial-and-error techniques for inferring surface-temperature history from borehole-temperature observations. When control methods are coupled with an adequate physical model of ice-sheet thermal evolution, the inferred climate history is more detailed and more faithful to observed borehole temperatures than is the case for trial-and-error. The enhancement of detail in the inferred climate history depends sensitively on very small improvements in the fit between calculated and observed borehole temperatures. This sensitivity suggests that it is unlikely for a trial-and-error method, which is hampered by human subjectivity, to out-perform the computational algorithm provided by control theory.
As demonstrated here, control methods are easily implemented in borehole-temperature analysis, and involve the introduction of relatively few additional mathematical concepts beyond those describing the forward problem. The principal additional feature introduced by the control method is the adjoint trajectory, which is interpreted as the measure of the sensitivity of the current borehole temperatures to the past climate. Control methods are also easily programmed using software utilities available for matrix manipulation and inversion; they also run efficiently on modern computers. Computer time required for the implementation of control methods is comparable to that required to solve the forward problem during two cycles of the trial-and-error method. We thus recommend control methods as an alternative to trial-and-error for future borehole-temperature analysis.
Our analysis of the Dye-3 borehole temperatures was intended as a tutorial. The results of this analysis, however, point to the benefits of using control methods in future borehole-temperature analysis: temperature-depth profiles can provide an independent check on surface-temperature histories derived from isotopic analysis. Whether our inferred Dye-3 temperature history supports the possibility of a Younger Dryas event in Greenland remains to be established. Our history displays a cold period between 10 000 and 7500 years BP; but this event does not possess all of the characteristics implied by the isotopic record. This disagreement suggests that Reference FairbanksFairbanks (1989) may be correct in suggesting that the oxygen-isotope profile at Dye-3 is more a reflection of the glacial meltwater input to the ocean than the local surface temperature of Greenland. Thus, to establish the validity of ice-core oxygen-isotope interpretation and to determine the areal influence of the Younger Dryas event, we recommend paleothermometry by control methods.
Acknowledgements
V. Barcilon taught us the mathematical methods used in this study and provided valuable criticism of this manuscript. We thank C. Wunsch for introducing us to control theory and for bringing to our attention its océanographie applications. S. Kakuta provided insight on the sources of uncertainty. D. Dahl-Jensen is grate-fully acknowledged for support in conducting the numerical calculations, for providing data and for pointing out the uncertainty which arises from our assumptions about the accumulation rate. B. L. Hansen and N. Gundcstrup provided the Dye-3 borehole-temperature data. II. Alley and R. Frölich are thanked for reviewing the manuscript and for suggesting improvements in the explanations. We thank G. Clow for suggesting the model-resolution matrix as a means to evaluate the uncertainty associated with our inversion procedure. This research was supported by U.S. National Science Foundation grant DPP89–14938 to the University of Chicago and U.S. NSF grant DPP86–13935 to the University of Washington. Computer support was provided by grants from the U.S. National Center for Atmospheric Research and the U.S. National Center for Supercomputing Applications in Champaign-Urbana.
APPENDIX Unconstrained Surface-Temperature History
If ∈ is chosen to be zero in the performance index 
of Equation (34), the consequence is that Equation (38), the relation between the adjoint trajectory and the surface-temperature schedule, becomes
This changes the algorithm for determining T s(n) as follows.
An expression for T(N) in terms of the initial condition and surface-temperature schedule can be derived by recursively applying Equation (27) N-1 times.
This expression is substituted into Equation (37) to provide an initial condition on the adjoint trajectory. Equation (36) is then used to determine expressions for λ(n) for n = 2, …, N These expressions, when substituted into Equation (A1), yield a series of N-1 linear equations for Ts (n), n = 2, ... , N.
where
and
and where
The (N-1) × (N-1) matrix K in Equation (A3) is entirely different from the 25 × 25 matrix K in Equation (40) because it is dimensioned by the number of time steps. Equation (A3) can be solved for the surface-temperature schedule without any further references to T(N) or λ(n).
Problems can arise in attempting to solve Equation (A3). When the number of time steps is large, or the time-step size is large, the rows of K corresponding to values of T s(t) in the distant past become arbitrarily close to zero. Physically, this represents the loss of thermal memory. The effect of this problem is to make the matrix equation too ill-conditioned to be solved by normal inversion procedures. Singular-value decomposition (SVD) should be used in this circumstance. In effect, SVD places a constraint on the surface-temperature schedule of the distant past; thus, SVD is equivalent to the control method described in the main body of this paper.
