# Distribution and intensity of water ice signature in South Xanadu and Tui
Regio

Maélie Coutelier11affiliation: Groupe de Spectrométrie Moléculaire et
Atmosphérique - UMR CNRS 7331 Campus Moulin de la Housse - BP 1039 Université
de Reims Champagne-Ardenne 51687 REIMS – France , Daniel Cordier11affiliation:
Groupe de Spectrométrie Moléculaire et Atmosphérique - UMR CNRS 7331 Campus
Moulin de la Housse - BP 1039 Université de Reims Champagne-Ardenne 51687
REIMS – France , Benoît Seignovert22affiliation: Université de Nantes, Univ
Angers, CNRS, LPG, UMR 6112, F-44000 Nantes, France , Pascal
Rannou11affiliation: Groupe de Spectrométrie Moléculaire et Atmosphérique -
UMR CNRS 7331 Campus Moulin de la Housse - BP 1039 Université de Reims
Champagne-Ardenne 51687 REIMS – France , Alice Le Gall33affiliation:
LATMOS/IPSL, UVSQ Université Paris-Saclay, Sorbonne Université, CNRS, Paris,
France ,44affiliation: Institut Universitaire de France (IUF), Paris, France ,
Thibaud Cours11affiliation: Groupe de Spectrométrie Moléculaire et
Atmosphérique - UMR CNRS 7331 Campus Moulin de la Housse - BP 1039 Université
de Reims Champagne-Ardenne 51687 REIMS – France , Luca
Maltagliati55affiliation: Nature Publishing Group, LONDON, UNITED-KINGDOM ,
Sébastien Rodriguez66affiliation: Université de Paris, Institut de physique du
globe de Paris, CNRS, F-75005 Paris, France<EMAIL_ADDRESS>
###### Abstract

Titan’s surface was revealed by Cassini’s instruments, showing the presence of
liquid hydrocarbons in lakes, and features like dry riverbed. In order to
study the sediment transport in Titan’s channels and to map distribution of
the water-ice signature in these terrains, we use a radiative transfer model
to retrieve the surface albedo, after we estimated VIMS error via an original
method. We also establish a criteria related to the intensity of the water ice
signature. The tuning of the radiative transfer model shows that the fractal
dimension of Titan’s aerosols is higher than previously thought, around 2.3 -
2.4. We find spots of increasing signal of water ice downstream, at the
margins of Tui Regio, that could correspond to alluvial fans, deltas or crater
rims. We also observe a very low water ice signal on Tui Regio, with a
positive gradient between the central region and the boundary of the area,
possibly due to the thickness variation of an evaporitic layer. The riverbeds
show within the error bars a decreasing grain size from the top to the bottom
of the channels.

††slugcomment:

##  1\. Introduction

Titan, the biggest moon of Saturn is one of the most complex moon of the Solar
System, with dense atmosphere, hydrocarbon cycle, and water-ice crust. Titan
has been targeted by two major space missions: Voyager I (1980) and Cassini-
Huygens (2004-2017), and is still the subject of many studies. Lakes and seas
of liquid hydrocarbons were discovered by Cassini’s ISS and RADAR in Titan’s
polar regions (Stofan et al., 2007; Turtle et al., 2009). The RADAR instrument
also revealed in equatorial regions geomorphological structures related to the
presence of liquid, like fluvial valleys incised in the _bedrock_ , and
alluvial fans (Le Gall et al., 2010; Birch et al., 2016). The existence of
evaporitic terrains where also suggested (Barnes et al., 2011; Cordier et al.,
2013; MacKenzie et al., 2014), often in place of paleo-sea (Moore & Howard,
2010). Water ice signature is not ubiquitous on Titan, contrary to most of
Saturn’s and Jupiter’s moons, like Encelade and Europa. It was detected with
Cassini’s Visual and Infrared Mapping Spectrometer (VIMS) in Titan’s infrared
dark region, often at the transition between dark and bright unit, and mixed
with a dark material (McCord et al., 2006; Soderblom et al., 2007; Brossier et
al., 2018). Griffith et al. (2019) also highlighted an equatorial corridor of
possibly exposed water-ice using a principal component analysis (PCA)
technique applied on VIMS data, showing a large scale terrains with low and
high water ice content. While PCA is extremely useful for large data sets, and
to study weak feature, it is not as detailed as a radiative transfer modeling
and can not be used to retrieve quantitative information.

Barnes et al. (2007) have studied the channels of Mabon Macula using a
combination of Synthethic Aperture Radar Image (SAR) and VIMS measurements.
They find that the spectral signals of the channels share characteristics with
Titan’s “dark blue” terrains according to the adopted color composition. The
considered terrains show an enrichment in water ice (Brossier et al., 2018;
Rodriguez et al., 2006) and the presence of benzene. These dark blue areas
should result from the transport of material eroded from the mountains by the
channels.This has been confirmed by Langhans et al. (2012) and Brossier et al.
(2018). The latter studied the “dark blue” terrains found with VIMS on Titan.
They used an accelerated version of radiative transfer model with
precalculated tables to compare the spectral variations of different terrains.
The IR-blue units were found to be enriched in water-ice (20 to 50 %), even if
it is not the main component. Brossier et al. (2018) also described different
types of IR-blue units. First, they have identified outwashed plains that are
located between IR-bright terrains and dune field, created by sediments
transported from mountain IR-bright units. Secondly, they found patches of IR-
blue units located inside IR-bright terrains, produced by the erosion of these
mountainous terrains. The last type is the ejecta blanket from impact craters.
They also found that rivers are IR-blue, probably enriched in water-ice
pebbles as well.

We focused our study on Titan’s river bed, where the erosion due to liquid
hydrocarbon could have revealed the bedrock, and where we could have clues on
the sediment transport. Here we propose a study of the water ice signal
distribution in South Xanadu and Tui Regio, where there are river channels (Le
Gall et al., 2010), paleo-seas of hydrocarbons (Hofgartner et al., 2020), and
terrains classified as IR-blue patches by Brossier et al. (2018). We employ a
photometric analysis combined with an original method based on indicators to
characterize the relative strength of the water-ice signature on the chosen
VIMS cubes.

We focused on two radar bright channels investigated by Le Gall et al. (2010),
larger than those reported by (Barnes et al., 2007) in Mabon Macula. These
channels are located south of Xanadu, emerging directly from a mountain range.
They are at the north of the evaporitic region Tui Regio (MacKenzie & Barnes,
2016). They form dry dendritic fluvial valleys (Langhans et al., 2012), with a
width up to $8$ km. Some areas studied in this work have also been
investigated by Griffith et al. (2019), where evaporitic terrains were
classified as ice-poor.

##  2\. Radar and Infrared Observations

To investigate our region of interest, we used two complementary instruments
of the Cassini’s mission: the RADAR, and VIMS. Cassini’s RADAR had four
different modes: it could operate as a SAR, scatterometer, altimeter and
radiometer at a wavelength of 2.2 cm. In SAR mode, it provided images of the
surface that are sensitive to local slope and terrains roughness: smooth
surfaces appear dark while rough terrains are radar-bright. About 70% of
Titan’s surface has been mapped with this mode with a resolution ranging
between 300 m/pixel and 4 km/pixel, revealing channels, lakes, mountains and
dune fields.

VIMS Infra-Red instrument onboard Cassini scanned in a spectral range from
0.88 $\mu$m to 5.12 $\mu$m. In Titan’s case, due to the methane absorption,
VIMS could only observe the surface in $7$ atmospheric windows centered at
$0.93$, $1.08$, $1.27$, $1.59$, $2.01$, $2.78$ and $5$ $\mu$m (Sotin et al.,
2005; Barnes et al., 2005). Parts of the spectrum not located in an
atmospheric window are called _bands_.As the RADAR spatial resolution was
higher than VIMS, the SAR map allows a better characterization of
geomorphologic structures in VIMS cubes.

Figure 1.—  Grey-scale: (a) RADAR swath containing the T48 and T44 flyby (Le
Gall et al., 2010). (b) Annotated RADAR map, with footprints of VIMS cube
1590648776_1 (C15, blue) and 1809727868_1 (C18, orange). Purple lines: the
rivers bright channels on which we focus our study. Tui Regio is roughly
delimited in yellow, Xanadu mountains are similarly delimited in red, and the
units defined by Lopes et al. (2019) and Moore & Howard (2010) are annotated
in cyan. Green dots: channels that end in fan (Radebaugh et al., 2018). Blues
lines: areas of high water-ice signal. The row of dark spots cutting the map
vertically in half comes from HiSAR data, where the surface resolution is
lower than SAR data. As an other useful central portion of this map comes from
HiSAR data, we have chosen to keep the complete map instead of working only
with SAR flyby.

The radar-bright river channels from Le Gall et al. (2010) have been imaged
during the T44 flyby (2008-28-05) using SAR mode. We employed the T44 SAR
image to select VIMS images overlapping these channels. Among them, we chose
the cubes with the highest resolution in this area. Only two cubes have a
pixel scale smaller than 18 km/pixel (see Tab. 1 and Fig. 1), which is twice
the width of the bright channels. A weaker resolution may not contain enough
channel signal to be distinguished. The first cube 1590648776_1 (labeled C15
hereafter) was acquired in 2008 during the same flyby as our reference radar
swath. For this cube, VIMS was configured in High Resolution mode which
increases the spatial resolution in the samples direction. The second cube
(C18 hereafter) was taken in normal mode 7 years later in 2015 during the T111
flyby. For C18 we do not work with the entire cube, only a square made of 38
rows of samples and lines located in our region of interest. These cubes were
calibrated with the RC19 calibration published by Clark et al.
(2018)111https://pds-imaging.jpl.nasa.gov/data/cassini/cassini_orbiter/vims-
calibration-files/, using ISIS3222https://isis.astrogeology.usgs.gov/ USGS
tools. They are available at the Cassini VIMS Data Portal333https://vims.univ-
nantes.fr (Le Mouélic et al., 2019). The navigation of the VIMS-IR cubes were
computed using the latest C-smithed kernels (December 2017) and the python
SpiceyPy toolkit (Annex et al., 2020). We took into account the offseted
pointing boresight and the corrected initial scan mirror sample position for
cubes acquired in HighRes mode, as recently reported by Nicholson et al.
(2020). Where changes exist, they are about 1 normal VIMS pixel compared to
the original ISIS navigation and within the expected accuracy of the Cassini
kernels ($\sim$0.12 mrad).

Table 1VIMS cubes details. C15 And C18 are used in the study of the bright channels, while C14 is employed to test the validity of the model at Huygens landing site. Cube number | Name in | Flyby | Date | Exposure | Mean ground | Size | Purpose
---|---|---|---|---|---|---|---
| the text |  | (DD-MM-YYYY) | (ms) | resolution | (Samples x Lines) |
1590648776_1 | C15 | T44 | 28-05-2008 | 325 | 15 km/pxl | (28 x 14) | Bright channels
1809727868_1 | C18 | T111 | 07-05-2015 | 122 | 14 km/pxl | (64 x 64) | Bright channels
1481624349_1 | C14 | TB | 13-12-2004 | 160 | 16 km/pxl | (24 x 12) | Model test
Figure 2.—  Equirectangular projection of (a) C15 and (b) C18 using Barnes et
al. (2007) colors: 5.0 $\mu$m as red, 2.00 $\mu$m as green, and 1.29 $\mu$m as
blue. This composition highlights evaporites in bright peach, sand in dark
brown, enhanced water ice signal in dark and light blue.

