© P. Rannou et al. 2022

1 Introduction

Saturn’s moon, Titan, is a complex world with active cycles, lakes, and seas of liquid hydrocarbons (e.g. Stofan et al. 2007; Mitri et al. 2007; Mastrogiuseppe et al. 2014), clouds (e.g. Griffith et al. 2005, 2006; Brown et al. 2010; Rodriguez et al. 2009, 2011), and rains (e.g. Turtle et al. 2011; Dhingra et al. 2019, 2020). Methane and nitrogen, the two main gases, are photolysed in atmosphere layers above 1000 km and yield a complex photochemistry, a large variety of molecules (e.g. Vuitton et al. 2009), and a planetary-scale photochemical haze. Haze hides the lower atmosphere and the surface from view. Thanks to their optical instruments probing in near-infrared and thanks to the RAdio Detection And Ranging (RADAR), Cassini and Huygens revealed an incredibly rich and various surface: mountains, dry river beds, ancient seas or lakes, dunes, craters, and a wide variety of other surface features (e.g. Tomasko et al. 2005; Radebaugh et al. 2008; Turtle et al. 2011; Lopes et al. 2020). This is a signature of an active and wet climate. During the time of the Cassini mission, some observations also demonstrated that the surface is still actively modelled by rains and winds, with obvious changes with time in some places (Turtle et al. 2011; MacKenzie et al. 2019).

A first motivation for detailed studies of atmospheric properties is to better understand and predict Titan’s climate. Gases and particles in the atmosphere play a major role in Titan’s climate through their interactions with the radiation field, which further affects the thermodynamical equilibrium and atmosphere circulation. The haze layer produces a strong anti-greenhouse effect, heating the stratosphere while it cools the troposphere (McKay et al. 1991; Bézard et al. 2018). Gases produce a greenhouse effect through collision-induced absorption, and selective spectroscopic absorption in the stratosphere produces a cooling to space. The main effect of clouds is to control the vertical distribution of gases, especially methane, and to produce haze scavenging (Rannou et al. 2006; Barth & Toon 2006). They probably also participate in the distribution of radiative energy in polar regions, where they have a long lifetime (Le Mouélic et al. 2012, 2018; West et al. 2018), but this has not been studied yet. In return, atmospheric circulation also acts on the spatial distribution of these components and on their formation processes. It was shown, for instance, that the spatial structure of haze and clouds is a good tracer of atmospheric circulation (e.g. Rannou et al. 2002, 2006; Larson et al. 2015; Rodriguez et al. 2011; West et al. 2018; Turtle et al. 2018).

Efforts to characterise the surface from photometric observations depend on our understanding of atmosphere properties: that is, either optical or spectral properties and spatial distributions. In many previous works, surface properties have been categorised with their colour as seen through atmosphere in methane windows or with synthetic colour based on intensity ratios (e.g. Barnes et al. 2007, 2011; MacKenzie et al. 2019), or they have been analysed from a statistical perspective as the principal analysis component (e.g. Solomonidou et al. 2014; Karkoschka & Schröder 2016; Griffith et al. 2019). In another category of studies, radiative transfer modelling is used in order to remove the influence of the atmosphere scattering and absorption to derive the surface reflectivity (e.g. Griffith et al. 2012; Hirtzig et al. 2013; Solomonidou et al. 2016; Brossier et al. 2018; Coutelier et al. 2021). In these works, the validity of the results strongly depends on the quality of the models used to account for the atmosphere properties. The constrast in reflectivity is large enough between different regions or areas of Titan’s surface and between different parts of the spectra to retrieve useful surface properties with radiative transfer modelling. However, because of the lack of knowledge concerning atmosphere properties, we still lack the fine details of surface spectra within methane windows (especially the windows at 2 and 2.8 µm), and even the absolute surface reflectivity for each window could also be biased. A better characterisation of haze and gas absorption and their scattering properties are also compulsory for producing valuable analysis of Titan surface.

In this work, we analysed four observations acquired by the infrared part of the Visual Infrared Mapping Spectrometer (VIMS), hereafter VIMS-IR, onboard the Cassini spacecraft, in occultation mode (Maltagliati et al. 2015). In this observation mode, the instrument monitors the light emitted by the sun or a distant star and transmitted through the atmosphere while the spacecraft moves along its orbit. This is a very efficient way to probe different layers of an atmosphere without disturbance due to multiple scattering produced by the surface or the deeper atmosphere. VIMS-IR is able to measure Titan spectra between about 0.884 and 5.122 µm, and it is sensitive to haze and gas extinction. At short wavelengths, both haze and gas have strong absorption features; in this case, it is possible to probe up to 400 km. At long wavelengths, the haze is optically thin, and in the 5 µm methane window, only refraction by the atmosphere can limit retrievals at low altitude. In this spectral interval, the extinction continuum due to haze and gaseous absorption bands enables us to characterise the morphology of the aerosols, the vertical distribution of haze extinction, and of some gas mixing ratios.

In the next parts of this article, we display the data analysed in this work and information related to these observations (Sect. 2), and we discuss the different contributions to opacity in Titan’s atmosphere (Sect. 3). We then discuss the effect of atmosphere refraction and of the possible impact of multiple scattering on the observed intensities (Sect. 4). We carried out an initial analysis of haze extinction in the methane windows (the continuum) between 0.884 and 2.5 µm to constrain the haze’s spectral behaviour, and we analysed aerosol morphology in the altitude range of 200–300 km (Sect. 5). In this paper, we describe the procedure used to retrieve haze properties and gas mixing ratios. We then explain the main principles of the Bayesian inference complemented with a Tikhonov regularisation to remove spurious oscillations (Sect. 6). In the following, we split the observed spectra in several intervals, and we study haze, methane, deuterated methane, and carbon monoxide separately. For each retrieved species, we dedicate a sub-section to displaying and discussing the results (Sect. 7). The last section features our conclusions.

2 VIMS-IR occultation data

In occultation observation mode, VIMS-IR observes the light emitted by a distant source (the Sun or a star) along lines of sight through the atmosphere at the limb. As the spacecraft moves around Titan, it probes different layers. We used the data (Luca Maltagliati, personal communication) presented in the paper of Maltagliati et al. (2015), where data processing is explained in detail. The 2006 dataset (from the T10 flyby) was also analysed by Bellucci et al. (2009). In theory, such observations are self-calibrated, because observations along lines of sight passing outside the atmosphere give the reference flux Φ₀, and the quantity of interest is the transmission along the line of sight through the atmosphere, calculated from the attenuated flux Φ, as T = Φ/Φ₀. In this work, we use the term ‘flux’ for a power surface density and the term ‘intensity’ for a power surface density per solid angle unit.

In practice, there are many effects that must be accounted for. Bellucci et al. (2009) corrected for a drift in the observed flux before the beginning of the occultation. Maltagliati et al. (2015) considered ten occultations, but for some observations the instrument pointing was not stable enough to ensure that the image of the source was located on the same group of pixels on the VIMS-IR detector. When the image moved on the detector, this produced flux variations due to different pixel responses on the detector array that cannot be corrected well enough. Then, only four observations with good pointing stability were retained for our analysis (Table 1).

Stray light also appeared in three out of four observations. It showed up as a strong peak superimposed on the real signal in the transmission curves. The stray light was not uniform on the VIMS-IR detector and changed with time. Maltagliati et al. (2015) modelled the observed signal as the true signal from the observed source with a 2D Gaussian plus a flat signal, and the contaminating stray light as a quadratic surface. The flat signal and the stray light were then removed from the observed signal. This sophisticated procedure finally yielded the transmission curves that were used in our analysis.

The final data are a set of transmissions T through the atmosphere along a tangential path at the limb (Fig. 1). This also gives direct access to the tangential total opacity of the atmosphere from the relation τ = − log(T).

3 Opacities in the atmosphere

Light transmission through a tangential path at limb is affected by the absorption due to the haze and several absorbing gases. In this work, we used the gas absorption lines published in the HIgh-resolution TRANsmission molecular absorption database (HITRAN; Gordon et al. 2017), except for the methane lines, which are those published by Rey et al. (2016, 2018). Gas absorption lines are converted into absorption coefficients with a treatment including the correction of the line strength with temperature, the line convolution by a Voigt profile with a collision line broadening half-widths as determined by Menard-Bourcin et al. (2007), and a cut-off of the Lorentzian far wings at ∆_v = 26 cm⁻¹ from the line centre and a sub-lorentzian decay σ = 120 cm⁻¹ (e.g. de Bergh et al. 2012; Hirtzig et al. 2013; Rannou et al. 2018). Finally, we cast gas absorptions in correlated-k tables (Goody et al. 1989) and used 16 terms.

In the spectral range probed by VIMS-IR, Titan spectra are essentially dominated by haze and methane properties between 0.88 and 4.25 µm, with some contributions of minor gases in the methane windows (Fig. 2). Beyond 4.25 µm, absorption is due to CH₃D, CO, and haze. This allows us to split the dataset into spectral intervals where haze and gaseous species can be retrieved separately, as displayed in Fig. 2.

The gas cross-sections are calcutated at different pressures between 1.46 × 10⁵ and 5 × 10⁻³ Pa. For each pressure level, we used the corresponding temperature T_h(P) retrieved by the Huygens Atmospheric Structure Instrument (HASI; Fulchignoni et al. 2005) at this pressure level and two other temperatures, $T_{\text{h}}^{\pm }\left(P\right)=T_{\text{h}}\left(P\right)\pm \text{Δ}T\left(P\right)$ (Fig. 3). We used ∆T(P) = 5 K in the troposphere and 15 K in the stratosphere. Then, for any temperature and pressure profiles, we interpolated through the table of computed cross-sections with a linear interpolation in log(P) and with the following second-order interpolation in log(T). This represents the observed variations of k(P, T) with P and T well:

$$\begin{array}{lll}\log k\left(P,T\right)\hfill & =\hfill & \log k\left(P_0,T_0\right)+\frac{\partial \log k\left(P,T\right)}{\partial \log P}\times \left(\log P-\log P_0\right)\hfill \\ \hfill & \hfill & +\frac{\partial \log k\left(P,T\right)}{\partial \log T}\times \left(\log T-\log T_0\right)\hfill \\ \hfill & \hfill & +\frac{1}{2}\frac{{\partial }^2\log k\left(P,T\right)}{\partial \log T^2}\times {\left(\log T-\log T_0\right)}^2.\hfill \end{array}$$(1)

