Inverse radiation modeling of Titan's atmosphere to assimilate solar aureole imager data of the Huygens probe

Available online at www.sciencedirect.com

Planetary and Space Science 51 (2003) 147 – 158 www.elsevier.com/locate/pss

Inverse radiation modeling of Titan’s atmosphere to assimilate solar aureole imager data of the Huygens probe B. Griegera;∗ , M.T. Lemmonb , W.J. Markiewicza , H.U. Kellera a Max-Planck-Institut b Physics

fur Aeronomie, Max Planck Strasse 2, Katlenburg-Lindau D-37191, Germany Department, Texas A & M University, College Station, TX, USA

Abstract During the descent of the Huygens probe through Titan’s atmosphere in January 2005, the Descent Imager/Spectral Radiometer (DISR) will perform upward and downward looking measurements at various spectral ranges and spatial resolutions. This internal radiation density could be estimated by radiative transfer calculations for Titan’s atmosphere. However, to do this, the optical properties—i.e. volume extinction coe8cient, single scattering albedo and scattering phase function—have to be prescribed at every altitude, and these are apriori not known. Herein, an inverse approach is investigated, which retrieves the single scattering albedo and the phase function of the aerosols from DISR observations. The method uses data from a DISR subinstrument, the Solar Aureole imager (SA), to estimate the optical properties of the atmospheric layer between two successive observation altitudes. A unique solution for one layer can in principle be calculated directly from a linear system of equations, but due to the sparseness of the data and the unavoidable noise in the measurements, the inverse problem is ill-posed. The problem is stabilized by the regularization method requiring smoothness of the resultant solution. A consistent set of solutions for all layers is obtained by iterating several times downward and upward through the layers. The method is tested in a simulated radiation density scenario for Titan, which is based on a microphysical aerosol model for the haze layer. Within this scenario, the expected coverage of SA data allows a reconstruction of the angular dependence of the scattering phase function with an explained variance better than 90%. ? 2003 Elsevier Science Ltd. All rights reserved. PACS: 96.30.M; 96.30; 95.55.P; 42.68.A Keywords: Titan; Huygens

1. Introduction Saturn’s satellite Titan is the only moon in the solar system with a considerable atmosphere. At the surface, the pressure is 1.5 times the pressure of Earth’s atmosphere, with the main constituent also being nitrogen. Early ground-based observations had already revealed the presence of an atmosphere containing substantial amounts of methane (Kuiper, 1944). The photolysis of methane was expected to produce solid organic aerosols (Khare and Sagan, 1973). This was conArmed by the Voyager observations in 1980 which showed an optically thick haze completely obscuring the ∗

Corresponding author. Tel.: +49-2561-9595-6137. E-mail addresses: [email protected] (B. Grieger), [email protected] (M.T. Lemmon), [email protected] (W.J. Markiewicz), [email protected] (H.U. Keller). URL: http://www.huygens.info

surface (Rages et al., 1983; Hunton et al., 1984). The known properties of Titan in general and of the haze in the atmosphere have recently been reviewed by Taylor and Coustenis (1998) and McKay et al. (2001), respectively. In the wake of the Voyager 1 Gyby of Titan, several initiatives were started to promote a mission to Titan. In 1997, the Cassini/Huygens mission was launched, which will reach the Saturnian system in July 2004. The Cassini spacecraft will deliver the Huygens probe to Titan and then explore the Saturnian system for 4 –10 years. The descent of Huygens through Titan’s atmosphere was originally scheduled for November 2004, but due to a communication link problem, Cassini’s trajectory for the probe delivery has been changed, resulting in a delay of Huygens’ descent to January 2005. The Huygens payload comprises six instruments (Lebreton and Matson, 1997), Ave of which will carry out measurements during the descent through the atmosphere (starting at about 160 km altitude). One of these instruments

0032-0633/03/$ - see front matter ? 2003 Elsevier Science Ltd. All rights reserved. PII: S 0 0 3 2 - 0 6 3 3 ( 0 2 ) 0 0 1 4 2 - 3

148

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

is the Descent Imager/Spectral Radiometer (DISR, Tomasko et al., 1997), which will take images and spectral measurements in diLerent directions at various spectral ranges and spatial resolutions. The detectors of the DISR are a 512 × 254 CCD for the visible, two linear arrays of 132 photodiodes each for the infrared and two photodiodes for the violet spectral range. Light from the foreoptics of three surface imagers (high resolution, medium resolution and side looking), two visible spectrometers (upward and downward looking) and the solar aureole imager (SA) is conducted by Aber optics to the CCD. In the present paper, we investigate an inverse approach to reconstruct the scattering phase functions of Titan’s haze from the SA data. Optical properties of aerosols have been derived from observations inside a planetary atmosphere for Mars (e.g. Markiewicz et al., 1999). But as only measurements from one altitude (the surface) were available, no altitude dependence could be retrieved. Here we present a method to reconstruct the phase functions of scattering particles at diLerent altitudes—i.e. in the layers between successive observations—from SA data. For one layer, this is an ill-posed linear problem. Such problems frequently occur in image restoration but also in many other astronomical applications, cf. Craig and Brown (1986), e.g. Grieger et al. (1991). The SA will take images of the solar aureole in two wavelengths—500 ± 25 and 935 ± 35 nm—and two polarizations. We only consider the 500 nm SA images, because in this wavelength range the optical properties of the atmosphere are dominated by the haze. Scattering and absorption by methane is not very well known for Titan’s pressure and temperature conditions, but in the visible, gas opacities are negligible (McKay et al., 2001). Therefore, in our radiative transfer calculations only the scattering properties of the aerosols are taken into account. While this is a reasonable approximation for the blue SA channel at 500 nm, the red channel at 935 nm should be studied with a more comprehensive radiative transfer model. The phase functions of the scattering particles obtained from Huygens DISR data can be used to constrain the size and shape of the aerosols present in Titan’s atmosphere. For this purpose, it would be desirable to extend the method presented herein to take polarization into account and to make use of the SA measurements at diLerent polarizations. Together with the phase functions, the radiation density Aeld in the atmosphere is estimated. The diLuse radiation can inGuence measurements of the spectrophotometry of the surface (as carried out by the downward-looking visible spectrometer of the DISR) by signiAcantly altering the apparent spectrum of the surface. This is an important eLect on Mars (Thomas et al., 1999) and can be expected to be even more important below the optically much thicker haze layer of Titan. Accurate models of the sky brightness are required to separate the atmospheric illumination from the direct one, and hence to comprehend the composition of the surface.

