Vertical distribution of dust in the martian atmosphere:
OMEGA/MEx limb observations

Emiliano d'Aversa, Fabrizio Oliva, Francesca Altieri, Giuseppe Sindoni,
Filippo Giacomo Carrozzo, Giancarlo Bellucci, François Forget, Anna
Geminale, Arnaud Mahieux, Shohei Aoki, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Emiliano d'Aversa, Fabrizio Oliva, Francesca Altieri, Giuseppe Sindoni, Filippo Giacomo Carrozzo, et al.. Vertical distribution of dust in the martian atmosphere: OMEGA/MEx limb observations. Icarus, 2022, 371, pp.114702. 10.1016/j.icarus.2021.114702 . hal-03457069

HAL Id: hal-03457069

https://hal.science/hal-03457069v1

Submitted on 25 Jul 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

# Vertical distribution of dust in the Martian atmosphere: OMEGA/MEx limb observations

Emiliano D'Aversa(1), Fabrizio Oliva(1), Francesca Altieri(1), Giuseppe Sindoni(2), Filippo

Giacomo Carrozzo(1), Giancarlo Bellucci(1), François Forget(5), Anna Geminale(1), Arnaud

Mahieux(3,4), Shohei Aoki(4), Marilena Amoroso(2)

(1) IAPS-INAF, Via del Fosso del Cavaliere 100, 00133 Rome, Italy.

(2) Italian Space Agency ASI, Rome, Italy.

(3) The University of Texas, Austin, TX, USA.

(4) IASB-BIRA, Bruxelles, Belgium.

(5) LMD, CNRS, Paris, France.

## Abstract

The OMEGA spectrometer on board Mars Express acquired several observations of the Martian limb that are still largely unexploited. Here, we explore the information content of these data in terms of abundance and size of the Martian atmospheric dust. We approach the dust retrievals by applying a full spherical multiple scattering Monte Carlo 1D radiative transfer code to model the 0.5-2.5 µm spectral range (VNIR and SWIR OMEGA channels), in order to derive dust effective radius and number density variation with altitude, between approx. 8 and 50 km. This is the first time this approach is applied to OMEGA limb data; therefore, we only present three study cases where water ice is below the detectability level, in order to focus on the methodological issues, assumptions, and performances before future more extensive applications. The model fully includes multiple scattering effects, which are known to be responsible for the coupling between the limb spectra taken at different altitudes and the surface. An extended three-dimensional modeling of the surface reflectance, acting as surface-shine for limb spectra, is developed. The VNIR channel is found useful in reducing the degeneracy of the radiative transfer solutions. The derived dust vertical distributions, yielding dust effective radii of the order of 0.85 ± 0.15 µm (corresponding to a modal radius $r_m$ of ~0.3 µm) between 15 and 30 km are quite in agreement with Global Circulation Models (GCMs), but the abundances lie about one order of magnitude above model predictions, as previously reported from comparison between models and Mars Climate Sounder (MCS) data. An overall agreement is in fact found with MCS data, and in one case the OMEGA-retrieved dust is compatible with a local storm rising from the Hellas basin. Such storms are not well represented in the Mars Climate Database, which provides monthly mean statistics. Our results demonstrate the capability of OMEGA limb data to quantitatively contribute to Martian dust studies, despite the complex and slow radiative transfer computation scheme required to accurately model multiple scattering in the probed spectral range. The desirable application of the retrieval method to the whole OMEGA limb dataset, also expected to help assess the occurrence of local dust storms, requires further work aimed to include water ice aerosols and possibly thermal emission. Being the first attempt to use a Monte Carlo modeling approach to OMEGA limb data, this work represents a useful benchmark for validating faster, although less accurate, radiative transfer models.
# 1. Introduction

Dust is one of the main components of the Martian environment, widespread throughout the planet's atmosphere and surface, with changing abundance and distribution in both space and time. Beyond their role as a passive tracer of atmospheric motion, Martian dust grains, similar as on Earth, influence the atmospheric energy balance through radiative feedback. General circulation models suggest that both abundance and size of dust grains strongly affect atmospheric temperatures and dynamics (Kahre et al. 2008; Madeleine et al. 2011). Moreover, the presence of suspended dust heavily influences radiative transfer processes in the atmosphere both in the visible and near infrared (IR) spectral ranges. Consequently, dust controls the measurement of gases abundance through remote sensing, which often requires accurate modelling of dust light scattering.

Although the majority of remote sensing data at Mars are obtained in nadir-looking geometry, limb observations were exploited over the years to gain insights about the vertical structure of Martian aerosols, since the beginning of the exploration of the planet (e.g. (Anderson and Leovy 1978; Chassefière et al. 1992). In fact, the thin Martian atmosphere makes the nadir observations usually poorly sensitive to the vertical location of dust grains, with the exception of optically thick conditions during dust storms, when the top altitude of the bulk dust load can be locally retrieved (e.g. Oliva et al. 2018), or when specific observing schemes are aimed to this purpose (for example using the Emission Phase Function, Clancy and Lee 1991; Clancy, Wolff, and Christensen 2003). On the contrary, limb observations are naturally more sensitive to the vertical atmospheric structure, at the expense of loss in horizontal resolution, allowing the investigation of both global-scale dust vertical distributions and local features such as dust layering, storms, and clouds (e.g. Jaquin, Gierasch, and Kahn 1986; Montmessin et al. 2006; Sánchez-Lavega et al. 2015, 2018).

Several instruments on spacecraft orbiting Mars are devoted to the spectroscopy of Martian limb, or acquired limb data as a byproduct of their main activities. In the thermal IR range, limb scans produced by the Mars Climate Sounder (MCS) and Thermal Emission Spectrometer (TES) instruments on board Mars Reconnaissance Orbiter (MRO) and Mars Global Surveyor (MGS), respectively, are used to systematically retrieve vertical profiles of dust density and grain size, useful from a climatological point of view (Kleinböhl et al. 2009; Schofield et al. 2018). Dust grain size distribution centered around 1.5-2.0 µm, in effective radius, $r_{\text{eff}}$, are often inferred at lower altitudes, in general agreement with the local measurements taken on the surface (Lemmon et al. 2004). Smaller sizes are usually detected at higher altitudes, with distributions centered around 0.5-1.0 µm $r_{\text{eff}}$.

The analysis of limb observations acquired in the UV range by SPICAM/MEx spectrometer demonstrated the presence of a smaller dust population, peaking between 10 and 100 nm (Rannou et al. 2006). Using a wider spectral range from UV to near IR, accessed through both SPICAM channels, a bimodal distribution of aerosol size was also inferred, with
a coarser mode populated by dust and water ice grains with average radii of 0.7 and 1.2 µm respectively, and a finer, Aitken particles, mode ($r < 0.1$ µm) (Fedorova et al. 2014).

The CRISM/MRO limb spectra, covering the visible and the near IR spectral range, were also studied in details, suggesting a global size distribution for the dust always around 1 µm $r_{eff}$ without significant vertical variations below 40 km (Guzewich, Smith, and Wolff 2014; Clancy et al. 2019).

It is worth noting that limb analyses always suffer from a much higher complexity level both in terms of modeling and quantitative retrievals. In fact, in general they cannot take advantage of the plane-parallel approximation usually adopted in the atmospheric radiative transfer calculations, and the inclusion of full light scattering processes in a spherical geometry is more complex and computationally demanding than in the plane-parallel approximation. This enhanced complexity explains the smaller amount of quantitative outcomes from limb compared to the nadir observations reported in the past decades, and the under-exploitation of significant datasets, such as the solar channels of MCS and TES instruments (Wolff and Pankine 2016) as well as the data by OMEGA/MEx (e.g. Gondet and Bibring 2018; Vincendon et al. 2011). More comprehensive analyses are emerging in recent years, stimulated by the development for Earth studies of more efficient tools and computational power (e.g. Eluszkiewicz et al. 2017; Rozanov et al. 2017).

In this paper we assess the potential of OMEGA/MEx limb observations in quantitative retrieval of the vertical distribution of the dust density and grain size. OMEGA hyperspectral data cover a wide spectral range, 0.4-5 µm, very sensitive to light scattering by dust and clouds, with a spatial sampling enabling vertical resolution at the limb between 2 and 12 km in most cases. We will describe in detail the OMEGA data content, the processing and the data selection in Section 2, while in Section 3 we discuss the modelling efforts. The adopted retrieval technique, applied to the solar range (0.5 to 2.5 µm) of the OMEGA limb spectra, is based on a Monte Carlo radiative transfer algorithm in a spherical shell atmospheric model. The model, aimed at retrieving the dust properties in the absence of water ice, was applied to selected OMEGA study cases only. We describe the procedures used for tuning all the parameters and evaluate the assumptions and the performances before wider applications.

## 2. The OMEGA observations of the Martian limb

### 2.1. Data description

OMEGA observations consist of hyperspectral images in the visible and near IR spectral ranges, covered by three separate spectrometers, spatially co-aligned and spectrally partially overlapped: the VNIR channel, covering wavelengths from 0.35 to 1.0 µm with a ~7 nm spectral resolution, the SWIR and LWIR channels, covering respectively from 0.93 to 2.69 µm and from 2.53 to 5.1 µm, with a ~20 nm spectral resolution. A detailed description of the instrument and its calibration can be found in Bibring et al. (2004). The raw digital number (DN) data, obtained from the ESA Planetary Science Archive (PSA) along with associated image navigation data, were calibrated to radiance factor (I/F) applying a standard
pipeline (Carrozzo et al. 2012; Bibring et al. 2005), after a despiking process. Anomalous signals can be found in all three OMEGA channels more and more frequently over time. These signal failures, as well as bad pixels, are manually identified and removed from the analysis.

While most OMEGA data are acquired in a nadir-looking geometry, the Martian limb appears in several cubes taken with the Mars Express spacecraft in inertial pointing mode. Limb geometry parameters provided in the PSA archive refer to the limb tangent points of each pixel, i.e. the point of the pixel's line-of-sight closest to the Martian surface. As explained below (Section 3.4), some aspects of our radiative transfer analysis (mainly related with the surface treatment) required a wider geometric reconstruction, in order to model the illumination to and from the surface region below each tangent point. In these cases *ad hoc* geometrical calculations were implemented by taking advantage of topography data from MOLA/MGS (Smith et al., 1998) and of JPL's NAIF SPICE kernels and library (Acton et al. 2018).

## 2.2. Data selection

The features of our radiative transfer approach (Section 3) strongly constrain the OMEGA data to which the retrieval is applicable. First of all, since thermal emissions are excluded from our analysis, we are only interested in the diurnal limb of the planet. Further constraints come from the requirements that the OMEGA cube contains meaningful data in all three channels and that water ice is not detectable over the whole atmospheric column (see Section 2.4).

These criteria strongly narrow down the dataset suitable for a robust retrieval of dust properties. In particular, the requirement of good-quality limb spectra in all three channels greatly favors the first (nominal, from 11/2004 till 06/2006) phase of the Mars Express mission, in which the instrumental performances were optimal. The absence of water ice requirement appears also to be strongly restrictive, since spectral features ascribable to water ice, with different degrees of purity, are very widespread in the OMEGA dataset (more details on this kind of selection are given in Section 2.4).

The cubes selected for the present investigation are listed in Table 1, along with their main parameters. Two of them (ORB2008_0 and ORB2130_0) were acquired at the beginning of the summer season (at mid and circumpolar southern latitudes during MY27), and the other one (ORB0072_2) in the middle of the Northern winter (northern mid-latitudes during MY26). The location of the limb tangent regions is shown in Figure 1. The cases ORB2008_0 and ORB2130_0 are quite close both in time (about one month apart), in latitude, and spatial resolution (6.5-7.6 km/pixel). The case ORB0072_2 is instead located on the opposite side of the planet and is characterized by a much better spatial resolution (2.3 km/pixel).

<table><caption>Table 1: Summary of the main parameters for the OMEGA observations selected for limb dust retrievals. All geometric parameters refer to the limb tangent points.</caption><thead><tr><td>obs id</td><td>date</td><td>Ls (°)</td><td>location</td><td>latit (°)</td><td>longit (°)</td><td>solar zenith</td><td>solar phase</td><td>local time</td><td>spatial resol</td></tr></thead></table>
<table><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th>angle (°)</th><th>angle (°)</th><th>(h)</th><th>(km)</th></tr></thead><tbody><tr><td>ORB0072_2</td><td>2004-Feb-01 (MY26)</td><td>343.2</td><td>Amazonis Planitia</td><td>32N</td><td>202E</td><td>63.9</td><td>30.6</td><td>15.5</td><td>2.3</td></tr><tr><td>ORB2008_0</td><td>2005-Aug-07 (MY27)</td><td>264.6</td><td>Hellas Planitia</td><td>40S</td><td>71E</td><td>25.2</td><td>96.8</td><td>13.8</td><td>6.5</td></tr><tr><td>ORB2130_0</td><td>2005-Sep-11 (MY27)</td><td>286.0</td><td>Noachis Terra</td><td>64S</td><td>14E</td><td>47.5</td><td>117.3</td><td>14.7</td><td>7.5</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 1 - Location of the limb tangent regions for the selected OMEGA limb scans listed in Table 1.</figcaption></figure>

## 2.3. Limb data processing

As described in Section 3, our analysis of limb spectra takes advantage of both VNIR and SWIR channels, which are known to be significantly spatially misaligned with respect to each other in the limb observation configuration. Three effects can be invoked as main sources to explain such misalignment: a) a physical tilt of the optical axes; b) a slight difference in pixel instantaneous field of view (IFOV) size; c) the effective large difference in exposure times, due to the different acquisition schemes (push-broom for VNIR, whisk-broom for IR). Both b) and c) effects yield to a wider point spread function for the VNIR channel, with a loss of spatial information with respect to SWIR. The co-registration method developed to obtain coherent VNIR-SWIR spectra of the Martian surface (Carrozzo et al. 2012) cannot work in the limb case, because the large variation of observing distance when approaching the planetary limb magnifies the misalignment and results in large discrepancies in tangent altitude between the two channels. In order to compensate for this effect, we
implement an *ad hoc* two-step method for limb observations. First, we maximize the correlation between a VNIR and SWIR image, at a similar wavelength and encompassing contrasted surface features, by applying a 2-dimensional shift and smearing of the SWIR pixels. Then, the best shift and smearing values obtained along the image axis closest to the limb-parallel direction are used to optimize the shift in the direction perpendicular to the limb (which corresponds to the most critical value for tangent altitudes) by matching the drops of reflectance while crossing the limb between the two channels. The aim of this method is to link each SWIR limb spectrum to a visible counterpart taken from the VNIR pixel whose footprint best matches the SWIR one. Spatial interpolations (smearing and shifting) are only used for allowing meaningful comparisons between VNIR and SWIR data, but are not applied to the final spectra, in which the IR part and its best-matching visible counterpart are kept unchanged.