For the four observed occultations, we used the temperature profiles derived from the work published by Vinatier et al. (2015) and Mathé et al. (2020). However, analysis of the Composite InfraRed Spectrometer (CIRS) observations is based on an a priori hypothesis about the methane mixing ratio in the stratosphere, assumed to be $x_{{\text{CH}}_4}=1.48\times {10}^{-2}$ in these works. We remark that the temperature retrieved with CIRS at a latitude close to the equinox in 2009 (CIRS profile chosen for occultation T53, Table 1) is ≃5 K lower than the profile observed by HASI (Fulchignoni et al. 2005). This could be explained by the negative correlation between the assumed methane mole fraction and the retrieved temperature with CIRS observation (e.g. Lellouch et al. 2014). We also consider that low vertical resolution of CIRS retrieval may decrease the temperature constrast, especially near the stratopause where the temperature reaches a maximum. Finally, it could also be due to the seasonal cycle since these two observations are taken four years apart. However, the CIRS temperature for 2009 near equator is in very good agreement (within ≃2 K) with the temperature profile retrieved with the radio-occultation experiment (Schinder et al. 2020). Both differ in the same way from the Huygens HASI profile. We then decided to use CIRS temperatures without any correction. However, a sensitivity test including both the effect of temperature on the absorption properties and on the amount of molecules for a given molar fraction through the induced change in the hydrostatic structure shows that a temperature difference of ≃5 K does not change the final result in a significant way.

In order to analyse occultations, only extinction properties of aerosols are needed. We used a model of scattering and absorption by fractal aggregates of identical spheres able to account for the fractal dimension (Botet et al. 1997; Rannou et al. 1997). This model has long been validated against the Discrete Dipole Scattering (DDSCAT) model and the T-Matrix model (e.g. Tazaki & Tanaka 2018) and a comparison with the model used by Tomasko et al. (2008) shows an excellent agreement for aggregates with a fractal dimension of 2 and with 3000 monomers of 50 nm in radius (Fig. 3).

Fig. 1 Transmission curves as function of wavelength observed by VIMS-IR during the four flybys. For each flyby, each curve corresponds to transmissions at a given altitude between the surface and 400 km. The transmission increases with altitudes; the transmissions at 400 km are close to 1, and those near the surface are close to 0.

Fig. 2 Vertical opacity in our model in the middle stratosphere (at 200 km altitude) as a function of the wavelength for five gases and for the haze with a constant refractive index (m = 1.6 + 0.01i), The grey horizontal bars show where ethane has weak (thin bar), medium (medium bar), or strong (thick bar) absorption features (Maltagliati et al. 2015). The vertical bars delimit intervals where different components can be studied separately: (1) haze and CH4 between 0.88 and 2.8 µm; (2) haze and CH3D between 2.8 and 2.95 µm; (3) haze, CH4, C2H2, and HCN between 3.10 and 4.20 µm; (4) haze, CH3D, and CO beyond 4.2 µm. In the absence of any reliable line list for C2H6 in the spectral range of VIMS-IR, the spectral ranges where ethane has strong features cannot be studied and spectral ranges where ethane has medium strength should be studied with caution.

Fig. 3 Preliminary control of model background properties: aerosol optical cross-sections and atmosphere temperature profiles Left: extinction (upper curves) and scattering (lower curves) cross-sections for an aggregate of 3000 monomers of radius 50 nm with a fractal dimension Df = 2, computed with the mean field approximation model (Botet et al. 1997; Rannou et al. 1997) and with the model used by Tomasko et al. (2008). Right: temperature profile as a function of pressure observed with HASI (Fulchignoni et al. 2005) and profiles based on HASI with a temperature shift ±∆T (red dots). We also plotted the temperature profiles after CIRS observations (Vinatier et al. 2015; Mathé et al. 2020) that best correspond to the dates and the latitudes of the occul-tations treated in this work (Table 1). The thick black dashed line shows the temperature profile retrieved by Schinder et al. (2020) with radio-occultation measurements during the flyby T52E (4 April 2009 at lat = 12°S).

4 Correction of the data due to physical effects

Before performing complete haze and gas retrievals, we have to consider two specific effects that may affect observations at low altitude. During occultation measurements, VIMS-IR makes observations of the solar flux transmitted through Titan’s atmosphere. This transmitted flux, compared to the source flux, is used to determine the transmission spectra through the atmosphere. If effects other than extinction occur along the line of sight, this could be mistakenly considered as being due to extinction only. In this part, we consider two effects that may have an impact on the observed fluxes.

4.1 Impact of the haze forward scattering on transmitted flux

Maltagliati et al. (2015) accounted for the haze forward scattering along the line of sight. Their work included the single and multiple scattering effect. They did not treat the problem explicitly, but they incorporated a change of 10% in the retrieved methane mixing ratio following de Kok & Stam (2012). We have two remarks about this. First, de Kok & Stam (2012) studied this effect by integrating the light on the total field of view of VIMS, while only the pixels where the image of the Sun appears were actually used to produce the transmitted flux. For our study, as shown in the following, we find the forward scattering to be negligible compared to the transmitted light. Second, in their data processing Maltagliati et al. (2015) substracted the background intensity from the image of the Sun through the atmosphere. In doing so, we assume that they also removed all the contributions not from the direct Sun, including multiple scattering due to the surrounding atmosphere. Thus, bringing a supplementary correction to these data in order to remove contributions due to scattering is of no use. However, here we show that the single scattering contribution, which dominates multiple scattering in high layers, can be neglected in our observations. The contribution to the scattered intensity by a segment of the line of sight, with an optical length dτ, can be written as follows:

$$\text{d}I={\varphi }_0\frac{\overline{\omega }-P\left(0^{\text{o}}\right)}{4\pi }\exp \left(-\tau \right)\exp \left(-\left({\tau }_{1\text{os}}-\tau \right)\right)\text{d}\tau ,$$(2)

where ϕ₀ is the flux of the source before the entry in the atmosphere, $\overline{\omega }$ the average single scattering albedo, P(0°) the haze phase function in the forward direction, and τ_los the total opacity along the line of sight. The opacities before and after the scattering are given, respectively, by τ and τ_los − τ, where τ is the opacity counted from the entry point. To compare the flux transmitted through the atmosphere and the flux scattered by haze along the line of sight, we integrated dI along the line of sight and then over the field of view of VIMS-IR. At a given altitude, the scattered intensity integrated along the line of sight in the forward direction yields

$$I_{\text{scat}}={\varphi }_0\frac{\overline{\omega }P\left(0^{\text{o}}\right)}{4\pi }\exp \left(-{\tau }_{\text{los}}\right)\times {\tau }_{\text{los}}.$$(3)

We note the pixel field of view Ω_VIMS and the field of view of the image of the sun is given as N_pxl × Ω_VIMS, where N_pxl is the number of pixels composing the image of the sun on the detector. The width of a pixel field of view (5 × 10⁻⁴ radian) is much smaller than the angular width of the phase function forward peak for Titan aerosols. The width of the first diffraction lobe is ∆θ_D ≃ 1.22 × λ/d, where d is twice the aerosols radius, and we estimate that the phase function is flat, within 10%, for scattering angles between θ = 0 and up to ∆θ_D/5. For aggregates of several thousand monomers of 50 nm, ∆θ_D/5 is several degrees (°) (>2 × 10⁻² radian), and, in the most constraining case considered in this work (fractal aggregate of 10⁶ monomers of 50 nm with D_f = 2), ∆θ_D/5 is about 0.1° (≃1.7 × 10⁻³ radian). With these conditions, we then always assume that P(θ) = P(0°) inside the field of view. We also assume that the atmosphere properties do not vary inside the field of view. The projection of a pixel field of view on the plane of Titan’s limb has a size of 3–5 km, depending on the distance to the planet (Table 1). At this scale, we assume no significant spatial variation of the atmosphere properties. We also estimate that the vertical variation of the tangential opacity along the line of sight due to the variation of densities with altitude, and that of the geometrical length of the line of sight with altitude, cause relative variations of less than 4 × 10⁻⁴ on the integral over a pixel field of view, which is negligible. We then conclude that the integration over the field of view can be performed as a simple product of P(0°) and Ω_VIMS. The same rule applies if the image of the sun is spread over several pixels (N_pxl). In this case, we used the solid angle N_pxl × Ω_VIMS and the width of the total field of view, $N_{\text{pxl}}^{1/2}\times {5.10}^{-4}$, which is still smaller than ∆θ_D/5. The equation that we finally use for the scattered flux is valid only under these specific conditions. In the case of Titan aerosols and with the VIMS spectral range, these conditions are always met, even in the most constraining cases.

The scattered flux mixed with the image of the sun can be written ϕ_scat = I_scat × Ω_vimsN_pxl. The quantity exp(−τ_los) × τ_los in the description of I_scat is bounded between 0 (at τ_los = 0 or τ_los → ∞) and e⁻¹ (at τ_los = 1), and for one pixel Ω_vims is 0.5 × 0.5 mrad², (i.e. 2.5 × 10⁻⁷ sr). The single scattering albedo $\overline{\omega }$ is physically bounded by 1, and P(0°) is several tens sr⁻¹. For instance, the largest value for the spectral range of VIMS-IR given by Tomasko et al. (2008) is P(0°) = 228 sr⁻¹ at 0.8224 µm. The maximum ratio between the flux of scattered light and the solar flux $\left(r_{\text{scat}}^{\max }={\varphi }_{\text{scat}}/{\varphi }_0\text{\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}}{\tau }_{\text{los}}=1\right)$ takes the following form:

$$r_{\text{scat}}^{\max }\left(\frac{I_{\text{scat}}\left({\tau }_{1\text{os}}=1\right){\Omega }_{\text{VIMS}}{N}_{\text{px}1}}{{\varphi }_0}\right)=\frac{\overline{\omega }P\left(0^{\text{o}}\right)}{4\pi }{e}^{-1}{\Omega }_{\text{VIMS}}{N}_{\text{px}1}.$$(4)

According to Maltagliati et al. (2015), N_pxl is generally between 2 and 4, and in any case lower than 9. Thus, in the worst case (N_pxl = 9), we still find that the scattered flux is 10⁻⁴ times smaller that ϕ₀. In our analysis, we never treat transmission T = ϕ_obs/ϕ₀ smaller than 10⁻², that is, 100 times larger than the threshold beyond which forward scattering would matter. In reality, the ratio $r_{\text{scat}}^{\max }={10}^{-4}$, taken at τ_los = 1, should be compared to the transmission for the same value of τ_los, which is T = e⁻¹. We then conclude that single scattering in the forward direction is negligible compared to the transmitted light in our data.

4.2 Flux attenuation due to refractive divergence