C15 and C18 are re-projected in a equirectangular view for better viewing and
comparison purpose in Fig. 2, with a color composition similar to that of
Barnes et al. (2007).

##  3\. The Radiative Transfer Model

In order to retrieve and study the surface albedo in different areas of
interest from VIMS-IR spectra, we used a model of photometry inherited from
Rannou et al. (2010, 2016). We first defined the average scattering and
absorption features of the atmosphere, derived from the properties of its
different components (particles, gases), as a function of wavelength and
altitude. These optical properties were cast under appropriate forms to be
used by radiative transfer (RT) solvers. This model depends on free parameters
including surface albedo.

###  3.1. Atmospheric properties and composition

Most of Titan’s atmospheric properties are known from Cassini and Huygens
observations. The atmosphere temperature and pressure vertical profiles have
been measured in situ by the Huygens Atmosphere Structure Instrument (HASI)
(Fulchignoni et al., 2005) during Huygens probe descent into Titan’s
atmosphere. The methane mixing ratio has been measured by the Gas
Chromatograph Mass Spectrometer (GCMS) onboard Huygens (Niemann et al., 2005,
2010). Isotopic ratios of CH3D and 13CH4 are those published by Nixon et al.
(2012) (see table 2). For CO, HCN and C2H2 mixing ratios and vertical profiles
(see table 2), we used results derived from CIRS observations (Vinatier et
al., 2015, also available in VESPA
database444http://vespa.obspm.fr/planetary/data/), extrapolated for the low
stratosphere and troposphere, and limited by the saturation vapor pressure.

###  3.2. Gas spectroscopic absorption

Molecular spectroscopic data (table 3) are mainly provided by
HITRAN555https://hitran.org database (Rothman et al., 2013), with a notable
exception for CH4 and its isotopes CH3D and 13CH4 for which we employed the
up-to-date theoretical line lists computed by Rey et al. (2017, 2018). These
data are publicly available at the Theoretical Reims-Tomsk Spectral
database666https://theorets.tsu.ru and https://theorets.univ-reims.fr (Rey et
al., 2016). Following de Bergh et al. (2012) for the line broadening, we
represented the line profiles by a Voigt function with a cut-off at 26 cm-1
from the line center. They were then modified with a sublorentzian decay with
a wavenumber characteristic length of 120 cm-1 for all species except CO. For
the latter, we used the specific cut-off recommended by Lellouch et al.
(2003). As VIMS channels wavelengths change through the mission, we followed
Clark et al. (2018) recommendations and calculated the spectroscopic
absorption consequently.

Table 2List of CH4 isotopes ratio taken into account, from Nixon et al. (2012), and other absorbing gases mixing ratio derived from CIRS observations (Vinatier et al., 2015). Methane | Isotope | Absorbing | Mixing
---|---|---|---
| ratio | gaz | ratio
D/H |  1,59$\times$10-4 | CO |  3.3$\times$10-4
12C/13C | 86,5 |  C2H2 |  4$\times$10-4
CH3D/CH4 |  6,36$\times$10-4 | HCN |  3$\times$10-7
13CH4/CH4 |  1,156$\times$10-2 |  |

###  3.3. Aerosol model

Considering the previous works of Tomasko et al. (2008) and Doose et al.
(2016), we divided the aerosols altitude distribution into two layers (see
Fig. 3). The upper haze layer extends from the tropopause at 55 km to the top
of the atmosphere. It is composed of photochemical aerosols assumed to be
fractal aggregates. The lower layer, hereafter called _mist_ , is an undefined
mixture of condensed material and photochemical haze.

To model the haze optical properties such as the single scattering albedos and
the extinction cross-sections, we adopted the scattering model for fractal
aggregate developed by Rannou et al. (1997). The aerosol phase functions that
we use are adapted for the backscattering. The preliminary analysis of a
sequence of VIMS observations of Titan acquired in November 2012 (flyby T88)
led to propose a modification of the aerosol phase function initially
recommended by Doose et al. (2016). During this 5-min long sequence, VIMS
acquired 26 consecutive datacubes over the same location, but with varying
emission and phase angles. In order to keep the haze opacity constant over the
26 datacubes, one had to change the aerosol phase functions by adding 20% more
of backscattering at all wavelengths. We used these modified aerosol phase
functions in the present study. A comprehensive analysis of these observations
is out of the scope of the present article, but will be soon subject to a
separate dedicated article.

This model allows us to study different morphology of aggregates, whereas
Tomasko et al. (2008) phase functions are only valid for a restricted set of
parameters (i.e. a fractal dimension of 2 and 3000 monomers of 50 nm). In our
model, we keep the fractal dimension as a free parameter. The aerosols
refractive indices were derived from those published by Rannou et al. (2010).
They were adapted to fit the haze single scattering albedo in agreement with
the prescription proposed by Doose et al. (2016). Since mist particles remain
poorly known, we made no differences between haze and mist particles phase
functions. The haze and mist extinction $k(z,\lambda)$ spectral slopes were
obtained with this scattering model, because Doose’s results based on the
Descent Imager / Spectral Radiometer (DISR) data from Hyugens probe are only
available between $\sim$ 0.4 and 1.6 $\mu$m.

Consistently, our mist scattering albedo $\omega_{\rm mist}$ follows Doose et
al. (2016) prescription, i.e. $\omega_{\rm mist}=(0.565+\omega_{\rm
haze})/1.5$. Similarly, the transition altitude between ”mist scattering
regime” and ”haze scattering regime” and the vertical extinction profile
$k(z,\lambda)$ came from the same source.

The haze and mist optical depths were scaled with respect to the optical
thickness obtained by Doose et al. (2016) at 1 $\mu$m. Those profiles were
adjusted thanks to a minimization procedure with two free parameters: $F_{h}$
for the haze, and $F_{m}$ for the mist, allowing for different ratios for the
haze and mist. As a consequence, the haze and mist opacity at an altitude $z$
and a wavelength $\lambda$ are:

$\Delta\tau_{h}(z,\lambda)=F_{h}\Delta\tau_{h_{D}}(z)\frac{\sigma(\lambda)}{\sigma(\lambda_{0})}$
(1)
$\Delta\tau_{m}(z,\lambda)=F_{m}\Delta\tau_{m_{D}}(z)\frac{\sigma(\lambda)}{\sigma(\lambda_{0})}$
(2)

with $\sigma(\lambda)$ the extinction cross-section of fractal aggregates,
$\Delta\tau_{h_{D}}(z)$ and $\Delta\tau_{m_{D}}(z)$ the normalized haze or
mist opacity of Doose at 1 $\mu$m for an atmospheric layer at the altitude $z$
times the total opacity of the atmospheric column at $\lambda_{0}=$ 1 $\mu$m.

###  3.4. Other components and technical aspects

The scattering coefficients due to Rayleigh scattering were derived from
equations (2.30) and (2.31) of Hansen & Travis (1974). They adapted to Titan’s
atmosphere using the CH4 column density. Our model also includes N2-N2 and
N2-H2 collision-induced absorption (CIA) (Lafferty et al., 1996; McKellar,
1989). Titan’s surface is assumed to scatter as a lambertian surface.

We discretized the atmosphere into $70$ layers uniformly distributed from
ground level up to an altitude of $700$ km, since it was found as an
acceptable compromise between accuracy and computation time (McKay et al.,
1989; Cours et al., 2020). In addition, for each layer we computed the average
optical depths, phase functions and single scattering albedos for all VIMS-IR
channels. We have the choice between two solvers to solve numerically the
radiative transfer equations: SHDOMPP developed by K.
Evans777http://coloradolinux.com/shdom/shdomppda/ (Evans, 2007) that runs in a
plan-parallel approximation, and SPSDISORT from the library
libRadtran888http://www.libradtran.org/doku.php (Mayer & Kylling, 2005; Emde
et al., 2016) that includes a pseudo-spherical geometry. The first solver was
found to compute ten times faster thanthe second, which has a wider validity
range since it includes spherical corrections. The optical depths of the
different gases were calculated according to the correlated-k approximation
(Goody et al., 1989), including four terms for each channel to represent the
opacity of the gases.

In radiative transfer models such as SHDOMPP or SPSDISORT, the phase functions
have to be expanded on a basis of Legendre polynomials, and we use 100 terms.
In the model, radiance is calculated for a given number of directions (that
is, streams). The resolution and computation time increase with the number of
streams. The section 5.1 treats in details the optimal number of streams to
use in the model.

The synthetic spectra were fitted to observational data in the bands using the
Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963). For that
purpose we employed the implementation of this method available under the form
of the LMDIF subroutine, part of the MINPACK library (Moré et al., 1980)
written in FORTRAN.

Table 3List of gaseous species accounted for in our radiative transfer model. Species | Instrument | Spectrum
---|---|---
CH4 | GCMS | TheoReTS
CH3D | CIRS/GCMS | TheoReTS
13CH4 | CIRS/GCMS | TheoReTS
CO | CIRS | HITRAN
C2H2 | CIRS | HITRAN
HCN | CIRS | HITRAN
Figure 3.— Sketch of the atmospheric profile we assumed in this work. Fractals
are from Botet et al. (1995), and this illustration if from Seignovert &
Rannou (2019). According to the adopted global picture, a low altitude layer,
extending from the ground to $55$ km, is filled by a misty air containing more
or less methane-wet particles. This first layer is topped by a less dense but
more extended region of haze settled by fractal particles. In our model, these
particles are assumed to be made of $3000$ monomers, each with an individual
radius of $50$ nm, the adopted fractal dimension is $\mathcal{D}_{f}=2.3$. The
single-scattering albedo $\omega_{\rm mist}$ of mist _particles-droplets_ has
been scaled to the single-scattering albedo $\omega_{\rm haze}$ according to
the law: $\omega_{\rm mist}=(0.565+\omega_{\rm haze})/1.5$ (Doose et al.,
2016).

##  4\. On the error of VIMS

The photometric uncertainties in VIMS calibrated cubes are not well known, and
depend on the considered channel (Brown et al., 2004). Sromovsky & Fry (2010)
established a global average error of 2.5% of the $I/F$ values, with an $I/F$
offset error of 5$\times$10-4. In order to provide a more accurate description
of the instrument noise, we developed a new method based on the analysis of
inter-pixels variability. With this method we can estimate combined errors due
to the instrument, data calibration and post treatment. Assuming a slow pixel
to pixel variation, we perform an estimate of the intrinsic variability of
adjacent pixels.

First, for a given data cube, we chose about 10 pixels evenly distributed over
the image. We took into account the 8 first neighbors around them, forming a
3$\times$3 box around the selected pixel. For each set of 9 pixels, we
calculated at each wavelength the mean $\overline{I/F}$ and standard deviation
$\sigma$. Finally we derived an estimator of the relative standard deviation
$\Sigma$ of each 9-pixels set:

$\Sigma=\text{std}(I/F)/\text{mean}(I/F)$ (3)