A further issue affecting VNIR data is the presence of straylight, i.e. a spurious signal due to light scattered inside the spectrometer. Since this kind of contamination is additive in nature, its effect is mostly evident at the edges of the spectral range, where the signal level is lower, but it actually affects all VNIR wavelengths. A spectrally flat straylight level is obtained, after applying the spatial mapping between the two channels described above, as the spectral offset that has to be subtracted to each VNIR spectrum to match the SWIR level. The matching is calculated by a second-order extrapolation of the closest reliable parts of the IR spectrum to the overlapping range.

OMEGA signal-to-noise ratio (SNR) is of the order of 100 (Bibring et al. 2004), but uncertainties associated with the spectra are also affected by the absolute radiometric calibration accuracy. Final relative errors have been estimated to be of the order of 10% (Carrozzo et al., 2012), even though comparisons with TES suggested errors close to 3% at worse (Vincendon et al. 2015). In our case, intrinsic instrumental errors must be combined with other sources of uncertainties, such as tangent altitude binning, co-registration and straylight removal (the latter two affecting the VNIR error level). Hence, the actual uncertainties associated with limb spectra have to be evaluated case-by-case. Generally speaking, uncertainties are found to be non-linearly dependent on the signal level, such that the relative error still depends on the reflectance level itself. This seems particularly true for the low signal levels usually characterizing limb spectra. For SWIR, the uncertainties span the range 2-5% for most of the considered limb scans, reaching 10% value only for very low reflectance levels (<0.04). The VNIR case is a bit worse, with errors between 5% and 10% (or even up to 20% in some case) for most points of the scans.

In order to quantitatively analyse the OMEGA limb spectra, we extracted from each cube under investigation one or more series of spectra at different tangent altitudes, forming a so-called "spectral scan" of the limb. Given the very low reflectance values ($10^{-2}-10^{-1}$) typically pertaining to limb signal, we improved the SNR of the spectra by averaging individual spectra over small limb regions. The vertical boundaries of these regions are determined by the pixel resolution. The horizontal boundaries are chosen on the basis of the signal variability along the limb, constrained to be less than the measurement uncertainty, and in any case with a limited variation in latitude and longitude (5°), ensuring that the final set of spectra is actually a vertical scan of the same atmospheric column (Figure 2).
<figure><img src="image_4.png" /><figcaption>Figure 2: Illustration of a limb scan mapping from a portion of an OMEGA cube. Surface coordinates grid (thin black lines), planetary limb position (thick black curve), tangent altitude isolines (white curves, labelled in km) and directions of vertical atmospheric columns (dashed lines) are overlapped to the 1.1 µm OMEGA-SWIR reflectance image (grey-scale). Individual pixels assigned to a given limb scan are filled with colors associated with the tangent altitude bins detailed in the legend. The line-filled region encompasses those OMEGA pixels, excluded from the scan, whose footprint is likely crossing the limb. The bottom-right panel shows a larger view of the cube where the limb position is highlighted by the red rectangle.</figcaption></figure>

A typical result of this processing is a spectral scan consisting of 10 to 20 reflectance spectra covering tangent altitudes below about 50 km. The spatial resolution is variable, with most of the data ranging within 2 and 13 km/pixel, depending on the spacecraft altitude over the planetary surface (resolution at limb is typically 30% worse than in the nadir-pointing case, due to the greater distance of the limb from the observer).

## 2.4. Spectral characteristics of the limb observations

Typical limb reflectance spectra (some examples are shown in Figure 3) show a peak value at wavelengths around 0.65-0.75 µm, within the VNIR range, a positive slope shortward and a negative slope longward up to about 2.5 µm. The shape in the LWIR range is more variable, and it is driven by thermal emission (by gases and aerosols) longward of 4 µm. Several gaseous bands are noticeable in the OMEGA spectra, whose identification is highlighted in the same Figure 3. The strongest CO₂ absorption bands are evident around 4.3 µm (saturated band often at all altitudes), 2.7 µm (saturated below ~30 km), and 2 µm (unsaturated). Minor CO₂ bands are spread in the whole range, and show weak absorptions in the OMEGA spectra. A CO band is clearly seen in the 4.6-4.7 µm range, whereas a small H₂O band is often noticeable at 2.6 µm (Melchiorri et al. 2007; Maltagliati et al. 2011), very
262
263

close to the 2.7 µm CO₂ band shoulder. The non-LTE O₂ emission band at 1.27µm appears strongly variable in intensity among the data (Altieri et al. 2009).

<figure><img src="image_5.png" /></figure>

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289

Figure 3: Examples of the variability of the spectral shape among OMEGA limb reflectance measurements (colored lines with error bars), taken from different data cubes at different altitudes. For the sake of comparison, all spectra are normalized to their maxima in the visible range. In the upper part of the plot the main spectral features of several atmospheric components are shown as black lines: the CO₂ strong (thick line) and weak (thin line, signal multiplied by 10) absorption bands, the CO bands around 4.6 µm, the locations of emission/absorption bands by O₂, CO₂ and H₂O are indicated by arrows at 1.27 µm, 4.3 µm, and 2.6 µm respectively. Moreover, the dashed and dot-dashed lines illustrate the (normalized) single-scattering reflectance of 1.0 µm effective radius grains of dust and water ice respectively. Dashed areas indicate the overlapping ranges among OMEGA channels. The red arrow at the bottom-right indicates the spectral range dominated by thermal emission, whereas the black arrow at the upper-left denotes the spectral range used for the dust retrievals in this work.

The overall limb reflectance in the 0.4-4.0 µm range is mainly due to light scattering by aerosols present in the Martian atmosphere. The observed reflectance values are strictly related to the aerosol optical depth in the probed atmospheric column, but the retrieval of the aerosol number density profile implies the knowledge of the aerosols scattering efficiency, i.e. its microphysics. Both dust and water ice, whose abundance and grain sizes are greatly variable in space and time, are expected to significantly contribute to the limb scattering. The single scattering properties of the aerosols, evaluated through Mie theory, suggest that the spectral slope in the 0.7-2.5 µm range becomes more and more negative with decreasing particle radius. No clear discrimination is possible in that region between dusty and icy composition, leaving this spectral range almost bereft of diagnostic features (the water ice bands at 1.5 and 2.0 µm become noticeable only for uncommonly large grain sizes). On the
contrary, the short-wavelength side 0.4-0.7 µm is affected by a dust extinction much larger
than water ice, but the same behavior is also common to the spectral albedo of many Martian
terrains, which can significantly affect limb spectra through secondary illumination (see
Section 3.4.3 for more details). These considerations prevent the definition of general spectral
parameters in the VNIR+SWIR range for a quick classification of limb spectra in terms of
their physical properties.

On the other hand, both dust and water ice show diagnostic features between 2.7 and
4.0 µm, within the OMEGA LWIR channel. Longward of the 2.7 µm band, the spectral slope
is positive for both dust and grains of water ice comparable in size with the wavelength (as in
both the orange and blue spectra in Figure 3). However, in the case of water ice, a double-
peak feature may appear, composed of the narrow Fresnel peak near 3.1 µm and a wider peak
in the 3.4-3.6 µm range. This double-peak merges in a unique peaked structure for ice grains
much smaller than a given wavelength and, as pointed out by several authors (e.g. Vincendon
et al. 2011), in such a case, the peak position (variable between 3.2 and 3.6 µm) is a
diagnostic feature of the grain size in high altitude clouds (this is actually what is observed in
the grey and green spectra in Figure 3). In the general case, the interplay of dust and ice
scattering appears as a “mixed” spectral shape (as for the red spectrum in Figure 3) and
explains the enhanced variability of the spectral shape in this wavelength range. We used this
spectral range for evaluating the presence of water ice in the sounded atmospheric column. In
particular, in the selection process of OMEGA cubes suitable for applying the retrieval
algorithm, which is described in Section 3, we rule out all the limb scans that, at any tangent
altitude, either include a significant change of spectral slope in the 3-4 µm range or show a
Fresnel peak. The combination of these conditions ensures that water ice is not present in the
complete atmospheric column sounded by the limb scan, at least at the OMEGA sensitivity
level.