A second optical effect acts on the transmission of light. As a beam crosses the atmosphere tangentially, it can be bent according to the Snell-Descartes law because it travels through layers with varying real refractive indices. The bending globally increases with decreasing altitude and increasing atmosphere density. When observing the solar or a stellar flux through an atmosphere, this differential bending causes the atmosphere to produce divergence of an incident plane wave. For an instrument observing occultation, it appears as an attenuation of the beam. If one is interested in extinction due to gases or haze in the atmosphere, this refractive attenuation effect should be corrected for.

In this work, we used a first-order correction. We integrated the deviation angle θ along the line of sight from the source to the instrument. We first differentiate the Snell-Descartes law n(λ) sin i = c₀, where n(λ) is the refractive index of the atmosphere at the wavelength λ, i is the incident angle relative the local vertical direction and c₀ a constant. n(λ) has the form 1 + α_E(λ)/2 × ρ_n(z), where α_E(λ) is the molecular polarisabil-ity and ρ_n(z) is the molecule number density in the atmosphere. The differentation of the Snell-Descartes equation along the line of sight (axis x) yields

$$\frac{1}{\sin i}\frac{\partial \sin i}{\partial x}=-\frac{1}{n(\lambda )}\frac{\partial n(\lambda )}{\partial x},$$(5)

where x is the distance taken along the line of sight, with x = 0 being the abscissa of the point at the lowest altitude z₀ reached by the line of sight in the atmosphere (that is, the impact parameter). To compute the deviation angle, we changed ∂n(λ)/∂x into ∂n(λ)/∂z × ∂z/∂x, ρ_n(z) into P(z)/k_BT(z), where k_B is the Boltzmann constant, and we used $z=\sqrt{{\left(R_{\text{T}}+z_0\right)}^2+x^2}-R_{\text{T}}$. Here, R_T is Titan’s radius. We assume, at first order, that the inflexion of the path of the light is negligible. It is then considered almost as a straight line, and the value of the angle i can be directly related to the location x and R_T + z₀ along the line of sight. We write

$$\frac{1}{\sin i}\frac{\partial \sin i}{\partial x}=\frac{1}{\text{sin\hspace{0.17em}}i}\frac{\partial \sin i}{\partial i}\times \frac{\partial i}{\partial x},$$(6)

and here we assimilate di to the deviation dθ along the path, while the remaining part of the term, (tan i)⁻¹, only depends on x and R_T + z₀. After a straightforward calculation, we finally obtain

$$\text{d}\theta =-\frac{1}{n\left(z\right)}\frac{R_{\text{T}}+z_0}{\sqrt{{\left(R_{\text{T}}+z_0\right)}^2+x^2}}\times \frac{{\alpha }_{\text{E}}(\lambda )}{2}\frac{\partial {\rho }_n(z)}{\partial z}\text{d}x;$$(7)

dθ has to be integrated from x = −∞ to x = +∞ to give the total deviation θ(z₀) along the line of sight as a funtion of the impact parameter z₀. This requires a temperature and pressure profile and information about the polarisability α_E. To fix α_E(λ), we used the results of Washburn (1930) for pure nitrogen, as reported by Sicardy et al. (2006). We found this to be in agreement with Ciddor (1996), relevant for Earth’s atmosphere, within 2%. We also found that α_E(λ) is almost constant in the VIMS-IR spectral range within 0.9% relatively to the value at 5 µm. The deviation angle as a function of the altitude (the impact parameter) is shown in Fig. 4 (left).

The attenuation produced by atmosphere refraction is due to the divergence of the incident light plane wave. If two parallel rays are separated by an altitude difference, ∆z, before entering the atmosphere, their deviation at emergences differ by ∆θ = ∂θ/∂z × ∆z. Due to energy conservation, in a purely refractive case we know that all the energy concentrated between these two rays in the initial plane wave has to be conserved. This conservation imposes the following relation:

$$A=\frac{{\varphi }_{\text{out}}}{{\varphi }_0}=\frac{1}{1+D\times \text{Δ}\theta /\text{Δ}z},$$(8)

where D is the distance between the planet and the spacecraft and ϕ_out/ϕ₀ is the refractive attenuation that could be mistakenly treated as extinction. Figure 4 (middle) shows examples of refractive attenuation found near the equator of Titan for three representative values of D. Figure 4 (right) shows the altitude actually probed as a function of the impact parameter of the observation. We also compare our results with previous similar computations made by Bellucci et al. (2009) (Fig. 5). To perform our analysis, we first corrected the probed altitude by accounting for the bending of the light rays, and then we corrected the observed fluxes from the refractive attenuation. This correction essentially affects the altitude below 100 km and has an impact on the retrieval of CO and haze extinction in the 5 µm methane window.

Fig. 4 Effect of atmosphere refractivity for the observation T53E, calculated with the pressure and the temperature profiles retrieved by CIRS (flyby T51, 27 March 2009 at 1°S, Distance to the spacecraft : 6300 km - Table 1). Left: vertical profile of the deviation angle caused by the refraction in the atmosphere. The angular resolution of VIMS-IR (0.5 mrad) is indicated. Middle: vertical profile of refractive attenuation for three values of the distance from Cassini to Titan. Right: altitude actually probed by the light (continuous line) and the apparent altitude assuming that the light follows straight paths (dashed lines). This altitude assuming straight light is the impact parameter of the line of sight.

Fig. 5 Altitude where transmitted flux is attenuated by 1 and 50% as a function of the distance D of Cassini to Titan. The distance for each observation is reported in Table 1. The shapes (triangles and dots) correspond to the evaluation made by Bellucci et al. (2009). The blue and red lines show our calculations using the pressure and temperature profiles retrieved by CIRS (flyby T51, 27 March 2009 at 1° S) and by HASI onboard Huygens Fulchignoni et al. (2005), respectively (figure based on Bellucci et al. 2009).

5 Aerosol morphology

In this section, we first determine the morphological properties of the haze particles (aerosols) using transmission only due to aerosols in windows. In transparency windows, log(τ) as a function of log(λ) appears as straight lines between 0.8 and about 2.8 µm. This means that haze opacities at different altitudes follow power laws of the type τ_H(λ) = τ_H(λ₀) × (λ/λ₀)^α. The spectral slope α can then be used to characterise the spectral opacity of aerosols.

5.1 Tholin optical constant between 1 and 2.8 µm

Before considering aerosol scattering, we first discuss the aerosol material optical properties. A large sample of real and imaginary parts of the refractive indices of several tholins made in laboratory experiments are reported in Brassé et al. (2015). Using this review, we consider that the refractive indices appear fairly constant in the spectral range between 1 and 2.8 µm. This is consistent with the general observation that experimental tholins have no absorption peak in this spectral range. Generally, absorption peaks are found beyond ≃3 µm and have an impact on the refractive index only beyond about 2.8 µm. More specifically, it is shown that tholins have a constant value for the refractive index up to wavelength of 2.8 or 3.3 µm depending on the strength of the absorption due to the various types of N-H bonds. These are around ≃3 µm (Mahjoub et al. 2012; Sciamma-O’Brien et al. 2017).

In the first part of this work, we thus assumed that the complex refractive index is constant between 1and 2.8 µm. As consequence, the spectral slope of transmission in CH₄ windows in this wavelength range is only due to aerosol morphology, not to spectral variations of the aerosol optical constants.

5.2 Constraints on aerosol morphology with VIMS-IR

In this study, we assume that the real part of the refractive index is n = 1.6, which is representative of the tholin produced in laboratory (Brassé et al. 2015). We set the imaginary part of the refractive index κ and few morphological characteristics of the aerosols, assumed to be aggregates of spherical grains distributed with a fractal structure (West 1991; Cabane et al. 1993), as free parameters. The two free parameters linked to aerosol morphology are the fractal dimension of the aggregates, D_f and the number of grains per aggregates N. The radius of each grain (monomer) is set at r_m = 0.05 µm, following the evaluation made with photometric observations by Tomasko et al. (2008). Below, we also report tests carried out with the real part of the refractive index n = 1.5 and n = 1.7, and with a monomer radius of r_m = 0.04 µm, following Tomasko et al. (2009).

Considering n = 1.6 and r_m = 0.05 µm, we scanned the aerosol properties in the parameter space (N,κ,D_f) and calculated a goodness parameter to evaluate the best fit of the scaled tangential opacities. The imaginary refractive index κ varies between 10⁻³ and 1 , the number of grains N between 100 and 10¹⁰, and the fractal dimension D_f between 1 .8 and 2.5. We normalised all the tangential opacities at a given wavelength close to the centre of the spectral range of interest (λ₀ ≃ 1.2 µm) in such a way that all opacity points are aligned along a straight line in a log – log plot. For this step, we selected data in the altitude range between 200 and 300 km. This demonstrates that a simple power law allows to describe the spectral behaviour of haze extinction between 1 and 2 µm.

A systematic study allows us to draw, for each observation and for each fractal dimension D_f, a χ² map indicating which aerosol parameters, κ and N, yield the best agreement (Fig. 6). A wide range of fractal dimensions potentially allow for fitting data with the same goodness, and the corresponding fits are indiscernible from each other (Fig. 7). The parameters N corresponding to the best fits are shown as a function of the fractal dimension for each occultation data (Fig. 8a). For a given occultation, there is a broad set of couples N and D_f that fit the data, and their χ² are statistically equivalent. Therefore, the extinction spectrum alone does not allow us to discriminate between these different cases. We need a supplementary constraint.

Fig. 6 Maps of χ2 for fits of the haze continuum in windows as a function of the imaginary refractive index κ and of the number of grains per aggregate N for the four observations. The values of χ2 are adjusted for the minimum value of χ2 to be equal to the number of data minus the number of free parameters. This was necessary to be able to draw the 1σ error envelope, while there are differences in aerosol structure between different observations, as we later show. For clarity, we only show the contours of the 1σ error bars for fractal dimensions Df = 2, 2.1, 2.2, 2.3, 2.4, and 2.5.

5.3 Constraints on aerosol morphology with VIMS-IR and DISR

The reference case for Titan’s aerosols (Tomasko et al. 2008), featuring fractal aggregates of ≃3000 grains of radius 0.05 µm and a fractal dimension of D_f = 2, does not fit our data (Fig. 8a). This aggregate morphology is thus not consistent with observations by VIMS-IR made in the occultation mode. However, this reference case is consistent with observations of the solar aureole made by the Descent Imager Spectral Radiometer (DISR). Of course, the analysis of DISR is extremely complicated, and, because multiple scattering plays an important role, the solution is not unique and the choice of the fractal dimension has an influence on the final result. From Tomasko et al. (2008), we retain that the forward peaks of the phase functions for the reference aerosol actually fit DISR data and are representative of the real haze phase functions in the forward direction. Therefore, this provides a good way to be more selective in our analysis. We calculated the relative least-square values between a phase function P(θ) used by Tomasko et al. (2008) in the near-infrared (wavelength 1.75 µm), at scattering angles between 10° and 90°, and with phase functions computed with various grain numbers and fractal dimensions. Before discussing the results, it is important to remind the reader that for fractal aggregates, the shape of the phase function bears direct information about the fractal dimension (e.g. Sorensen 2001; Tazaki et al. 2016). Therefore, in this comparison between phase functions, we expect a much better agreement for the case D_f = 2 because this fractal dimension was initially assumed by Tomasko et al. (2008) to generate their phase functions. We then expect the agreement to worsen with increasing D_f. Our purpose is then to define, separately for each fractal dimension, the best value for the number of monomers N, regardless of the absolute goodness of the fit, which is of no importance here (Fig. 8 b).