The structure of the paper is as follows: In Section 2 we describe the microphysical model of the aerosols composing the haze layer. The results of this model provide the scenario for testing the inverse method, thus based on these results we compute the internal radiative Aeld (Section 3.1). This Aeld is than sampled to simulate the expected SA data (Section 3.2). From the simulated data the phase functions are reconstructed with the inverse method (Section 4). The results are presented and compared with the original scenario in Section 5. The Anal Section 6 contains discussion and conclusions with respect to the validity of the used procedure. 2. The aerosol model 2.1. Light scattering by fractal aggregates West and Smith (1991) showed that aggregates of smaller spheres could be both highly polarizing and very forward scattering, as required by spacecraft observations of Titan (Tomasko and Smith, 1982; Rages et al., 1983). Microphysical models by Cabane et al. (1993) showed that such aggregates would be expected to form in Titan’s upper atmosphere. The aggregates are frequently termed “fractal aggregates” because their internal structure is characterized by a fractal dimension, D. Cabane et al. (1993) assumed a loose aggregation (D = 2, corresponding to cluster–cluster aggregates) based on microphysical arguments, while West and Smith (1991) modeled a more tightly bound diLusion-limited aggregate (D = 2:5). A sphere, like a close-packed aggregation, has D = 3. Rannou et al. (1995, 1999) showed that aggregates with a fractal dimension near D = 2 provided the best At to Titan’s geometric albedo and to the intensity proAle near the limb in high-phase angle images. To model the optical properties of Titan’s aerosols, we use the aggregate scattering model of Lemmon (1994). In this model, light scattering by individual particles is modeled by replacing the aggregate with a set of dipoles; this is known as the discrete dipole approximation (DDA, Draine, 1988), and is used in a way very similar to West (1991) and West and Smith (1991). Lemmon (1994) modeled Titan’s polarization and geometric albedo to show that the bulk of Titan’s haze could be made up of aggregates of ≈ 16 monomers with a monomer radius between 0.060 and 0:070 m, assuming a fractal dimension of 2. A more heterogeneous distribution with large particles above smaller particles was not considered. For the purpose of investigating the vertical distribution of the aerosols, it is assumed that the monomer size is Rmon = 0:06 m. DDA scattering calculations were performed for aggregates comprising up to 128 monomers (individual particles making up the aggregate) and for wavelengths from 0.2 to 2:2 m. Properties of larger aggregates were estimated based on a simple extrapolation from the set of calculations, and only made up

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

2.2. Microphysical modeling of the distribution The vertical distribution of the aerosols that is used here, taken from Lemmon (1994), was determined with a one-dimensional cloud physics model. A microphysical model following the algorithm of Toon et al. (1980) was solved for an aerosol created at high altitude that coagulates and falls with time. The model is solved for a series of 47 altitude bins from the surface to 400 km and for a series of size bins. Each size bin contains particles with double the mass of the previous size bin. Size bins below the monomer radius (0:06 m) are treated as spheres; size bins larger are treated as fractals with N = 2; 4; 8; 16; : : : monomers. Important changes to the microphysics come from the assumption of a fractal aggregate structure. For particles smaller than the monomer radius, the microphysics is unchanged. For aggregate aerosols, there are several changes. The fall velocity is reduced due to the larger cross-section of the diLuse aggregates (following Meakin et al., 1989). A key impact of this is that the fall velocity increases only very slowly with number of monomers, as the cross-sectional area of the aggregate increases nearly as fast as the mass of the aggregate. The rate of coagulation is also signiAcantly inGuenced (Meakin et al., 1985). A more thorough review of the aggregate microphysics as applied here appears in Cabane et al. (1993). The one-dimensional model is a clear limitation. Titan’s haze varies between the north and south, and the asymmetry changes seasonally (Sromovsky et al., 1981; Lorenz et al., 1999). However, for the purpose of determining a test model, this approach is satisfactory and it has the advantage over a more complicated model that only a few free parameters need to be constrained. Attempts to simulate the seasonality of the haze distribution within a general circulation model have only partly been successful (Tokano et al., 1999). 2.3. How the model was constrained The Lemmon (1994) model was determined by iterating a microphysical model run, a radiative transfer model run, and comparing the results to certain data sets. The primary data constraining the model were Titan’s high-phase angle scattering (Rages et al., 1983); Titan’s polarization at low and medium-phase angles in the ultraviolet blue, red, and near-infrared (Tomasko and Smith, 1982; West et al., 1983); and Titan’s geometric albedo (NeL et al., 1984; Courtin et al., 1991). All radiative transfer calculations were done accounting for the eLect of polarization: Titan’s albedo is

208-218 km 168-176 km 128-136 km 88-96 km

0.003

Mass density [mic-g/m^3]

a small amount of the total scattering at high levels (where the phase function is constrained by data). The aggregates were assumed to be made up of Titan tholin, with optical properties as determined by Khare et al. (1984), except that the imaginary index of refraction was multiplied by 4=3 at all wavelengths.