A more comprehensive evaluation of the information content of limb spectra, as well
as any quantitative retrieval of aerosols properties and their vertical variation, needs full
radiative transfer considerations. In fact, the coupling among different atmospheric layers
introduced by multiple scattering enhances the correlation of a given limb spectrum with all
the other spectra above or below the same limb scan. It induces a the mixing of the
information contained at different altitudes dependent on the total abundance of scattering
particles in the atmospheric column. This is clearly more significant in case of enhanced
aerosol optical depth (e.g. storms).

3. Modelling

**3.1. Predictions from Global Climate Models**

Global Climate models, like the LMD Mars GCM (MCD, Forget et al. 1999;
Madeleine et al. 2011), offer an estimate of the amount of dust we can expect in the
atmosphere during the selected OMEGA observations. In this paper we access the GCM
prediction through the Mars Climate Database version 5.2 (Millour et al. 2018) which
provides monthly mean statistics of the fields simulated by the GCM. In all three cases a
332

moderate dust content is predicted by the MCD, with columnar visible (700 nm wavelength) optical depths ranging between $\tau_{700}=0.2$ and $\tau_{700}=0.6$, coherently with the climatological scenarios given by Montabone et al. (2015). It is worth stressing that these optical depth values encompass the full content of the vertical atmospheric column, which is not the quantity the limb observations are directly sensitive to. The atmospheric curvature is such that the line of sight of a limb observation at a given tangent altitude crosses the overhead atmospheric layers at locations even hundreds of kilometers away from the tangent point. Moreover, the path lengths of the line of sight in each atmospheric layer are much longer in the case of a limb observation than for a nadir slant case, especially for low altitudes limb observations. These issues must be taken into account when comparing quantities retrieved from limb observations with those derived from atmospheric models. In particular, horizontal spatial average, i.e. in latitude and longitude, has to be included in order to simulate the limb viewing conditions. In our case, we obtain model-predicted vertical profiles of dust number density and grain effective radius by extracting from the MCD the vertical profiles of dust mass mixing ratio and dust effective radius over a wide atmospheric section across tangent points, in the given observing conditions of each study case. These 3D dust distributions are then applied to a spherical stacking of layers and averaged over the lines-of-sight defined for every tangent altitude involved in the observations. The results of this process, shown in Figure 4, are 1D profiles representing the model-based predictions closest to the dust-related quantities observable when scanning the planetary limb. A pure feldspar ($2.6 \text{ g/cm}^3$ mass density) composition is assumed for converting dust mass densities to number densities. The variability associated with the profiles (filled areas in the same Figure 4) is a representation of the horizontal variations of the modelled atmosphere along the averaging lines-of-sight.

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

From this analysis, we see that the number density values are expected to span several orders of magnitude in the sounded altitude range. For all three cases, the model predicts an amount of dust between 1 and $10 \text{ particles/cm}^3$ at lowest altitudes (<10-15 km), comprising 1.2-1.7 µm effective radius grains. A decrease of dust effective size with altitude is moreover expected (more pronounced for the ORB2130_0 case), reaching effective radii as small as 0.1-0.2 µm at the top of the probed atmosphere.

361

362

363

364

365

These quantities define the overall parameters space where OMEGA data analysis is expected to take place, since they represent our knowledge about dust properties prior to the OMEGA limb retrievals. They are therefore useful for defining the *a priori* values needed for initializing the inversion algorithm, as will be described in Section 3.4.1.
<figure><img src="image_6.png" /></figure>

Figure 4: Predicted vertical profiles of dust number density (left panel) and dust grain effective radius (right panel) derived from the Mars Climate Database (MCD v5.2) for the observing conditions of the selected OMEGA cubes (using MCD standard dust scenarios based on Martian year). Dashed curves represent the values after a spatial average process along the lines-of-sight, aimed to approximate the limb viewing conditions. Filled areas indicate the variability range associated with the variability of atmospheric parameters along the averaging lines-of-sight.

## 3.2. The Radiative Transfer (RT) scheme and forward model

The treatment of radiative transfer in limb geometry turns out to be complicated by the fact that light scattering through the atmosphere has to be treated in a spherical geometry. Since multiple scattering is routinely managed within most of the plane parallel radiative transfer (RT) codes, pseudo-spherical approximations were developed over the years for enabling the usage of these codes in a limb geometry setup. This can be achieved by considering the line-of-sight as a curved path in a plane parallel atmosphere, and consequently adapting the path lengths within each atmospheric layer. These methods were tested for specific datasets and used for Mars applications (see e.g. Smith et al. 2013 for the CRISM data). However, in order to apply this approach to a new dataset, a validation procedure against a more robust method is desirable.

The most reliable approach for this kind of calculation is the Monte Carlo technique, which enables the description of the physical processes involved in the RT with an accuracy only limited by computational cost. We take advantage of the backward Monte Carlo RT solver MYSTIC (Mayer 2009), developed and validated for Earth studies, for tracing photons from the observer back to the source and calculating the RT forward model. The solver is incorporated in the MITRA RT tool, already exploited to study Mars (Oliva et al. 2018), Saturn (Oliva et al. 2016) and Jupiter (Sindoni et al. 2017) spectroscopic data. We use MYSTIC in a 1D configuration, i.e. discretizing the atmosphere in a stack of horizontally homogeneous spherical shells. This assumption implies some approximation, since some degree of aerosol horizontal variability is quite common in the Martian atmosphere. Moreover, limb observations come along with an intrinsic mixing of information in the vertical and horizontal directions in the atmosphere, due to both spherical geometry and
multiple scattering, and the choice of the best modelling configuration is not straightforward. It is also interesting to note that simultaneous nadir and limb measurements at the concurrent location, although very hard to implement, could in principle help disentangling these kinds of inhomogeneities. In any case, also the retrieval of the same atmospheric component from limb and nadir views can still diverge since their radiative transfer weighting functions along the lines-of-sight might be very different. We consider the 1D RT modelling configuration as the best trade-off in most cases, given the low level of environmental information. Further details on such mixing of information, levels of approximations, and their consequences are given in Sections 3.3, 3.4.4, and 3.4.5.

The main drawback of using the MYSTIC solver is that it is not yet fully validated for thermally emitted photons (C. Emde, personal communication). Therefore, although the OMEGA LWIR channel would add valuable information on dust properties through its thermal emission, we keep this channel and the associated wavelength range out of the analysis. Since we also discard the UV side of the VNIR channel because of strong straylight contamination, our retrievals are limited to the 0.5-2.5 µm spectral range. This limitation makes the retrieval insensitive to the atmospheric temperature. Consequently, thermal profiles are assumed fixed to the values obtained from the global climate model via the MCD for the observing condition of each data cube. These profiles are used for calculating both atmospheric densities and gas absorption, the latter being included in RT modelling and testing, but without actually contributing to the results shown in this paper. Gaseous absorption coefficients for CO₂ and H₂O molecules and single scattering properties for dust and water ice are calculated using the look-up tables and interpolation method described in Oliva et al. (2019). They are based on line-by-line technique for gases, using spectroscopic data from HITRAN 2012, and the Mie model for scattering by spherical particles, using the dust refraction index from Wolff et al. (2009) and the water ice refraction index from Warren and Brandt (2008). The vertical profiles of molecule abundances are also derived from MCD and kept fixed during a retrieval. Gas scattering is included by means of Rayleigh cross sections linked to the gas density profile, but atmospheric refraction is neglected given the low total density level encountered on Mars. The surface is treated in the Monte Carlo code as a Lambertian reflector, with a fixed spectral albedo. Since the true surface reflectance of Mars can significantly change in the large region probed by the limb observations, due to both composition and topography, we adopted an albedo definition based on an average evaluation of the effect of the diffuse “surface shine”, which will be detailed in Section 3.5.

### 3.3. Multiple scattering effects

The use of a Monte Carlo RT model enables an easy quantification of the relevance of multiple scattering in limb geometry, by evaluating the fraction of the synthetic radiance field associated with photons which undergo a single interaction in the atmosphere (without including interactions with the surface). As shown in Figure 5, even in a low dust abundance case ($\tau_{700} \sim 0.05$) significant fractions of the observed signal (e.g. from 10% at $\lambda \gtrsim 2$ µm above 33 km, to more than 50% at 0.7 µm below 10 km) are carried by multiply scattered photons. Significant amount of multiple scattering is found at all altitudes in the investigated spectral
442
443
444
445

range and this percentage grows even larger in dustier conditions. This fact strongly discourages the use of retrieval schemes based on independent layers, such as the onion peeling approach, which cannot deal with the large correlation length of the scattering processes.

<figure><img src="image_7.png" /></figure>

446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466

<caption>Figure 5 - Evaluation of the effect of multiple scattering in limb spectra. The contour plot illustrates the fraction in percent of the emerging simulated radiance field due to photons undergoing more than one interaction with aerosols (multiply scattered) between 0.5 and 2.5 µm. The estimation is obtained under low-dust atmospheric condition, not accounting for scattering by the planetary surface.</caption>

Solar photons scattered by the planet surface can give a significant contribution to the radiance observed in limb scans. This is also a kind of multiple scattering, since solar photons scattered by the surface must have scattered again in the atmosphere toward the detector to contribute to the limb observation. In principle, a wide surface region underlying the atmospheric path of the line-of-sight is expected to contribute to a measurement at a given tangent altitude. Before discussing the approach for defining a surface albedo spectrum effective for limb observations (Section 3.5), we can evaluate the significance of this kind of secondary illumination by comparing synthetic limb spectra obtained with and without a surface albedo. Figure 6, corresponding to a model atmosphere with low dust content (columnar $\tau<1$), shows that the percentage of the radiance associated with photons undergoing at least one surface interaction can be as much as 20-40% at almost all wavelengths and altitudes. This fraction begins to decrease when the dust optical thickness increases beyond some threshold value, and in very dusty conditions the surface contribution gradually vanishes. The exact threshold value above which the dust completely masks the
surface contribution is dependent on wavelength as well as on the geometry of the observation, and must be established on a case-by-case basis. In any way, we cannot know a priori the dust content contributing to a limb spectrum; hence, a comprehensive evaluation of the surface spectral albedo is needed. Further discussions on the accuracy required in treating the surface can be found in Section 3.4.3, and the final approach followed in this work is detailed in Section 3.5.

<figure><img src="image_8.png" /><figcaption>Figure 6: Weight in percent of the surface albedo on the limb reflectance, in the low-dust condition ($\tau \le 1$). Contour lines represent the percentage of the simulated radiance due to photons interacting at least once with the planetary surface.</figcaption></figure>

## 3.4. Retrieval method and performances

The retrieval of aerosol-related quantities is performed by coupling the forward RT model described above with a Bayesian inversion algorithm, based on the iterative Gauss-Newton method with $\chi^2$ minimization (Rodgers, 2000). The correlation induced by multiple scattering among different tangent altitudes forces the retrieval simultaneously on all spectra of a given limb scan. The same retrieval algorithm was already validated and successfully applied to several datasets on Mars and Saturn, and further descriptions can be found in Oliva et al. (2018; 2016).