We find that whatever the fractal dimension, the phase function derived from the analysis of DISR indicates that aggregate should have around 3000 ± 1000 grains. Such an analysis cannot lead to a strict evaluation of the best fit with error bars. This is because we do not use the initial data but a final outcome of Tomasko et al. (2008). Instead, we established a lower and an upper limit for the number of grains, as the numbers of grains with a relative least-square value that is twice the minimum least-square value. Thus, for each fractal dimension, the best value and the limits of the interval give indicative constraints. Yet, the combination of this constraint from DISR and the constraints previously obtained with VIMS-IR occultations are indeed extremely selective. It totally discards the case with D_f = 2 (Figs. 6 and 8b). It also indicates that, to be consistent simultaneously with DISR and with VIMS-IR occultations (this work), the fractal dimension of aerosols should be around 2.3 ± 0.1. The corollary consequence of this result is that fitting DISR data should be possible with another fractal dimension, by adjusting other physical quantities such as the vertical profile of the haze and the haze single scattering albedo differently.

We performed the same analysis assuming a monomer radius r_m = 0.04 µm and the reference refractive index n = 1.6 on one hand, and the reference monomer radius r_m = 0.05 µm and the real part of refractive index n = 1.5 or n = 1.7 on the other hand. A change in the monomer radius causes all our results to be shifted towards larger numbers of monomers by a factor of (5/4)². This factor corresponds to the conservation of the aggregate apparent area. Therefore, the best result is obtained for aggregates of ≃4900 monomers with the same range for the fractal dimension as before, that is, 2.3 ± 0.1. If the real part of refractive index is changed, this essentially causes the fractal dimension to be shifted according to ∂D_f /∂n ≃ 0.30, without any significant change in the number of monomers per aggregate. The change of ±0.03 in D_f induced by a change of ±0.1 in n is quite moderate and does not change our conclusion.

In order to account for the constraints given by DISR and VIMS-IR occultations, we considered the fractal aggregate characteristics D_f = 2.3, grain radius r_m = 50 nm, and number of grains per aggregate N ≃ 3160 as a reference. This is given assuming the real part of the refractive index n = 1.6. Notably, this choice is of little importance for the retrieval of the haze extinction or for the gas absorption, provided that we choose an aerosol morphology that gives the correct extinction spectral slope. The corresponding reference spectral slope, which we used in the following steps, is given by α₀ = −2.18. On the other hand, the choices for the reference aggregate structure and real part of the refractive index would have an impact on the retrieval of the imaginary part of refractive index, but this aspect is not discussed in this article.

Fig. 7 Best agreements obtained between data (scaled opacities between 200 and 300 km) and model for each observation. Good fits (blue line) can be obtained for any fractal dimension between 1.8 and 2.5, albeit with different values of N and κ. All these fits are superimposed (only cases between Df = 2 and 2.5 are reported here). The green line shows the opacity slope for aggregates with N = 3000 and Df = 2.0, as in Tomasko et al. (2008). This plot and Fig. 6 show that this structure of aggregates can be discarded.

Fig. 8 Parameters producing best fits for the four observations and for different fractal dimensions Df between 2.0 and 2.5. (a) Values of N obtained for the best fits of the spectral opacity (Fig. 7) as a function of Df. The constraints on N retrieved from the fit of the phase functions from DISR analysis are shown with black dots and dashed curves (reported from the graph (b)). (b) Relative least-square values for comparisons between phase functions derived with DISR (Tomasko et al. 2008) and phase functions computed for aggregates of various fractal dimensions and number of grains.

6 Retrieval procedure for the haze and gas vertical profiles

To go further in the analysis, we developed a retrieval model following the technique of the Bayesian inference (Rodgers 2000). Using the Bayes theorem, this technique allows for retrieving a parameter set x along with an estimated error σ_x from a set of observations y with observation errors σ_y. Each set is represented by a vector. The strength of this approach essentially relies on the fact that we only need knowledge of the direct relationship between y and x. This relation is formally given by y = F(x), where in our case F stands for the radiative transfer model that links observed transmissions (y) and the atmosphere parameters (x) such as the haze extinction or gas mixing ratios.

To proceed, the function F has to be linearised around a given solution, ${y\left(x\right)=y_0+\partial F\left(x\right)/\partial x|}_{x_0}\times \left(x-x_0\right)$, and where y₀ = F(x₀). In this way, the set of partial derivatives ${\partial F\left(x\right)/\partial x|}_{x_0}$ yields a matrix K, and the estimation for the best value of $\overline{x}$ can be written as follows (Rodgers 2000):

$$\overline{x}=x_0+{\left(K^T{S}_{ϵ}^{-1}K\right)}^{-1}{K}^T{S}_{ϵ}^{-1}\left(y-y\left(x_0\right)\right),$$(9)

where S_є is the diagonal error matrix $\left(S_{ϵ,\left(ij\right)}={\delta }_{i,j}/{\sigma }_{y_i}^2\right)$, $S={\left(K^T{S}_{ϵ}^{-1}K\right)}^{-1}$ is the covariance matrix, and y is the observation vector. If x₀ is far from the best solution, the first estimated value $\overline{x}$ will not be the best possible solution, but it will be better than x₀. So, we can start an iterative process to converge towards a most likely solution. In general, 20–40 iterations are needed to obtain an acceptable solution. It depends on the size of the vector x. Notably, although we start with an initial guess for x₀, we do not include any a priori value in the retrieval matrix. The errors on the retrieved parameter x are given by the square root of the diagonal terms in the covariance matrix ${\sigma }_{xi}=\sqrt{S_{ii.}}$.

In order to avoid instabilities in the converged solution, we used a Tikhonov regularisation. This consists of adding a process analogous to a diffusion in the x-space to attenuate spurious oscillations. It takes the form of a matrix, H, inserted in the covariance matrix as

$$S^{†}={\left(K^T{S}_{ϵ}^{-1}K+\beta H\right)}^{-1}={\left(S^{-1}+\beta H\right)}^{-1}$$(10)

where H is defined from the second-order derivative matrix L as H = L² (e.g. Quémerais et al. 2006; Koskinen et al. 2011) and β is analogous to a diffusion coefficient. The order of magnitude of the matrix terms in L is around (∆Z)⁻², with ∆Z being the vertical distance between two consecutive observations. Setting β is a matter of tuning since there is no related physical process behind. We scaled β against the larger term of S⁻¹ and defined a free factor γ, giving us $\beta =\gamma \text{Δ}{Z}^4\times \max \left(S_{i,j}^{-1}\right)$. We tested several values of γ around 1 in order to remove small-scale oscillations but keep the variations comparable to the atmosphere scale height (≃40 km). We found that γ ≃ 0.316 gives the desired effect for retrieval performed with VIMS-IR data at wavelengths shorter than 3 µm. Beyond this wavelength, data are much noisier, and as a direct consequence we needed to use a larger value such as γ ≃ 3.16 to obtain smooth results at length consistent with the atmospheric scale height.

7 Results and discussions

7.1 Haze vertical profile

Occultation data contain information about haze in methane windows. We used the VIMS-IR channels corresponding to methane windows from 0.88 up to 2.8 µm. The observation vector y is formed by the observed radiance factors I/F recorded in the N_Z levels where we expect to have information. Each level contains a certain number of channels N_D(i = 1, N_Z) from which we use the I/F value. Then, the size of the observation vector y is $N_y={\displaystyle {\sum }_{i=1,N_Z}{N}_D\left(i\right)}$. For haze, we decided to use two parameters for each level, the extinction at a reference wavelength λ₀ = 1.4 µm, and the spectral slope ∆α = ∂ ln(k_ext)/∂ ln(λ) − α₀. The slope α₀ = −2.18 is the slope determined from the previous step about aerosol morphology (aggregate of ≃3160 monomers of fractal dimension D_f = 2.3). Thus, we only seek variations around this reference value. The parameter vector x therefore has 2 × N_Z values.

We used the iterative Eq. (9) and obtained convergence on the parameter x with about 40 iterations. It should be noted that the retrieved values in the few uppermost layers depends on the assumption made for atmosphere above the first level. We assumed that the haze extinction decreases with a simple scale height H = 50 km at high altitudes.

The haze extinction profiles retrieved from the four observation sets are shown in Fig. 9. They present some noticeable structures. The observation T10E was taken three years before equinox in the southern circumpolar region and during a period of steady state for haze, and T53E was taken a few months before the equinox when circulation had just started its turnover (e.g. Seignovert et al. 2021). For the observation T10E, the main haze layer has an almost constant vertical scale height extending up to 400 km. For T53E, the top of the main layer seems to be located at ≃320 km, with a sharp cut-off of the main haze layer. The observations T78I and T78E were taken in 2011, two years after equinox and during the collapse of the detached haze layer (West et al. 2011, 2018; Seignovert et al. 2021). The detached haze layer clearly appears at 335 ± 10 km, consistently with a previous observation made by the Imaging Science Subsystem (ISS) onboard Cassini and almost at the same altitude as when it was observed by ISS onboard Voyager 2 on Titan one year before (e.g. Rages & Pollack 1983). A secondary feature can be observed around 200 km in observations T53E, T78E, and T78I. This feature does not appear, or it is very subtle, in the T10E observation. With only these four profiles, we cannot draw a firm conclusion about its origin. We think that it could be a relic of the previous detached haze, which settled after the previous equinox in 1995.