149

0.002

0.001

0.000

0.1

1.0

Effective Radius [mic-m] Fig. 1. Mass distribution of particles for four diLerent altitudes. The dots on the curves correspond to the particle size bins for which mass densities have been modeled.

increased by ≈ 10% due to the large single scattering polarization of the aerosols. Parameters constrained by data from high altitudes were constrained Arst. These parameters were the pressure level where the aerosols form, the eddy diLusion multiplier (eddy diLusion is a reduced version of the Yung et al., 1984, diffusion proAle), and the aerosol charging. The high-phase angle scattering and the short-wavelength polarization and albedo are sensitive primarily to the top few tenths of an optical depth of the haze. A broad family of solutions was determined; albedo within near-infrared methane bands was then calculated. These bands probe more deeply into the atmosphere, but do not probe to the surface. With the additional constraints, the aerosol mass production rate was also constrained. The weaker methane bands are too dark to have large amounts of haze in the lower stratosphere. Because of this a “rain-out” was implemented at 88 km altitude, with a much-reduced scattering opacity below this level. 2.4. The used atmosphere scenario The result of the model is a distribution of the number of aerosols in each size and altitude bin. The distribution of mass is shown for some altitudes in Fig. 1. The size bins shown correspond to aggregates of Nmon = 2; 4; 8; : : : ; 215 monomers. The eLective radius of these aggregates is ReL = 1=3 Rmon Nmon , with Rmon = 0:06 m being the radius of one monomer. At 300 km and above, spherical particles with ReL ¡ 0:06 m dominate, but the total amount of mass is small. The penetration of sunlight into the atmosphere, based on this model, is shown in Tomasko et al. (1997) and McKay et al. (2001). The cumulative optical depth as well as the single scattering albedo corresponding to the modeled aerosols are

150

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

function

300

P(’) = fPg1 (’) + (1 − f)Pg2 (’);

Cum. opt. depth Single scatt. alb.

250

(1)

Altitude [km]

where 200

Pg (’) =

150

1 − g2 (1 + g2 − 2g cos ’)3=2

(2)

with f=0:763; g1 =0:62; g2 =−0:294, cf. Figs. 2 and 3. The location and thickness of any clouds is not well constrained by the existing observations from outside the atmosphere. Moist convection can occur under conditions plausible for the near surface of Titan, but such plumes are expected to occupy only a small fraction of the surface area (Awal and Lunine, 1994).

100 50 0 0.0

0.5

1.0

1.5

2.0

Albedo / Optical depth Fig. 2. The cumulative optical depth and single scattering albedo of the model atmosphere as a function of altitude ( = 500 nm).

100.0 Phase function at 156-160 km Phase function below 88 km

Phase function

10.0

1.0

0.1 0

30

60 90 120 Scattering angle [deg]

150

180

Fig. 3. Two examples of phase functions in the model atmosphere. Above the assumed rain-out altitude of 88 km, the phase functions are determined by the modeled aerosols and vary with altitude, but they are all quite similar to the one shown for 156 –160 km altitude. Below the rain-out altitude, condensation clouds are represented by a Henyey–Greenstein phase function.

plotted versus altitude in Fig. 2 for a wavelength =500 nm. Fig. 3 shows the scattering phase function near 160 km altitude, where the probe will begin taking data, and below the rain-out, where the atmosphere is assumed to be free of aerosols. Observations indicate the presence of reGective methane condensation clouds in the lower troposphere (Gri8th et al., 1998). Such clouds are represented in our model by an optical depth of 0.4 uniformly distributed between 88 km altitude and the surface with a single scattering albedo of one and a double Henyey–Greenstein phase

3. Simulating solar aureole imager observations 3.1. Modeling the radiation density
B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

151

The downward radiation at the bottom of layer i is given by I (i) (j↓ ) = I (i−1) (j↓ )(1 − ai (|j↓ |)) (i) F ai (|j↓ |)wi P (i) (|j↓ −  |) 4 1
(i−1) I (k↓ ) + ai (|j↓ |)wi N=2

+

k↓

× P (i) (|j↓ − k↓ |) + ai (|j↓ |) × wi

1
(i) I (k↑ )P (i) (|j↓ − k↑ |); N=2

(3)

k↑

where ai () = 1 − e

−i Rzi =cos 

(4)

is the opacity of layer i. The four terms of the sum on the right-hand side of Eq. (3) represent: (1) the downward radiation at the top of the layer attenuated by passing through the layer, (2) the direct solar beam scattered in the layer, (3) the downward radiation at the top of the layer scattered in the layer and (4) the upward radiation at the bottom of the layer scattered in the layer, respectively. Analogously, the upward radiation at the top of layer i is given by I (i−1) (j↑ ) = I (i) (j↑ )(1 − ai (|j↑ |)) (i) F ai (|j↑ |)wi P (i) (|j↑ −  |) 4 1
(i) I (k↑ ) + ai (|j↑ |)wi N=2

+

k↑

× P (i) (|j↑ − k↑ |) + ai (|j↑ |) × wi

1
(i−1) I (k↓ )P (i) (|j↑ − k↓ |): N=2 k↓

(5) This formulation of the radiative transfer problem takes pattern from the matrix operator theory (Plass et al., 1973). Starting from the known downward radiation I (0) (j↓ )=0 at the top of the atmosphere (TOA) and in a Arst iteration downward through the layers assuming all upward radiations I (i) (j↑ ) = 0, all downward radiations I (i) (j↓ ) can be calculated from Eq. (3). The downward radiation at the surface is reGected by multiplication with a surface albedo of 0.1 to give the Arst estimate of the upward radiation at