For C15 and C18, we calculated $\Sigma$ for the whole spectrum as we need it
to estimate the surface albedo error bars. The results are represented by dots
in Fig. 4-a. The method is only strictly valid in the bands because in the
windows we get a variability of the signal related to the natural variability
of the surface. We calculated $\Sigma$ in the atmospheric windows for
comparison purpose. The geometry changes between each pixels and can influence
the retrieved $\Sigma$. This will be corrected by another step of the method.

By analyzing different set of pixels, over a few images, we can first notice
in Fig. 4-a that $\Sigma$ depends on the cube and the wavelength. We
determined specific profiles of $\Sigma$ for each image of interest. They were
chosen manually, following the general shape of the results ($\Sigma_{1}$ and
$\Sigma_{2}$ uncorrected in Fig. 4-a for C15 et C18 respectively). We chose to
cap $\Sigma$ below 2.5 $\mu$m and at 5 $\mu$m by a high value so we do not
underestimate the error bars. To obtain the uncorrected $\Sigma$ in
atmospheric windows, we interpolate the uncorrected $\Sigma$ of adjacent
bands.

We then determined the error due to variations in viewing geometry from one
pixel to an other. We obtained it with the same method, but instead of the
cube spectra, we used synthetic spectra with fixed parameters for the surface
albedo and the aerosols. For these synthetic spectra, only incidence and
emergence angles change, similarly as those of the cube. In that way, the only
differences in $I/F$ between pixels of the same cube are due to angle
variations. With the differences in incidence and emergence angles between C15
and C18, the thickness of the atmosphere layer crossed by the beams is greater
for C18 than for C15. As a result, contrary to C18 where incident angles are
higher, the error due to angle variations is lower in the bands for C15.

We subtracted from our previous profiles the $\Sigma$ due to angle variation.
Final profiles of $\Sigma$ are shown in Fig 4-b. Here $\Sigma_{1}$ is adapted
for C15, whereas the profile $\Sigma_{2}$ is adapted for C18. Theoretically,
in the bands we should have a negligible influence from the surface. In Fig.
4-a for C15 and C18 (red and blue points) we tend to have a lower $\Sigma$ in
the windows. VIMS uncertainties depend on the number of photons. The greater
$I/F$, the smaller uncertainties. Because of the interpolation of
$\Sigma_{1,2}$ uncorrected for the windows, we may overestimate the error in
the atmospheric windows where we have a greater $I/F$ than in the bands.

The corrected $\Sigma$ for C15 and C18 in Fig. 4-b show the influence of
incident and emergent angles. The corrected $\Sigma$ for C15 is higher in the
windows than in the bands, under 2 $\mu$m. For C18 it is slightly lower,
except at 5 $\mu$m. Generally, there is an enhanced variability and error
around 2.35 $\mu$m, and beyond 3 $\mu$m. Our results there are quite different
from the 2.5% estimated by Sromovsky & Fry (2010), but below 2$\mu$m the
values are of the same order in the windows. The influence of the temperature
inside VIMS, causing IR radiation, is greater at large wavelengths. That could
partly explain the increasing uncertainties with the wavelength. The low
number of photons at greater wavelength is also another factor.

Figure 4.— (a) Chosen spectral profiles of $\Sigma$ (lines), and retrieved
values of $\Sigma$ (points) for C15 (red) and C18 (blue). Titan’s atmospheric
windows are in gray rectangles. (b) Final spectral profiles of $\Sigma$
(lines) for the C15 ($\Sigma_{1}$) and C18 ($\Sigma_{2}$), where we subtracted
from our previous uncorrected profiles the $\Sigma$ due to the angle variation
(dots).

These $\Sigma$ are significant in the bands up to 2.5 $\mu$m for the
calculation of the $\chi^{2}$ and the minimization of the function
$\mathcal{Y}$ used to retrieve atmospheric parameters Fh and Fm, the scaling
factors for the haze and the mist:

$\mathcal{Y}_{n}=\dfrac{|(I/F)_{n}^{obs}-(I/F)_{n}^{mod}|}{\Sigma_{n}\times(I/F)_{n}^{obs}}$
(4) $\chi^{2}=\frac{1}{N}\sum_{n=1}^{N}\mathcal{Y}_{n}^{2}$ (5)

with $N$ the number of channels we used to determine the parameters,
$I/F_{obs}$ the observed $I/F$ and $I/F_{mod}$ the modeled $I/F$.
$\Sigma_{\lambda}$ is also used later to estimate the albedo error bars.

Once surface albedos are determined for each channel, we can finally evaluate
the uncertainties of the surface albedo for a wavelength $\lambda$. We need
the values of $\Sigma_{\lambda}$ previously calculated and the derivative of
the albedo $A_{\lambda}$ with respect to $I/F_{\lambda}$:

$\sigma_{A_{\lambda}}=\left|\dfrac{dA_{\lambda}}{dI/F_{\lambda}}\right|\Sigma_{\lambda}$
(6)

These errors will be used to calculate our results uncertainties.

##  5\. Test and preparation of the model for surface analysis

In this section we present sensitivity tests performed with the radiative
transfer model, using SPSDISORT and SHDOMPP. The scope of this section is to
define the number of stream used to discretize the intensity field in the
model and to define the limit incident and emergent angles for using the
parallel plane model beyond which we must account for the atmosphere
sphericity. Then, we also describe how we set the atmosphere free parameters
and our observed data fitting method.

###  5.1. Tests of the model against the number of streams

A solver for the radiative transfer equations calculates the intensity field
in a finite number of directions for each atmospheric layer. This is what we
call the number of streams. In order to optimize the model execution time, we
tested the solvers accuracy with different number of streams. A high number of
streams increases the accuracy, but also the execution time.

We define a testing case with fixed haze and mist parameters $F_{h}=F_{m}=1$.
The surface albedo is set at a uniform value of 0.05. We use a near-nadir
viewing with an incidence angle of 3∘, an emergence angle of 5∘ and a phase
angle of 7∘. Those are extreme conditions for the SHDOMPP because the nadir
configuration requires more accuracy, i.e. a higher number of streams than
other geometries (Evans, 2007).

We set a reference with a simulation made with 100 streams, as we would never
use a number so high due to the runtime, even if the accuracy of the results
is better. We then compare the model performance with different number of
streams. In order to estimate the minimum number of streams required to have
less than 0.5% difference in $I/F$ at each wavelength with the reference case,
we start with 16 and we increase numbers by steps of 4 streams (Fig. 5). As
the differences between the reference case and the tests decrease, the
differences between the tests also decrease, allowing us to choose an optimum
number of stream.

For the solver SHDOMPP, we have less than 1% difference when using more than
28 streams. With 32 streams the difference falls below the 0.5% limit. The
solver SPSDISORT requires a minimum of 36 streams to have less than 0.5%
difference compared to the 100-streams test case. In view of the results in
Fig. 5, and taking into account the geometry and the increase of computation
time with the number of streams, we decided to use for SHDOMPP a number of 32
streams up to 1.7 $\mu$m, and 24 streams for the rest of VIMS-IR wavelengths.
For the SPSDISORT, we work with 40 streams up to 1.7 $\mu$m, and 24 streams
for the reminder of the spectrum. Doing so, we have less than 0.5 % of
difference with the reference test at each wavelength, with a smaller
computation time. For example, for a cube made of 4096 pixels, the RT model
using SPSDISORT has a runtime of less than 24 hours when using 500 processors.

Figure 5.— Differences in % between the $I/F$ calculated with SHDOMPP (top
panel) and SPSDISORT (bottom panel) with different number of streams compared
to the results with 100 streams. Incidence: 3∘, emergence: 5∘, phase: 7∘. The
relative error is: $\Delta r=|\frac{I/F_{Nstr}-I/F_{100}}{I/F_{100}}|$ .

###  5.2. Tests of the model in a spherical viewing

Barnes et al. (2018) consider as acceptable the use of parallel-plane
approximation up to 60∘ for incidence or emergence angles. In this work, we
compared the results obtained respectively with SHDOMPP in a parallel-plane
approximation and SPSDISORT in a pseudo-spherical geometry at different
incident angles (see Fig. 6). For this testing case, we use 32 streams for
both models, and a emergence angle of 10∘. We change the incidence angle from
10∘ to 80∘ with a 10∘ increment.

Fig. 6 shows that for incidence angles lower than 30∘, we have differences
between solvers smaller than 1%, which is in agreement with Barnes et al.
(2018) result with a Monte-Carlo RT model. Errors reach more than 5% in
Titan’s atmospheric windows at 60∘. Giving these results, we decided for our
study to employ SHDOMPP with C15 (Fig. 13-a) when the incident angles are
below 50∘. For C18 (Fig. 13-b), we chose to work with SPSDISORT, as in our
region of interest incidence angles go beyond 55∘.

Figure 6.— Differences (%) between $I/F$ calculated with SHDOMPP and SPSDISORT
for different incidence angles. Emergent angle: 10∘. Relative error: $\Delta
r=|\frac{I/F_{SHD}-I/F_{SPS}}{I/F_{SPS}}|$

###  5.3. Haze and mist properties for surface analysis

The last step of the radiative transfer model setting consisted in defining
the haze and mist properties that best reproduce the spectra in the bands,
where surface has no influence. With our model, we retrieved haze and mist
properties by fitting a spectrum taken in a dark tropical area of Titan in
order to minimize the influence of surface. More specifically, we chose 16
test pixels where surface albedo is known to be extremely dark (Griffith et
al., 2012), from the cube v1567239055, acquired on August, 31, 2007. In that
way, we have set the model in such a way a negative value of the surface
albedo is not possible even where surface albedo is especially low. VIMS
uncertainties defined in section 4 were taken into account in the algorithm.

We used the standard case for aerosols morphology described in section 3.3. In
Fig. 7 with a fractal dimension of $\mathcal{D}_{f}$=2 we obtain a good fit
except for the band around 1 $\mu$m. This is a recurring problem for each VIMS
cube. To solve this problem we first performed several tests, using other
tholins refractive index, varying the single scattering albedo with a constant
factor, or using the condition Fh=Fm. However this did not improve the fit
around 1 $\mu$m. We finally tested other values of fractal dimension.

Varying $\mathcal{D}_{f}$ has an impact on the spectral slopes of the haze and
mist optical depths (see Fig. 8), and on the phase functions. We kept the
number of monomers at 3000 since it was defined from the solar aurea with
DISR, which depends in first place on the amount of grains, regardless the
fractal dimension.

Then, in order to retrieve the surface albedo and the $\chi^{2}$ (see Eq. 5)
corresponding to each $\mathcal{D}_{f}$, the parameters Fh and Fm are set free
to fit the data in the bands. We found that the fit is greatly improved with
$\mathcal{D}_{f}>$ 2.0 in each of the selected pixels at short wavelengths
(see Fig. 7 for one of the test cases), especially for the band around 1
$\mu$m, and between 1.8 and 2 $\mu$m. The $\chi^{2}$ is generally reduced by a
factor of 2 between $\mathcal{D}_{f}$=2 ($\chi^{2}=250$ in the example) and
$\mathcal{D}_{f}$=2.3 ($\chi^{2}=120$ in the example). We have obtained a
better fit at 1 $\mu$m, at the short wavelength of the 1.6 $\mu$m window, and
around 1.8-2.0 $\mu$m with $\mathcal{D}_{f}=2.3-2.4$, than with
$\mathcal{D}_{f}=2$ as proposed by Tomasko et al. (2008). The surface albedo
is also impacted with the change of the fractal dimension. The slope of the
surface albedo at 1.28 $\mu$m for $\mathcal{D}_{f}>2.1$ is reversed and the
shape at 1.55 $\mu$m is flatter even if the graphic scale tends to smooth the
variations. Lastly, the surface albedo globally decreases when increasing
fractal dimensions, except at 0.93$\mu$m (Fig. 7-bottom).

