Wirtanen

Pergamon www.elsevier.nl/locate/asr

Adv. Space Rex Vol. 23, No. 7, pp. 1333-1336, 1999 0 1999COSPAR. Published by Elsevier Science Ltd. All rights reserved Printed in Great Britain 0273-I 177/99$20.00 + 0.00 PII: SO273-1177(99)ooo46-0

A METHOD TO INVERT MUPUS TEMPERATURE RECORDINGS FOR THE SUBSURFACE TEMPERATURE FIELD OF PNVIRTANEN A. Hagermann and T. Spohn lnstitut fiir Planetologie, Wesrfiilische Wilhelms-Universitdt Wilhelm-Klemm-Stz IO, D-48149 Miinster; Germany

Miinstel;

ABSTRACT The thermal sensors on the penetrator of the MUPUS experiment package selected for the ESA Rosetta mission will enable us to determine the near-surface energy balance of the nucleus of comet P/Wirtanen by measuring the subsurface temperature profile and the thermal conductivity of the near-surface layers. Model calculations suggest that the penetrator itself will perturb the ambient temperature field such that the temperature profile will be smoothed if the thermal diffisivity of the nucleus is significantly smaller than that of the penetrator tube. It is possible, however, to calculate the undisturbed temperature profile from the data using a method based on a solution of the transient inverse heat conduction problem. Our model calculations show that a satisfactory estimate of the undisturbed temperature field can be obtained in comparatively little computing time by calculating the temperature distribution in the model volume from temperature histories at discrete points representing the penetrator temperature sensors. 01999 COSPAR.Publishedby Elsevier Science Ltd.

INTRODUCTION In situ heat-flow measurements enable us to determine the energy balance of solar system bodies and draw conclusions with regard to their origin and evolution. The only in situ heat-flow experiments carried out on a body other than Earth are the Apollo lunar heat-flow measurements in the 1970s (e.g. Langseth et al., 1970). Recently, another lunar heatflow measurement is being planned for the Japanese LUNAR-A mission (Mizutani, 1995). On a comet nucleus heat-flow measurements serve a different purpose. In most cases, heat will flow form the surface into the interior which is colder than the insolated surface. The heat flow into the interior is an important element of the surface energy balance. The energy balance determines the rate of coma gas production by sublimation. Moreover, heat flowing into the interior will cause modifications of the pristine matter of the nucleus through sublimation, recondensation and sintering. The pristine matter of the nucleus is the prime target of the Rosetta mission. A heatflow experiment such as MUPUS proposed by Spohn et al.( 1995) will therefore help cometary physicists understand the formation of the coma as well as it will help cosmochemists to assess the question how pristine the samples gathered from near- surface layers by the Rosetta lander will be. In the case of the MUPUS heat-flow experiment a penetrator (PEN) is inserted into the nucleus surface of comet P/Wirtanen at some distance from the Rosetta lander in order to avoid the thermal influence of the lander (Hagermann & Spohn, 1998). Along the penetrator tube, which consists of a fibre compound, temperature sensors are located at different depths. In passive mode these sensors measure the temperature, in active mode they are heated and supply information about thermal diffusivity and thermal conductivity. The PEN itself, however, has a significant influence on the subsurface temperature field. If its thermal diffusivity is larger than that of the ambient material, the isotherms are no longer parallel to the surface (provided the surrounding material is isotropic and only horizontally layered), but are bent such that the vertical temperature gradient is smoothed (Hagermann & Spohn, 1998). The undisturbed temperature field (represented by the temperature gradient at a considerable distance from the penetrator) can be obtained by solving the inverse heat conduction problem (IHCP), as we will show below. 1333

1334

A. Hagermann and T. Spohn

THE INVERSE APPROACH

l

In contrast to the direct problem in which the temperature field within a volume is calculated from known boundary conditions the treatment of the inverse problem involves the extrapolation of temperature measurements within the volume onto the boundaries. Hills and Mulholland (1979) give a summary of methods used, including graphical, exact analytical and numerical approaches. A problem that arises in all types of methods was already pointed out by Stolz (1960) in one of the first publications on the IHCP: “Whatever method is used, the solution will suffer from an inherent uncertainty. In a heat conduction system the effect of boundary conditions is always damped at interior points, and the inverse problem involves basically the extrapolation of the damped datum to the surface...” The method used in the present work is based on a solution of the IHCP given by Kurpisz (1991). For numerical implementation we use the method of finite control volumes. Mathematical

formulation