Retrieved quantities include aerosol particles number density and effective radius for each spectrum. Particle size distributions are assumed lognormal in shape, with an effective variance $v_{eff}$ constant with altitude. After testing $v_{eff}$ values between 0.1 and 0.5 (at a 0.1 step) a $v_{eff} = 0.5$ value has been chosen as yielding smaller $\chi^2$ values when fitting the spectral scans (the same value is also adopted in the MCD).

The retrieval process for a single limb scan is quite time consuming, since the Jacobian matrices needed for inversion are computed numerically and need to be updated at
each iteration step. Computational time is found strongly dependent on (1) the number of photons $N_γ$ initializing the Monte Carlo code, (2) the number $N_b$ of the measured spectral points, (3) the number of phase function expansion coefficients and (4) the number of iterations needed for reaching convergence. On the other hand, the calculation time for gas opacity is by far negligible being performed only once per limb scan, since the gas column is not fitted in our retrievals. The current release of MYSTIC takes advantage of the “Variance Reduction Optimal Options Method” (VROOM, Buras and Mayer 2011) for speeding up simulations but does not implement truncation methods for the phase function (e.g. $δ$-M, Wiscombe 1977), forcing to compute a large number of expansion coefficients to correctly describe the aerosols phase functions. This is especially true for materials with enhanced opposition peaks, like water ice, which are found much more computationally demanding than dust. The number of simulated photons $N_γ$ affects the amplitude of the fluctuations of the simulated radiant energy in a given interval $Δλ$, i.e. for a given spectral point. Increasing the number of adjacent spectral points $N_b$ is useful to accurately estimate the statistics of these fluctuations, hence ensuring keeping them below the uncertainty level of the measurements. We found values of $10^6$ photons in 158 spectral points as a good trade-off between calculation time and accuracy of the simulations. The spectral selection does not imply a change in spectral resolution, but simply the removal from each OMEGA spectrum of some spectral points or some specific ranges (e.g. the range across the 2 µm CO$_2$ band).

Our tests, executed in a parallel calculation scheme, on highly performant 24-core machines, take an average time of 2-3 hours for a single iteration for a limb scan composed of 8 spectra. Up to 38-hour calculation time is needed in some cases for a single limb scan fitting. Calculations involving both dust and water ice profiles as free parameters are at least 10 times more time expensive making the treatment of water ice hardly feasible in the current implementation. Further work is needed to find optimal configurations for simultaneous retrievals of both dust and water ice, out of the aim of the present work.

### 3.4.1. A priori vertical profiles

Since the whole retrieval process unavoidably relies on several assumptions, we performed extensive testing of the sensitivity of the results to model and algorithm parameters. First, we validated the retrieval module in reproducing the same aerosol profiles used by the forward model to produce a set of synthetic spectra. These closed-loop-tests demonstrated a rather high level of sensitivity of the retrievals to the *a priori* profiles, for both density and size of aerosol particles. Far from being a mere feature of the algorithm, this sensitivity is indeed intrinsic to the atmospheric physical model itself: as above mentioned, the level of degeneracy in the radiative transfer problem in the investigated spectral range is high enough to let different solutions be compatible with the same limb spectra within the OMEGA accuracy. Hence, after testing several families of *a priori* profiles, we found the best results in terms of absolute $\chi^2$ when adopting vertically constant *a priori* profiles for both dust mixing ratio and grain size.

The same approach was therefore followed in application to OMEGA data. We took advantage of dust profiles predicted through MCD, as explained in Section 2.2 and shown as an example in Figure 2, to define vertically constant values to be assigned as *a priori* profiles.
539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

This choice implies that the retrievals do not rely upon the shape of MCD profiles and cannot be therefore biased in that direction. Only the expected orders of magnitude of both dust density and size are introduced in the process. In principle, one cannot completely exclude some bias of the results towards profiles more flattened than the true ones, but this represents the minimal information content that can be introduced in the process. The confidence level associated with the *a priori* has been derived from the horizontal fluctuations of the MCD values in the region of the atmosphere crossed by the limb lines-of-sight. The effectiveness of this kind of *a priori* has been verified in closed-loop tests, where we fit a set of synthetic limb spectra, corresponding to predefined dust profiles, by means of the retrieval algorithm, and the accuracy of retrieved profiles can easily be judged by comparing them with the original ones.

### 3.4.2. Dust phase function and VNIR channel

In calculating the scattering properties of Martian dust, we used the Mie theory applied to the refraction indices by Wolff et al. (2009). Since the actual phase function of dust is not yet clearly assessed (at least for the backscattering peak at visible wavelengths, see Tomasko et al. 1999), we tested our sensitivity to this parameter as well. As expected, assuming that the whole phase function shape is probed by multiple scattered photons only, its change affects the retrievals mainly at lower altitudes. In fact, using an empirically truncated phase function (as suggested e.g. by Vincendon et al. 2015) resulted in significant changes of radius and abundances (up to 1 order of magnitude) at the lowest altitudes of the OMEGA limb scans under investigation. However, a clear worsening of the quality of the spectral fits in the VNIR range is resulting, being fit residuals at VNIR wavelengths significantly larger than measurement uncertainties. On the contrary, both kinds of phase function are able to fit SWIR data equally well. It is worth stressing that if we only considered the SWIR channel, we could not discriminate between the two solutions. Hence, including the VNIR range to the analysis does not only have the effect of reducing final uncertainties on the retrieved parameters, but it also helps resolving some degrees of degeneracy of the solution.

Figure 7 illustrates this fact in a simplified case: a given individual limb OMEGA SWIR spectrum can be fitted equally well, within the given uncertainties, by two synthetic spectra with different dust effective radius (with lognormal distribution with the same effective variance, and hence different number density). However, only one of the two solutions is effective in fitting also the VNIR counterpart of the same spectrum, and this fact is clearly demonstrating the added value of the VNIR spectrum in resolving some kinds of degeneracies. This specific test is performed on a high altitude limb spectrum (at 41 km from the ORB0072_2 scan) in order to minimize multiple scattering effects and the influence of other spectra in both the measurement and the model.
<figure><img src="image_9.png" /><figcaption>Figure 7: An OMEGA limb spectrum (black line with error bars) is here used to illustrate how the VNIR channel helps in better constraining the dust retrievals. The green line represents a synthetic spectrum obtained as a best fit of the whole VNIR+SWIR spectrum, yielding r<sub>eff</sub>=0.62 µm. The red line is the same result using the SWIR channel only, which yields r<sub>eff</sub>=0.85 µm. Both models (corresponding to different dust effective radii and hence number densities but to the same kind of lognormal distribution with the same effective variance) are able to reproduce the IR part equally well within the given uncertainties, but only one is also well fitting the visible range.</figcaption></figure>

3.4.3. Sensitivity to surface albedo

As mentioned in the previous section, the reflectivity of the planetary surface represents a significant secondary source of illumination in limb geometry. The region of the surface actually contributing as a source can be rather wide and heterogeneous, depending on the local conditions of illumination and the altimetry. From the spectral point of view, eventual surface spectral features (e.g. the reddening slope in the visible or the 1 µm mafic band in the IR) propagates towards atmospheric aerosols and is observed in the limb spectra as well. Thus, an improper definition of the surface spectral albedo can in principle bias the retrievals, forcing to ascribe all the spectral features to aerosols scattering.

We found that different assumptions on the surface albedo can allow equally good spectral fits with quite different profiles of retrieved parameters, especially at low altitudes. As mentioned above (Section 3.3), the impact of this surface contribution is dependent on the wavelength, on the overall dust content of the atmospheric column, and on the geometry of the observation. A first order sensitivity to the surface albedo in a given observing condition can be inferred by applying the retrieval process to synthetic spectra obtained from predefined dust profiles. Since the “true” dust profiles are known in this case, it is straightforward to evaluate how much of the retrieved dust profile diverges from the true one
607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

while changing the surface albedo. For the geometry of the 0072_2 observation, we found that, in low dust condition (columnar $\tau_{700}$ ~ 0.1), at a wavelength-averaged albedo percentage variation $\delta a$ corresponds a similar percentage variation of retrieved single-layer optical depth $\delta\tau$ ~ $\delta a$, while variations of effective radius up to $\delta r_{eff}$ ~ 0.5 $\delta a$ and of number density up to $\delta n$ ~ 1.5 $\delta a$ can be found. These results, although attenuated for higher columnar dust contents, which is not known a priori in an observation, are sufficient to highlight that an independent and accurate definition of the surface spectrum is needed to keep the final uncertainties under control. The method we applied for defining a surface spectral albedo suitable for limb observations will be described in Section 3.5.

627

628

629

630

631

632

633

634

635

Every parameter used to construct the 1D atmospheric structure, and hence the forward models, comes with its own uncertainty. Anyway, it is worth noting that the vertical location of the limb scans exclusively relies on geometric calculations. The accuracy of the vertical scale is constrained by the OMEGA angular resolution (1.2 mrad/pixel, worsened by about 50% in the vertical binning process) and by the overall spacecraft pointing accuracy (estimated ~0.7 mrad in inertial attitude, ESOC personal communication). This is also true for the GCM, which, even if internally working on pressure-based grids, is always enquired at fixed scales of altitudes.

Another point worth stressing is the fact that the aerosols number density is always managed in absolute units (number of particles per unit volume), making the retrieved dust profiles largely insensitive to the uncertainties on atmospheric pressure and temperature. Anyway, the P/T profile is involved in the conversion of dust abundance from number density to mass mixing ratio and vice versa, used for allowing the comparison of the results with other data (see Section 4). In this case, a 1D gas density profile has been derived by averaging MCD values within the atmospheric plane interested by the limb lines-of-sight, while the horizontal fluctuations of the parameters are used to estimate the associated uncertainties (analogously to the dust averaging process illustrated in Figure 4).

636

637

638

639

640

641

642

643

644

The evaluation of the final errors on the retrieved quantities is made complicated by the significant correlation among atmospheric layers demonstrated by the forward RT model. This fact makes significant the off-diagonal elements of the error covariance matrix, and requires a suitable mixing of the errors to obtain a single value for each altitude. In other words, the errors are not uniquely defined by the Bayesian approach. After the detailed analyses on sensitivities, we are confident that the error budget is globally captured by the model, but anyway the uncertainties associated with individual retrieved points can vary as a consequence of specific optimization algorithms and weights assigned in the minimization of cost functions.

645

646

647

Moreover, we noted that the residuals between the measured spectra and the best fit appear log-normally distributed, hence we preferred to treat separately positive and negative deviations, resulting in final error bars not always symmetrical on a linear scale.

648

649

650

### 3.4.4. Errors

#### 3.4.5. Horizontal inhomogeneity and intrinsic limitations
A basic assumption of the adopted 1D version of the RT solver is that the spherical atmospheric layers are horizontally homogeneous in all their parameters. As already mentioned, this could represent in principle a serious limitation of the model, since the lines of sight in a limb geometry can encompass significant variations of dust optical depth to the same extent in both horizontal and vertical directions. A full 3D treatment of the problem could in principle shed light on the horizontal distribution of the aerosols by linking the variations of gaseous absorption bands in the limb spectra with changes in optical depth (analogously to the nadir-looking case, when a dust vertical profile can be inferred by mapping optically thick gas bands). Anyway, besides bursting computational times, such 3D treatment in the lack of measurements of the gaseous atmospheric structure and its horizontal variability would introduce a strong dependence of the retrievals from the GCM fields of pressure and temperature, needed to establish the gaseous opacity along the limb lines-of-sight.