To compare the total haze opacity with Doose et al. (2016) analytical model,
we have chosen a test case with $F_{h}=F_{m}=1$. The optical depths obtained
in Fig. 8 for different fractal dimensions are adjusted for a better viewing.
The differences of $I/F$ at 1 $\mu$m in Fig. 7 are correlated with the
variations of haze opacity $\tau$ as a function of wavelength in Fig. 8: the
spectral slopes are a close equivalent to the one recommended by Doose et al.
(2016) beyond $\sim$ 1.3 $\mu$m. Below this threshold, the smaller is
$\mathcal{D}_{f}$, the steeper is the spectral slope. While opacity remains
the same above 1.3 $\mu$m, fractal aggregates lead to a too large opacity at 1
$\mu$m for $\mathcal{D}_{f}$ = 2.0, and yield an opacity comparable to the
expected opacity at $\mathcal{D}_{f}$ $\sim$ 2.3 - 2.4. As Fig.8 shows,
$\mathcal{D}_{f}$ = 2.4 is a good candidate, but below 0.8 $\mu$m the result
is less consistent with Doose et al. (2016) than for $\mathcal{D}_{f}$ = 2.3.

These results imply that aerosols are more compact than expected at least, at
this location and date. They also imply that the fractal dimension value is
between 2.3 and 2.4. This is the case for C15 and C18 as well, as the fit is
also improved for the pixels we have chosen to do the comparison with.

Figure 7.— Fits of the pixel (38,19) of the cube v1567239055, using different
fractal dimension (top), and the corresponding surface albedo (bottom). Dotted
lines: uncertainties for $\mathcal{D}_{f}$ = 2 and $\mathcal{D}_{f}$ = 2.4.

Based on these results, we decided for our study to employ a fractal dimension
of 2.3, even if we obtain almost similar results with $\mathcal{D}_{f}$ = 2.4.
As seen in Fig. 7, with $\mathcal{D}_{f}$ = 2.4 the fit is the best at 1
$\mu$m, but $\chi^{2}=150$ which is larger than for $\mathcal{D}_{f}$ = 2.3
($\chi^{2}=120$), meaning that the fit is not as great at other wavelengths.

Figure 8.— Total haze opacity with Fh and Fm equal to 1, and angles of
incidence and emergence at 10∘, with different values of the fractal
dimension. The curves follows the microphysics aerosol model from Rannou et
al. (2004). Results are adjusted to fit Doose et al. (2016) curve at 2.88
$\mu$m

###  5.4. Tests of the final model on Huygens landing site

Thanks to Huygens probe, the aerosol vertical profile is known at Huygens
landing site i.e., at 192.3°W 10.5°S Kazeminejad et al. (2011) on January, 14,
2005. To confirm the model validity, we ran a test on a cube containing
Huygens landing site, and at a date close to Huygens landing day. For cubes
respecting these criteria, the retrieved values of Fh and Fm should be close
to 1 as we use the vertical extinction profile retrieved from Huygens given by
Doose et al. (2016). The tests where made with our final model settings,
$\mathcal{D}_{f}=2.3$ and the stream number chosen in section 5.1. Incidence
and emergence angles are in the validity domain of SHDOMPP.

The closest cube in time from Huygens landing date is the cube 1481624349_1
(labeled C14 hereafter), acquired with the High Resolution mode on December
13th, 2004, one month before the landing. We performed the RT analysis on the
entire cube. Results for the haze and mist parameters are presented in Fig. 9.
The histogram scales are the same, for better comparison of the results. They
show a distribution of the parameters $F_{m}$ and Fh close to 1.1 for the case
$F_{h}=F_{m}$ (Fig. 9-a) where the normalized vertical profile is equally
multiplied. This result is good considering that we adjust the parameter using
a larger wavelength range with VIMS-IR than for DISR, and that we do not have
the same incident and emergent angles. The pixel (10,2) overlaps Huygens
landing site, where the $\chi^{2}$ from Eq. 5 is 90 for $F_{m}=F_{h}=1.09$,
and 75 for $F_{h}\neq F_{m}$ with $F_{h}=1.08$ and $F_{m}=1.24$. Having
different haze and mist parameters thus leads to a better fit.

Figure 9.—  (a) Histogram of retrieved haze and mist parameters for
$F_{h}=F_{m}$ for the cube C14. (b,c) Histogram of retrieved haze and mist
parameters respectively for the cube C14 for the case $F_{h}\neq F_{m}$.
Results at Huygens landing site are indicated in red. The relative standard
deviation is noted $\Delta F_{h}$* in order not to be confused with $\Sigma$.
(b)

There is no reason for having the same haze and mist parameters, so for the
second test case, their values are allowed to be different. The resulting
best-fit mean values are 1.22 for Fm (Fig. 9-b) and 1.09 for $F_{h}$ (Fig.
9-c). The mist parameter is more dependent on the lower atmosphere, and is
more variable than the haze parameter. This is why the result differs from the
first test case. $F_{h}$ is close to the previous test value while $F_{m}$
differs, showing that the mist only plays a second order role in the model.
Using two parameters instead of one thus results in a gain of information on
the mist. The 10 % differences between Doose et al. (2016) values and ours are
acceptable given the different initial data.

Since DISR instrument on board Huygens made the unique in situ measurement of
a Titan’s terrain albedo, we propose here a comparison between these
observations and surface reflection coefficient obtained using our method on
C14 cube. DISR wavelengths range is less extended than VIMS-IR, so the
comparison is only possible below 1.6 $\mu$m (Schröder & Keller, 2008;
Karkoschka et al., 2012).

In Fig. 10, we plotted our results derived form VIMS cubes against two DISR
observations. We can not compare the absolute albedos directly because of the
different conditions like the pixel resolution (14 km for our cube versus a
few meters at most for DISR after landing), and the geometry of incidence and
emergence angles. In our model we assumed a lambertian surface with isotropic
reflectivity, but Huygens landing site reflectivity is more complex and angle-
dependent (Karkoschka et al., 2012; Karkoschka, 2016). We do not expect to
obtain exactly the same surface albedo as the previous studies for these
reasons. However, the global decreasing shape is similar, except at 1.08
$\mu$m. We also note at 1.55 $\mu$m a similar bowl shape of the albedo curve.
The slight difference at 1.08 $\mu$m could mean that this pixel also contains
absorbing material around 1.08 $\mu$m. It may also be due to a missing
absorption in this window in our model, or a difference in the atmosphere
properties, for a small change in the single scattering albedo of the mist has
a greater effect at short wavelength ($<1.3\mu$m) (see the light blue results
in Fig. 10). With an observed field much larger for VIMS, the chemical
composition and spectra of considered terrains may be heterogenous for fields
as large as $\sim 100$ km2; explaining the differences that we have noticed.
It could also be caused by the slight VIMS wavelength shift that is more
noticeable in the firsts windows since they are much narrower. The differences
occurs for other cubes taken in this area and close to the landing date, but
not with cubes taken years later, maybe because of the shift.

Figure 10.—  Huygens landing site surface albedo, derived from DISR and VIMS
instruments. The surface albedo of pixel (10,2) in the cube C14 is retrieved
for the cases $F_{h}=F_{m}$ (purple), $F_{h}\neq F_{m}$ (green) and for a case
with the single scattering albedo (SSD) of the mist decreased by -0.03 (light
blue). Results from Schröder & Keller (2008) (red) and Karkoschka et al.
(2012) (dark blue), based on DISR dataset are indicated. Literature data are
for a null phase angle. We approximately scaled them following Karkoschka et
al. (2012) for a phase angle of 18∘.

Globally, this VIMS/DISR comparison is reassuring even if, for above mentioned
reasons, it cannot be considered as a formal validation of our approach.

##  6\. Surface albedos of Xanadu and Tui regio

###  6.1. Method

To retrieve the surface albedo from VIMS pixel data, we fit the observed
spectra in the methane bands by adjusting both scaling factor Fh and Fm, using
the minimization method described in Sect. 4. Once the scaling factors are
adjusted, we retrieve the surface albedo in the methane windows. Fig. 4-a
shows the detailed windows. We routinely ran three simulations for each
window, assuming a surface albedo of 0.05, 0.15 and 0.25, independently of the
wavelength. For each VIMS channel sensitive to the surface, the surface albedo
is finally determined by linear interpolation.

###  6.2. Criteria for water ice presence

Using VIMS spectra collected on Enceladus and shown in Fig.11-a, we determined
the water ice absorptions. They are centered at 1.04, 1.27, 1.51, 2.03 and 3
$\mu$m. Titan’s atmospheric windows are not exactly located where these
absorptions occur, but they are influenced by them even when water ice
absorptions are located at the edges of atmospheric windows.

Figure 11.—  (a) Spectra of Enceladus, from 5 different VIMS cubes, (b)
Experimental water-ice spectra with different grain size, from the GHOSST
database. Grey columns represent the channels used to calculate the ratios.
Their widths depends on the number of channels averaged.

Due to Titan’s atmosphere strong absorptions, and because sedimented aerosols
partly cover the surface, we do not expect a pure water-ice spectra on Titan.
This is why we looked for traces and influence of the water-ice in Titan’s
surface spectra.We have then defined our criteria to measure the water ice
signature on Titan. These criteria are based on a combination of conditions on
the ratios between the water-ice absorption bands (Barnes et al., 2007;
Robidel et al., 2020). Following the shape of the water ice spectrum, those
mean and ratios values likely indicate the presence of water ice signature if:

  * •

$r_{1}$ : $1.58/1.28$ $<1$, as water ice absorbs more at 1.58 $\mu$m

  * •

$r_{2}$ : $2.03/1.28$ $<1$, as water ice absorbs more at 2.03 $\mu$m

  * •

$r_{3}$ : $2.03/1.58$ $<1$, as water ice absorbs more at 2.03 $\mu$m

  * •

$r_{4}$ : $2.7/2.8$ $>1$, as water ice absorbs more at 2.8 $\mu$m

  * •

$\overline{A_{5}}$ : Mean of the $5$ $\mu$m window low as water ice spectrum
is low at $5$ $\mu$m compared to the rest of the spectrum.

We chose not to use the ratio $1.28/1.08$ (Barnes et al., 2007), which,
according to the considered spectra of Enceladus (Fig. 11-a), can either be
greater or lower than 1 - possibly because of impurities - even if in the pure
water ice case from GHOSST999https://ghosst.osug.fr/ database (Fig. 11-b), the
value of the ratio is strictly greater than 1. Not all ratios have the same
impact: we tend to have greater error bars at short wavelengths for the
albedo, and there are known missing absorptions (Cours et al., 2020), like for
ethane at 2.7-2.8 $\mu$m (Maltagliati et al., 2015). Furthermore, we work with
the VIMS data calibration RC-19 which is more recent, but less reliable for
the 2.7-2.8 windows compared to the RC-17 as it tends to change the
intensities values in this window by about 5 % (Clark et al., 2018). With the
assumption that the effect of the missing ethane absorptions and the
calibration is the same for each pixel, working in relative and not in
absolute is acceptable. Even so, we chose to weight the ratios and mean by
their known calculated uncertainties. For the mean albedo $\overline{A_{i}}$
of a window of $N$ values, the uncertainties $\Delta\overline{A_{i}}$ is:

$\Delta\overline{A_{i}}=\frac{1}{N}\sqrt{\sigma_{A_{i_{1}}}^{2}+...+\sigma_{A_{i_{N}}}^{2}}$
(7)

with $\sigma_{A_{i}}$ the uncertainty of the surface albedo at a given
wavelength (see Eq. 6). Then for each ratio
$r=\overline{A}_{i}/\overline{A}_{j}$, the relative uncertainty of the ratio
is:

$\sigma_{r}=\dfrac{\Delta
r}{r}=\dfrac{\Delta\overline{A_{i}}}{\overline{A_{i}}}+\dfrac{\Delta\overline{A_{j}}}{\overline{A_{j}}}$
(8)

We defined a multi-criteria analysis parameter $\delta$ to estimate the
presence of water ice. This parameter combines the relevant individual $I/F$
ratios and the 5$\mu$m mean previously mentioned. For the multi-criteria
analysis, these ratios and 5$\mu$m-mean are normalized between 0 and 1 with
the following method. For a ratio $r_{1}$, the normalized ratio $r_{1_{n}}$
is:

$r_{1_{n}}=\frac{r_{1}-r_{1_{min}}}{r_{1_{max}}-r_{1_{min}}}$ (9)

With $r_{1_{min}}$ and $r_{1_{max}}$ the minimum and maximum values of the
ratio for the considered cube. The normalized values are then added if the
expected value of the non-normalized ratio or mean in case of pure water ice
signature is greater than 1 (i.e. $r_{2.7/2.8}$). They are subtracted if the
expected value is lower than 1 (i.e. all other ratios and 5$\mu$m-mean), as
shown in equation 10. In this way, theoretically the highest possible value
for $\delta$ is 1 indicating a prominent water ice-rich signature, and the
lowest value is $-4$, for an extremely poor water-ice signal. In addition, the
normalized ratios and mean at 5 $\mu$m are weighted by (1- $\sigma_{r}$), with
$\sigma_{r}$ the relative error of each VIMS channels ratio $r$. We also set a
lower limit for (1- $\sigma_{r}$) which cannot be lower than 0, even if it
does not happen in our results. The parameter $\delta$ is explicitly:

$\delta=-(1-\sigma_{A_{5}})\overline{A_{5_{n}}}+\sum_{k=1}^{4}\pm(1-\sigma_{r_{k}})r_{n_{k}}$
(10)

With $\overline{A_{5_{n}}}$ the normalized mean surface albedo at 5 $\mu$m,
and $k$ the number of the normalized ratio $r_{n}$. $\delta$ increases when
the probability of water ice presence increases. With Eq. 7 and 8, we are able
to establish our ratios’ error, and determine which results are relevant. With
ratios weighted by their uncertainties, -4 and 1 are ideal values that can not
be obtained. Furthermore, $\delta$ must not be interpreted alone since it is
obtained using normalized values: there is no specific threshold value of
$\delta$ which define is indeed water ice or not. It shows subtle changes in
the surface composition. The possible results consistency between the cubes
reinforces the interpretation. We also need VIMS colormaps and the RADAR map
to scout the terrains beforehand and to better interpret the results. The
criteria $\delta$ should be considered as an indicator of the strength of the
water ice signature and not as a direct estimator of water ice abundance. It
should be regarded as a tool for terrain classification.

###  6.3. Haze and mist results

Figure 12.—  Maps of the haze and mist coefficients for C15 (a) and (c)
respectively, and for C18 (b) and (d) respectively, overlaid to the RADAR map.

Fig. 12 shows Fh and Fm regional variation in C15 and C18, overlaid on the
RADAR map. Each point corresponds to a cube pixel. The haze factor seems to be
more stable than the mist factor, which is expected since the mist factor
should be more sensitive to local atmospheric variations on the surface. Both
factors are not correlated, but we note a correlation with the incident and
emergent angles in Fig 13. C15 coefficient Fh (Fig. 12-a) is uncorrelated, but
Fm in C18 (Fig. 12-d) seems to be correlated with the emergent angles, while
Fm in C15 and Fh in C18 (Fig. 12 c and b) appear to be anti-correlated with
the incident angles. Yet, as no regular pattern stands out, there is no
argument for an actual correlation between the angles and the haze and mist
parameters values.

Figure 13.— (a) Pixel footprints of the VIMS cube 1590648776_1. Colors are
mapped with VIMS spectral channels ratios: $1.58/1.28$ as red, $2.03/1.28$ as
blue, $1.28/1.08$ as green, with contourlines corresponding to the incident
angles (orange lines) and emergence angles (blue lines). (b) Similar figure
for the VIMS cube 1809727868_1. This color code highlights water ice in dark
blue.

We acknowledge a bias in our model: the atmosphere is assumed to be laterally
homogeneous and infinite with no change of the aerosols vertical profile,
whereas in our results the parameters Fh and Fm change for each pixels.
However the bias is small, as at most rays go through 5 different columns of
atmosphere before reaching the ground pixel for C15, and 8 for C18, and the
values of the haze and mist parameters do not change in a significant manner
within 5 or 8 adjacent pixels. The model can be tuned with free parameters but
the haze and mist vertical profiles, the particle scattering phase functions
and the gas vertical profiles that we used may be different from the actual
profiles on Titan. It would modify the angular description of the modeled
intensities and of the retrieved surface albedo. Therefore, it could induce
another bias because the two images are not taken with the same observation
geometry.

###  6.4. Water-ice distribution results

The maps of different non-normalized ratio values for both VIMS cubes are in
Fig. 14, and the $\delta$ values are presented in Fig. 15. Given our estimated
uncertainties, the extreme values of the ratios (in red and blue) are
significant. Furthermore, similar structures stand out on both cubes, although
we cannot compare the absolute values as the illumination and viewing
geometries are different between C15 and C18 that, in addition, were acquired
7 years apart. It is also important to remind that the ice crust on Titan is
covered by sedimented aerosols (Janssen et al., 2016; Brossier et al., 2018).
The sediment layer thickness has an influence on VIMS spectral channels ratios
(Rannou et al., 2016). In both Fig. 14 and 15, units exhibiting different
strength of the water ice signal according to the multi-criteria analysis of
C18 (Fig. 15-d) are outlined. In C15, we do observe the highest water ice
signature on Xanadu mountains, and a medium-high signature at 24∘S 120∘W, next
to the dune unit. Surprisingly, in C18, the strongest signals are not in the
same place as in C15. They are at the ends of the channels in blue dotted
areas, which only displays a medium water ice signature in C15. There is a
clear correlation between the structures on the RADAR map, and the signature
of water ice. For example the blue dotted areas tend to correspond to dark and
smooth gray structures in the lowest resolution part of the RADAR map. The
area with low water-ice signature around 23.5∘S 116∘W corresponds to the dune
unit in the RADAR map, and its surrounding dark gray halo. The halo has a
lower water-ice signature than its center. A dark gray terrain at 22∘S 122∘W
on the RADAR map is also well delimited by a low $\delta$ on C15, and a
medium-low one in C18. This is also the case for the medium-low $\delta$ on
the top corners of C18. The largest terrain with low $\delta$ is located in
Tui Regio. Ratios are of medium value in the channel beds, suggesting a
slightly higher water ice signature than in the surrounding terrains (see Fig.
14-d-e for the most flagrant). Even if the error on these regions reduces the
significance of these variations, we do observe the same trend in both cubes.
This support the hypothesis that the channels descend from the mountains, and
contain water ice cobbles, probably covered by organic sediments. Above all,
we notice the presence of an area of high $\delta$ at the end of each channel.

The color code in Fig. 13 is set to highlight in dark blue terrains with high
water ice signature. However, our results (Fig. 14 and 15) reveal areas with
strong water ice signal in dark brown for C15 (Fig. 13-a). For C18 the highest
water ice signatures are in dark brown, and dark blue areas. This confirms the
need of a radiative transfer analysis to extract the surface signal from the
atmospheric signal, and to characterize the nuances of the water-ice signal
intensity.

Figure 14.—  Values of the ratios 2.7/2.8 (top panels), 2.0/1.28 (middle
panels), and 1.59/1.28 (bottom panels) for the VIMS cubes C15 (a,c,e), and C18
(b,d,f) respectively, superimposed on the RADAR map. Each point corresponds to
the position of a pixel of the cubes. Channels are outlined by purple lines,
and areas of low and high water-ice signal are circled in yellow and blue,
respectively. The absolute error for each ratio is: (a) $\pm 0.05$, (b) $\pm
0.15$, (c) $\pm 0.05$, (d) $\pm 0.05$, (e) $\pm 0.07$, (f)$\pm 0.06$. Figure
15.—  Values of the multi-criteria analysis parameter $\delta$ for the VIMS
cubes C15 (a), and C18 (b) respectively, superimposed to the RADAR map. Each
point corresponds to the position of a pixel from the cubes. Channels are
outlined by purple lines. Regions of high and low $\delta$ are annotated in
blue and yellow lines respectively, and fans highlighted by Radebaugh et al.
(2018) are represented by green dots. Figure 16.—  $I/F$ spectra of pixels
(34,16) and (15,25) of the cube C18 with low (red) and high (blue) values of
our water-ice signal criteria $\delta$ (top panel), and their respective
albedo (bottom panel).

Fig. 16 displays the spectra and surface albedo of two pixels taken in zones
of low and high water-ice signature for the cube C18 (respectively,
$\delta=0.5$, and $\delta=-3.4$). The presence of water-ice darkens the albedo
in most of Titan’s atmospheric windows (see Fig. 11 water-ice absorptions),
consequently lowering the $I/F$, whereas where the water-ice signal is low,
the albedo is increased. This effect is especially significant at 1.5, 2, 2.8
and 5 $\mu$m. However the shape of the albedo for the 2.8 $\mu$m window
changes the most.

##  7\. Discussion on the water ice signature and on geological structures

With the color maps, we can discuss the meaning of the ice index and the
connection with surface features and processes. At first glance, the two
dataset C15 and C18 yield, for each color ratios and for the multicriteria
($\delta$), similar and consistent maps. However, we note that the absolute
values reported on the maps systematically differ between C15 and C18. We also
remark that the nature of some areas clearly differs, which reveals changes
between the two flybys. We describe in this section the factors that help
understanding our maps and the parameters that could influence and modulate
the retrieved surface albedos, ratios and multicriteria index in the maps.

###  7.1. Grain size and sediment cover