In the following, we model the axially symmetric experimental setup by a solid cylinder representing the PEN of radius rpelz, length .+.. and thermal diffusivity ~~~~ p laced in a homogeneous cylinder of radius rma+ > rpen and depth zrnaz > max(.+,,, zhw), where aw is the penetration depth of the diurnal heat wave in the nucleus. The heat conduction equation in cylindrical coordinates can be written as

using the dimensionless quantities p = r/rpen and C = z/rpen for radius and depth, T = t/At for time, Q = for thermal diffusivity and 0 = T/< T(T) > lz=o for temperature, respectively. The normalizer At (~AW$l Our e temperature sensors are located at 2, = (pi, &) with is the time between two temperature measurements. temperature responses Ri satisfying the conditions R(r)

=

WC,

R(k’(T)

=

dkfi(T) < k

z

i=l...C

(P, C) -+ (Pi> CJ,

r),

k=

o.

dr

vr

(2)

1,2,...

Let the model be divided into m finite control volumes. The nodes at which temperature histories are measured are identified by si E [l, 2,. . . , m] , i = 1, . .if, and temperature is discretized as O(p, I, T) := 0, (T), j = 1, . . m. For ease of formulation we define a sensor matrix S = Sij, i = 1, . . , I; j = 1,. . . , m with Sij = 0 for J’ # si and S, = 1 for j = si and a model matrix M including D = Dij, i = 1, . . . , d; j = 1,. . . , m, which is the discretized FU-IS of eq. (1) and known or estimated discretized boundary conditions B = Bij, i = 1,. . . , r; j = 1, . . . , m. M andSarebothsubmatricesofH=Hi,,i = l,..., d+r+l,j = l,..., m,i.e. D

M=

(4)

B, [

1

Note that D can be a submatrix of the corresponding direct problem if direct and inverse models have the same dimensions. The approximate solution of the inverse problem can now be written as (Kurpisz, 1991) 0,(r)

=

5

&;)

k=O

i=l

Ri’k’(~).

(5)

The q(k) are the solutions of Y

(6) where pri

= 1 for n = d + T + i and pfj

= 0 else. For k > 0 the &)

(k) = $;-” are defined recursively as pni

nd+r. If Rk(H) < m, ze solution of eq. (6) becomes problematic or impossible. by adding a priori information such that H is sufficiently well+onditioned.

for

In our case, however, B can be extended

A 9

20

Inverse

Methodto Invert MUPUS Temperature

solution,

1335

k=O

0

20 2 4 Nermallzed

6 8 radius p

2

4

6 P

8

2

4

6 P

0

2

4

6

6

P

Fig. 1: Direct (“true”) temperature field (left, PEN marked with dotted lines and sensor locations marked by small rectangles) and “evolution” of the inverted temperature field with increasing order k of the solution in normalized coordinates. The dashed vertical line in the last plot indicates the location of the estimate of the undisturbed profile shown in figure 2. Note that these results represent a subset of those shown in figure 2. Time derivatives and errors The high order derivatives in eq. (6) are extremely difficult to obtain for measurements with realistic errors taken into account. Small errors in the temperature measurements can result in substantial errors in the calculated derivatives. Higher order time derivatives can be calculated by solving the corresponding Volterra’s integral equation of the first kind with the aid of Tikhonov regularization of order 0 and 1. Although a choice of K 2 m would be ideal, we obtained the most stable results by setting K = 2, i.e. by neglecting time derivatives of higher order. Of course, the choice of a small value for K suppresses short wavelength components of the temperature field. The regularization parameters X0 and X1 acting as filters have to be chosen such that the solution is neither oversmoothed nor errors are overamplified. In the first case, the temperature profile obtained would not be significantly better than if we assumed that temperature at the PEN were undisturbed, in the second case measurement errors would dominate the results. In our test cases, a parameter ratio Xl/X0 in the order of l/100 had the desired effect. A choice of X0 << Tr(H) prevents oversmoothing, X1 (< X0 warranties that the result is continuous and errors are small. As the surface temperature signal is damped with depth we increased X0/X1 linearly with depth. To rule out very short-period oscillations due to measurement errors, an additional five-point least-squares smoothing of the time series is applied before processing.

RESULTS In order to test the feasibility of the inversion a solution of the forward problem has been calculated for several upper surface temperature histories, bottom surface temperature values and thermal diffusivities of the ice. Note that the model cylinder in the forward calculations is larger than that used in the inverse solution, which enables us to test the use of boundary condition estimates rather than exact boundary conditions. The number of nodes of the inverted field is limited by K. In our tests of the inversion algorithm we have calculated direct solutions for a number of scenarios with varying surface temperature history, initial temperature profiles and thermal parameters. For a given timestep pi we performed 100 numerical Monte Carlo experiments where temperature histories of normalized duration 50 at given nodes were recorded with Gauss-distributed errors with a standard deviation of 0.01. These randomly perturbed temperature histories were then used to perform the inversion of the temperature field. A normalized parameter value commmon to all models is C = 20 x p. We present two cases where he temperature gradient is assumed to be negative because temperature decreases with depth as the comet approaches the sun. The results presented in figures 1 and 2 correspond to a case where a triangular temperature signal with period 100 and amplitude 0.5 was applied at the surface. The timestep of the solution was chosen as pi = 70. Inside the PEN we chose (Y = 10 and outside (Y = 0.5. The number of temperature sensors was assumed to be -!J= 10. Figure 1 shows the forward temperature field on the left as a function of radius and depth calculated (PEN marked with dotted lines) and the evolution of the inverse solution with increasing k. Figure 2 shows the temperature profile from the forward solution and the results of 100 inversions at p = 0.5 and p = 8.5 which corresponds to the C-axis and the dashed line in figure 1 (right-most panel) respectively. A second model presented in figure 3 shows the results of 100 inversions with a triangular surface temperature signal of period 50, e = 5 and (Y = 1 outside the PEN. The timestep of the inversion is 7i = 100. As

D @-Profile,

=

0.5

A. Hagermann and T. Spohn

8

-P,rqfi!e?

CI-

O0

5.8::

cl

f +

50

e 0 m p 100 g 0

5

k,-10

z I3 b 15, z 0 20 I3 0.5 0.6 0.7 0.8 0.9 1.0

15

I , 1 20 I. * .P 0.5 0.6 0.7 0.8 0.9 1.0

Fig. 2: Temperiture as a function of de\th at ~0.5 (left) and p=8.5 (right). Solid lines indicate the direct (“true”) temperature profile. The results of 100 inversions are marked by dots. Sensor depths are indicated by squares on the left.

2oP. Y. I 1 0.6 0.7 0.8 0.9 1.0 1.1

20 0.6 0.7 0.8 0.9 1.0 1.1

Fig. 3: Temperatuye as a function of depth afp=O.5 (left) and p=8.5 (right) for 5 sensors, and a shorter wavelength surface temperature signal for model 2.

becomes obvious from the results, the undisturbed profile can be reconstructed errors corresponding to a standard deviation of 1%.

despite the presence of measurement

CONCLUSIONS It has been shown that the inevitable thermal influence of a penetrator in a half space can be compensated by mathematical means, namely by solving an inverse problem. Apart from the reduction of computing time, which is a significant advantage compared to forward simulations, the method also smoothes the influence of measurement errors on the inverted field. The solution, which basically consists of a series of temperature signals of increasing frequency (i.e. higher derivatives) must be truncated early to prevent divergence of the solution, which results in a loss of accuracy. However, our models suggest that this effect is only of minor importance. The results provide a theoretical background for data processing of the MUPUS heat-flow experiment on comet P/Wirtanen.

ACKNOWLEDGMENTS This research was supported by Deutsche Forschungsgemeinschaft Raumfahrt e.V. (DLR grant no. WE 150 OH 9503 7-ZA).

(DFG) and Deutsches

Zentrum

fiir Luft- und

REFERENCES Hagermannn, A. and T. Spohn, Evaluation of MUPUS data and the inverse heat conduction problem, Annal. Geophys. 16 Suppl., Cl012 (1998) Hills, G. and G. P. Mulholland, The accuracy and resolving power of one dimensional transient inverse heat conduction theory as applied to discrete and inaccurate measurements, Int. J. Heat Mass Transfer 22, pp. 1221-1229 (1979) Kurpisz, K., Numerical Solution of One Case of Inverse Heat Conduction Problems, J. Heat Transfer 113, pp. 280286 (1991) Langseth, M.G., A.E. Wechsler, E.M. Drake, G. Simmons, S.P. Clark and J.Chute, , Apollo 13 Lunar Heat Flow Experiment. Science 168, 10 (1970) Mizutani, H., Lunar Interior Exploration by Japanese Lunar Penetrator Mission, LUNAR-A, J. Phys. Earth 43, pp. 657-670 ( 1995) Spohn, T, M. Banaszkiewicz, J. Benkhoff, W.H. Ip, C. Jaupart et at., MUPUS, Proposal for a Suite of Instruments to Make Extended in situ Physical Measurements at the Surface of an Evolving Comet Nucleus (1995) Stolz, G., Numerical Solutions to an Inverse Problem of Heat Conduction for Simple Shapes, J. Heat Transfer 82, pp. 20-26 ( 1960)