In the 1D approximation all the horizontal variations are collapsed in a unique vertical profile. In this sense, the retrieved dust properties may not represent true vertical profiles in the atmosphere, since in principle the portions of the atmosphere mostly contributing to the signal measured at a given tangent altitude may actually pertain to different locations along the lines of sight, hence to different true altitude above the surface. An extreme example of such mixing of information is the case of a localized, very optically thick, aerosol clump (a cloud-like structure) which prevents probing the atmosphere along the line of sight behind it. A similar situation is also encountered at wavelengths where strong absorptions of atmospheric gases are encountered: in such a case, the penetration capability of the limb observation is largely reduced, and only the upper atmospheric layers can be probed even if the tangent altitude is much lower. This information problem results in difficulties to fit the CO₂ absorption bands in a 1D-limb model, which explains the exclusion of the 2 µm CO₂ band from the retrieval process.

As a consequence, even if the dependence of the dust retrievals on the GCM is largely reduced in the 1D case with respect to the 3D case, the 1D approximation can introduce a methodological issue when comparing the results with the GCM profiles extracted vertically at the altitudes given by the limb tangent points: while the former are retrieved as if the aerosol would be perfectly uniform horizontally distributed, the latter are picked from a variable 3D field. In our approach, we minimize such a bias by averaging the MCD values over the lines of sight of the limb observations (for both dust and gas components), and by considering the horizontal spatial fluctuations as model uncertainties.

### 3.5. Surface-shine

As mentioned above, an independent and comprehensive evaluation of the surface spectral albedo is needed to reduce the degrees of degeneracy of the limb RT problem. The radiative geometric relationship between a limb observation and the surface situation is sketched in Figure 8, where the reflection of illuminated surface terrains towards a single point along the limb line-of-sight is drawn. Obviously, similar contributions take place for every atmospheric element along the line-of-sight, since the measured limb spectrum is the
695
696

result of the integration along the full line-of-sight. The overall contribution, acting as an
extended secondary illumination source for limb spectra, can be called surface-shine.

<figure><img src="image_10.png" /></figure>

697
698
699

<caption>Figure 8: Sketch of the geometrical configuration of the surface-shine toward a single point along a limb viewing line-of-sight in a plane parallel configuration.</caption>

700
701
702
703
704
705
706
707
708

The variability of the surface albedo underneath the limb line-of-sight can be ascribed
to both the intrinsic spectral change of the involved terrains and the actual condition of
illumination of the surface. The first component can be investigated by using the OMEGA
dataset itself, after removing gaseous atmospheric features (to this purpose we apply the
method by Langevin et al. 2005). An example is provided in Figure 9, where a variability as
large as 50% of the surface albedo underneath the track of the limb line-of-sight makes very
hard to define a unique surface spectrum for an effective use in radiative transfer calculations.

<figure><img src="image_11.png" /></figure>

709
710
711
712
713
714
715
716
717
718
719

<caption>Figure 9: Example of albedo variability of the terrains underneath a limb line-of-sight. Left panel: The surface track of a given limb observation (actually pertaining to cube ORB0072_2) is drawn as a dashed line over a surface albedo map taken from another nadir-looking OMEGA cube (ORB0351_4 reflectance factor at 1.85 µm, corrected for atmosphere with the method of Langevin et al. 2005). The central cross represents the tangent point (the lowest point of the line-of-sight), whereas the labelled ticks at both sides indicate how the altitude of the line-of-sight changes while going away from the tangent point. Right panel: The three curves with symbols represent surface spectra from different regions close to the line-of-sight (numbered circles in the left hand-side map). The dashed curve with error bars represents the weighted average surface spectrum, described in the text as surface-shine.</caption>
720

721

722

723

724

725

726

727

728

729

730

731

732

733

734

735

736

737

The second component is in principle well known since it only depends on the 3D geometry of the observation at the time of the limb observation. Once the needed geometric parameters are known (altimetry and slope of a set of surface elements, along with their local solar incidence angles, emission and zenith angles relative to several limb points along any line-of-sight, which we calculated using MOLA data and NAIF SPICE kernels and libraries), this calculation yields to the distribution of the efficiency of illumination from the surface. The final efficiency is obtained by averaging the efficiencies for the individual path points, weighted by the inverse distance to the surface itself. An example of such results is shown in the left panel of Figure 10. Since the adopted RT scheme is unidimensional, we need to collapse all the spatial information on the surface albedo underneath the limb observation in a single averaged spectrum. To achieve this purpose, we select the best OMEGA cubes covering the surface region (e.g. in Figure 10, middle panel), then the final surface spectral albedo effective for a given limb observation is obtained by spatially integrating the OMEGA nadir atmosphere-corrected reflectance factors weighted by the surface-shine efficiency. An example of the result is also shown in the right panel of Figure 10.

The described procedure for the surface-shine spectral albedo $S_S(\lambda)$ can be described analytically as:

738

$$S_S(\lambda) = \frac{h_{top}}{A N} \sum_{i=1}^{N} C_i \iiint dxdydl \rho(x, y) \varepsilon_i(x, y, l) / h_i(l)$$

739

740

741

742

743

744

745

746

747

748

749

750

751

where the integrations are performed over the surface coordinates x, y and on the limb line-of-sight coordinate l; A is the involved surface area; N is the number of limb spectra; $h_i$ is the altitude of the line-of-sight for the i-th limb spectrum; $h_{top}$ is the maximum considered altitude of the limb scan; $\rho$ is the surface reflectance at normal incidence obtained from OMEGA surface cubes corrected for their own incidence angles; $\varepsilon_i = \mu_0^t v_e / (\mu_0^t + v_t)$ is a weight function dependent on the cosines of the incidence and emission angles with respect to the local topography, $\mu_0^t(x, y)$ and $v_t(x, y, l)$ respectively, and on the emission with respect to the planetary ellipsoid, $v_e(x, y, l)$, where the emission angles are referred to individual points of the i-th line-of-sight; and $C_i$ is a normalization constant ensuring that the overall weights numerically evaluated for every limb line-of-sight sum up to unity.

We point out that other formulations can work equally as well, depending on the adopted surface reflectance model and the availability of a suitable surface spectral reflectance field in the spectral range of interest.
<figure><img src="image_12.png" /><figcaption>Figure 10: Illustration of the surface averaging method for obtaining the surface-shine spectrum, in the case of ORB2008_0. Left panel: The contour lines represent an example of the illumination efficiency from the surface for a line-of-sight at 30 km tangent altitude (black dashed line). Middle panel: The filled polygons indicate the footprints of other OMEGA nadir-looking cubes covering the region underneath the line-of-sight and suitable for measuring the surface albedo. Right panel: Example of the application of the efficiency field to a single surface cube (ORB01209_6), with the same contour lines as in the left panel. The grayscale images represent MOLA altimetry data in left and middle panels, while the OMEGA atmospheric-corrected reflectance factor at 1.2 µm in the right panel.</figcaption></figure>

## 4. Results and discussion

The OMEGA limb spectral scans for the three study cases here investigated are shown in Figure 11 (black curves). For each scan, the individual 0.5-2.5 µm spectra are presented side-by-side in order to better illustrate the concept of the retrieval process, which attempts to model all the spectra simultaneously. The top altitude (40 km in the case ORB0072_2, 56 km in the case ORB2008_0, and 39 km in the case ORB2130_0) and the number of limb spectra included in each scan (8, 8, and 7 respectively for ORB0072_2, ORB2008_0, and ORB2130_0) are basically determined by the noise level of the measurements: all the limb spectra above those altitudes, although not completely noisy, include more and more spectral points with signal-to-noise ratio very close or below unity, yielding retrieval uncertainties larger than 100%.

The best fits obtained for these limb spectral scans are also shown in Figure 11, as red curves. The overall fitting quality is good, with reduced $\chi^2$ values of the order of unity for all three scans. We can see that the model correctly reproduces all the spectral shapes and intensities. The different concavity of the spectral shapes (the ORB0072_2 case exhibits enhanced visible peaks at almost all altitudes) can be ascribed to the change of dust scattering efficiencies with phase angle (which is rather variable among the observations, see Table 1). The poorest fit results are found at low altitudes in the ORB2008_0 case, where the synthetic
spectra are not able to fully capture the spectral slope. Minor misfits also appear in the long-wavelength side (2.2-2.5 µm) of some spectra (e.g. in the ORB2130_0 case at 16 and 22 km). The source of this kind of misfits can be primarily ascribed to horizontal inhomogeneities which likely affect the measured spectra and that cannot be taken into account in the model (see Section 3.4.5). The ORB2008_0 case, overlying the Hellas basin, is also that one where stronger horizontal variations are expected (from the MCD) to happen along the limb line-of-sight in both temperature, winds, and aerosol distribution. More complex size distributions and/or variations in grain shape may yield selective changes of the spectral shape. Anyway, letting all these quantities vary in the retrievals means exploring a much wider region of the parameters space, while the available data and radiative transfer assumptions cannot provide enough information for constraining all of them. It is worth stressing again that the limb spectra are highly correlated among each other (10% to more than 50% of the signal is carried by multiply scattered photons), and hence the RT solution is to be found simultaneously for all of them: a change of a synthetic spectrum at a given altitude always affects the spectra at the other altitudes (mainly the adjacent ones, but not only). Therefore, the fit quality must be judged for all limb scans considered as a whole. The CO₂ band at 2.0 µm is not shown in the synthetic spectra since it is not included in the fitting process and χ² calculation (see Section 3.4.5).

<figure><img src="image_13.png" /><figcaption>Figure 11: Best fits of the limb scans under investigation. The OMEGA 0.5-2.5 µm spectra composing the scan are shown side-by-side in black, as long as the associated 1-σ</figcaption></figure>
uncertainties indicated by the shaded grey regions. Boundaries of individual spectra are drawn as vertical lines, whereas the mean tangent altitude value for each spectrum is reported at the top. Synthetic best fitting spectra, computed from the parameters displayed in Figure 12, are provided as red lines.

The vertical profiles of the dust parameters corresponding to such best fitting models are shown in Figure 12 (thick lines and dash-filled areas). For sake of comparison, the corresponding MCD-based profiles already discussed (Section 3.1 and Figure 4) are reported here as well.

As far as number density is concerned, it is evident that the values retrieved at the base of the atmospheric column are larger than the GCM ones by about one order of magnitude, as an average for all three cases. It is worth remembering that the MCD profiles shown here depend on the assumed gas density profile, and therefore suffer from the uncertainties on GCM pressure and temperature fields (see Section 3.4.4). All the profiles exhibit a depletion at a different altitude level, decreasing by about 1 order of magnitude at about 18, 38, and 25 km respectively for the ORB0072_2, ORB2008_0, and ORB2130_0 cases. As mentioned above, limb spectra above the top altitudes are not suitable for inversion due to their overall noise level (uncertainties on retrieved parameters grow to more than 100%). This is especially true for the ORB2130_0 case, where the very low reflectance measured for the uppermost spectrum (<0.03) makes the retrieved number density value (~0.1 cm⁻³) very close to the OMEGA dust lower detection limit.