Fig. 4. Location of the two frames of the SA taken during one observation cycle, one towards the Sun (black) and one in the opposite direction (white), illustrated looking upward—i.e. with the zenith being in the center and a 180◦ Aeld of view—in the hemisphere of downward radiation, as simulated for 112 km altitude. The units correspond to a solar Gux at TOA normalized to 4.

the surface. This implies the assumption of Lambertian reGection. The true albedo and reGection characteristic of Titan’s surface will be known from the downward looking instruments of the DISR and can be taken into account in the inverse procedure described in Section 4.3. Then in an upward iteration through the layers all upward radiations I (i) (j↑ ) can be calculated from Eq. (5). By iterating in this way downward and upward through the layers, convergence is achieved after several times (typically 10 iterations), yielding a consistent solution for the radiation density Aeld throughout the atmosphere. The results of these radiative transfer calculations for a solar zenith angle | |= 50◦ —which is the nominal value for Huygens’ descent— are illustrated showing the intensities of downward radiation for one altitude in Fig. 4. Further results are shown below in the left hand images of Figs. 10–14. The units correspond to a solar Gux at TOA normalized to 4. With the solar Gux being F (500 nm)=13:8 W m−2 m−1 , the normalized units used in the images correspond to 1:1 W m−2 sr −1 m−1 . 3.2. Sampling simulated observations One frame of the SA (Tomasko et al., 1997) occupies a 6 × 50 pixel area of the DISR’s CCD. It nominally maps to a 6◦ × 50◦ Aeld of view (FOV) in azimuth and zenith angle, respectively, which is vertically centered around the nominal solar zenith angle of 50◦ . This nominal geometry is used here. The actual geometry of the onboard instrument has been measured and its slight deviation from the nominal one can easily be taken into account. In an observation cycle, two frames are taken during one rotation of the probe, one towards the Sun (but slightly oL

152

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

160 km altitude 112 km altitude 0 km altitude

log intensity

1e+01

1e+00

Fig. 6. Illustration of the downward and upward radiation Aeld at the altitudes of two successive observation cycles. The scattering properties—volume extinction coe8cient, single scattering albedo and phase function—are assumed to be constant in the atmospheric layer between the measurements.

1e-01

0

30

60 90 120 150 Angular solar distance [deg]

180

Fig. 5. Simulated SA measurements for three altitudes. Observations will start at 160 km. The maximum intensity modeled is at 112 km. The units correspond to a solar Gux at TOA normalized to 4.

the Sun to avoid the direct solar beam) and the other in the direction opposite to the Sun. These FOVs are illustrated in Fig. 4. It is planned to transmit the SA frame taken towards the Sun completely—i.e. all 6 × 50 pixels—while the frame opposite to the Sun is averaged over six pixels before transmission. Herein, we assume that both frames are averaged in this way, i.e. together they provide intensity measurements at 100 locations of the sky. To simulate these measurements, the modeled radiation intensities are mapped to the pixels of the SA frames, a random noise error is added, and then the frames are average over six pixels. The noise error of the DISR’s CCD is lower than RN =30 photons. With E =hc==3:98×10−19 J being the energy of one photon with wavelength =500 nm, A=3:91×10−10 m2 being the area of one pixel of the CCD, R=3:05×10−4 sr being the 1◦ × 1◦ FOV of the pixel, R = 0:05 m being the bandwidth of the Alter and Rt = 0:14 s being the integration time, the intensity error is RI =

RNE W : = 0:014 2 ARRRt m sr m

(6)

Although the integration time that will actually be used may be considerably longer for the frames taken in the direction opposite to the Sun, here a random noise error corresponding to Eq. (6) is added to all the simulated intensity values. The resultant simulated observations are illustrated for some altitudes in Fig. 5. 4. The inverse approach 4.1. Radiation intensity measurements and phase function In order to reconstruct the phase function of the scattering particles from DISR measurements, we assume the optical

properties to be constant between the altitudes of two successive observation cycles, cf. Fig. 6. The thickness of such a layer will typically be about 4 km. The downward radiation at the bottom of the layer coming from a direction j↓ , I (i) (j↓ ), can be expressed by Eq. (3). While in Section 3.1 this equation referred to a very thin layer, we are now considering much thicker layers between two observations. We still consider only single scattering in this layer. This will introduce some error in our inverse reconstruction. One purpose of the work presented herein is to test if it is nevertheless possible to accurately reconstruct the phase function of the scattering particles. From the SA we obtain discrete sparse measurements of the downward radiation, Iˆl(i) , with i indexing the observation altitudes and l = 1; : : : ; 100 indexing the viewing directions l of the SA, cf. Section 3.2. Introducing RIˆl(i) := Iˆl(i) − Iˆl(i−1)

(1 − ai (|l |));

P˜ (i) := wi P (i)

(7) (8)

and normalizing the solar Gux at TOA F to 4, from Eq. (3) we see that RIˆl(i) = ai (|l |)P˜ (i) (|l −  |) + ai (|l |)

1
(i−1) I (k↓ )P˜ (i) (|l − k↓ |) N=2 k↓

+ ai (|l |)

1
(i) I (k↑ )P˜ (i) (|l − k↑ |): N=2

(9)

k↑

If the optical properties in the layer between two altitudes of SA observation were known and the downward radiation at the top and the upward radiation at the bottom of the layer were given, the intensity diLerences of the two altitudes as measured by the SA could be calculated from Eq. (9). However, we are here concerned with the inverse problem of reconstructing the phase function from SA measurements.

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

4.2. The inverse problem for one layer

matrixes,