The retrievals of the spectral slopes are displayed in Fig. 10. They have noticeable variations around the reference value. We may associate changes in spectral slopes with variations in aerosol properties. The reference spectral slope is α₀ = −2.18, and variations of this slope with D_f and N are, respectively, ∂α/∂D_f ≃ 0.71 and ∂α/∂ log(N) = 0.106. A change of ±0.2 in the spectral slope, as plotted in Fig. 10 for T10E and T78E, can be produced by a change of ±0.28 in D_f or by multiplying or dividing N by a factor ≃6.5, or a combination of the two effects. For T78I, this change is even more pronounced, and ∆α = −0.4 would mean a change of −0.56 in the fractal dimension D_f, a decrease by a factor 42 in N or a combination of the two. Such variations in the fractal dimension cannot be due to a change in aerosol growth mode in the main haze, because aerosols are produced by long term processes. For instance, the growth mode for Titan aerosols is thought to be comparable to the theoretical ballistic cluster-cluster aggregation mode. This produces aggregates with a theoretical fractal dimension ≃2 (Cabane et al. 1992), and the production of aerosols of several thousands of grains should exceed a Titan year.

The observation T10 was made in the summer southern hemisphere, where the atmosphere was submitted to strong ascending winds (Rannou et al. 2004; Lebonnois et al. 2012). A variation in the fractal dimension or in aerosol size in all the stratospheric column may be explained by a dynamical segregation in ascending motions since the settle speed of aerosols strongly depends on their compacity (i.e. fractal dimension) and on their size, as proposed by Larson et al. (2015). This would promote, in upwelling motions, aerosols with low fractal dimensions rather than aerosols with large fractal dimensions and small aggregates rather than a large one. We may be witnessing such an effect acting on both the fractal dimension and the number of grains simultaneously. Larson et al. (2015) explicitly accounted for the change in fractal dimension of aerosol that they correlated with aerosol size. Unfortunately, they did not explore the related effect on the aerosol distribution during different seasons in detail. However, this could be the key to understanding why haze vertical profiles in global climate models (GCMs) are so different from observations. We also note that a strong variation in the spectral slope, but not in the extinction coefficient, occurs in the low stratosphere below 150 km, which indicates a specific process strongly altering the properties of the haze layer in some way. This feature appears correlated with the humid layer reported by Rannou et al. (2021) and by us in the next section; it is attributed to an intrusion of air from the tropical region. We note that these results for the haze layer obtained with T10E can be compared to those retrieved by Bellucci et al. (2009). They found comparable vertical profiles, and they also reported a change in the aerosol spectral slope but did not gave explanation for it.

The three other observations were taken at moderate latitudes around the equinoctial transition. With T53E, observed before the equinoctial turnover, the aerosol spectral slope is fairly constant over all the altitude range up to 350 km. The sharp increase in spectral slope above this level is related to the sharp drop in extinction that marks the transition between the main haze and the detached haze layer. For T78E and T78I observations, two years after the northern spring equinox, the spectral slope appears fairly constant with altitude up to about 250 or 270 km and then has strong variations apparently related to the transition between the main haze and the detached haze layers. These variations appear quite complex. They are consistent with a smaller value of α inside the detached haze (indicating smaller aerosols) relatively to the main layer. The peak at 300 km in the two profiles, T78E and T78I, corresponds to the zone that is depleted in aerosol and located between the main and detached haze. This peak, considered along with the peak at 350 km retrieved with T53E, shows that aerosols in the depleted layer may have a very specific population compared to the main haze below and the detached haze above. The detached haze layer is a zone where aerosols of different origins co-exist with size distributions different to those in the main haze layer (Rannou et al. 2004; Cours et al. 2011; Larson et al. 2015) and where rapid dynamical processes occur (Porco et al. 2005; Seignovert et al. 2021). The sharp changes in the aerosol spectral slopes indicate sharp transitions in the aerosol distributions. Both their sizes and their fractal dimensions play a role in this case.

Fig. 9 Haze extinction as a function of altitude, retrieved for the four observations, at wavelengths 0.884 µm (channel #97), 1.540 µm (channel #137), and 2.199 µm (channel #177). The extinction profiles retrieved by Seignovert et al. (2021) with ISS onboard Cassini at wavelength 338 nm (CL-UV3 filters) are shown with green lines (labelled S2020). Those from Vinatier et al. (2010) or Vinatier et al. (2015), scaled at the wavelength 1µm, are shown with black lines (V2010 or V2015). The profiles in cyan (RP83) are the extinctions retrieved by Rages & Pollack (1983) at 30°N in August 1981 (wavelength 0.5 µm). The differences in the detached haze altitudes between VIMS-IR (Ls = 2б°), ISS (Ls = 14.8°) and ISS onboard Voyager 2 (Ls = 18°) are their dates while the detached haze layer is falling down (West et al. 2018; Seignovert et al. 2021). The grey line shows the haze profile by Doose et al. (2016) with DISR in 2005 at 10° S (labelled D2016).

Fig. 10 Difference in spectral slope of haze extinction, ∆α = α − α0, where α = d ln(kext)/d ln(λ)) and the reference slope is α0 = −2.18, as a function of altitude retrieved with the four observations.

7.2 Methane vertical profile

In a second step, we analysed the methane abundance by fitting the methane bands in the range between 0.88 and 2.5 µm. We used the methane line list published by Rey et al. (2016, 2018), and at this step we only set the methane mixing ratio as free parameter. The haze properties were fixed according to our previous analysis. The parameter vector x has now the size N_Z and describes the vertical profile of the methane mixing ratio. This time, the observation vector y includes all the channels in the studied spectral range, except the channels between 1.60 and 1.65 µm, which give spurious measurements (The VIMS Wavelength and Radiometric Calibration 19; final Report, NASA Planetary Data System, hereafter RC19¹). We used Eq. (9) to converge towards a solution, and we used the same rule to define the parameter γ for the Tikhonov regularisation.

The wavelength associated with VIMS-IR channels has changed with time (Sromovsky et al. 2012). This is obvious thanks to the sharp signatures of gaseous absorptions. A complete study of the VIMS-IR calibration was performed, and the last version of the calibration, RC19, gives the wavelength shift as a function of time. The methane lines are now extremely well defined thanks to recent works about methane spectrometry (Campargue et al. 2013, 2016; Rey et al. 2016, 2018), and we are able to define the wavelength shift with errors of about 0.1 to 0.2 nm depending on the quality of the observation (Table 1). We evaluated the shift using the two bands that we modelled in full at 1.18 and 1.4 µm. We did not model the core of the 1.60 µm band, and the 2.3 µm has an unknown absorption. However, we observe no significant difference (that is, within the error bars) when using these two bands separately. The wavelength shift changes with time and our estimations are very close to those already found by Sromovsky et al. (2012) and Maltagliati et al. (2015), and they significantly differ from RC19 values by a factor of 1.5 to 2 depending on dates. Although this wavelength shift mainly affects the retrieval of gas abundances, we also accounted for it in all the retrieval procedures concerning the haze presented in the previous sections.

Figure 11 shows the retrieved vertical profiles for the four observations in the two wavelength intervals: between 0.88 to 2 µm (top) and between 2.0 and 2.8 µm (bottom). Figure 12 shows the corresponding fits of the data. From our results, we immediately remark that the 2.3 µm band cannot allow us to safely retrieve methane abundance, because another species strongly interacts with methane absorption. This was already well known for the C-H fundamental band at 3.4 µm, where methane does not absorb enough to explain the data. As suggested by Maltagliati et al. (2015), ethane could be the absorbing gas, although it would appear surprising that a trace gas could compete with methane in one of its bands. If ethane is not responsible for the extra absorption in the 2.3 µm band, it would mean that something is not yet understood in Titan’s atmosphere opacity. For now, we consider that retrieval of methane with the 2.3 µm band is not reliable, and it is only shown for information purposes.

With the methane bands between 0.88 and 2.00 µm, we find four valuable methane mixing ratio profiles from about 100 km up to 350 or 400 km. The profiles derived before the northern spring equinox (that is, T10E and T53E) were already analysed in Rannou et al. (2021), and we simply remind the reader here of the main conclusions of their studies. We also discuss the retrieval made with T78E and T78I observations.

With the T10E observation, the methane profile found at 70^°S has a sharp, methane-rich layer at 165 km, with a mixing ratio around 1 .45 ± 0.05% in a relatively dry background atmosphere. We note that this differs slightly from the previous evaluation by Rannou et al. (2021), that is, 1.65 ± 0.05%, after a spurious value in data was corrected (Figs. 11 and 12). The mixing ratio is quite constant (around ≃1.0%) in the stratosphere above about 190 km. Below the humid layer, it drops to a minimum value lower than 0.8% at around 100 km. This is consistent with the results obtained by Lellouch et al. (2014), which found 1.1 ± 0.09% at 60°S, 1.54 ± 0.19% at 80°S, considered as an average over the altitude range of 125 ± 50 km.

We considered the profile inferred with T53E (near the equator and before the equinox) along with results retrieved with Huygens and CIRS. Our methane profiles are consistent with retrievals made by Lellouch et al. (2014) with CIRS at 125 and 225 km, showing a low value of methane abundance. We note that results found with CIRS are associated with a broad altitude range (about ±50 km) due to a low spatial resolution. Two independent retrievals of the methane mixing ratio are also made with the Gas Chromatograph and Mass Spectrometer (GCMS; Niemann et al. 2010) and with DISR (Bézard 2014; Rey et al. 2018), both onboard Huygens. They both find a methane mixing ratio of approximately 1.4% around the tropopause, and they diverge above it. The first analysis with DISR gives about 1.4–1.5% up to to 100 km and a higher and uncertain value around 1.9% at 110 km (Bézard 2014). A quite constant value of 1.4% is found up to 140 km with GCMS (Niemann et al. 2010). The most recent analysis of DISR data, performed with the latest methane spectroscopic data, yields a decreasing methane mixing ratio above the tropopause, which reaches a low value around 1.0% at 110 km (Rey et al. 2018). At this step, this seems to imply that the methane mixing ratio decreases from 5 to 1.4% due to cloud activity and then further decreases from 1.4% to reach a low value of about 1.15% around 200-220 km.

Rannou et al. (2021) compared the atmospheric circulation pattern, the vertical profiles of the methane mixing ratio obtained with T10E and T53E and mixing ratios in the stratosphere and at the cold trap derived from CIRS observations (Lellouch et al. 2014). They came to the conclusion that the methane-rich layer at 70°S and at 165 km (T10E) has a mixing ratio much larger than the that at the cold trap just below. On the other hand, the circulation patterns show that this layer is directly connected (via a circulation stream) to the summer tropical convergence zone, where convective clouds are frequently observed (Griffith et al. 2005; Rodriguez et al. 2009, 2011; Turtle et al. 2011) and where clouds are also modelled by GCMs (e.g. Rannou et al. 2006; Lora et al. 2015). These ascending motions at the summer convergence zone continue vertically in the stratosphere. This stream diverges in the stratosphere, but a part of the stream goes up to the polar region where the methane-rich layer is observed. Eventually, convective clouds could also inject methane ice directly above the cold trap or directly into the stratosphere (Barth 2010; Barth & Rafkin 2010), and, therefore, they could send the methane into the stratosphere with a mixing ratio that exceeds the value expected from the cold-trap limitation.