First of all, the intensity of water ice signatures depends on the size of ice
grains. Data displayed in Fig. 11-b clearly illustrates this dependence. This
stresses that a change in the strength of water ice signature could be due to
a change in the average size of ice grains composing the uppermost layer of
Titan’s ground, or to a thicker sedimentary layer covering or mixed with the
ice layer (Brossier et al., 2018). Based on this figure, the ratio $2.7/2.8$
could be (with caution) used as an indicator of the grain size due to the
strong absorption of water ice at 2.9 to 3.25 $\mu$m. The higher the ratio is,
the smaller the grain size is. It could be directly used to classify different
terrains if the sedimentary layer is homogeneous. If the sedimentary layer is
not homogeneous, disentangling the two effects is not possible because the
effect of a coverage by sediments is not known. A specific model is needed to
interpret the details of the color ratio in term of physical parameters, like
what was done in Brossier et al. (2018), where they noted in the IR-blue
terrains a grain size generally larger than 100 $\mu$m, and half as many
solutions with a grain size smaller than 0.1 $\mu$m. Another difficulty arises
because the 2.7-2.8 $\mu$m window is incompletely treated with ethane and
possibly other molecules absorptions missing in the model.

Therefore, in this work, we can compare and classify the different terrains
with the color ratio and the multicriteria index. It means that we discuss the
apparent presence of water ice, but we have to keep in mind that more or less
prominent ice signature may have multiple origin and can not be directly
linked to the real amount of ice at the surface.
Using a ratio instead of absolute band values allow us to still be more
sensitive to water ice absorption. Fig. 14-a and -b show different patterns
between C18 and C15. Similarity still arises: on both cubes, at the beginning
of the channels, the $2.7/2.78$ ratio is lower than at the bottom. However the
differences are within the estimated error bars. Taking this into account, as
noted previously, there is either a variation of the grain size in the cube
indicating a larger grain size in Xanadu, a variation of abundance of water
ice, and/or a disparity in the dust layer thickness. In the case of sediment
transport, smaller and lighter grains travel generally further as less energy
is required to move the particles (Gyr & Hoyer, 2006). We then expect larger
grains on Xanadu at the channels origin, and that further along the river the
grain size decreases, as inferred here.

###  7.2. Seasonal cycle and dust storm

As noted previously, we observe differences between C15 and C18, either for
the haze and mist scale factors, or for $\delta$. Surface albedo are retrieved
with a specific setup of absorption and scattering properties in the
atmosphere, given in detail in previous sections. But, this should be
considered as an operational description, with some parameters determined by
previous studies, predicted by models or set by free parameters. Constraining
the atmosphere properties for retrieving the same surface albedo with C15 and
C18 would be possible. However, the atmosphere characteristics (particle
vertical and gas profiles, phase functions,…) may also change with time. The
two observations were taken 7 years apart, over the north spring equinox. Some
atmospheric properties probably have changed between these two dates,
especially the haze and mist vertical profiles (Vinatier et al., 2010;
Seignovert et al., 2020). This explains well why the amount of particles
needed to fit data changes between these two dates. But, the actual change in
haze and mist profiles can not be accounted in detail by simple parameters
such as the haze and mist factors ($F_{h}$ and $F_{m}$). The global shifts in
color ratio and in the multicriteria index between the two observation sets
can then be also explain by a real change in atmosphere properties.

The maps derived from observations also differs in some areas with obvious
changes in water ice signatures compared to the surrounding. For instance, a
strip of several pixels width at the northern edge of the maps retrieved with
C15 systematically has a stronger ice signature than the pixels a bit south.
This strip does not exist in C18. Notably, this strip corresponds to the
southern edge of Xanadu mountains, delineated in Figure 1 and 15. These
changes in color ratios and in $\delta$ indicate changes in the composition of
the uppermost surface layer.

Two specific processes, which are not mutually exclusive, may explain these
regional changes. First, the amount of deposited material rate of the aerosols
is low on Titan, but dust storms probably occurred in between 2008 and 2010
due to the equinox (Charnay2015; Rodriguez et al., 2018), and one particularly
at the East side of Xanadu. It could have influenced the amount of deposited
material in this period, and the number of large particles in the lower
atmosphere issued from surface erosion. Such an event may have deposited a
sufficiently thick layer of particles to change the layer’s color relative to
the surrounding. Secondly, Turtle et al. (2011) report an episode of heavy
rain during the equinoctial period in 2010, in the equatorial region. This
occurred at 20oS, but at an another longitude than our region of interest. The
rain flooded a wide region and produced an obvious and persistent change in
surface reflectivity. A similar or less intense flooding may have occurred
elsewhere in the same period of time. Rodriguez et al. (2011) report a
significant cloud coverage in the south tropical band between 2007 and 2010.
They report a cloud event in 2008 above 120o W and 20oS, comparable in term of
apparent area to the event reported by Turtle et al. (2011). Such an event can
probably alter or remove light sedimentary layers, eventually ice gravels or
small pebbles and change locally both surface reflectivity and color ratios.
Our maps, indeed, clearly show that an event changed the surface
characteristics of some area compared to their surrounding.

###  7.3. Dunes

The water-ice signature variations can help identifying geological structures,
by the shapes and values of those variations. Based on morphological
arguments, Lopes et al. (2019) classified the terrains in C15 and C18 as
”plains”, ”hummocky” for Xanadu mountains, ”dunes” for the yellowish structure
(at the bottom of Fig. 13-a and annotated in Fig. 15-a), and ”dry lakebed”
(Moore & Howard, 2010) for the RADAR-gray structure at the bottom left of the
dunes (around 25∘S,118.5∘W) located in Tui Regio. Hummocky terrains are
expected to be areas where the icy crust of Titan is more exposed than
elsewhere on Titan (i.e., covered by a thinner layer of organic material)
(Janssen et al., 2016), which is supported by our results on Xanadu. The so-
called ”dune unit” (Fig. 15-a), on the contrary, should be poor in water ice.
Yet, both the albedo (especially at 5 $\mu$m) and inferred delta parameter in
this area suggest the presence of a water ice component. The SAR image does
not resolve any individual dunes but there are most likely present at a scale
smaller than the RADAR resolution here and with interdune corridors enriched
in water ice (Barnes et al., 2008; Bonnefoy et al., 2016).

###  7.4. Precipitations and Alluvial fans

We find an accumulation of high water-ice signature at the end of the
channels. With arguments based on radar-brightness, Le Gall et al. (2010) have
suggested that the bed of these channels could be ice-rich. Unfortunately,
they are too narrow to be characterized in detail at VIMS resolution. Brossier
et al. (2018) interpret ’IR-blue’ VIMS patches like our region of interest as
possible accumulated erosional products from the icebed, that cover the layer
of organic materials, or eroded icy terrains where rainfall washed out the
organic layer. The revealed ice-rich signature features could be alluvial fans
where sediments have been transported and deposited from the mountains of
Xanadu. Such alluvial fan results from the abrupt decrease of the terrain’s
slope, inducing a slower stream flow and sediments deposit in the plains.
There are generally triangular-shaped, similar to the features discussed here,
with coarse sediments close to the mouth of the stream and smaller ones
further away. The channels that carry materials to the fans can be dry most of
the time, as this is most likely the case here and in desert regions on Earth,
and occasionally be fed occasionally by violent rainfalls or storms.

As discussed in a previous section, rain showers and cloud activity occur
during equinox at these latitudes (Turtle et al., 2011; Rodriguez et al.,
2011; Turtle et al., 2018). As further evidences, other VIMS cubes with lower
pixel resolution from the flyby T44 show a large and discontinuous cloud, with
an orientation W-E. Strikingly, Griffith et al. (2005) detected a rapid
evolution of mid-latitude clouds, that eventually dissipated as a rain fall
within a few hours. Such hydrologic events may be the source of liquid
alimenting sporadically the channels considered here. A high discharge of
liquid and a non-porous soil is needed to transport sediments over 300 km in a
dry region. This aspect is also discussed by Faulk et al. (2017) and Birch et
al. (2016), who have considered occasional intense rainfalls at 20∘S of
latitude, that can create alluvial fans.

The extrapolated topography from Corlies et al. (2017) does not show abrupt
slope decrease at this particular site (maybe because of a low resolution),
but Birch et al. (2016) do find equatorial fans in Xanadu region. However they
do not classify these blue dotted areas as fans. Radebaugh et al. (2018)
mention fans in our area of study (green dots in Fig. 1 and 15 ), but these
patterns appear to be much smaller than what we observe. Although, when we
look closely at our results in Fig. 15-b, we notice that in most cases the
green dots’ signature of water-ice is a bit higher than the area surrounding
the dots. Another origin may also explains the presence of these enhanced
water-ice signature at the foot of channels, between 17 and 19∘S.

###  7.5. Impact craters

To investigate the nature of the spots of high water-ice signature, it may be
relevant to put our zone of study in a broader context. Few craters impact
exist on Titan, compared to other moons (Wood et al., 2010; Neish et al.,
2013). This feature can be explained by the presence of a dense atmosphere,
and can be the consequence of and deposition by eolian and fluvial activities.
Resurfacing by tectonic or cryovolcanism may also play a role. Impact craters
on Titan, particularly if they are young, are windows to the crust and
therefore have ice-rich inner and outer rim (Janssen et al., 2016; Griffith et
al., 2019). Compared with the rest of Titan, their rims are radar-bright with
low emissivity at 2.2cm, like that of the hummock terrains (Lopes et al.,
2019). On Titan, low emissivity appears to be associated with water-ice
(Malaska et al., 2016; Janssen et al., 2016).

Xanadu is a complex zone with multiple ancient impact crater candidates (Wood
et al., 2010). The global mosaic map of VIMS-ISS (Seignovert et al., 2019) in
Fig. 17-a, reveals 2 giant circular shapes on Xanadu (crater #1 and #2 in this
Fig. 17). These features could be artifacts due to the enhanced resolution of
VIMS dataset in these areas compared to the rest of Xanadu. However, crater #1
is also visible on the emissivity map (Fig. 17-b) as a large-scale low
emissivity region. Brown et al. (2011) argue in favor of a giant crater called
XCF. The spots of enhanced water-ice downstream are located at the limit of a
circular feature that could correspond to the rim of crater #2 but there is no
noticeable ring-like feature on the emissivity map, probably due to the low
resolution of the bottom half part of the map for this potential crater.
Instead, they could also correspond to the inner ring of crater #1. Finally,
crater #4 may be the inner ring of crater #1.

Though, giant impacts are not required to create these putative rims. Smaller
potential craters exist in Xanadu. They are numbered form 3 to 7 in our region
of interest. We select potential craters when they check at least 2 of these
following criteria: a circular shape in Fig. 17-c, a low emissivity spot on
Fig. 17-b and a corresponding geological crater-like feature in Fig. 17-a.
Putative craters #6 and #7 may be associated with water-ice enhanced signal on
the rims. Putative craters #3, #4 and #5 are too far from the area of enhanced
water ice signal, but it would be interesting to investigate those features in
more details. Remarkably, high $\delta$ spots and channels often correspond to
value of medium emissivity.

Despite these correlations, we do not find them convincing enough to ascertain
that impact craters were the cause of the enhanced water-ice signatures at the
channels’ ends.

Figure 17.—  (a) Zoom of Titan’s ISS_VIMS map from Seignovert et al. (2019) on
Xanadu region,(b) Titan’s 2.2cm wavelength emissivity map from Janssen et al.
(2016), and (c) Titan RADAR map, with outlined potential crater areas. The
dashed rectangle is our zone of study, and the thin black line represent the
channels.

###  7.6. Delta and evaporites