Let the phase function (or rather the phase function multiplied by the single scattering albedo) P˜ (i) in layer i be represented by a polygonian (’k ; P˜ k(i) )k=1; :::; m , thus for an arbitrary scattering angle ’ with ’k−1 ¡ ’ ¡ ’k it is

(10)

Then the three terms on the right-hand side of Eq. (9) can all be expressed by multiplication of the vector P˜ k(i) representing the phase function with a corresponding matrix, i.e. RIˆl(i) =

m 
k=1

 Klk + Klk↓ (I (i−1) ) + Klk↑ (I (i) ) P˜ k(i) ;

(11)

where matrixes K ↓ and K ↑ depend on the downward radiation at the top and the upward radiation at the bottom, respectively, of layer i. For convenience deAning the sum of the three matrices K (i) := K  + K ↓ (I (i−1) ) + K ↑ (I (i) ):

(12)

Eq. (11) can be written RIˆ(i) = K (i) P˜ (i) :

(13)

This is a linear equation system for the vector P˜ (i) representing the unknown phase function, but it is ill posed. The available observations RIˆ(i) do not contain all the information necessary to uniquely determine the unknown vector P˜ (i) , implying that the matrix K (i) does not have full rank. To stabilize the solution we make use of the regularization method (Franklin, 1970; Turchin et al., 1971; Craig and Brown, 1986), introducing a priori constraints to minimize the norm of the Arst and the second derivative—i.e. the noise—in the reconstructed phase function. We then obtain a minimization problem RIˆ(i) − K (i) P˜ (i) 2 + 1

m
k2

(i) 2 (P˜ k(i) − P˜ k−1 )

m−1
(i) (i) 2 + 2 (P˜ k−1 − 2P˜ k(i) + P˜ k+1 ) k=2

(i) (i) 2 + (2P˜ m−1 − 2P˜ m )

(15)

This linear equation system has full rank and can uniquely be solved for P˜ (i) .

To solve Eq. (15) for P˜ (i) , the matrix K (i) has to be known But it depends on the downward radiation at the top and the upward radiation at the bottom of layer i, I (i−1) (j↓ ) and I (i) (j↑ ), which are also (at Arst) unknown. A consistent solution for the phase functions and the complete radiation density Aelds can be obtained iteratively. Similar to the radiation modeling with known phase function described in Section 3.1, we start assuming all upward radiations I (i) (j↑ ) = 0. From the known downward radiation at TOA, I (0) (i↓ ) = 0, we proceed in the following way: (1) Propagate downward for layer i = 1 → n: (a) Inversely estimate the vector P˜ (i) from Eq. (15), (b) calculate radiation I (i) (j↓ ) for all j↓ at bottom of layer from Eq. (3). (2) ReGect at surface to get upward radiation I (n) (j↑ ). (3) Propagate upward for layer i = n → 1: (a) Calculate radiation I (i−1) (j↑ ) for all j↑ at top of layer from Eq. (5). (4) Repeat from (1) until convergence (which is usually achieved after about 10 iterations). The total number of layers considered for the inverse reconstruction of phase functions (which equals the number of assumed observation altitudes) is n = 23. The uppermost layer extends from TOA at 300 km down to the altitude of the Arst observation at 160 km. Between 160 and 88 km, the altitude diLerence between two successive observations was assumed to be 4 km. Below 88 km, where rain-out yields lower optical thickness, less observations were used corresponding to a layer thickness of 22 km. This increases the diLerence RIˆl(i) , cf. Eq. (7), and thus its signal-to-noise ratio, which enters into the calculations. 5. Results

 = min;

(K (i)T K (i) + 1 H1 + 2 H2 )P˜ (i) = K (i)T RIˆ(i) :

4.3. Iterative assimilation for all layers

’ − ’k−1 (i) (i) P˜ (i) (’) = P˜ k−1 + (P˜ k(i) − P˜ k−1 ) ’k − ’k−1 ’ − ’k−1 ˜ (i) ’k − ’ ˜ (i) P + P : = ’k − ’k−1 k ’k − ’k−1 k−1

153

(14)

where 1 weights the norm of the Arst derivative and 2 the (i) norm of the second derivative. The additional term (2P˜ m−1 − (i) 2 ˜ 2P m ) ensures smoothness at the total back scattering direction of ’m = 180◦ . Derivation of Eq. (14) with respect to P˜ k(i) for all k yields, with H1 , H2 being the respective smoothing

For the results presented herein, the smoothing parameters, cf. Eq. (15), were chosen according to 1 =2 =5×10−4 . The only exception to this was the very thick top layer between TOA at 300 km and the altitude of the Arst observation at 160 km where 1 = 2 = 5 × 10−3 . We also tested the L-curve method (Hansen and O’Leary, 1993; Reginska, 1996) to choose automatically an individual smoothing parameter for each layer, but the very noisy resultant phase functions indicate that the method picks much too small values. The reason is probably that the L-curve method— like all approaches to the objective choice of the smoothing

154

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

50

10

45

True at 156-160 km Reconstructed at 156-160 km Reconstructed at 88-92 km

40

True at 0-88 km Reconstructed at 66-88 km Reconstructed at 0-22 km

8

35 6

25

w*P

w*P

30

20 15

4 2

10 5

0

0 (a) -5

-2

0

30

60 90 120 Scattering angle

150

180

8 Fig. 8. Comparison of true and reconstructed phase function multiplied by single scattering albedo for two layers below the haze, where a Henyey–Greenstein phase function was assumed, cf. Eqs. (1) and (2).

w*P

6 4 2

140 120

(b)

0

30

60 90 120 Scattering angle [deg]

150

180

Fig. 7. Comparison of true and reconstructed fractal aggregate phase functions within the haze layer multiplied by single scattering albedo, shown completely (upper panel) and on an expanded scale (lower panel). The true curve for 88–92 km would be almost indistinguishable from the one for 156 –160 km.