On the other hand, near the equator (T53e), the larger mixing ratio retrieved with DISR (Bézard 2014; Rey et al. 2018) and with the GCMS (Niemann et al. 2010) at the tropopause than in the mid-stratosphere shows that the methane in the stratosphere is not directly linked to its value at the tropopause. This was also noticed by Lellouch et al. (2014). This shows that methane is not mainly transported at a global scale through the cold trap and then through the tropopause. Instead, there are some locations, probably around tropics and maybe elsewhere depending on the season, where methane is injected massively and then dispersed by the global circulation. The relative flux of methane reaching the stratosphere by a global slow transport through the tropopause and by intrusion at the tropics cannot be known from existing data.

The GCM by Rannou et al. (2006) accounts for the haze and cloud microphysics and their radiative feedbacks, but not for the thermodynamical feedback or for the corresponding convective processes. Although this model is able to reproduce the main features of the cloud cover, especially the mid-latitude clouds at the tropical convergence zone (TCZ), it does not yield any sharp constrast in the stratospheric methane mixing ratio. Instead, it yields a uniform methane mixing ratio in the stratosphere with a value very close to the mixing ratio at the tropopause between tropics and at mid-latitudes. We then infer that cloud convection, which is lacking in this model, plays a crucial role in the methane cycle and is an important player for the stratosphere humidification.

In this work, we analysed two more profiles, the T78E and T78I observations obtained in 2011, two terrestrial years after the north spring equinox. At that time, models predicted that the seasonal circulation turnover had started (Lebonnois et al. 2012), and several post-equinoctial changes can be noticed (e.g. West et al. 2018; Teanby et al. 2019; Seignovert et al. 2021; Mathé et al. 2020; Turtle et al. 2018). The methane vertical profiles are fairly constant with altitude, and they also show a low methane mixing ratio. This reinforces the idea of a low methane mixing, around 1.1% at the global scale in the stratosphere. However, there are noticeable differences between the two profiles taken in 2011. First, the T78E observation is taken at 27°N, and the T78I is taken at 40°N. Between these two latitudes, the stratospheric circulation may differ. The profile obtained from T78E is almost constant from 100 to 360 km, with small oscillations. With T78I, the methane mixing ratio have a slight increase with altitude and has marked oscillations between 150 and 250 km. The marked structure below 250 km may be the signature of a complex circulation pattern where the interaction between different streams occurs. Similar differences between profiles of gas trace species were observed at these two latitudes and at the same point in the season (Vinatier et al. 2015).

The vertical profiles retrieved by CIRS in the northern hemisphere (Lellouch et al. 2014), reported in Fig. 11, also increases with the altitude, while those retrieved in the southern hemisphere do not. They found a methane mixing ratio as large as 1.88 ± 0.44% at 30°N (observation 28 January 2010), and with VIMS-IR we find a much smaller and constant vertical profile at about the same latitude in September 2011. This could be the signature of a strong injection by an equinoctial intrusion of methane in the stratosphere due to a strong cloud activity (similar to the one reported by Turtle et al. 2011) transported northward and then dissipated by the circulation. It cannot be excluded that cloud events and intrusions occur preferentially at certain longitudes (Rodriguez et al. 2009) or that major cloud events near the equator during the circulation turnover are sufficiently rare that an excess in methane mixing ratio is not zonally homogeneous.

Methane and aerosols are tracers of the circulation. They obviously have different sources and behaviours in the atmosphere. Thus, to fully understand the methane cycle, these methane vertical profiles in the stratosphere must be analysed with a climate model accounting for the cloud convection and the transport of species. The profiles retrieved in this work, along with other retrieved values, will be strong constraints for climate models.

Fig. 11 Methane mixing ratio retrieved with the four observation sets, with data between 0.88 and 2 µm (top) and between 2 and 2.8 µm (bottom). We also plot the methane mole fraction retrieved with the GCMS onboard Huygens (Niemann et al. 2010) and with DISR (Bézard 2014) and CIRS (Lellouch et al. 2014). The green dashed profile in the upper left graph shows the evaluation made by Rannou et al. (2021).

Fig. 12 Comparison between observed transmissions (red dots) and the model (blue dots) in the spectral range used to retrieve methane mixing ratio. The spectra between 1.6 and 1.7 µm were not used in this study due to unreliable spectels (RC19). The green dashed lines for the transmission of T10E show the spurious point that was modified.

7.3 CH₃D vertical profile

In this step, we analysed the deuterated methane abundance by fitting the methane bands in the ranges between 2.77 and 2.91 µm (range R1) and 4.22 and 4.49 µm (range R2) (Fig. 2). In these spectral ranges, haze also participates in extinction. The two wavelength ranges R1 and R2 are far from each other, and haze absorption peaks are expected in the range around 3 µm to 4 µm, that is, between the intervals R1 and R2 (e.g. Sciamma-O’Brien et al. 2017). Thus, we cannot assume that the haze extinction spectra follow simple power laws from λ₀ = 1.4 µm up to 4 µm. Instead, we assume that the haze spectral properties already retrieved can be extrapolated up to the range R1. However, for the range R2, we account for a ratio of r_4.4, retrieved for each altitude, between haze extinction at 4.4 µm and 2.8 µm (i.e. r_4.4(z) = k_ext(λ = 4.4 µm,z)/k_ext(λ = 2.8 µm,z)). This allowed us to fix the haze in R1 and to the retrieve haze extinction in R2 separately.

We then retrieved the ratio r_4.4 and the CH₃D mixing ratio simultaneously. The parameter vector x thus has the size 2 × N_Z and contains the ratio r_4.4 and the vertical profile of the CH₃D methane mixing ratio. The vector y includes the observed radiance factor data at the selected wavelengths except the channels between 1.60 and 1.65 µm. We used Eq. (9) to converge towards a solution, and the parameter for the Tikhonov regularisation is set to γ = 3.16 because data at wavelengths larger than 3 µm are much noisier than at shorter wavelengths (RC19). A stronger smoothing function is needed to erase spurious structures on small scales. The vertical profiles of the CH₃D mixing ratio are given in Fig. 13, and the corresponding spectra are displayed in Figs. 14 and 15.

The spectral features related to CH3D are quite moderate at 2.8 at 4.4 µm. A single marked peak emerges from haze extinction at 4.52 µm. However, this structure is strongly affected by VIMS-IR instrumental noise beyond 3 µm. Thus, the error bars for this species are quite large and our analysis does not bring more information than that already known from previous works. We find a mixing ratio in the range of 1–10 ppm, with large error bars. This is consistent with previous estimates that yielded between 4 and 8 ×10⁻⁴ times the methane abundance (e.g. Coustenis et al. 2003; Penteado et al. 2005; Penteado & Griffith 2010; Bézard et al. 2007; Nixon et al. 2012).

Fig. 13 Deuterated methane mixing ratio retrieved with the four observation sets from spectra between 2.78 and 2.90 µm and 4.23 and 4.50 µm. The haze extinction is fixed thanks to the spectra in methane windows between 0.88 and 2.03 µm and adjusted between 4.23 and 4.50 µm.

Fig. 14 Comparison between observed transmissions (red dots) and the model (blue dots) in the spectral range used to retrieve the deuterated methane mixing ratio for the first spectral feature at 2.85 µm. The haze extinction is set with spectra between 0.88 and 2 µm (Figs. 7 and 12).

7.4 CO vertical profile

As for CH₃D, we simultaneously analysed the haze extinction and the CO gaseous component. Between 4.55 and 4.95 µm, transmission is sensitive to extinction by CO and haze, and beyond 4.95 µm only haze produces extinction. This latter part of the spectral range allows us to retrieve haze extinction at 5 µm independently from CO. The parameter vector x has the size 2 × N_Z and describes the vertical profile of haze extinction and the CO mixing ratio. The observation vector y includes the observed radiance factor data at the selected wavelengths. We set the parameter γ = 3.16 for the Tikhonov regularisation because of the noisy data.

The retrieved vertical profiles and the corresponding fits are displayed in Figs. 16 and 17. These results are consistent with a constant values of CO in the stratosphere with a value of 30 to 50 ppm. From the lifetime of CO (500 Myr, Lellouch et al. 2003) and because it does not condense in Titan’s atmosphere, a constant vertical mixing ratio is expected. Thus, the structures of the vertical profile (vertical slope for observation T78I and oscillations) that we find are very uncertain because the spectra are noisy. However, these values are comparable to previous retrievals made by other groups with different observations (32 ± 10 ppm by Lellouch et al. 2003, 2004, 60 ppm by Lopez-Valverde et al. 2005, 52 ± 12 by Gurwell & Muhleman 2000, and 47 ± 8 ppm by de Kok et al. 2007). With the same dataset as ours, Maltagliati et al. (2015) also found similar, although not identical, vertical profiles, because the treatment of the haze and the gas absorption are not identical.

Fig. 15 Comparison between observed transmissions (red dots) and the model (blue dots) in the spectral range used to retrieve the deuterated methane mixing ratio for the second spectral feature, between 4.25 and 4.45 µm. The haze extinction is set with spectra between 0.88 and 2 µm and adjusted with a factor that accounts for the ratio between the haze opacity of 4.4 µm and 2.8 µm (Figs. 7 and 12).

Fig. 16 Vertical profiles of CO mixing ratio retrieved with the four observation sets and using data between 4.55 and 5.2 µm. The haze extinction is fixed by fitting the data with the spectral range between 4.85 and 5.0 µm (Fig. 17).

7.5 Haze extinction spectra at 5 µm

Analysis of the CO band and the 5 µm window offers an opportunity to probe the haze layer down to the low stratosphere, just above the tropopause. Figure 18 shows the haze profiles at a short wavelength, along with the haze retrieved at 5 µm. Observations at long wavelengths are much noisier and the corresponding errors are large (Fig. 17); however, the structures of the vertical profile found at 5 µm are consistent, within the error bars, with the profiles found at shorter wavelengths. The extinction profile for T53E is quite different from the profiles retrieved with DISR (Tomasko et al. 2008; Doose et al. 2016). Of course, the T53 observation was taken in 2009, four years after the Huygens descent in July 2005. In 2009, we find a transition in the haze profile around 80 km, as with DISR, but occuring in a very different way. This change could be linked to interactions between the haze and condensing species. At these altitudes, methane is largely undersaturated, and only minor species can reach saturation above the tropopause (e.g. Barth & Toon 2006).

