## Info

ig T

Fig. 24. Cumulative probability distributions for annual mean temperature record for Northern Hemisphere presented in Figure 3 (solid line) and for Gaussian process (dashed line). T is expressed in terms of standard units: (T-mean)/s.d. (mean = — 0.260K, s.d. = 0.081 K).

Actually such extreme events occurred eight times in the near 1800-year long Mann and Jones (2003) temperature time series, thus, with much shorter average return period of near 220 years.

It should be mentioned that the soft limits result from exact calculations, and not simply from physical plausibility arguments; thus, they seem to be more appropriate than the hard limits. The simplest way to incorporate soft limits in the solution is the examination of the extreme solutions (Jackson, 1979). The maximum (or minimum) value of bTV0(t) is given by

where C = ATS—1A, V0 = CATS—1T, and ^ = [Qmax - Q0)/(bTCb)1/2. In the last expression Q0 = s TS—ls, where s = AV0 — T, and Qmax is the maximum allowable residual criterion for a model that fits the data in a satisfactory way. The value of Qmax should be chosen so that Tmax,min make values of bTV0 to fall as close as possible into the given soft limits.

Moment vector of the solution b is usually used with the unit coefficients. Some kind of weighting that is reasoned by the well-known progressive smoothing down into the past of the reconstructed climatic fluctuations (see examples presented in the next section) is also possible.

Different kinds of additional information are not equally important in all cases. Thus, they may also be weighted for specific problems.

2.3.5 Least squares inversion in functional space (FSI)

This algorithm is based on the theory by Shen and Beck (1983, 1991) and modified by Shen et al. (1995). The problem is conventionally formulated as the 1-D pure heat conduction in a laterally homogeneous subsurface with depth-dependent thermal properties (Eq. (4)), where transient component arises from the surface temperature variations. Initial and boundary conditions are the same as described above for SVD and ramp/step methods. It should be mentioned that the use of the 1-D approach is not inspired by the wish to simplify the inversion problem only. An essence of the fact that all inversion techniques are based on the 1-D equation, which obviously neglects the influences of numerous terrain effects such as medium heterogeneity, topography, and groundwater circulation, is that the 2-D inversions will not necessarily improve the GST histories. The 2-D approach will significantly raise the number of degrees of freedom of the inverse problem (underground structure, thermophysical parameters, and pattern of the steady-state temperature field), which implies the corresponding limitation of a priori parameter range treated.

As previously, the subsurface temperatures are divided into steady-state and transient components according to Eq. (8). The model is denoted as a set of parameters m = [^(z), pc(z), U0, Qm, V0(t)]. If necessary, the radioactive heat generation can be also included as the parameter into the model m. However, it is justified only for deep holes in the case of the reconstruction of remote GST changes (see Section 2.2 of this chapter). As can be seen, the thermal properties of the medium are formulated as unknown parameters and, in contrast with two previous methods, should be estimated together with the quantities describing the initial thermal state of the medium and the past climate history. In this way the uncertainties of the knowledge of the physical properties of the medium can be also taken into account. It is assumed that the deviation of the true GST history V0(t) from the a priori version V(t) is a stationary Gaussian process with an assumed exponential auto-covariogram

<[V0(t + t) - V (t + T)],[V0(t + t) - V (t + t)]> = a 2exp

where the standard deviation a and the correlation time tc characterize a priori constraint imposed on the GST history. An exact shape of the autocovariance function is not a crucial factor of inversion. Thus, in the study by Shen and Beck (1991) the authors used the Hanning window. Wang (1992) applied function a2exp(-T2/i;;) at the right side of Eq. (27). The only requirement is that this function should be sufficiently tapered at large t. The really important parameters of the autocovariogram are a and tc. The former v Tc 7

is responsible for the amplitude of the GST changes, and the latter controls the smoothness of the quantity [V0(t)_V(t)]. As in the case of SVD method, the correct choice of values for a and tc should be based on the use of the existing climatologic records.

In its most recent version the FSI can be accomplished in six general steps:

(1) Set the model m and define a priori model m0. The values and variances for K(z) and pc(z) have been given. It is assumed that they are uncorrelated. To complete m0 one needs the values for remaining parameters U0, Qm, and V0(t). Generally, these quantities are not well known. It is accepted that V0(t) takes small values close to zero, is bounded, and is rather smooth function of time t. Very large a priori standard deviations are assigned to U0 and Qm values to avoid undue bias.