parameter—has to assume that the convolution matrix K (i) , cf. Eq. (13), is known exactly. But in our case this matrix depends on the estimated radiation density Aeld at the top and the bottom of layer i, which may deviate to some extent from the true radiation. Neglecting this additional source of errors besides the noise in the observations RIˆ(i) , the L-curve method chooses the smoothing parameter too small in order to reduce the residual RIˆ(i) − K (i) P˜ (i) 2 in Eq. (14) more than achievable with a reasonable result for P˜ (i) . Examples of reconstructed phase functions are presented in Figs. 7 and 8. The layers shown in Fig. 7 correspond to the layer between the two highest altitudes of simulated observations (156 and 160 km) and to the layer between two observations near the bottom of the haze layer (88 and 92 km altitude). The layers shown in Fig. 8 are located below the haze, where the used true phase function is the same Henyey–Greenstein function from 88 km altitude down to the surface, cf. Section 2.4, Eqs. (1) and (2). The accuracy of the reconstructions for all layers is shown in Fig. 9 in terms of the explained variance. If P˜ (i)∗ is the

Altitude [km]

0 -2

Phase function Radiation density field

100 80 60 40 20 0 0.0

0.2

0.4 0.6 Explained variance

0.8

1.0

Fig. 9. Accuracy of the reconstructed phase functions and the corresponding radiation density Aeld, given in terms of the explained variance, cf. Eq. (16).

vector representing the true phase function (multiplied with the single scattering albedo) that was used in the simulation of the radiation density Aeld and P˜ (i) is its reconstruction, the explained variance is m (P˜ (i) − P˜ k(i)∗ )2 ; (16) vi = 1 − k=1 k2 #P˜(i)∗ thus it is just 1 minus the relative mean squared error between true and reconstructed phase function. Also shown in Fig. 9 is the accuracy of the radiation density Aeld that has been obtained by the assimilation according to Section 4.3. Some examples of these reconstructed radiation intensities are presented and compared with the

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

155

Fig. 10. Simulated radiation intensities (left), reconstructed intensities (center), and the diLerences (right), side looking with 180◦ Aeld of view, at an altitude of 300 km (TOA). The white dot marks the position of the Sun. The units correspond to a solar Gux at TOA normalized to 4.

Fig. 11. As Fig. 10, but at 160 km altitude.

Fig. 12. As Fig. 10, but at 112 km altitude.

modeled (true) radiation that has been used to simulate the SA measurements in Figs. 10–14. 6. Discussion and conclusions The results show that—within the limitations discussed below—the phase functions of the scattering particles and the radiation density Aeld inside Titan’s atmosphere can be reconstructed from Huygens’ SA observations with an explained variance larger than about 90%, cf. Fig. 9. The most

pronounced errors occur in the tail of the phase function, cf. Figs. 7 and 8. This corresponds to the SA frames taken in the opposite direction of the Sun, where the radiation intensity is much lower, cf. Fig. 5, and thus the signal-to-noise ratio is lower because noise of the same absolute magnitude has been added to all simulated observations. The SA frames opposite to the Sun will probably be taken with a much higher integration time, which increases the signal-to-noise ratio. A varying absolute accuracy of the individual measurements can be taken into account using the so-called statistical regularization (Craig and Brown, 1986).

156

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

Fig. 13. As Fig. 10, but at 88 km altitude.

Fig. 14. As Fig. 10, but at the surface.

If the integration time is not negligible compared to the rotation period of the probe, the 1◦ × 1◦ FOV of one SA frame pixel will be smeared out. This eLect can be incorporated in the matrixes in Eq. (11). For the frame taken towards the Sun, the integration time will be limited by the demand not to exceed a 3:5◦ rotation of the probe. Thus the loss of resolution will be smaller than that imposed by the averaging over six pixels, which was assumed in the calculations but will probably not be done for the frame towards the Sun, cf. Section 3.2. For the frame opposite to the Sun, the integration time will not exceed a 20◦ rotation. Because of the small gradients of the sky brightness in this direction, this will not imply a considerable loss of information. Of course the presented results are only valid for the Titan scenario that was used to test the approach, i.e. the aerosol model described in Section 2. If the optical properties of Titan’s atmosphere are very diLerent from this scenario, the accuracy of reconstructed quantities can be diLerent. A sensitivity experiment in which the optical thickness of the haze layer was increased by a factor of four shows that in the optically thickest region near the bottom of the haze layer the inGuence of multiple scattering becomes signiAcant. The reconstructed phase function—which represents the scattering properties of the complete layer—is then very diLerent from the true single scattering phase function. Even in this case it could be possible to retrieve the single scattering phase

function by subtraction of the multiply scattered intensity (Nakajima et al., 1983), but some information would be lost and thus the accuracy of reconstructions would decrease. In order to prevent the layers between two successive measurements to become optically thick, observation cycles should be carried out with respective frequency. More frequent observations would also increase the vertical resolution for reconstructed optical properties. On the other hand, if the optical thickness of the layer between two considered measurements is lower than about 0.1, it turns out that the diLerence of the two measurements, cf. Eq. (7), is so small that the low signal-to-noise ratio increases the errors in the reconstructed phase functions. Therefore, less frequent observation cycles where used in the optically thinner region below the haze layer, cf. Section 4.3. The solar zenith angle has been assumed to be constant during Huygens’ descent, but in fact it may vary due to considerable drift caused by a global zonal wind system (Hubbard et al., 1993; Flasar et al., 1997; Flasar, 1998) by up to about 10◦ . This does not pose a principle problem and could be taken into account by assimilating the SA data simultaneously into a set of radiation models running for various solar zenith angles. But of course this would considerably increase the computation expense. All radiative transfer calculations herein are based on the plane parallel approximation, the inverse reconstructions as

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