Observations T10E, T78I, and T78E also suggest a change in haze profile at around 80 km, but not in a significant way regarding the error bars. The T78I and T78E observations are taken in the mid-northern hemisphere almost at the end of the equinoctial circulation turnover. Cloud activity is quite weak at that time (Turtle et al. 2018), and interaction with minor species may be sensitive to the changing circulation. Unfortunately, the region just above the tropopause cannot be probed with VIMS-IR at that distance due to the refractive effect. Occultations at shorter distances, if they exist, may improve the retrieval in this region.

Fig. 17 Comparison between observed transmissions and the model in the spectral range used to retrieve the CO mixing ratio (between 4.55 and 4.85 µm) and the haze extinction (between 4.85 and 5.0 µm). For clarity, we only plotted selected spectra and the corresponding fits.

Fig. 18 Same plots as in Fig. 9, but for short wavelengths (SW, thin red and blue lines and dots. From the right to the left, the wavelengths are 0.884 µm (channel #97), 1.540 µm (channel #137), and 2.199 µm (channel #177) along with the haze vertical profile at 5 µm (LW, thick red and blue lines and red squares). The red dots and squares are the retrieved values, and the blue lines show the 1σ envelope.

8 Retrievals in nadir viewing with a low methane abundance

As stated in the introduction, several works have used radiative transfer modelling to retrieve surface reflectivities or information about the aerosol layer. For these works, the methane mixing ratio in the stratosphere retrieved with Huygens instruments (e.g. Niemann et al. 2010) generally prevails. However, it now appears that the methane mixing ratio gradually decreases from about 1.4% at the tropopause (Niemann et al. 2010; Bézard 2014) to about 1–1.2% in the stratosphere, as found in this work and in previous works (e.g. Lellouch et al. 2014; Maltagliati et al. 2015).

In order to assess the consequence of a different amount of methane on the retrieval of the aerosol opacity and on the surface reflectivity, we carried out two simple tests. We selected a bright and a dark pixel from the region already studied by Coutelier et al. (2021), and we performed a retrieval of the column opacity of aerosols (photochemical haze above 55 km and condensate mist below 55 km) and of the surface reflectivity in methane windows between 0.88 and 2.6 µm. The model of Coutelier et al. (2021) represents, for our test, the reference case. With this reference case, the methane mixing ratio in the stratosphere is 1.4%. In the first test, we simply decreased the methane mixing ratio in the stratosphere to 80% of its initial value. This is approximately the methane mixing ratio retrieved with the T53E observations, that is, 1.12%. In the second test, we wanted to account for the extra absorption detected in the 2.3 µm methane band. We then used the reduced methane mixing ratio, 1.12%, in the first bands up to 2.0 µm. However, we also increased the methane absorption in the 2.3 µm band to mimic the absorption pattern observed for T53E (Fig. 11 – bottom). To do so, we simply multiplied the methane absorption by a Gaussian function f(z) = 1.25 + 1.10 × exp (−(z − 190)²/(2 × 65²)), where z is the altitude in kilometres. In the third test, we used the reference case again, but with the temperature profile related to the observation T53E, as reported in Fig. 3, instead of the temperature profile retrieved with Huygens (Fulchignoni et al. 2005), as used by Coutelier et al. (2021).

We first note that the fits of the data with different gas absorption set-ups cannot be distinguished from each other (Fig. 19 – top). The amount of photochemical haze and condensated mist needed to match the data are given relatively to the haze and mist column opacity given by Doose et al. (2016) at 1 µm with two parameters, F_h and F_c (Table 2). To reach the same fits, the amount of haze has to be decreased, while the amount of mist in the lower layers must be increased. This is because outgoing photons in methane bands come from the low stratosphere and are insensitive to the amount of mist (F_c). Thus, less stratospheric methane implies less haze to reach the same radiance factor in the cores of the methane bands. Outgoing photons in the wings of the windows come from much deeper layers, typically in the troposphere, and radiance factors are then sensitive to both F_h and F_c. For these probed depths, changing the methane mixing ratio in the stratosphere without changing its value in the troposphere weakly affects the integrated column opacity of methane because lower layers are much more massive. Thus, to compensate for the decrease in haze opacity at wavelengths where methane absorption has barely changed, the opacity of the mist has to be increased. We find that dimming the methane mixing ratio to 80% of its initial value requires a decrease in the haze opacity and an increase in the mist opacities by several percent. The significance of these results can be estimated relatively to the uncertainties on the retrieved values of F_h and F_c, around 2.5 and 3%, respectively, as evaluated by Coutelier et al. (2021). We note that the uncertainties on the relative variations are about twice the uncertainties on the values themselves. In windows, methane is weakly absorbed, and thus retrieval of the surface reflectivity only depends on haze and mist opacity. Increasing haze and decreasing mist opacity turn out to be quite neutral on the retrieved values of the surface albedo (Fig. 19 – bottom).

For the second test, adding opacity in the 2.3 µm band (relatively to the first test), logically compensates for the need to increase the haze opacity band. The adjustment for the mist opacity is the opposite of what it was in test #1, but only on a virtually negligible percentage level. For the third test, the corrections are much more moderate than for tests #1 and #2. The changes in F_h and F_c can be considered negligible. The exact correction to bring to the haze and the mist opacities will depend on the exact nature of the missing absorption, but we find here that these corrections should be quite moderate, although not negligible. However, one has to keep in mind that the actual methane profiles are more complex than in this test. The methane mixing ratio is not uniform in the stratosphere; its vertical profiles have quite marked structures that probably change with time.

Fig. 19 Sensitivity of model fits and retrieved surface albedos to methane-mixing ratio and temperature vertical profiles. Top: radiance factor observed in the region of a possible tropical lake (Griffith et al. 2012; Coutelier et al. 2021) over a bright and a dark area. The lines are for the results of the retrieval model, used with the reference set-up and with the modified gas absorption for tests #1 and #2 and with the modified temperature for test #3. The different fits of each spectrum, obtained with different model setups, cannot be distinguished from one another. Bottom: retrieved surface albedo for the bright and the dark pixels for the reference model (black dots) and with the modified gas absorption for tests #1 (red dots), #2 (blue dots), and #3 (pink dots). For clarity, we display the error bars for the reference model only. The error bars obtained with the other models are very similar.

9 Conclusion

In this study, we analysed four observations made by VIMS-IR in occultation mode. Although this is a small set of observations that were already studied in the past (Bellucci et al. 2009; Maltagliati et al. 2015; Cours et al. 2020), they nevertheless yield new and interesting results about the haze and gases within the framework of this study. We clearly find that the spectral shift of VIMS-IR channel changes with time differently to what was published in the official calibration made for VIMS-IR. We find a shift of 2.64 nm in January 2006 instead of 1.1 nm, 9.40 nm in April 2009 instead of 6.0 nm, and 10.44 and 10.63 nm in September 2011 instead of 7.2 nm (Table 1). This is consistent with previous evaluation by Sromovsky et al. (2012) and Maltagliati et al. (2015).

We found that aerosol morphology is more compact than previously thought, with a fractal dimension D_f = 2.3 ± 0.1 instead of ≃2.0 (Cabane et al. 1993; Tomasko et al. 2008) and with a number of monomers per aggregate, around 3000, which is comparable to previous estimates. This has an impact on the photometric properties (Coutelier et al. 2021), and we also expect significant impacts on the mechanical properties of the aerosols (Cabane et al. 1993; Larson et al. 2015). More compact aggregates may be segregated in size by large-scale ascending circulation, mainly in the stratosphere. This probably explains differences between GCMs and observations concerning the vertical profile of haze and some aspects of the detached haze layer, as noted previously (e.g. West et al. 2011, 2018; Seignovert et al. 2021).

We find various haze profiles depending on latitude and time. These vertical profiles are consistent with other works (e.g. Rages & Pollack 1983; Seignovert et al. 2021; Vinatier et al. 2010, 2015). We clearly see the detached haze layer around an altitude of 320 km, consistently with the findings of West et al. (2018) and Seignovert et al. (2021). With our analysis, we are also able to characterise the spectral properties of the haze layer. We find a strong variation of the haze properties in ascending motions in the polar region (70°S) during the 2006 (T10) observation. The spectral slope decreases with altitude, consistently with a decrease in aerosol size or a decrease in fractal dimension. This is fully consistent with a segregation of the aerosols based on their settle speeds, which allows smaller and fluffier aerosols to reach higher levels. This would not be possible with aggregates of a unique fractal dimension D_f = 2, because they all have the same settle speed in the stratosphere (Cabane et al. 1993). The spectral slope of the haze also changes significantly below and inside the detached hazes observed in 2011 (T78E and T78I) and in 2009 (T53E). For the latter observation, we essentially probed the main haze layer, and the very top of the layer appears to be related to the transition between the main and the detached haze layers. The aerosol distribution in the detached haze is known to differ from the distribution in the main haze (Rannou et al. 2004; Cours et al. 2011; Larson et al. 2015). The changes that we observe here are consistent with this statement.

We were able to retrieve vertical profiles of methane. In theory, we could use the four methane bands between 0.8 and 2.8 µm. The absorption band around 3.4 µm is known to be far too thick to be explained by methane only, and we have long known that one or more gaseous or solid compounds should absorb light in this band (Bellucci et al. 2009; Maltagliati et al. 2015; Cours et al. 2020). We found that methane bands between 0.8 and 2.0 µm can be used to determine a methane profile, but the 2.3 µm band is also too optically thick to be explained by methane alone. The excess of opacity cannot be attributed to aerosols, and thus another unknown extinction or absorption must exist. The three bands that can be used for methane retrieval give a very accurate estimation. We find that the methane mixing ratios in the stratosphere are significantly smaller than the canonical value for the low stratosphere given by the GCMS (Niemann et al. 2010) and by DISR (Bézard 2014), both onboard Huygens. On the other hand, our values are consistent with those retrieved by Lellouch et al. (2014) and Rey et al. (2018). All these values and our results are consistent if we consider that the methane mixing ratio around the equator decreases from ≃1.4% above the tropopause to about 1.1–1.2% around 200 km. The observation T10, made in 2006 at 70°S, also shows that methane mixing ratio can be strongly modulated by circulation. The excess in methane mixing ratio at ≃165 km is explained by the stratosphere circulation driving warm and methane-rich air from a tropical region to dryer and cold regions (Rannou et al. 2021). The methane profiles retrieved with T78E and T78I, after the northern spring equinox, also have features in the stratosphere that clearly indicate an interaction with the atmosphere circulation.

The mixing ratios of CH₃D and CO that we retrieve are consistent with previous works. However, due to significant noise in VIMS-IR data beyond 3 µm, we find large error bars in our results. Consequently, our results reveal nothing new compared to previous results.

Two important improvements are foreseen related to this analysis. First, there are other valuable occultation sets taken by VIMS-IR that would complete this work. They will allow us to follow, with a good spatial resolution but scarce sampling, the evolution of the methane mixing ratio and the haze layer in the stratosphere. Our work shows that the occultation mode can yield extremely accurate information not available with other observation modes of VIMS-IR and with other instruments. The second expected improvement comes from the use of future line lists for C₂H₆. Many signatures not accounted for in our analysis are known to be produced by ethane (e.g. Maltagliati et al. 2015); in our analysis, this lack of information prevented us from studying the full spectra in depth. However, ethane probably cannot be responsible for the excess of absorption in the 2.3 µm band or in the 3.4 µm band. In any case, modelling ethane will allow us to safely retrieve the C2H2 and possibly the HCN mixing ratios. Moreover, it would give access to the haze extinction at around 3 µm, which is a key feature used to constrain haze chemical composition. Of course, ethane spectroscopy line lists are also of utmost importance when it comes to analysing nadir observations. This is done in combination with occultations in order to characterise Titan’s surface reflectivity in the 2 µm and in the 2.8 µm windows. Finally, it may allow us to determine the ethane mixing ratio in the troposphere, which is an important piece of information that is currently lacking.

References

All Tables

All Figures

	Fig. 1 Transmission curves as function of wavelength observed by VIMS-IR during the four flybys. For each flyby, each curve corresponds to transmissions at a given altitude between the surface and 400 km. The transmission increases with altitudes; the transmissions at 400 km are close to 1, and those near the surface are close to 0.
In the text

Fig. 2 Vertical opacity in our model in the middle stratosphere (at 200 km altitude) as a function of the wavelength for five gases and for the haze with a constant refractive index (m = 1.6 + 0.01i), The grey horizontal bars show where ethane has weak (thin bar), medium (medium bar), or strong (thick bar) absorption features (Maltagliati et al. 2015). The vertical bars delimit intervals where different components can be studied separately: (1) haze and CH4 between 0.88 and 2.8 µm; (2) haze and CH3D between 2.8 and 2.95 µm; (3) haze, CH4, C2H2, and HCN between 3.10 and 4.20 µm; (4) haze, CH3D, and CO beyond 4.2 µm. In the absence of any reliable line list for C2H6 in the spectral range of VIMS-IR, the spectral ranges where ethane has strong features cannot be studied and spectral ranges where ethane has medium strength should be studied with caution.

In the text

Fig. 3 Preliminary control of model background properties: aerosol optical cross-sections and atmosphere temperature profiles Left: extinction (upper curves) and scattering (lower curves) cross-sections for an aggregate of 3000 monomers of radius 50 nm with a fractal dimension Df = 2, computed with the mean field approximation model (Botet et al. 1997; Rannou et al. 1997) and with the model used by Tomasko et al. (2008). Right: temperature profile as a function of pressure observed with HASI (Fulchignoni et al. 2005) and profiles based on HASI with a temperature shift ±∆T (red dots). We also plotted the temperature profiles after CIRS observations (Vinatier et al. 2015; Mathé et al. 2020) that best correspond to the dates and the latitudes of the occul-tations treated in this work (Table 1). The thick black dashed line shows the temperature profile retrieved by Schinder et al. (2020) with radio-occultation measurements during the flyby T52E (4 April 2009 at lat = 12°S).

In the text

Fig. 4 Effect of atmosphere refractivity for the observation T53E, calculated with the pressure and the temperature profiles retrieved by CIRS (flyby T51, 27 March 2009 at 1°S, Distance to the spacecraft : 6300 km - Table 1). Left: vertical profile of the deviation angle caused by the refraction in the atmosphere. The angular resolution of VIMS-IR (0.5 mrad) is indicated. Middle: vertical profile of refractive attenuation for three values of the distance from Cassini to Titan. Right: altitude actually probed by the light (continuous line) and the apparent altitude assuming that the light follows straight paths (dashed lines). This altitude assuming straight light is the impact parameter of the line of sight.

In the text

Fig. 5 Altitude where transmitted flux is attenuated by 1 and 50% as a function of the distance D of Cassini to Titan. The distance for each observation is reported in Table 1. The shapes (triangles and dots) correspond to the evaluation made by Bellucci et al. (2009). The blue and red lines show our calculations using the pressure and temperature profiles retrieved by CIRS (flyby T51, 27 March 2009 at 1° S) and by HASI onboard Huygens Fulchignoni et al. (2005), respectively (figure based on Bellucci et al. 2009).

In the text

Fig. 6 Maps of χ2 for fits of the haze continuum in windows as a function of the imaginary refractive index κ and of the number of grains per aggregate N for the four observations. The values of χ2 are adjusted for the minimum value of χ2 to be equal to the number of data minus the number of free parameters. This was necessary to be able to draw the 1σ error envelope, while there are differences in aerosol structure between different observations, as we later show. For clarity, we only show the contours of the 1σ error bars for fractal dimensions Df = 2, 2.1, 2.2, 2.3, 2.4, and 2.5.

In the text

Fig. 7 Best agreements obtained between data (scaled opacities between 200 and 300 km) and model for each observation. Good fits (blue line) can be obtained for any fractal dimension between 1.8 and 2.5, albeit with different values of N and κ. All these fits are superimposed (only cases between Df = 2 and 2.5 are reported here). The green line shows the opacity slope for aggregates with N = 3000 and Df = 2.0, as in Tomasko et al. (2008). This plot and Fig. 6 show that this structure of aggregates can be discarded.

In the text

Fig. 8 Parameters producing best fits for the four observations and for different fractal dimensions Df between 2.0 and 2.5. (a) Values of N obtained for the best fits of the spectral opacity (Fig. 7) as a function of Df. The constraints on N retrieved from the fit of the phase functions from DISR analysis are shown with black dots and dashed curves (reported from the graph (b)). (b) Relative least-square values for comparisons between phase functions derived with DISR (Tomasko et al. 2008) and phase functions computed for aggregates of various fractal dimensions and number of grains.

In the text

Fig. 9 Haze extinction as a function of altitude, retrieved for the four observations, at wavelengths 0.884 µm (channel #97), 1.540 µm (channel #137), and 2.199 µm (channel #177). The extinction profiles retrieved by Seignovert et al. (2021) with ISS onboard Cassini at wavelength 338 nm (CL-UV3 filters) are shown with green lines (labelled S2020). Those from Vinatier et al. (2010) or Vinatier et al. (2015), scaled at the wavelength 1µm, are shown with black lines (V2010 or V2015). The profiles in cyan (RP83) are the extinctions retrieved by Rages & Pollack (1983) at 30°N in August 1981 (wavelength 0.5 µm). The differences in the detached haze altitudes between VIMS-IR (Ls = 2б°), ISS (Ls = 14.8°) and ISS onboard Voyager 2 (Ls = 18°) are their dates while the detached haze layer is falling down (West et al. 2018; Seignovert et al. 2021). The grey line shows the haze profile by Doose et al. (2016) with DISR in 2005 at 10° S (labelled D2016).

In the text

	Fig. 10 Difference in spectral slope of haze extinction, ∆α = α − α0, where α = d ln(kext)/d ln(λ)) and the reference slope is α0 = −2.18, as a function of altitude retrieved with the four observations.
In the text

Fig. 11 Methane mixing ratio retrieved with the four observation sets, with data between 0.88 and 2 µm (top) and between 2 and 2.8 µm (bottom). We also plot the methane mole fraction retrieved with the GCMS onboard Huygens (Niemann et al. 2010) and with DISR (Bézard 2014) and CIRS (Lellouch et al. 2014). The green dashed profile in the upper left graph shows the evaluation made by Rannou et al. (2021).

In the text

	Fig. 12 Comparison between observed transmissions (red dots) and the model (blue dots) in the spectral range used to retrieve methane mixing ratio. The spectra between 1.6 and 1.7 µm were not used in this study due to unreliable spectels (RC19). The green dashed lines for the transmission of T10E show the spurious point that was modified.
In the text

	Fig. 13 Deuterated methane mixing ratio retrieved with the four observation sets from spectra between 2.78 and 2.90 µm and 4.23 and 4.50 µm. The haze extinction is fixed thanks to the spectra in methane windows between 0.88 and 2.03 µm and adjusted between 4.23 and 4.50 µm.
In the text

	Fig. 14 Comparison between observed transmissions (red dots) and the model (blue dots) in the spectral range used to retrieve the deuterated methane mixing ratio for the first spectral feature at 2.85 µm. The haze extinction is set with spectra between 0.88 and 2 µm (Figs. 7 and 12).
In the text

Fig. 15 Comparison between observed transmissions (red dots) and the model (blue dots) in the spectral range used to retrieve the deuterated methane mixing ratio for the second spectral feature, between 4.25 and 4.45 µm. The haze extinction is set with spectra between 0.88 and 2 µm and adjusted with a factor that accounts for the ratio between the haze opacity of 4.4 µm and 2.8 µm (Figs. 7 and 12).

In the text

	Fig. 16 Vertical profiles of CO mixing ratio retrieved with the four observation sets and using data between 4.55 and 5.2 µm. The haze extinction is fixed by fitting the data with the spectral range between 4.85 and 5.0 µm (Fig. 17).
In the text

	Fig. 17 Comparison between observed transmissions and the model in the spectral range used to retrieve the CO mixing ratio (between 4.55 and 4.85 µm) and the haze extinction (between 4.85 and 5.0 µm). For clarity, we only plotted selected spectra and the corresponding fits.
In the text

Fig. 18 Same plots as in Fig. 9, but for short wavelengths (SW, thin red and blue lines and dots. From the right to the left, the wavelengths are 0.884 µm (channel #97), 1.540 µm (channel #137), and 2.199 µm (channel #177) along with the haze vertical profile at 5 µm (LW, thick red and blue lines and red squares). The red dots and squares are the retrieved values, and the blue lines show the 1σ envelope.

In the text

Fig. 19 Sensitivity of model fits and retrieved surface albedos to methane-mixing ratio and temperature vertical profiles. Top: radiance factor observed in the region of a possible tropical lake (Griffith et al. 2012; Coutelier et al. 2021) over a bright and a dark area. The lines are for the results of the retrieval model, used with the reference set-up and with the modified gas absorption for tests #1 and #2 and with the modified temperature for test #3. The different fits of each spectrum, obtained with different model setups, cannot be distinguished from one another. Bottom: retrieved surface albedo for the bright and the dark pixels for the reference model (black dots) and with the modified gas absorption for tests #1 (red dots), #2 (blue dots), and #3 (pink dots). For clarity, we display the error bars for the reference model only. The error bars obtained with the other models are very similar.

In the text