(2) Solve the problems posed by the heat conduction Eq. (4) and by its steady-state analog for the steady-state U(z) and the transient components V(z,t) (Eq. (8)) of the temperature field, respectively. Calculated data can be represented as T(z,t) = Pu[U(z)] + P[V(z,t)]. This equation describes the "theory" that projects the data on the model, where P and Pu are corresponding projectors. For the discrete data this equation, similar to the SVD, can be written in the algebraic form. However, the fundamental difference between the FSI formulation and the SVD concept is that the former works with operators, while the latter deals with matrices.

(3) Calculate the data residual Sm = 1/2[T(z,t)_T(zk)]T[T(z,t)_T(zk)].

The problem of parameter estimation then becomes one of minimizing the weighted least squares misfit function (Tarantola and Valette, 1982a, b; Menke, 1989)

5(m) = 1(T _ T,)TC_j(T _ T0) + 1(m _ m0)TCm'(m_ m0) (28)

Cd and Cm are a priori covariances describing uncertainties in T(zk) and m0, respectively. The latter equation highlights the fact that for given T0 and Cd, the deviation of the output model from the a priori model depends critically on Cm. It can be illustrated with the simple example. At large Cm (no confidence in the validity of a priori model m0) the second term in the right hand of Eq. (28) will participate only insignificantly in the minimization of S(m), so that the initial model become irrelevant. In other words, the model is well resolved by the data. For a small Cm (strong validity of m0), the data will become irrelevant, so the model appears to be poorly resolved by the data. Thus, the application of the FSI strongly depends on the proper assessment of Cm.

(4) The quasi-Newton method is applied for minimizing 5(m) (Shen and Beck, 1991; Wang, 1992). This method can directly compute (approximate) a posteriori model covariance. Its computationally convenient form for jth iteration is mj+1 = m0 _ CmMj[MjCmMj]_1 X [T _ T _ Mj (mj _ m,)] (29)

where M is the derivative operator that maps the model space into the data space, defined by 5T = M5m. In a functional space formulation the operators M and its transpose Mr are not explicitly computed and stored. Only the results of their mapping are calculated. The process of M mapping is presented in detail in the Sections 6 and 8.2 by Shen and Beck (1991) and the mapping of Mr in Sections 8.4-8.6 of that work.

During the step (4) one should compute the symmetric, positive definite matrix E;. = MjCmM j + Cd and solve the vector algebraic equation E;.3(d;.) = 3(m.) for the weighted data residual ¿(d.).

(5) Compute the gradient ¿(g.) = M;.r3(d) and model correction 3(mcorr) = Cm3(g;).

The quantification of the uncertainties in the estimated model parameters can be done with the a posteriori covariance operator Cm that is given by (Tarantola and Valette, 1982a; Shen and Beck, 1991)

where denotes the derivative operator M evaluated at m = Eq. (29) corresponds to one iteration of the quasi-Newton expression (28). Generally, the diagonal of Cm is interpreted as the variance of and off-diagonal as its covariance. Because of significant computation time needed to examine the whole a posteriori covariance operator, it is common to compute only diagonal entries (the a posteriori variance of mm). They are usually called the SD ratios. The smaller the SD ratio is, the better the estimated parameter is resolved by the measured data. Physically, it means that the data contain sufficient information about this parameter. On the contrary, close to unity SD ratio hints that the data cannot resolve this parameter.

### 2.3.6 General comparison of the methods

Table 3 summarizes the main characteristics of the parametrization schemes and the data inversion for main three methods of the GST reconstruction described above. All methods are based on the 1-D theory of heat conduction and employ the least squares inversion theory. As mentioned above, the use of the 1-D approach is not certainly shortcoming of the inversion problem. The 1-D inversion techniques obviously neglect the influences of numerous 2- and/or 3-D terrain effects such as lateral medium heterogeneity, topography, and groundwater circulation. On the other hand, the development of the 2-D techniques will not necessarily improve the GST histories. The 2-D approach will significantly raise the number of degrees of freedom of the inverse problem (underground structure, thermophysical parameters, and pattern of the steady-state temperature field), while we may only handle finite amount of the measured data. The application of a 2-D approach means that we use more parameters to describe the unknowns than could be uniquely determined by the data. In practice the use of a 2-D approach implies severe limitation of a priori parameter range treated. In all inversion problems some optimal relation between resolution and variance should be established. In other words, one should

 Property Ramp/step method