well as the forward simulation of observations. If the plane parallel approximation is not valid for the optical conditions in Titans atmosphere for the solar zenith angle during Huygens’ descent and the viewing directions of the SA, the accuracy of reconstructions based on this approximation will decrease correspondingly. To investigate this potential problem, radiative transfer calculations in full spherical geometry will be carried out. A radiative transfer model for a spherical planetary atmosphere is available (Rozanov et al., 2001) but has not yet been applied to Titan. In the reconstructions it was assumed that the opacity, Eq. (4), of the considered layer was known, cf. Eqs. (9) and (11), as it enters into the matrix K (i) of the linear equation system in Eq. (13). The opacity of a layer between the altitudes of two DISR observation cycles can in principle independently be determined from another subinstrument, the upward-looking visible spectrometer (ULVS, Tomasko et al., 1997). However, any error in the estimation of this optical thickness might inGuence the reconstruction of phase functions, which has not been considered herein. If the angular dependence of the true radiation Aeld is not known, the opacity estimated from ULVS data can have an error of about 10%. However, due to the fact that the radiation Aeld varies very smoothly with altitude, the error imposed by this on phase function reconstructions is much smaller. The radiation density Aelds also estimated by SA data assimilation can in turn be used to correct the opacity estimations. The best approach may be to reconstruct phase functions and opacities simultaneously combining the use of SA and ULVS data, but then the problem would be nonlinear even for one layer. The inverse method described herein can also be useful for other planetary missions where a solar aureole is observed. This addresses instruments carried by landers or other Gying craft (i.e. planes, helicopters, balloons, airships) that perform upward looking measurements at diLerent altitudes. Besides Titan, the planets Venus and Mars are targets worth such investigations. Acknowledgements This work was supported by the German Aerospace Center (DLR), contract number 50 OH 98044. The program used for the icosahedron-based pixelization of the sphere has been made available on the World Wide Web by Max Tegmark. References Awal, M., Lunine, J.I., 1994. Moist convective clouds in Titan’s atmosphere. Geophys. Res. Lett. 21, 2491–2494. Cabane, M., Rannou, P., ChasseAUere, E., Israel, G., 1993. Fractal aggregates in Titan’s atmosphere. Planet. Space Sci. 41, 257–267. Courtin, R., Wagener, R., McKay, C.P., Caldwell, J., Fricke, K.-H., Raulin, F., Bruston, P., 1991. UV spectroscopy of Titan’s atmosphere, planetary

157

organic chemistry, and prebiological synthesis. II. Interpretation of new IUE observations in the 220 –335 nm range. Icarus 90, 43–56. Craig, J.J.D., Brown, J.C., 1986. Inverse Problems in Astronomy. Adam Hilger Ltd, Bristol. Draine, B.T., 1988. The discrete-dipole approximation and its application to interstellar graphite grains. Astrophys. J. 333, 848–872. Flasar, F.M., 1998. The dynamic meteorology of Titan. Planet. Space Sci. 46, 1125–1147. Flasar, F.M., Allison, M.D., Lunine, J.I., 1997. Titan zonal wind model. In: Wilson, A. (Ed.), Huygens—Science, Payload and Mission, Vol. SP-1177. ESA Publications Division, ESTEC, Noordwijk, The Netherlands, pp. 287–298. Franklin, J.L., 1970. Well-posed stochastic extensions of ill-posed linear problems. J. Math. Anal. Appl. 31, 682–716. Grieger, B., Kayser, R., Schramm, T., 1991. The deconvolution of the quasar structure from microlensing light curves. Astron. Astrophys. 252, 508–512. Gri8th, C.A., Owen, T., Miller, G.A., Geballe, T., 1998. Transient clouds in Titan’s lower atmosphere. Nature 395, 575–578. Hansen, P.C., O’Leary, D.P., 1993. The use of the L-Curve in the regularization of discrete ill-posed problems. SIAM J. Sci. Comput. 14, 1487–1503. Hubbard, W.B., et al., 1993. The occultation of 28 Sgr by Titan. Astron. Astrophys. 269, 541–563. Hunton, D.M., Tomasko, M.G., Flasar, F.M., Samuelson, R.E., Strobel, D.F., Stevenson, D.J., 1984. Titan. In: Gehrels, T., Matthews, M.S. (Eds.), Saturn. University of Arizona Press, Tucson, pp. 671–759. Khare, B.N., Sagan, C., 1973. Red clouds in reducing atmospheres. Icarus 20, 311–321. Khare, B.N., Sagan, C., Arakawa, E.T., Suits, F., Callcott, T.A., Williams, M.W., 1984. Optical constants of organic tholins produced in a simulated Titanian atmosphere: from soft X-ray to microwave frequencies. Icarus 60, 127–137. Kuiper, G.P., 1944. Titan: a satellite with an atmosphere. Astrophys. J. 100, 378–383. Lebreton, J.-P., Matson, D.L., 1997. The Huygens probes: Science, payload and mission overview. In: Wilson, A. (Ed.), Huygens— Science, Payload and Mission, Vol. SP-1177. ESA Publications Division, ESTEC, Noordwijk, The Netherlands, pp. 5 –24. Lemmon, M.T., 1994. Properties of Titan’s haze and surface. Ph.D. Thesis, University of Arizona. Lorenz, R.D., Lemmon, M.T., Smith, P.H., Lockwood, G.W., 1999. Rapid and complex seasonal change on Titan observed with the Hubble Space Telescope. Icarus 142, 391–401. Markiewicz, W.J., Sablotny, R.M., Keller, H.U., Thomas, N., Titov, D., Smith, P.H., 1999. Optical properties of the Martian aerosols as derived from Imager for Mars PathAnder midday sky brightness data. J. Geophys. Res. 104, 9009–9017. McKay, C.P., Coustenis, A., Samuelson, R.E., Lemmon, M.T., Lorenz, R.D., Cabane, M., Rannou, P., Drossart, P., 2001. Physical properties of the organic aerosols and clouds on Titan. Planet. Space Sci. 49, 79–99. Meakin, P., Chen, Z.-Y., Deutch, J.M., 1985. The translational friction coe8cient and time dependent cluster size distribution of three-dimensional cluster–cluster aggregation. J. Chem. Phys. 82, 3786–3789. Meakin, P., Donn, B., Mulholland, G.W., 1989. Collisions between point masses and fractal aggregates. Langmuir 5, 510–518. Nakajima, T., Tanaka, M., Yamauchi, T., 1983. Retrieval of the optical properties of aerosols from aureole and extinction data. Appl. Opt. 22 (19), 2951–2959. NeL, S.J., Humm, D.C., Bergstralh, J.T., Cochran, A.L., Cochran, W.C., Barker, E.S., Tull, R.G., 1984. Absolute spectrophotometry of Titan, X Icarus 60, 221–235. Uranus, and Neptune: 3500 –10500 A. Plass, G.N, Kattawar, G.W., Catchings, F.E., 1973. Matrix operator theory of radiative transfer. 1: Rayleigh scattering. Appl. Opt. 12 (2), 314–329.