MacKenzie & Barnes (2016) classify Tui Regio as an evaporitic terrain. This is
in good agreement with our results showing a very low water-ice signature
inside this area (i.e. the yellow dashed area in Fig. 15). This region is
suspected to be a paleo sea (Moore & Howard, 2010; Hofgartner et al., 2020).
While it is not incompatible with an old impact crater formation, as seas may
form in craters, it casts a doubt on the origin of the water-ice signature
enhancement at the border of Tui Regio, where the channels end. If the
channels where active at the same time of the sea, these features could be
delta-like formations, occurring when the stream is slowed down at the contact
of the sea, resulting in an accumulation of fine deposits at the margins of
Tui Regio. If these structures are deltas, it could imply that the channels
supplied the paleo-sea. One ice-rich (high $\delta$) area
(130$\pm$1oW,16$\pm$1.5oS) does not seem to end a channel, but it may be
because it was too old or too narrow to be detected on the RADAR map. If
precipitations are sufficient, rainstorms also could make small or subsidiary
channels not resolved with the RADAR resolution overflows and deposit washed
off sediments at the margins of the paleo-sea.

Another ice-rich area (119$\pm$2oW, 24$\pm$1.oS) seems to be connected to a
remote channel through a broad and diffuse river bed, not different but much
longer that the structures that connects the other channels to the ice-rich
zones. This is a possible old channel, highlighted in Fig. 15 (in purple
dotted line).

Our results also show a gradual increase of $\delta$ from the center towards
the borders of Tui Regio. This result is consistent with the pattern predicted
by Cordier et al. (2016), in their Fig. 10 showing a decrease of the
evaporitic layer’s thickness with the distance from the center of the sea.
However VIMS only probes the surface firsts tens of microns, so if there is a
thickness variation, it is on a micrometric scale. A more probable
interpretation would be a change of composition in the evaporitic layer due to
the different saturation values of the components in the liquid (Cordier et
al., 2013, 2016).

North to the deltas, depending on the observation, color on Xanadu mountains
are contrasted. In C15, Xanadu mountains appear with a high ice-index (forming
the ice-rich strip discussed previously). The plain at the foot of the
mountains has a low ice signal and it is then probably covered by sediment. In
C18, Xanadu mountains and the plains below appear to be quite uniform,
regardless the topography. Differences between the two observations, as
explained above, may be due to meteorological events or to the different
incident and emergent angles if the surface reflectivity is not isotropic as
assumed in our model.
The explanation of our ice index maps (ratios and $\delta$) with dry
(occasionnally active) river beds ending in delta formations at the edge of a
paleolake appears the most consistent with observations. These rivers seem to
flow from Xanadu mountain slopes. The site topology could have been formed by
impact craters but they do not seems to be the direct reason for ice rich
terrains to be apparent.

##  8\. Conclusion

In this work we used a model of photometry to account for atmosphere
scattering and absorption. The purpose is to retrieve the surface albedo and
to remove as much as possible the effects due to the atmosphere and the
observation geometries. The model setup is quite standard and derived from
previous models (Rannou et al., 2010). Among other specific characteristics,
we use a model of scattering by fractal aggregates (Rannou et al., 1997) to
calculate the optical properties of the haze. Our simulations show the need
for aerosols with a fractal dimension between 2.3 and 2.4, i.e., more compact
than initially assumed (2.0). This is obviously meaningful in term of aerosol
growth and in term of mechanical behavior in the atmosphere. This updated
fractal dimension reduces and flattens the surface albedo spectra. It also
changes its shape at 1.3 $\mu$m.

The haze opacity varies over large distances, and not always along the
latitude. The mist opacity factor changes over shorter distances than the haze
opacity factor, as expected since it is more sensitive to variations of the
lower atmosphere. Our study confirms that these two factors are not
correlated. Seasonal cycle may also increase the haze parameter, but further
studies are required to confirm the influence of a seasonal cycle on the haze
over several years.

The photometry model enabled us to retrieve the surface albedo in seven
methane windows, and for the pixels of two VIMS cubes probing the region
between Xanadu and Tui regio. This region is of some interest as it is assumed
to be an ancient sea. We established several criteria to detect water ice
signature (four ratios and a global index). The radiative transfer analysis
highlights areas of high water-ice signal on Xanadu mountains and also at the
far end of dry channels. Areas shapes of high water-ice signal at the stream
mouths are triangular as expected for alluvial fans, or deposits of deltas. In
the first case an abrupt slope change is needed, and in the second case the
channels supply the paleo-sea of Tui Regio. Impact craters with ice-rich rims
is a less likely explication for these features. The ice signal decreases
inside the channels, perhaps due to their sizes, smaller than the spatial
resolution of a VIMS pixel. Anyhow, our results support Barnes et al. (2007),
Langhans et al. (2012) and Brossier et al. (2018) hypothesis of sediment
transport from Xanadu mountains to the channels downstream, with a decreasing
grain size.

In addition, our analysis shows a gradient of low $\delta$ perfectly
superposed to Tui Regio. This gradient could correspond to the decreasing
thickness from the center to the border of the evaporite layer of Tui Regio,
as suggested by Cordier et al. (2016), or a variation of composition in the
evaporite layer. Here the overall shape of the albedo between a low and high
$\delta$ does not really change, except at 2.8 $\mu$m, i.e., in a domain where
the model is incomplete. Detailed analysis of the retrieved albedo could bring
further insight into the spectral composition of these terrains provided the
uncertainties are small enough. We used an original method to provide estimate
of the VIMS instrument intrinsic noise. This latter tends to increase with the
wavelength, with $\sim$ 1% at short wavelengths to 10-40% above 3.5$\mu$m. The
error also varies from one cube to another.

Comparing our results to previous global mapping and RGB composites (Fig. 2),
we see that the link between water ice and ’dark blue’, and ’bright blue’
areas is not as straightforward as it was previously thought (Barnes et al.,
2007). The ice signature can be modulated by grain size, dust cover that
darkens the albedo, or even the incident angles. Indeed, a light blue area in
Fig. 2 is actually associated with a low water ice signal in Fig. 15-a whereas
some black units corresponds to high values of $\delta$. Only the RT analysis
shows qualitatively the difference in intensity of water-ice signature between
those blue units.

A surface albedo map leads to a better understanding of the structures. It
allows removal of the atmospheric inputs and to correct effects of different
observation geometries - except for the surface photometric function since we
model surface reflectivity as a lambertian surface. These effects can deeply
alter the perception of the true surface features by inducing strong biases.
The different albedo ratios and the parameter $\delta$ that we consider in
this work also have a much better contrasts than a direct analysis of radiance
at the top of the atmosphere. They give a direct appreciation of the strength
of the water-ice signature. Nonetheless, we are conscious of the limits of
$\delta$. Since it is calculated using normalized ratios, meaning that for a
VIMS cube with absolutely no water-ice, $\delta$ could still reach 0.5. That
is why we must use it in addition to other methods of surface analysis. The
model can be improved to obtain identical absolute ratios for VIMS cubes in a
similar area, at the same time period and under different geometries. This is
possible, for example, with better aerosol phase functions and better
information about the vertical profile of haze and mist. More study is
required on the mist and haze to retrieve the absolute value of the contrasts
observed as water-ice signature. If our model is not perfect, it shows a
possible way to obtain more information about the surface in the future.

The terrains studied in this paper are of great interest in terms of
exobiology. Similar terrains on Earth, with sediment accumulation areas
(deltas, paleo-lakes) are places where bio-signatures are generally well
preserved (Meyers & Teranes, 2002). It would be particularly relevant for
Dragonfly to visit areas similar to the region studied in our work (Lorenz et
al., 2018). The mission will target zones with ice-rich material, which will
help alleviate the ambiguity of VIMS analysis by finding the composition of
the famous ’dark material’ often mixed with the ice (Barnes et al., 2020).

There are probably plenty of other areas like those identified but not
resolved by VIMS: the images from the descent of Huygens actually suggest the
existence of river networks that are undetectable by VIMS or the RADAR because
of their small dimensions (Owen, 2005).

Lastly, the approach proposed here can be applied to other regions in order to
study qualitatively the water-ice signature on Titan, like Shangri-La, future
landing site of Dragonfly and where Huygens landed. The averred flooded area
in Turtle et al. (2011) paper could also be of interest. The model needs no
specific adaptation for equatorial or tropical regions, but probably moderate
adaptation for mid-latitude regions (from $\pm$ 30o to 50o) because the haze
layer still has a vertical profile similar to the equatorial one. On the other
hand, to be able to study polar regions, we need to adapt the vertical profile
of aerosols and mist, and likely their spectral behavior.

## Data availability and final version

The RT analysis results on C14, C15 and C18 are available here:
https://doi.org/10.5281/zenodo.4675563. The final and peer-reviewed version of
this paper is available here: https://doi.org/10.1016/j.icarus.2021.114464.

## Acknowledgements

The present research was supported by the Programme National de Planétologie
(PNP) of CNRS-INSU co-funded by CNES, and also supported by the French HPC
Center ROMEO, the Région Grand-Est and the University of Reims Champagne-
Ardenne. Many thanks to Jérémie Burgalat for the IT help provided.