The disagreement between the MCD dust profile predictions and the observational fits reported in this work is not surprising. It corresponds to a well-identified problem in the Mars Global Climate Model prediction of the dust observations that was revealed by comparing GCM dust fields with Mars Climate Sounder Data (Heavens et al. 2014; Forget et al. 2014). During the seasons analyzed in this paper, it is likely that the enhanced dust concentration above 20 km may result from small scale processes not resolved in Global Climate Model such as the “rocket dust storms” that are thought to convectively lift dust to the middle atmosphere as a result of solar heating (Wang et al. 2018; Spiga et al. 2013). Improvements of the LMD GCM are currently developed by including parameterization of these processes (Wang et al. 2018), and future comparisons with our observations are expected to provide a better agreement.

As far as grain size is concerned, we can point out that retrieved effective radii are nearly constant in all the studied cases between 15 and 30 km of altitude, where a mean value of 0.85 ± 0.15 µm can be established. This value is also well in agreement with the model for the ORB0072_2 and ORB2008_0 cases, but is clearly larger than predicted for the ORB2130_0 case. The bottom part of the profiles (below ~10 km) shows somewhat larger discrepancies, but they are anyway in overall agreement with MCD values once uncertainties are taken into account. The most evident disagreement is found at the top of the ORB2008_0 case (the only one extending above 40 km), where the OMEGA retrieved value is nearly constant at 0.7 µm against a 0.3 µm MCD value.
<figure><img src="image_14.png" /><figcaption>Figure 12: Vertical profiles of dust parameters retrieved from the OMEGA limb scans, corresponding to the synthetic spectra in Figure 11. Left: Dust grains number density. Middle: Dust grains effective radius. Right: Dust mass mixing ratio. All: Thick curves and dash-filled areas show the OMEGA retrieved parameters and their associated uncertainties. Shade-filled areas indicate the corresponding quantities as predicted from MCD (same profiles as in Figure 2).</figcaption></figure>

Information on density and size can be combined to yield profiles of dust mass mixing ratio (i.e. the ratio of the total mass of the dust and of the gases in a unit volume of air), which are illustrated in the rightmost panel of Figure 12. Such conversion, if on one side will allow the comparison of our retrieval with other independent observations, on the other side needs further assumptions (namely about the bulk mass density of the mineral composing the dust and about the vertical profile of the gas number density, derived from the MCD pressure and temperature profiles), and is therefore subject to further uncertainties. If on one hand this quantity confirms the presence in the atmosphere of an overall mass of dust larger than predicted, on the other hand it can shed some more light on the physical interpretation of the results. E.g., the dust profiles in the ORB0072_2 and ORB2130_0 cases are rather different in terms of number density, but closer to each other in terms of mass. This is ascribable to the opposite trends of the grain size profiles and may indicate that the processes triggering the suspension of a similar amount of dust in the atmosphere were different in the two cases. In the ORB2008_0 case the behavior is clearly different, since the growing trends of both number density and radius yield a bulge in mixing ratio between 20 and 30 km of altitude.

Mass mixing ratio is also a useful quantity to be compared with measurements obtained by other instruments. A valuable reference dataset to check our results is that produced by MCS on board the Mars Reconnaissance Orbiter, provided as a climatological database of vertical profiles of dust abundance from Martian limb observations at thermal infrared wavelengths. Since MRO started data acquisition in 2006, there are no MCS profiles in the same Martian Year as the OMEGA observations here investigated (MY26 and 27). In Figure 13, our results are compared with MCS profiles retrieved at the same locations and Ls as OMEGA data but encompassing 4 Martian Years (between MY29 and MY33). MCS local
time is always near 15 h, not too far from the OMEGA values (see Table 1). This comparison shows that OMEGA and MCS results are in good agreement in the ORB2130_0 case, where both profiles lie well above the MCD predicted levels, and in relative agreement for the ORB0072_2 case as well. For the latter case, while the two measurements are coincident in the upper portion (30-40 km altitude), where MCS shows greater interannual variability, the OMEGA profile appears dust-depleted below ~25 km altitude (up to one order of magnitude). Anyway, the largest disagreement is found in the ORB2008_0 case: below 30 km, the shapes of the profiles are notably different but the uncertainties are large, hence the order of magnitude of the dust content can be assumed similar, and well above MCD levels in both cases. On the contrary, above 40 km, the MCS profile encompasses an abrupt drop witnessing the top of the dust bulk layer, while OMEGA retrievals show an amount of dust nearly constant up to the top of the retrieved profile, at ~55 km. Both profiles are dustier with respect to MCD predictions, but the OMEGA case is suggestive of enhanced dust activity. Since this case is located within Hellas basin, known as one of the largest dust reservoirs of the planet, a possible interpretation is that the ORB2008_0 observation captures a local dust storm rising from the basin in the early afternoon, and populating the complete atmospheric column with dust grains larger than usual.

A similar case is also reported in Oliva et al. (2018) where the presence of a local storm is made undoubtful by direct imaging. In that case, a nadir reflectance peaking at a value of about 0.35, close to the limb reflectance values here measured at the bottom of cases ORB2008_0 and ORB2130_0, yielded dust abundances within the bulk of the storm of the same order of the limb retrievals (about $10^2 \text{ cm}^{-3}$, see left panel in Figure 12 and Figure 8 of Oliva et al. 2018).

<figure><img src="image_15.png" /></figure>

900

901

902

<caption>Figure 13: Comparison of the OMEGA retrieved dust mass mixing ratio profiles (blue, green, and red solid lines with line-filled uncertainty areas, the same as in Figure 12)</caption>
with the analogous vertical profiles obtained by the MRO/MCS instrument (shaded polygons). The latter refers to the same Martian locations and Ls of the OMEGA ones, but to different Martian Years (MY30 to MY33 for the blue case, MY29 to MY33 excluding MY31 for the green case, and MY29 to MY32 for the red case). MCS local time is always 15h, not too far from the one of the OMEGA data.

As far as dust size is concerned, it is worth stressing that the effective radius is a useful parameter to simplify the treatment of radiative transfer but, being a cross-section averaged quantity, it must not be confused with the true radius of the dust grains. The physical distribution of dust sizes corresponding to the retrieved profiles can be obtained by fully taking into account the widths of the distributions. The comparison of $r_{\text{eff}}$ only is prone to misleading interpretation if $v_{\text{eff}}$ values are not properly taken into account. In the present case, since both our model and the MCD adopted the same kind of distribution (0.5 $v_{\text{eff}}$ lognormals), the comparison is simplified, and we can obtain the values of dust grains modal radius by multiplying the $r_{\text{eff}}$ values in Figure 12 by a constant factor of 0.36. Assuming a 0.1 uncertainty on $v_{\text{eff}}$, this yields an average value of the modal radius of $r_m=0.3 \pm 0.06$ µm between 15 and 30 km. This value is of course dependent on the adopted $v_{\text{eff}}$ value, which we selected from best fits of the spectral scans (see Section 3.4).

It is worth noting that the resulting widths of the 0.5-$v_{\text{eff}}$ distributions are quite large, encompassing significant amounts of nanosized particles as well as much coarser ones. The overall worsening of fitting performances when dealing with narrower distributions suggests that the measurements are somehow sensitive to this large range of sizes. In these situations, the effective parameters $r_{\text{eff}}$ and $v_{\text{eff}}$ do not always allow a direct interpretation of the physical processes driving the aerosol distribution, since the latter can be composed of aerosol populations of different origin. A slightly different view can be obtained by rearranging the number density and effective radius profiles in Figure 12 to give a series of vertical profiles of dust abundance with fixed particle sizes (or fixed size ranges). In other words, without adding nor cutting down any information with respect to Figure 12, we can calculate the retrieved aerosol distribution for every altitude and extract vertical sections for given classes of size, recovering a more physically comprehensive description of the dust content of the atmosphere. This is also useful for comparing the results with other kinds of measurement (e.g. *in situ*) or even with other OMEGA observations characterized by different $v_{\text{eff}}$ values.

In Figure 14, we can see the resulting density profiles for the three study cases, by using four classes of particle size: nanoparticles, with a radius $r < 0.1$ µm, fine particles with radii ranging $0.1 < r < 0.5$ µm, medium sized particles with radii ranging $0.5 < r < 1$ µm and coarser particles with radii ranging $1 < r < 10$ µm. We give in the following some considerations about the physical interpretation of the resulting profiles, based on the expected dynamical regime, with particular attention to the vertical winds field predicted by the MCD in the regions interested by the three observations. We do not attempt any quantitative analysis for supporting such interpretations, being far beyond the purpose of this paper, but they are just given to show the general physical feasibility of the retrieved profiles.

The result for ORB0072_2, showing a dust layer more and more flattened toward the surface for increasing grain size, is suggestive of size segregation, possibly related to a dust sedimentation process. Above this layer, the dust appears stably mixed (see the nearly
constant mass mixing ratio in Figure 13), with a small excess of Aitken class grains (nearly uniform number density, see left panel of Figure 14). This is also compatible with the weak subsidence that MCD predicts for the vertical wind field (<0.07 m/s downward velocities).

The ORB2008_0 is located just above Hellas Planitia and the dust distribution will be related to the local circulation pattern. The lowest part of this pattern, down to -7 km of altitude, is hardly accessible to OMEGA limb observations due to spatial resolution. An overall 0.1 m/s subsidence is expected within the crater for the whole investigated atmospheric column. However, no strong size segregation emerges from the retrieval, with the only possible exception of some depletion of large particles at its very bottom.

Also in the ORB2130_0 some evidence of size segregation appears: the grains growing trend with altitude suggested by the $r_{eff}$ profile in Figure 12 can be ascribed to slightly different scale heights in relation with the classes of particles which decrease with altitude. The greater complexity of the dust distribution in this case can reflect an enhanced complexity of the wind field, that MCD predicts to be characterized by alternate upwards (between 0-15 km and 30-60 km altitude) and downwards (between 15-30 km) vertical components. This wind pattern is also compatible with the dust depletion occurring in the mass mixing ratio profile at altitude above about 20 km (Figure 12, right panel), fully within the subsidence region, while a nearly homogenous dust mixing (below 18 km and above 27 km) is associated with the rising air at those altitudes.

<figure><img src="image_16.png" /><figcaption>Figure 14: Vertical profiles of retrieved dust number density for four size classes of grains. Nanometer particles (Aitken mode) are shown in the left panel. Sub-micrometer fine and medium ones are in the middle panels, while coarser micrometer particles are shown in the right panel.</figcaption></figure>

Although no meaningful statistics can be derived by the analysis of three cases only, the OMEGA retrievals seem to confirm a larger amount of dust in the Martian atmosphere than estimated by the GCM. Anyway, we may note that rather than sampling the most common Martian conditions, the criteria applied in the data selection may result in favoring
observations with enhanced dust content: in fact, since a large dust optical thickness increases the OMEGA detection limit for water ice, storm-like conditions can be statistically over-represented in the selected data sample.

## 5. Conclusion