158

B. Grieger et al. / Planetary and Space Science 51 (2003) 147 – 158

Rages, K., Pollack, J.B., Smith, P.H., 1983. Size estimates of Titan’s aerosols based on Voyager 1 high-phase-angle images. J. Geophys. Res. 88, 8721–8728. Rannou, P., Cabane, M., ChasseAere, E., Botet, R., McKay, C.P., Courtin, R., 1995. Titan’s geometric albedo: role of the fractal structure of the aerosols. Icarus 118, 355–372. Rannou, P., McKay, C.P., Botet, R., Cabane, M., 1999. Semi-empirical model of absorption and scattering by isotropic fractal aggregates of spheres. Planet. Space Sci. 47, 385–396. Reginska, T., 1996. A regularization parameter in discrete ill-posed problems. SIAM J. Sci. Comput. 17, 740–749. Rozanov, A., Rozanov, V., Burrows, J.P., 2001. A numerical radiative transfer model for a spherical planetary atmosphere: combined diLerential-integral approach involving the picard iterative approximation. J. Quant. Spectrosc. Radiat. Transfer 69, 491–512. Sromovsky, L.A., Suomi, V.E., Pollack, J.B., Kraus, R.J., Limaye, S.S., Owen, T., Revercomb, H.E., Sagan, C., 1981. Implications of Titan’s north–south brightness asymmetry. Nature 292, 698–702. Taylor, F.W., Coustenis, A., 1998. Titan in the Solar System. Planet. Space Sci. 46, 1085–1097. Tegmark, M., 1996. An icosahedron-based method for pixelizing the celestial sphere. Astrophys. J. Letts. 470, L81–L84. Thomas, N., Markiewicz, W.J., Sablotny, R.M., Wuttke, M.W., Keller, H.U., Johnson, J.R., Reid, R.J., Smith, P.H., 1999. The color of the Martian sky and its inGuence on the illumination of the Martian surface. J. Geophys. Res. 104, 8795–8808. Tokano, T., Neubauer, F.M., Laube, M., McKay, C.P., 1999. Seasonal-variation of Titan’s atmospheric structure simulated by a general circulation model. Planet. Space Sci. 47, 493–520.

Tomasko, M.G., Smith, P.H., 1982. Photometry and polarimetry of Titan: Pioneer 11 observations and their implications for aerosol properties. Icarus 51, 65–95. Tomasko, M.G., Doose, L.R., Smith, P.H., West, R.A., Soderblom, L.A., Combes, M., BZezard, B., Coustenis, A., deBergh, C., Lellouch, E., Rosenqvist, J., Saint-PZe, O., Schmitt, B., Keller, H.U., Thomas, N., Gliem, F., 1997. The descent imager/spectral radiometer (DISR) aboard Huygens. In: Wilson, A. (Ed.), Huygens—Science, Payload and Mission, Vol. SP-1177. ESA Publications Division, ESTEC, Noordwijk, The Netherlands, pp. 109 –138. Toon, O.B., McKay, C.P., Pollack, J.B., 1980. A physical model of Titan’s aerosols. Icarus 43, 260–282. Turchin, V.F., Kozlov, V.P., Malkevich, M.S., 1971. The use of mathematical-statistics methods in the solution of incorrectly posed problems. Sov. Phys. Usp. 13, 681–703. West, R.A., 1991. Optical properties of aggregate particles whose outer diameter is comparable to the wavelength. Appl. Opt. 30, 5316–5324. West, R.A., Smith, P.H., 1991. Evidence for aggregate particles in the atmospheres of Titan and Jupiter. Icarus 90, 330–333. West, R.A., Lane, A.L., Hart, H., Simmons, K.E., Hord, C.W., CoLen, D.L., Esposito, L.W., Sato, M., Pomphrey, R.B., 1983. Voyager 2 photopolarimeter observations of Titan. J. Geophys. Res. 88, 8699–8708. Yanovitskij, E.G., 1997. Light Scattering in Inhomogeneous Atmospheres. Springer, Berlin. Yung, Y.L., Allen, M., Pinto, J.P., 1984. Photochemistry of the atmosphere of Titan: comparison between model and observations. Astrophys. J. Suppl. Ser. 55, 465–506.