## References

  * Annex et al. (2020) Annex, A. M., Pearson, B., Seignovert, B., et al. 2020, Journal of Open Source Software, 5, 2050
  * Barnes et al. (2020) Barnes, J., Turtle, E., Trainer, M., et al. 2020, in American Astronomical Society Meeting Abstracts# 236, Vol. 236, 221–05
  * Barnes et al. (2005) Barnes, J. W., Brown, R. H., Turtle, E. P., et al. 2005, Science, 310, 92
  * Barnes et al. (2007) Barnes, J. W., Brown, R. H., Soderblom, L., et al. 2007, Icarus, 186, 242
  * Barnes et al. (2007) Barnes, J. W., Radebaugh, J., Brown, R. H., et al. 2007, Journal of Geophysical Research (Planets), 112, E11006
  * Barnes et al. (2008) Barnes, J. W., Brown, R. H., Soderblom, L., et al. 2008, Icarus, 195, 400
  * Barnes et al. (2011) Barnes, J. W., Bow, J., Schwartz, J., et al. 2011, Icarus, 216, 136
  * Barnes et al. (2018) Barnes, J. W., MacKenzie, S. M., Young, E. F., et al. 2018, The Astronomical Journal, 155, 264
  * Birch et al. (2016) Birch, S., Hayes, A., Howard, A., Moore, J., & Radebaugh, J. 2016, Icarus, 270, 238 , titan’s Surface and Atmosphere
  * Bonnefoy et al. (2016) Bonnefoy, L. E., Hayes, A. G., Hayne, P. O., et al. 2016, Icarus, 270, 222 , titan’s Surface and Atmosphere
  * Botet et al. (1995) Botet, R., Rannou, P., & Cabane, M. 1995, Journal of Physics A: Mathematical and General, 28, 297
  * Brossier et al. (2018) Brossier, J., Rodriguez, S., Cornet, T., et al. 2018, Journal of Geophysical Research: Planets, 123, 1089
  * Brown et al. (2011) Brown, R. H., Barnes, J. W., & Melosh, H. J. 2011, Icarus, 214, 556
  * Brown et al. (2004) Brown, R. H., Baines, K. H., Bellucci, G., et al. 2004, ssr, 115, 111
  * Charnay et al. (2015) Charnay, B., Barth, E., Rafkin, S., et al. 2015, Nature Geoscience, 8, 362
  * Clark et al. (2018) Clark, R. N., Brown, R. H., Lytle, D. M., & Hedman, M. 2018, NASA Planetary Data System, The Planetary Atmospheres Node
  * Cordier et al. (2013) Cordier, D., Barnes, J. W., & Ferreira, A. G. 2013, Icarus, 226, 1431
  * Cordier et al. (2016) Cordier, D., Cornet, T., Barnes, J. W., et al. 2016, Icarus, 270, 41–56
  * Corlies et al. (2017) Corlies, P., Hayes, A. G., Birch, S. P. D., et al. 2017, Geophys. Res. Lett., 44, 11
  * Cours et al. (2020) Cours, T., Cordier, D., Seignovert, B., Maltagliati, L., & Biennier, L. 2020, Icarus, 339, 113571
  * de Bergh et al. (2012) de Bergh, C., Courtin, R., Bézard, B., et al. 2012, Planetary and Space Science, 61, 85 , surfaces, atmospheres and magnetospheres of the outer planets and their satellites and ring systems: Part VII
  * Doose et al. (2016) Doose, L. R., Karkoschka, E., Tomasko, M. G., & Anderson, C. M. 2016, Icarus, 270, 355 , titan’s Surface and Atmosphere
  * Emde et al. (2016) Emde, C., Buras-Schnell, R., Kylling, A., et al. 2016, Geoscientific Model Development, 9, 1647
  * Evans (2007) Evans, K. F. 2007, J. Atmos. Sci., 64, 3854
  * Faulk et al. (2017) Faulk, S. P., Mitchell, J. L., Moon, S., & Lora, J. M. 2017, Nature Geoscience, 10, 827
  * Fulchignoni et al. (2005) Fulchignoni, M., Ferri, F., Angrilli, F., et al. 2005, Nature, 438, 785
  * Goody et al. (1989) Goody, R., West, R., Chen, L., & Crisp, D. 1989, J. Quant. Spectrosc. Radiat. Transfer, 42, 539
  * Griffith et al. (2012) Griffith, C. A., Lora, J. M., Turner, J., et al. 2012, Nature, 486, 237
  * Griffith et al. (2019) Griffith, C. A., Penteado, P. F., Turner, J. D., et al. 2019, Nat. Astron., 343
  * Griffith et al. (2005) Griffith, C. A., Penteado, P., Baines, K., et al. 2005, Sci, 310, 474
  * Gyr & Hoyer (2006) Gyr, A., & Hoyer, K. 2006, Sediment Transport, 1st edn. (Springer Netherlands)
  * Hansen & Travis (1974) Hansen, J. E., & Travis, L. D. 1974, Space Science Reviews, 16, 527
  * Hofgartner et al. (2020) Hofgartner, J. D., Hayes, A. G., Campbell, D. B., et al. 2020, Nat. Commun., 11, 2829
  * Janssen et al. (2016) Janssen, M., Gall, A. L., Lopes, R., et al. 2016, Icarus, 270, 443 , titan’s Surface and Atmosphere
  * Karkoschka (2016) Karkoschka, E. 2016, Icarus, 270, 339
  * Karkoschka et al. (2012) Karkoschka, E., Schröder, S. E., Tomasko, M. G., & Keller, H. U. 2012, Planetary and Space Science, 60, 342
  * Kazeminejad et al. (2011) Kazeminejad, B., Atkinson, D. H., & Lebreton, J.-P. 2011, Advances in Space Research, 47, 1622
  * Lafferty et al. (1996) Lafferty, W. J., Solodov, A. M., Weber, A., Olson, W. B., & Hartmann, J.-M. 1996, Appl. Opt., 35, 5911
  * Langhans et al. (2012) Langhans, M. H., Jaumann, R., Stephan, K., et al. 2012, Planet. Space Sci., 60, 34
  * Le Gall et al. (2010) Le Gall, A., Janssen, M. A., Paillou, P., et al. 2010, Icarus, 207, 948
  * Le Mouélic et al. (2019) Le Mouélic, S., Cornet, T., Rodriguez, S., et al. 2019, Icarus, 319, 121
  * Lellouch et al. (2003) Lellouch, E., Coustenis, A., Sebag, B., et al. 2003, Icarus, 162, 125
  * Levenberg (1944) Levenberg, K. 1944, Quarterly of Applied Mathematic, 2, 164
  * Lopes et al. (2019) Lopes, R., Malaska, M., Schoenfeld, A., et al. 2019, Nature Astronomy, 1
  * Lorenz et al. (2018) Lorenz, R. D., Turtle, E. P., Barnes, J. W., et al. 2018, Dragonfly: A Rotorcraft Lander Concept for Scientific Exploration at Titan, Tech. Rep. 3, Johns Hopkins APL
  * MacKenzie & Barnes (2016) MacKenzie, S., & Barnes, J. 2016, The Astrophysical Journal, 821, 0
  * MacKenzie et al. (2014) MacKenzie, S. M., Barnes, J. W., Sotin, C., et al. 2014, Icarus, 243, 191
  * Malaska et al. (2016) Malaska, M. J., Lopes, R. M., Williams, D. A., et al. 2016, Icarus, 270, 130 , titan’s Surface and Atmosphere
  * Maltagliati et al. (2015) Maltagliati, L., Bézard, B., Vinatier, S., et al. 2015, Icarus, 248, 1
  * Marquardt (1963) Marquardt, D. 1963, SIAM Journal on Applied Mathematics, 11, 431
  * Mayer & Kylling (2005) Mayer, B., & Kylling, A. 2005, Atmospheric Chemistry and Physics Discussions, 5, 1319
  * McCord et al. (2006) McCord, T. B., Hansen, G. B., Buratti, B. J., et al. 2006, Planet. Space Sci., 54, 1524
  * McKay et al. (1989) McKay, C. P., Pollack, J. B., & Courtin, R. 1989, Icarus, 80, 23
  * McKellar (1989) McKellar, A. 1989, Icarus, 80, 361
  * Meyers & Teranes (2002) Meyers, P. A., & Teranes, J. L. 2002, in Tracking environmental change using lake sediments (Springer), 239–269
  * Moore & Howard (2010) Moore, J. M., & Howard, A. D. 2010, Geophysical Research Letters, 37, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2010GL045234
  * Moré et al. (1980) Moré, J. J., Garbow, B. S., & Hillstrom, K. E. 1980, User Guide for MINPACK-1, Tech. Rep. ANL-80-74, pub-ANL, pub-ANL:adr
  * Neish et al. (2013) Neish, C. D., Kirk, R. L., Lorenz, R. D., et al. 2013, Icarus, 223, 82
  * Nicholson et al. (2020) Nicholson, P. D., Ansty, T., Hedman, M. M., et al. 2020, Icarus, 344, 113356, cassini Mission Science Results
  * Niemann et al. (2005) Niemann, H. B., Atreya, S. K., Bauer, S. J., et al. 2005, Nature, 438, 779
  * Niemann et al. (2010) Niemann, H. B., Atreya, S. K., Demick, J. E., et al. 2010, J. Geophys. Res., 115, E12006
  * Nixon et al. (2012) Nixon, C. A., Temelso, B., Vinatier, S., et al. 2012, The Astrophysical Journal, 749, 159
  * Owen (2005) Owen, T. 2005, Nature, 438, 756
  * Radebaugh et al. (2018) Radebaugh, J., Ventra, D., Lorenz, R. D., et al. 2018, Geological Society, London, Special Publications, 440, 281
  * Rannou et al. (1997) Rannou, P., Cabane, M., Botet, R., & Chassefière, E. 1997, Journal of Geophysical Research: Planets, 102, 10997
  * Rannou et al. (2010) Rannou, P., Cours, T., Le Mouélic, S., et al. 2010, Icarus, 208, 850
  * Rannou et al. (2004) Rannou, P., Hourdin, F., McKay, C. P., & Luz, D. 2004, Icarus, 170, 443
  * Rannou et al. (2016) Rannou, P., Toledo, D., Lavvas, P., et al. 2016, Icarus, 270, 291
  * Rey et al. (2018) Rey, M., Nikitin, A. V., Bézard, B., et al. 2018, Icarus, 303, 114
  * Rey et al. (2016) Rey, M., Nikitin, A. V., & Tyuterev, V. 2016, J. Mol. Spectrosc., 327, 138
  * Rey et al. (2017) Rey, M., Nikitin, A. V., & Tyuterev, V. G. 2017, ApJ, 847, aa8909
  * Robidel et al. (2020) Robidel, R., Mouélic], S. L., Tobie, G., et al. 2020, Icarus, 349, 113848
  * Rodriguez et al. (2006) Rodriguez, S., Le Mouélic, S., Sotin, C., et al. 2006, Planetary and Space Science, 54, 1510
  * Rodriguez et al. (2011) Rodriguez, S., Le Mouélic, S., Rannou, P., et al. 2011, Icarus, 216, 89
  * Rodriguez et al. (2018) Rodriguez, S., Le Mouélic, S., Barnes, J. W., et al. 2018, Nature Geoscience, 11, 727
  * Rothman et al. (2013) Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, J. Quant. Spectrosc. Radiat. Transfer, 130, 4
  * Schröder & Keller (2008) Schröder, S. E., & Keller, H. U. 2008, Planet. Space Sci., 56, 753
  * Seignovert et al. (2019) Seignovert, B., Le Mouélic, S., Brown, R. H., et al. 2019, Titan’s global map combining VIMS and ISS mosaics, doi:10.22002/D1.1173
  * Seignovert & Rannou (2019) Seignovert, B., & Rannou, P. 2019, PovRay model of fractal aerosols, doi:10.5281/zenodo.2559546
  * Seignovert et al. (2020) Seignovert, B., Rannou, P., West, R. A., & Vinatier, S. 2020, Icarus, In revision
  * Soderblom et al. (2007) Soderblom, L. A., Kirk, R. L., Lunine, J. I., et al. 2007, Planetary and Space Science, 55, 2025 , titan as seen from Huygens
  * Sotin et al. (2005) Sotin, C., Jaumann, R., Buratti, B. J., et al. 2005, Nature, 435, 786
  * Sromovsky & Fry (2010) Sromovsky, L., & Fry, P. 2010, Icarus, 210, 230
  * Stofan et al. (2007) Stofan, E. R., Elachi, C., Lunine, J. I., et al. 2007, Nature, 445, 61
  * Tomasko et al. (2008) Tomasko, M. G., Doose, L., Engel, S., et al. 2008, Planet. Space Sci., 56, 669
  * Turtle et al. (2018) Turtle, E., Perry, J., Barbara, J., et al. 2018, Geophysical Research Letters, 45, 5320
  * Turtle et al. (2009) Turtle, E. P., Perry, J. E., McEwen, A. S., et al. 2009, Geophys. Res. Lett., 36, 2204
  * Turtle et al. (2011) Turtle, E. P., Perry, J. E., Hayes, A. G., et al. 2011, Sci, 331, 1414
  * Vinatier et al. (2010) Vinatier, S., Bézard, B., Nixon, C. A., et al. 2010, Icarus, 205, 559
  * Vinatier et al. (2015) Vinatier, S., Bézard, B., Lebonnois, S., et al. 2015, Icarus, 250, 95
  * Wood et al. (2010) Wood, C. A., Lorenz, R., Kirk, R., et al. 2010, Icarus, 206, 334