In this work, we explore the information content of the OMEGA/MEx observations of the Martian limb, focusing on the vertical variations of dust. We approach the quantitative retrieval of dust properties from these data by applying a full spherical multiple scattering Monte Carlo 1D radiative transfer code (MYSTIC) as a forward model. The long-term purpose of this work is to derive a dataset of vertical profiles of dust density and grain size, which can contribute to our knowledge of the dust climatological cycles. Before extensive applications to the whole dataset, several methodological aspects must be well understood and optimized. Furthermore, being the first time this accurate modelling approach is applied to OMEGA limb data, it may represent a very useful benchmark for validating faster, and less accurate, radiative transfer models. For these reasons, in this paper we apply the retrieval method on three OMEGA study cases only, selected as not showing any spectral features of water ice, which is not yet included in the model. Spectral scans of the Martian limb with different vertical resolutions are extracted from VNIR and SWIR OMEGA channels (0.5 to 2.5 µm wavelength), hence ruling the thermal emission range out. We test the model assumptions needed to keep both computational time and accuracy of the results within feasible ranges. The sensitivity to the spectral range, to the single scattering phase function of the dust grains, to the *a priori* assumptions, as well as to the surface spectral albedo, are explored. For the latter, we implement a calculation based on a weighted-average of the surface reflectance below the limb tangent points, taken from other OMEGA nadir observations, in order to summarize with sufficient accuracy the overall surface-shine that acts as a secondary illumination source for the limb observations. The inclusion of the VNIR channel is found useful in reducing the degrees of degeneracy of the solutions, since the visible range covers the peak of the limb spectra reflectivity and provides slightly different constraints with respect to the near-infrared range.

In all the three study cases the retrieved dust number density appears enhanced with respect to the predictions of the LMD Global Climate Model (accessed through the MCD v5.2), by about one order of magnitude. An overall good agreement is found for the grain effective radius profiles, with an average value at mid altitudes (15-30 km) of $r_{eff} = 0.85 \pm 0.15$ µm, corresponding to a modal radius $r_m=0.3 \pm 0.05$ µm. Both smaller ($r_{eff} = 0.75$ µm; $r_m=0.27$ µm) and larger ($r_{eff} = 1.75$ µm; $r_m=0.63$ µm) values are found in the lowest layers (below 10 km), but associated with larger uncertainties ($\pm 0.25$ µm and $\pm 0.70$ µm respectively). When the OMEGA retrievals are compared to the dust mass mixing ratio profiles from MCS/MRO over the same Martian locations, we find an overall good agreement for two cases (ORB0072_2 and ORB2130_0). In the case of ORB2008_0 the OMEGA data overtakes MCS values as well, mainly at high altitude (40-55 km), suggesting the presence of a local storm rising from Hellas crater.

Obviously, the small sample of limb observations investigated here does not allow statistically meaningful conclusions on global behaviors. Moreover, the adopted selection
criteria probably favored the choice of data encompassing enhanced dust content, found in
local storms. Moreover, the fruitful comparison with MCS data demonstrates the possibility
of OMEGA data in contributing to quantitative studies of the Martian dust, although the
spectral range requires a more complex radiative transfer scheme with respect to the thermal
wavelength range. On the other hand, the enhanced sensitivity to multiple scattering makes
OMEGA limb data sensitive to very small amounts of aerosols. Therefore, a significant
return is expected from a wider application of the retrieval method to the whole OMEGA
limb dataset. In this sense, besides contributing to the general understanding of Martian dust
from a climatological point of view, OMEGA data may help assess the occurrence of local
dust storms and their relationship with dust circulation and global storm. However, these
studies require further work aimed to include water ice aerosols (and possibly thermal
emission) in the model, in order to unbias the retrieval of dust properties in every limb
observation. This approach is also expected to enable a global-scale study of suspended water
ice distribution and its relation with dust. We envisage that the computational burden required
by water ice inclusion is made more feasible not only by a trivial increase of computational
power, but also by some other optimization strategy, like e.g. the analytical evaluation of
Jacobian matrices.

In any case, being the first time a Monte Carlo modelling approach is applied to
OMEGA limb data, it represents a useful benchmark for validating faster, although less
accurate, radiative transfer models. Moreover, applications of the same methodology are
foreseen to the data currently acquired by the Exomars Trace Gas Orbiter NOMAD
spectrometer (Vandaele, 2018).

1041

1042 Acknowledgments

This study was performed within the UPWARDS project and co-funded in the context
of the European Union's Horizon 2020 Programme (H2020-Compet-08-2014), grant
agreement UPWARDS- 633127 and by the INAF PRIN 2014 project "Development of an
Atmospheric RT tool for the Exploitation of planetary limb data".

The authors wish to thank the Italian Space Agency for financial support through
contract ASI-INAF n. 2018-2-HH.0. FO and ED are also grateful to C.Emde for her support
in handling the MYSTIC solver, and to Sergio Fonte and Romolo Politi for preciously
supporting the management of the HPC machines within the INAF-IAPS DataWell facility.

We also acknowledge the OMEGA PI, J.P.Bibring (Institut d'Astrophysique Spatiale,  
France), the ESA Planetary Science Archive (PSA) and the whole MEx-OMEGA team for  
providing access to the OMEGA data and calibration pipeline. The MGS-MOLA data were  
obtained from the NASA Planetary Data System (PDS). We are also grateful to M.Vincendon  
for precious reviews and M.Lopez-Valverde for suggestions improving this work.

1056

1057

1058 References

1059 Acton, Charles, Nathaniel Bachman, Boris Semenov, and Edward Wright. 2018. "A Look towards the Future in the Handling of Space Science Mission Geometry." *Planetary and Space Science*.
1060 https://doi.org/10.1016/j.pss.2017.02.013.

1062 Altieri, F., L. Zasova, E. D'Aversa, G. Bellucci, F. G. Carrozzo, B. Gondet, and J. -P. Bibring. 2009.
1063 "O2 1.27µm Emission Maps as Derived from OMEGA/MEx Data." *Icarus*. https://doi.org/10.1016/j.icarus.2009.07.022.

1064 https://doi.org/10.1016/j.icarus.2009.07.022.

1065 Anderson, Ernest, and Conway Leovy. 1978. "Mariner 9 Television Limb Observations of Dust and

1066 Ice Hazes on Mars." *Journal of the Atmospheric Sciences*.

1067 https://doi.org/2.0.co;2">10.1175/1520-0469(1978)035<0723:mtlood>2.0.co;2.

1068 Bibring, Jean-Pierre, Yves Langevin, Aline Gendrin, Brigitte Gondet, François Poulet, Michel Berthé,

1069 Alain Soufflot, et al. 2005. "Mars Surface Diversity as Revealed by the OMEGA/Mars Express

1070 Observations." *Science* 307 (5715): 1576-81.

1071 Bibring, J-P, A. Soufflot, M. Berthé, Y. Langevin, B. Gondet, P. Drossart, M. Bouyé, et al. 2004.

1072 "OMEGA: Observatoire Pour La Minéralogie, l'Eau, Les Glaces et l'Activité." In *Mars Express*:

1073 *The Scientific Payload*, edited by Scientific Coordination: Agustin Chicarro Andrew Wilson, 37-

1074 49. ESA SP-1240, Noordwijk, Netherlands: ESA Publications Division.

1075 Buras, Robert, and Bernhard Mayer. 2011. "Efficient Unbiased Variance Reduction Techniques for

1076 Monte Carlo Simulations of Radiative Transfer in Cloudy Atmospheres: The Solution." *Journal*

1077 *of Quantitative Spectroscopy and Radiative Transfer*. https://doi.org/10.1016/j.jqsrt.2010.10.005.

1078 Carrozzo, F. G., F. Altieri, G. Bellucci, F. Poulet, E. D'Aversa, and J. -P. Bibring. 2012. "Iron

1079 Mineralogy of the Surface of Mars from the 1µm Band Spectral Properties." *Journal of*

1080 *Geophysical Research: Planets*. https://doi.org/10.1029/2012je004091.

1081 Chassefière, E., J. E. Blamont, V. A. Krasnopolsky, O. I. Korablev, S. K. Atreya, and R. A. West.

1082 1992. "Vertical Structure and Size Distributions of Martian Aerosols from Solar Occultation

1083 Measurements." *Icarus*. https://doi.org/10.1016/0019-1035(92)90056-d.

1084 Clancy, R. T., and S. W. Lee. 1991. "A New Look at Dust and Clouds in the Mars Atmosphere:

1085 Analysis of Emission-Phase-Function Sequences from Global Viking IRTM Observations."

1086 *Icarus* 93 (1): 135-58.

1087 Clancy, R. Todd, R. Todd Clancy, Michael J. Wolff, Michael D. Smith, Armin Kleinböhl, Bruce A.

1088 Cantor, Scott L. Murchie, et al. 2019. "The Distribution, Composition, and Particle Properties of

1089 Mars Mesospheric Aerosols: An Analysis of CRISM Visible/near-IR Limb Spectra with Context

1090 from near-Coincident MCS and MARCI Observations." *Icarus*.

1091 https://doi.org/10.1016/j.icarus.2019.03.025.

1092 Clancy, R. T., M. J. Wolff, and P. R. Christensen. 2003. "Mars Aerosol Studies with the MGS TES

1093 Emission Phase Function Observations: Optical Depths, Particle Sizes, and Ice Cloud Types

1094 versus Latitude and Solar Longitude." *Journal of Geophysical Research* 108 (E9).

1095 https://doi.org/10.1029/2003je002058.

1096 Eluszkiewicz, Janusz, Gennady Uymin, David Flittner, Karen Cady-Pereira, Eli Mlawer, John

1097 Henderson, Jean-Luc Moncet, Thomas Nehrkorn, and Michael Wolff. 2017. "A Fast Code for

1098 Channel Limb Radiances with Gas Absorption and Scattering in a Spherical Atmosphere."

1099 *Journal of Quantitative Spectroscopy and Radiative Transfer*.

1100 https://doi.org/10.1016/j.jqsrt.2017.02.010.

1101 Fedorova, A. A., F. Montmessin, A. V. Rodin, O. I. Korablev, A. Määttänen, L. Maltagliati, and J. -L.

1102 Bertaux. 2014. "Evidence for a Bimodal Size Distribution for the Suspended Aerosol Particles on

1103 Mars." *Icarus*. https://doi.org/10.1016/j.icarus.2013.12.015.

1104 Forget, F., E. Millour, A. Spiga, T. Navarro, J-B Madeleine, A. Pottier, L. Montabone, et al. 2014.

1105 "Simulating the Mars Climate with the LMD Mars Global Climate Model: Validation and

1106 Issues." In *The Fifth International Workshop on the Mars Atmosphere: Modelling and*

1107 *Observation, Held on January 13-16 2014, in Oxford, U.K.* Edited by F. Forget and M. Millour,

1108 edited by F Forget And, 1204.

1109 Forget, François, Frédéric Hourdin, Richard Fournier, Christophe Hourdin, Olivier Talagrand,

1110 Matthew Collins, Stephen R. Lewis, Peter L. Read, and Jean-Paul Huot. 1999. "Improved

1111 General Circulation Models of the Martian Atmosphere from the Surface to above 80 Km."

1112 *Journal of Geophysical Research: Planets*. https://doi.org/10.1029/1999je001025.

1113 Gondet, Brigitte, and Jean-Pierre Bibring. 2018. "More and Still Unexploited Atmospheric

1114 OMEGA/MEx Observations." In *European Planetary Science Congress 2018, Held 16-21*

1115 *September 2018 at TU Berlin, Berlin, Germany*, id.EPSC2018-2313.

1116 Guzewich, Scott D., Michael D. Smith, and Michael J. Wolff. 2014. "The Vertical Distribution of

1117 Martian Aerosol Particle Size." *Journal of Geophysical Research: Planets*.
1118 https://doi.org/10.1002/2014je004704.

1119 Heavens, Nicholas G., Morgan S. Johnson, Wedad A. Abdou, David M. Kass, Armin Kleinböhl,

1120 Daniel J. McCleese, James H. Shirley, and R. John Wilson. 2014. "Seasonal and Diurnal

