Comparative analysis of the past glacier surface temperatures based on various probe-functions

There have been carried out many reconstructions of the past surface temperatures based on the glacier borehole data and the grid functions. Also there has been derived that reconstructions are unique and stable if they based on the finite segments of the Fourier series. We compare results of glacier surface temperature reconstructions by both these methods.


Introduction
The studies of the past temperatures at the Earth surface is important problem for prediction of the climate changes. The systematic instrumental temperature measurements took place no more than two centuries. Thus, indirect estimations of the past temperatures present main information on the past climate. There are two different sources of information about paleotemperatures. These are the rock and ice boreholes [1,2,3], the tree rings [4], corals and sclerosponges [5] etc. The most reliable data are kept in the measured borehole temperatures that present response to the surface temperature history. Measured temperature profiles in glacier boreholes can keep information on the past temperature changes at the surface for last centuries and more. There are several methods to inverse the measured borehole temperatures in the past surface temperature. That is why it is needed to compare these methods to understand their advatages.
The known paleotemperature reconstructions are based on the control methods [6], the Monte-Carlo method [7] and the Tikhonov's regularization method [8]. In these approaches the past temperature at the surface was looking for the grid function. The uniqueness and stability properties of the solution of the inverse problem were not established. On the other hand it was found out that the solution of inverse problem in form of the finite series of the Fourier poses is unique and stable [9]. Therefore, these three methods can be used to derive the past surface temperature in form of the finite series of the Fourier. In this paper we compare solutions that derived these three methods.

Mathematical model
The mathematical statement of the inverse problem consists of the thermal conductivity equation that takes into account the vertical advection term, the initial condition, the boundary condition at the bottom of glacier and the re-determination condition. The measured-temperature-depth profile is used as the re-determination condition, χ(z), where z is vertical coordinate. Thus, the following one-dimensional inverse problem for the heat conduction equation has to be solved to find out the past surface temperature: where H is the glacier thickness, w(z) is the advection rate of glacier layers, a 2 is the thermal diffusivity, µ(t) is temperature variations on the surface in time from the moment t = 0 (µ(0) = 0) to the time of measurements of the borehole temperature profile is the steady-state surface temperature in the past: where U s is the steady-state surface temperature in the past, k is the thermal conductivity, q is the geothermal heat flux,

The control method
The control method is an analog of the method of the Lagrange undetermined multipliers. The initial functional minimized to reconstruct the surface temperature is a functional consisting of the difference and integral determining the standard deviation of the desired surface temperature from a certain function η(t) [6]: where R{η(t)} is the solution of the direct problem specified by equations (1) without the redetermination condition for temperature variations η(t) and β is the parameter matched with the accuracy of the input data (temperatures measured in a borehole). Additional conditions determining R{η(t)} include the equation describing the heat transfer process and the initial and boundary conditions. Taking into account these conditions, the problem of the reconstruction of the surface temperature is reduced to the minimization of the following functional (Lagrange function) with the Lagrange multiplier λ(z, t): If the functional variations are equated to zero ∂J ∂V (z,t) = 0, ∂J ∂Vz(0,t) = 0, ∂J ∂V (H,t) = 0, ∂J ∂V (z,t f ) = 0, ∂J ∂µ = 0, we arrive at the following system of equations for determining the Lagrange multiplier λ(z, t): Thus, the algorithm for the inverse problem solution by the control method is of the iteration procedure where η(t) equals to µ(t) from the previous iteration step. The past surface temperature at the first step can be given by any value.
Let us write the surface temperature in the form of finite set of the Fourier series: Then the Fourier series coefficients can be found out by the followinf iteration procedure:

The Monte-Carlo method
The Monte-Carlo method [7] tests randomly chosen variations in the surface temperature using them as the input data for the direct problem and taking into account the consistency degree between the calculated and measured temperature profiles. Any point of the model space can be taken as the start point for the random walk m 1 cur . The random walk can be described as follows: (1) the trial value of the surface temperature m n+1 trail is taken near the current surface temperature m n cur ; (2) m n+1 trail is tested with the probability P accept = min(1, L trail /L cur ); trail is accepted, then m n+1 cur = m n+1 trail , otherwise, m n+1 cur = m n cur . Thus, the necessary number of allowable model parameters the differences for which are smaller than a certain value S 0 (which is usually matched with the error of the borehole temperature measurements) can be collected. Using this statistics, histograms can be plotted for each parameter in the vector of model parameters. Such distributions usually have the maximum region, i.e., the most probable values of the corresponding parameters, which constitute the solution of the inverse problem.

The Tikhonov regularization method
The Tikhonov regularization method is the determination of the boundary temperature µ(t) minimizing a smoothing functional consisting of the difference and stabilizer [8]: where α is the regularization parameter matched with the accuracy of the input data. The functional Ω{µ(t)} is called the stabilizing functional or stabilizer where r is the stabilizer order, q j ≥ 0, and q r > 0. The procedure of the minimization of the smoothing functional Ψ can be performed by means of the gradient method and is an iteration procedure. The iteration procedure is carried out until the functional Ψ reaches the minimum with a given accuracy, which corresponds to the optimal solution of the inverse problem. We consider case when the surface temperature is in form of the finite set of the Fourier series (7). The initial Fourier coefficients are given in the first iteration step while the next n-th iterations are determinated by the following equations: Here, the profiles W a 0 (z), W am (z), and W bm (z) are the solutions of the direct problem with the boundary conditions on the surface µ(t) = 1/2, µ(t) = cos(2πt/T m ) and µ(t) = sin(2πt/T m ), respectively. It is easy to show that the term with the stabilizer in (9) when the boundary condition on the surface µ(t) has the form of the segment of the trigonometric Fourier series has the form: αΩ = α 0 Thus, to determine the Fourier coefficients of the boundary condition given by (7), the iteration procedure specified by (11) is performed with the use of (12) and (13).

Comparative analysis of the inverse problem solutions
We carry out comparison for the test problems for the following parameters of glacier: the thickness of glacier is H = 500 m, a 2 = 10 −6 m 2 /s, the advection rate is linear function decreasing from 0.3 m/year at the surface to 0 at the glacier bottom. The direct problems are solved for two probe functions: µ(t) = 2 sin(2πt/t f ) and µ(t) = 2 sin(πt/t f ). We consider the time interval from 0 to t f = 400 years. The solutions of these two direct problems are the corresponding temperature-depth profiles V (z, t f ). After that we solve the inverse problems (1) with the determined temperature-depth profiles θ(z)= V (z, t f ). Moreover, to check the stability of the methods we introduce the random errors in the input data with amplitude of 0.05 • C. The reconstructions are calculated by three methods using the surface temperature in the form of both the grid fuctions ( Fig. 1) and the finite series of the Fourier (Fig. 2). Comparison shows that all methods provide similar results in case of the exact input data. However, the errors in the input data result in the instability of the inverse problem solutions. The worst result takes place for the grid function solution obtained by the control method. The best results is provided by the Tikhonov regularization method (Fig. 1).
It is interesting to note that if we look for the inverse problem solution by the finite series of the Fourier then three considered methods are of similar solutions (Fig. 2).