1121 Variability of Detached Dust Layers in the Tropical Martian Atmosphere." *Journal of Geophysical Research: Planets*. https://doi.org/10.1002/2014je004619.

1122 Jaquin, Fred, Peter Gierasch, and Ralph Kahn. 1986. "The Vertical Structure of Limb Hazes in the Martian Atmosphere." *Icarus*. https://doi.org/10.1016/0019-1035(86)90050-3.

1123 Kahre, Melinda A., Jeffery L. Hollingsworth, Robert M. Haberle, and James R. Murphy. 2008.

1126 "Investigations of the Variability of Dust Particle Sizes in the Martian Atmosphere Using the

1127 NASA Ames General Circulation Model." *Icarus*. https://doi.org/10.1016/j.icarus.2008.01.023.

1128 Kleinböhl, Armin, John T. Schofield, David M. Kass, Wedad A. Abdou, Charles R. Backus, Bhaswar

1129 Sen, James H. Shirley, et al. 2009. "Mars Climate Sounder Limb Profile Retrieval of

1130 Atmospheric Temperature, Pressure, and Dust and Water Ice Opacity." *Journal of Geophysical Research: Planets*. https://doi.org/10.1029/2009je003358.

1131

1132 Langevin, Yves, François Poulet, Jean-Pierre Bibring, and Brigitte Gondet. 2005. "Sulfates in the

1133 North Polar Region of Mars Detected by OMEGA/Mars Express." *Science* 307 (5715): 1584–86.

1134 Lemmon, M. T., M. J. Wolff, M. D. Smith, R. T. Clancy, D. Banfield, G. A. Landis, A. Ghosh, et al.

1135 2004. "Atmospheric Imaging Results from the Mars Exploration Rovers: Spirit and

1136 Opportunity." *Science* 306 (5702): 1753–56.

1137 Madeleine, J-B, J.-B. Madeleine, F. Forget, E. Millour, L. Montabone, and M. J. Wolff. 2011.

1138 "Revisiting the Radiative Impact of Dust on Mars Using the LMD Global Climate Model."

1139 *Journal of Geophysical Research*. https://doi.org/10.1029/2011je003855.

1140 Maltagliati, L., F. Montmessin, A. Fedorova, O. Korablev, F. Forget, and J-L Bertaux. 2011.

1141 "Evidence of Water Vapor in Excess of Saturation in the Atmosphere of Mars." *Science* 333

1142 (6051): 1868–71.

1143 Mayer, B. 2009. "Radiative Transfer in the Cloudy Atmosphere." *The European Physical Journal*

1144 *Conferences* 1: 75–99.

1145 Melchiorri, R., T. Encrenaz, T. Fouchet, P. Drossart, E. Lellouch, B. Gondet, J.-P. Bibring, et al.

1146 2007. "Water Vapor Mapping on Mars Using OMEGA/Mars Express." *Planetary and Space*

1147 *Science*. https://doi.org/10.1016/j.pss.2006.05.040.

1148 Millour, E., F. Forget, A. Spiga, M. Vals, V. Zakharov, L. Montabone, and F. Lefèvre. 2018. "The

1149 Mars Climate Database (version 5.3)." In *Proceedings of the Mars Science Workshop “From*

1150 *Mars Express to ExoMars*, Held 27-28 February 2018 at ESAC, Madrid, Spain, 68.

1151 Montabone, L., F. Forget, E. Millour, R. J. Wilson, S. R. Lewis, B. Cantor, D. Kass, et al. 2015.

1152 "Eight-Year Climatology of Dust Optical Depth on Mars." *Icarus*.

1153 https://doi.org/10.1016/j.icarus.2014.12.034.

1154 Montmessin, Franck, Jean-Loup Bertaux, Eric Quémerais, Oleg Korablev, Pascal Rannou, François

1155 Forget, Séverine Perrier, et al. 2006. "Subvisible CO2 Ice Clouds Detected in the Mesosphere of

1156 Mars." *Icarus*. https://doi.org/10.1016/j.icarus.2006.03.015.

1157 Oliva, F., A. Adriani, M. L. Moriconi, G. L. Liberti, E. D'Aversa, and G. Filacchione. 2016. "Clouds

1158 and Hazes Vertical Structure of a Saturn's Giant Vortex from Cassini/VIMS-V Data Analysis."

1159 *Icarus*. https://doi.org/10.1016/j.icarus.2016.06.021.

1160 Oliva, F., E. D'Aversa, G. L. Liberti, G. Sindoni, F. Altieri, and E. Castelli. 2019. "A Database of

1161 Aerosols and Gases Coefficients for VIS-NIR Radiative Transfer in the Solar System Planetary

1162 Atmospheres." *Planetary and Space Science*. https://doi.org/10.1016/j.pss.2018.09.002.

1163 Oliva, F., A. Geminale, E. D'Aversa, F. Altieri, G. Bellucci, F. G. Carrozzo, G. Sindoni, and D.

1164 Grassi. 2018. "Properties of a Martian Local Dust Storm in Atlantis Chaos from OMEGA/MEX

1165 Data." *Icarus*. https://doi.org/10.1016/j.icarus.2017.07.034.

1166 Rannou, P., S. Perrier, J.-L. Bertaux, F. Montmessin, O. Korablev, and A. Rébérac. 2006. "Dust and

1167 Cloud Detection at the Mars Limb with UV Scattered Sunlight with SPICAM." *Journal of*

1168 *Geophysical Research*. https://doi.org/10.1029/2006je002693.

1169 Rodgers, Clive D. 2000. *Inverse Methods for Atmospheric Sounding: Theory and Practice*. Vol. 2.

1170 Series on Atmospheric, Oceanic and Planetary Physics. WORLD SCIENTIFIC.

1171 Rozanov, V. V., T. Dinter, A. V. Rozanov, A. Wolanin, A. Bracher, and J. P. Burrows. 2017.

1172 "Radiative Transfer Modeling through Terrestrial Atmosphere and Ocean Accounting for
1173 Inelastic Processes: Software Package SCIATRAN." *Journal of Quantitative Spectroscopy and Radiative Transfer*. https://doi.org/10.1016/j.jqsrt.2017.03.009.

1175 Sánchez-Lavega, A., H. Chen-Chen, I. Ordoñez-Etxeberria, R. Hueso, T. del Río-Gaztelurrutia, A.

1176 Garro, A. Cardesín-Moinelo, D. Titov, and S. Wood. 2018. "Limb Clouds and Dust on Mars

1177 from Images Obtained by the Visual Monitoring Camera (VMC) Onboard Mars Express."

1178 *Icarus*. https://doi.org/10.1016/j.icarus.2017.07.026.

1179 Sánchez-Lavega, A., A. García Muñoz, E. García-Melendo, S. Pérez-Hoyos, J. M. Gómez-Forrellad, C. Pellier, M. Delcroix, et al. 2015. "An Extremely High-Altitude Plume Seen at Mars' Morning Terminator." *Nature* 518 (7540): 525–28.

1182 Schofield, John, Armin Kleinboehl, David M. Kass, and Daniel McCleese. 2018. "The Mars Climate

1183 Sounder - Six Martian Years of Global Atmospheric Observations." In *42nd COSPAR Scientific*

1184 Assembly. Held 14-22 July 2018, in Pasadena, California, USA, B4.1–2 – 18.

1185 Sindoni, G., D. Grassi, A. Adriani, A. Mura, M. L. Moriconi, B. M. Dinelli, G. Filacchione, et al.

1186 2017. "Characterization of the White Ovals on Jupiter's Southern Hemisphere Using the First

1187 Data by the Juno/JIRAM Instrument." *Geophysical Research Letters*.

1188 https://doi.org/10.1002/2017gl072940.

1189 Smith, Michael D., Michael J. Wolff, R. Todd Clancy, Armin Kleinböhl, and Scott L. Murchie. 2013.

1190 "Vertical Distribution of Dust and Water Ice Aerosols from CRISM Limb-Geometry

1191 Observations." *Journal of Geophysical Research: Planets*. https://doi.org/10.1002/jgre.20047.

1192 Spiga, Aymeric, Julien Faure, Jean-Baptiste Madeleine, Anni Määttänen, and François Forget. 2013.

1193 "Rocket Dust Storms and Detached Dust Layers in the Martian Atmosphere." *Journal of*

1194 *Geophysical Research: Planets*. https://doi.org/10.1002/jgre.20046.

1195 Tomasko, M. G., L. R. Doose, M. Lemmon, P. H. Smith, and E. Wegryn. 1999. "Properties of Dust in

1196 the Martian Atmosphere from the Imager on Mars Pathfinder." *Journal of Geophysical*

1197 *Research: Planets*. https://doi.org/10.1029/1998je900016.

1198 Vandaele, A. C., the NOMAD Team, J-J Lopez-Moreno, M. R. Patel, G. Bellucci, F. Daerden, B.

1199 Ristic, et al. 2018. "NOMAD, an Integrated Suite of Three Spectrometers for the ExoMars Trace

1200 Gas Mission: Technical Description, Science Objectives and Expected Performance." *Space*

1201 Science Reviews. https://doi.org/10.1007/s11214-018-0517-2.

1202 Vincendon, M., J. Audouard, F. Altieri, and A. Ody. 2015. "Mars Express Measurements of Surface

1203 Albedo Changes over 2004–2010." *Icarus*. https://doi.org/10.1016/j.icarus.2014.10.029.

1204 Vincendon, M., B. Gondet, C. Pilorget, and Y. Langevin. 2011. "Spatial, Seasonal and Vertical

1205 Distribution of Water Ice Cloud Particle Size Inferred from OMEGA Limb Data." In *The Fourth*

1206 International Workshop on the Mars Atmosphere: Modelling and Observation, Held 8-11

1207 *February, 2011, in Paris, France.*, 207–8.

1208 Wang, Chao, François Forget, Tanguy Bertrand, Aymeric Spiga, Ehouarn Millour, and Thomas

1209 Navarro. 2018. "Parameterization of Rocket Dust Storms on Mars in the LMD Martian GCM:

1210 Modeling Details and Validation." *Journal of Geophysical Research: Planets*.

1211 https://doi.org/10.1002/2017je005255.

1212 Warren, Stephen G., and Richard E. Brandt. 2008. "Optical Constants of Ice from the Ultraviolet to

1213 the Microwave: A Revised Compilation." *Journal of Geophysical Research*.

1214 https://doi.org/10.1029/2007jd009744.

1215 Wiscombe, W. J. 1977. "The Delta–M Method: Rapid Yet Accurate Radiative Flux Calculations for

1216 Strongly Asymmetric Phase Functions." *Journal of the Atmospheric Sciences*.

1217 https://doi.org/2.0.co;2">10.1175/1520-0469(1977)034<1408:tdmrya>2.0.co;2.

1218 Wolff, M. J., and A. Pankine. 2016. "Limb Retrievals of TES solarband/IR Data (and MCS Solarband

1219 Data)." In *American Geophysical Union, Fall General Assembly 2016*, id.P21B – 2107.

1220 Wolff, M. J., M. D. Smith, R. T. Clancy, R. Arvidson, M. Kahre, F. Seelos, S. Murchie, and H.

1221 Savijärvi. 2009. "Wavelength Dependence of Dust Aerosol Single Scattering Albedo as

1222 Observed by the Compact Reconnaissance Imaging Spectrometer." *Journal of Geophysical*

1223 *Research*. https://doi.org/10.1029/2009je003350.

1224