# A spectral synergy method to retrieve martian water vapor column-abundance and vertical distribution applied to Mars Express SPICAM and PFS nadir measurements

Franck Montmessin, Stéphane Ferron

► To cite this version:

Franck Montmessin, Stéphane Ferron. A spectral synergy method to retrieve martian water vapor column-abundance and vertical distribution applied to Mars Express SPICAM and PFS nadir measurements. Icarus, 2019, 317, pp.549-569. 10.1016/j.icarus.2018.07.022. insu-01865343

HAL Id: insu-01865343

https://insu.hal.science/insu-01865343

Submitted on 28 Aug 2020

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A spectral synergy method to retrieve martian water vapor column-abundance and vertical distribution applied to Mars Express SPICAM and PFS nadir measurements

F. Montmessin<sup>a,*</sup>, S. Ferron<sup>b</sup>

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

## ABSTRACT

Up to now, all attempts to retrieve martian water vapor from nadir observations have focused their analysis on a single spectral domain and comparison between the results of various experiments showed the difficulty to reconcile the water vapor datasets together. Inspired by a methodology recently developed for the analysis of Earth observations, a spectral synergy approach has been tested on the water vapor extraction from the measurements of two instruments onboard Mars Express. These instruments cover near-infrared (NIR) and thermal infrared (TIR) domains within which water vapor possesses diagnostic absorption/emission signatures. Since the two instruments have operated concomitantly around Mars, co-located measurements could be selected and processed to create a dataset exploitable by synergy. The synergy relies on a Bayesian inference algorithm to find the best fitting values of water vapor and other parameters (e.g. atmospheric and surface temperature, dust opacity) simultaneously for both NIR and TIR intervals. Results demonstrate that synergy augments the information content about water vapor by up to >50% compared to traditional non-synergistic methods. Not only does the synergy provide an unbiased estimate of the total column abundance of H₂O but it also provides a correspondingly unbiased estimate of H₂O abundance in the first 5 km of the boundary layer, an atmospheric region long remained unexplored although being known for hosting key mechanisms controlling the fate of volatiles on Mars. Using synergy, the scientific return of nadir observations of current and future missions at Mars (such as the ExoMars Trace Gas Orbiter) can be fully optimized.

## 1. Introduction

On Earth, the global mapping of trace species has been accomplished thanks to measurements made by passive remote sensors in orbit onboard satellites (e.g. IASI, GOSAT, SCIAMACHY, etc.). In a nadir viewing geometry, these instruments are used to return a single information regarding the species: that is either the column-integrated abundances in the case of Near-Infrared measurements (NIR) or the middle atmosphere concentration in the case of Thermal Infrared measurements (TIR), leaving open the question of how these species are distributed throughout the vertical and in particular how they interact with the surface (sources and sinks). However, new methods for retrieving trace species in a nadir looking geometry have been introduced and have shown that it is possible to extract more information from these datasets than when each dataset is analyzed separately. This is the case for CO, as described in Worden et al. (2010) who confirmed the predictions made by Pan et al. (1995,1998). In this work, the demonstration is made that when CO is inferred from the combination of coincident TIR and NIR measurements, this opens an access to the concentration of CO in the near-surface layer, and thus brings knowledge about a key part of the atmosphere where this kind of gases are

emitted. Similar conclusions have been drawn for CO₂ (Christi and Stephens, 2004) as well as for CH₄ (Razavi et al., 2009) and O₃ (Costantino et al., 2017). The reason explaining why observations of the same species in distinct wavelength intervals provide constraints on the vertical distribution is that each spectral interval provides a distinct sensitivity along the vertical. While the NIR domain is sensitive to the full column, the TIR domain is mostly sensitive to layers located in the middle atmosphere.

On Mars, all the studies conducted to date aimed at exploring the space and time variability of atmospheric species have essentially relied on a single spectral domain retrieval approach (either UV, NIR or TIR). This has been done for a number of species like H₂O, CO, O₃, and aerosols (Fouchet et al., 2007; Perrier et al., 2006; Smith, 2002, 2009). Interestingly, the Mars Express (MEX) mission offers a unique opportunity to test the same NIR-TIR synergistic approach as the one used for Earth's observations. In particular, there is strong interest in testing this approach to study Mars' water vapor. Indeed, Mars possesses an active hydrological cycle characterized by significant seasonal activity of the atmospheric and surface reservoirs of water. The current properties of the martian atmosphere allow for the existence of water only in trace quantities (~100 ppmV) and only in the solid and gaseous phases.

<footnote><sup>a</sup> LATMOS/IPSL, UVSQ Université Paris-Saclay, UPMC University Paris 06, CNRS, Guyancourt, France</footnote>

<footnote><sup>b</sup> ACRI-ST, boulevard des Garennes, Guyancourt, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: franck.montmessin@latmos.ipsl.fr (F. Montmessin).</footnote>
Every year in summer, large amounts of water vapor sublime from the ice-covered surfaces of the polar regions. Water is then transported globally by winds and later condenses back at the poles to sublime again the following summer. This forms the classical view of the annual water cycle on Mars, which has exhibited a very high degree of repeatability since the first annual survey performed in the late 1970s. A detailed understanding of this cycle is a prerequisite to gain knowledge on present-day Mars climate, but also on its past configurations: a number of studies have predicted that the water cycle can explain the geomorphology carved during the Amazonian era, in particular all the glaciological features remnant of the changes in the water cycle subsequent to changes in Mars' orbital configuration (see e.g. Levrard et al., 2004; Forget et al., 2006; Madeleine et al., 2009). For these reasons, a relevant perception of Mars' water cycle relies on an acute description of the exchanges of water occurring between the lower troposphere and the surface / sub-surface, the latter being associated to adsorption-desorption of water molecules into/from the regolith pores.

The first spatial and temporal monitoring of water vapor abundance on Mars was obtained by Viking with the MAWD instrument (Mars Atmospheric Water Detector, an echelle-spectrometer which detected water vapor from the measurement of the solar flux reflected by the surface and the atmosphere around 1.38 µm, corresponding to a H₂O NIR vibrational-rotational absorption band). The MAWD climatology, covering more than one martian year (Jakosky and Farmer, 1982), made it possible to identify the northern polar ice cap as the main source of water vapor in the atmosphere of Mars when exposed to spring-summer insolation. For two decades, the MAWD data constituted the core knowledge of the water cycle, yet some intriguing features raised by these observations remained only tentatively explained, in particular the timing of release of water vapor in the atmosphere in the mid-latitudes (Jakosky, 1983; Haberle and Jakosky, 1990). In-depth analysis of this dataset has in fact highlighted the role played by the exchanges between the atmosphere and the surface. A climatology covering three martian years was then obtained by the TES (Thermal Emission Spectrometer) IR spectrometer on MGS (Mars Global Surveyor; Smith, 2004). TES measurement, based on the emission in the rotational band of water between 20 and 35 µm, has been considered to be more reliable than that of MAWD because of its lower sensitivity to suspended dust aerosols that were not accounted for in the first MAWD retrievals. However, the sensitivity of the measurement in the TIR domain is also dependent on the thermal contrast between the atmosphere and the surface, and thus decrease in cold regions. The MEX and the MRO (Mars Reconnaissance Orbiter) missions pursued the monitoring initiated by TES in the late 90s. MEX in particular has provided water vapor retrievals from several instruments and within a variety of spectral domains: the Planetary Fourier Spectrometer (PFS) is used for its TIR retrieval capability (Fouchet et al., 2007) while its NIR channel captures water vapor at 2.6 µm (Tschimmel et al., 2008), in the same spectral interval as the one used by the Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité (OMEGA) hyperspectral imager (Encrenaz et al., 2005). The band at 1.38 µm, previously sampled by MAWD, is covered by SPICAM (Spectroscopy for the Investigation of the Characteristics of the Atmosphere of Mars, Fedorova et al., 2006).

An important effort was conducted within the MEX community to perform a cross-comparison between all the MEX instruments capable of detecting water vapor. This cross-comparison revealed broadly dispersed retrievals. While the trends obtained by the MEX instruments were in agreement with MAWD and TES, they also indicated a drier atmosphere around the North Pole in summer. In this region, OMEGA and SPICAM gave comparable results while at mid-latitudes, OMEGA, in agreement with MAWD, gave 1.5 times more H₂O than SPICAM did. PFS TIR channel gave intermediate results between TES and SPICAM. These discrepancies were tentatively explained by the non-uniform mixing ratio of water vapor in the troposphere (Tschimmel et al., 2008). This non-uniform distribution is actually predicted by global climate models (GCM) as in Richardson et al. (2002), Montmessin et al. (2004)

and recently by Navarro et al. (2014) but this characteristic has remained poorly documented by measurements. Solar occultations by SPICAM on Mars Express have partially filled this observational gap (Fedorova et al., 2009; Maltagliati et al., 2011, 2013; Fedorova et al., 2018) but are however limited by their uneven time and space sampling as well as by the impossibility to make meaningful observations below a minimum of 5–15 km due to dust masking the solar disk.

In this study, we have used a combination of SPICAM and PFS nadir observations to elaborate a simultaneous retrieval of H₂O from two distant wavelength ranges (NIR and TIR). This is the first time this particular synergistic method is employed to explore the martian atmosphere. Eventually, this multi-spectral technique is intended to establish an extended multi-annual H₂O climatology that will contain both a continuous mapping of the column-integrated abundance of water vapor with concomitant information on its vertical distribution, thereby providing insights into the behavior of water vapor in the planetary boundary layer (PBL). This would address a part of the atmosphere that is otherwise inaccessible by non-synergistic methods performed from the martian orbit and that hosts key physical processes deciding of the fate of volatiles such as water.

The present study has been guided by the following ambitions: (i) provide an evaluation of the contribution of the NIR - TIR synergy to the exploration of water vapor within the context of the MEX mission, (ii) characterize, in terms of coverage, vertical resolution, and bias, the estimators of the water vapor column obtained with a synergistic retrieval, (iii) apply the synergistic technique to a selected set of MEX data representative of the range of variations encountered by Mars water vapor on an annual scale and at all latitudes.

The first part of the manuscript (Section 2) consists in a description of the MEX mission datasets that have been employed for this study. The selection and averaging processes used for the creation of a dataset compatible with a synergistic extraction of water vapor are described. In Section 3, the forward instrumental and radiative transfer models to be inserted in the data reduction algorithm are presented with a discussion on their known limitations. The Section 4 is dedicated to a description of the Bayesian and the maximum likelihood formalisms chosen for the study. This section also gives a description of the two approaches (parametric and non-parametric) designed to represent water vapor vertical distribution. The results of the synergistic method are detailed in Section 5. To evaluate the additional information content supplied by the spectral synergy, its results are compared to traditional methods analyzing a single spectral interval at a time. Finally, a discussion regarding the whole study is conducted in Section 6, while a conclusion opening on a variety of perspectives is provided in Section 7.

## 2. Mars express datasets

The dataset used in this study has been assembled out of spectra acquired during Martian Year (MY) 27 by the SPICAM and PFS spectrometers (described thereafter). It consists of an ensemble of calibrated spectra with their wavelength registration (also called Level 1, L1 products) that have been filtered out (see Section 2.3). Each L1 product corresponds to a series of collocated spectra that meet a number of selection criteria made to ensure (i) sufficient quality of every individual measurement, (ii) sufficient geographical and seasonal coverages so as to cover various dust or water ice cloud opacity configurations, (iii) a minimum error of radiative transfer modeling due to surface inhomogeneity (in terms of albedo, emissivity or topography) within the effective field of view covered by a L1 product. While OMEGA data could have been inserted in the synergistic study, this option has been discarded considering the results of preliminary performance tests that did not end up in a meaningful improvement of the synergy when adding OMEGA.
## 2.1. The SPICAM instrument

SPICAM is a dual UV-IR spectrometer (Bertaux et al., 2006). In the context of our study, only the IR part of SPICAM is considered since it is the one that covers the H₂O 1.38 µm absorption feature. A complete description of the SPICAM infrared channel (SPICAM-IR) can be found in Korablev et al. (2006). The instantaneous field of view (FOV) of the instrument is 1° (equivalent to a footprint of ~4 km at MEX orbit periapsis). This characteristic, in addition to a relatively good signal to noise ratio (SNR), has proven to be adequate for H₂O column abundance mapping (Fedorova et al., 2006 and Trokhimovskiy et al., 2015). In our analysis, we use the data from detector 1 between 1.34 and 1.43 µm (henceforth the NIR domain of our spectral synergy), which provides significantly higher performances compared to detector 2 (the SPICAM IR channel spatially separates the incoming flux into two polarization components as a consequence of the birefringent property of the TeO₂ crystal located in the active filter of the instrument). The radiometric noise is wavelength independent in this range (Korablev et al., 2006) and its absolute level is therefore fixed, as explained in Section 2.3.

We use a wavelength dependent parameterization of the instrument line shape function (ILSF) that was supplied to us by Anna Fedorova (IKI). Stray light has a huge effect on the retrieved water vapor column abundance. Its contribution has been evaluated based on a « grey » modeling of the ILSF far wings in a similar fashion as in Korablev et al. (2013). Details on the assessment procedure are given in S1. The final ILSF is shown in Fig. 1. The systematic error associated with the uncertainty of the stray light contribution is around ~10%, that is of the order of the statistical error.

## 2.2. The PFS instrument

A detailed description of the PFS (Planetary Fourier Spectrometer) Fourier transform spectrometer can be found in Formisano et al. (2005), only the most relevant characteristics are listed here. The full width at half maximum (FWHM) of the instantaneous FOV of the instrument is equal to 2.7°, corresponding to a 12 km-diameter footprint at periapsis. The unapodized spectral resolution is roughly 1.3 cm⁻¹. We use a sine cardinal parametrization of the instrument transfer function whose FWHM is adjusted from the data as explained in S2; together with a wavelength assignment parameter.

In this analysis, several overlapping spectral windows have been defined in the TIR domain. The TIR2 band extends from 20 to 35 µm and is used to retrieve the water vapor abundance. However, retrieving H₂O in the TIR domain requires a good and independent knowledge of the atmospheric temperature profile. To this end, we have considered the TIR1 band between 12 and 19 µm, which is characterized by the strong absorption of the CO₂ 15 µm vibrational transition. A third domain, TIR3, has been defined to constrain the surface temperature and

the dust model properties. It spans the regions between 8 and 10 µm µm and between 19 and 25 µm.

The TIR2 band has been used to map the water vapor column abundance in Fouchet et al. (2007). As explained in this paper, the signal-to-noise ratio (SNR) of a single spectrum is not high enough to infer a reliable value for the column abundance. Replicating what was done is Fouchet's study, we average 9 consecutive spectra upon which a retrieval of the water vapor is performed. Other parameters (i.e. temperature and dust optical depth) are retrieved this way as well whereas Fouchet et al. (2007) performed a temperature retrieval of individual spectrum. Furthermore, in order to reduce the influence of errors arising from inter and/or intra FOV observation conditions heterogeneities that are not represented in our modeling, we apply a set of selection rules namely: (i) the effective FOV corresponding to the 9 averaged spectra must span a distance less than 300 km on the martian geoid (that is 5° for a geoid radius of 3396 km), (ii) the airmass factor corresponding for each individual spectrum varies by less than 1%, (iii) the variation of topography within the 9 instantaneous FOV and between them is less than 1% of the mean scale height of the atmosphere (~0.1 km), (iv) the variation of the TES albedo (cf. Section 3.4) within the effective FOV is less than 3%, (v) at least one observation of SPICAM-IR lies within the effective FOV of PFS, (vi) no overlap occurs between effective FOV.

## 2.3. Synergy level 1

### 2.3.1. SYN spectra

The synergy method requires first to assemble a set of combined NIR and TIR spectra (SYN product) upon which the retrieval algorithm is applied. In order to ensure a common horizontal spatial resolution between NIR and TIR measurements, all SPICAM-IR spectra with a FOV inside the effective FOV of the SYN product (corresponding to 9 PFS FOV, cf. Section 2.2) are averaged together. An additional screening is performed that performs selection based on the following considerations: (i) the NIR airmass factor corresponding to an individual spectrum varies by less than 1%, (ii) the selected spectra are contiguous, (iii) the relative variation of the albedo in the NIR band stays within a 5σ dispersion of the radiometric noise (to do so, the albedo is estimated with a polynomial fit of the smoothed spectra).

The final dataset comprises 449 SYN spectra distributed among 133 orbits encompassing the entire MY27.

### 2.3.2. Noise consideration

The overall mean SNR for each channel is displayed in Fig. 2. In the TIR, the radiometric noise is obtained by smoothing the wavenumber dependent radiometric variance curve calculated from the 9 spectra. In the TIR1 band (temperature) the mean SNR lies typically between 70 and 140. In the TIR2 band (water vapor), SNR values vary between 90 and 220.

<figure><img src="image_3.png" /><figcaption>Fig. 1. The normalized instrument line shape function as it appears at the central wavenumber of the NIR band used throughout this study.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 2. SNR values for each spectral channel. NIR values obtained from individual spectrum are also shown.</figcaption></figure>

Inside a SYN spectrum, 10–20 individual SPICAM-IR spectra are typically averaged together. A variance analysis in the continuum range is used to determine the level of radiometric noise of the SPICAM-IR data (it is found to be roughly 60% of the noise value given in Korablev et al. 2006). For an individual spectrum, the SNR lies typically between 130 and 240. For an average spectrum, a covariance accounting for the FOV to FOV albedo variability is added to the radio-metric covariance. This accounts for the too permissive 5σ selection rule stated above and contributes to the high variability of the SNR. A comparison between water vapor column abundances inferred from individual SPICAM-IR spectra and the column abundances inferred from a SYN spectrum at its effective FOV is presented in S3. The observed fluctuations have a magnitude of the order of the estimated errors and a non-zero mean value.

It is to be noted that the selection method was deliberately restrictive to maximize the chances of drawing benefits from the synergy in the retrievals and to ensure synergy would not be evaluated against a biased or flawed dataset. Working with the best cases does not grant a universal application of the method but at least its potential can be established in a sound manner. In a subsequent part of the study, one might be interested to optimize the selection criteria so as to maximize the number of observations retained for the retrieval and thereby reach and extract the full “breadth” of the Mars Express dataset in the context of synergy. According to our preliminary evaluation, the Mars Express synergistic dataset potentially contains up to 200 000 co-located spectra matching the selection criteria, to be compared with the 449 selected and presented here.

## 2.4. OMEGA data

We did not use the OMEGA (Observatoire pour la Minéralogie, l'Eau, les Glaces et L'Activité) measurements at 2.56 µm because, as shown by a previous investigation based on synthetic simulated data, the information provided by the OMEGA and SPICAM spectrometers on the vertical distribution of water vapor are qualitatively redundant (in terms of degrees of freedom for the signal, DOFFS, and vertical coverage). Besides, the performances of OMEGA for retrieving water vapor are slightly degraded compared to SPICAM due to its lower spectral resolution. Indeed, at OMEGA spectral resolution, the 2.6 µm band of water vapor is more complex to analyze due to the nearby presence of a deep CO₂ absorption feature that is merged with that of H₂O vapor. As a consequence, no significant influence of the OMEGA data is expected to exist on the synergistic retrieval when SPICAM-IR data are used.

## 3. Forward model

A complete forward radiative transfer model was developed to reproduce the type of observations performed by the MEX instruments considered in this study. In that context, one has to account for a variety of processes involved in the radiative transfer configuration of each wavelength interval to be analyzed (the NIR, TIR1, TIR2 and TIR3 bands). This forward model is eventually inserted in the middle of the spectral inversion loop so as to perform the parameter adjustment, including those related to water vapor.

### 3.1. Molecular absorption

The HITRAN 2012 (Rothman et al., 2013) spectroscopic database is used as a baseline to compute the absorption coefficient by CO₂ and H₂O. CO₂ lines shape parameters come from precise measurements and/or from theoretical or semi-empirical predictions. When sources are available, they are in good agreement and the experimental information is supplied by HITRAN.

#### 3.1.1. Spectroscopic parameters

The collisional widths of H₂O lines have to be corrected to take into account the CO₂ atmosphere of Mars. In the NIR domain, multiplying the foreign collisional width parameters by a factor of 1.6 accounts for this. The sensitivity of the water vapor column retrieval to this parameter is investigated in S4. In the TIR domain, we rely on the calculations of Brown et al. (2007), which provide a modified value of the half width, its temperature dependence, and the pressure-shift parameters.

In the NIR spectral range, most of the molecular absorption is due to the H₂O vibrational band centered at 1.38 µm (the absorption by CO₂ is negligible). The interval between 1.4 and 1.42 µm is characterized by a low absorption of H₂O. This interval could help disentangle the aerosols and surface contribution to the continuum, which is an important issue since dust aerosols modulate the optical path length of photons and have therefore an impact on the retrieved water vapor content. However, the NIR domain does not allow for an unambiguous determination between dust optical depth and the water vapor column and a dust model has to be derived elsewhere.

#### 3.1.2. Line mixing

In the TIR1 and TIR2 ranges, the emission is due respectively to the 15 µm vibrational band of CO₂ and to the rotational band of H₂O
ranging from 20 and 40 µm. The line-mixing effect in the CO₂ 15 µm vibrational transition region leads to a slight decrease of the absorption coefficient in the far wings of the band (based on a calculation by Niro et al., 2004). Although it may have an impact on the derivation of the temperature profile, we decided to neglect it to limit the absorption coefficient calculation time. In the TIR3 band, the residual absorption by CO₂ between 9 and 11 µm (see Fig. 1 of Smith, 2004) is neglected and the H₂O contribution between 19 and 25 µm is removed using a low pass filter as explained in Section 4.1.

### 3.1.3. Isotope line intensities

In HITRAN, the line intensities are weighted by the Earth relative isotopic abundances, but no adaptation to the martian configuration exists. Partition function at local thermodynamic equilibrium is computed with the total internal partition sums (TIPS) program of Laraia et al. (2011). A parallelized version of a code (Wells, 1999) that implements the Humlicek algorithm (Humlíček, 1982) is used to calculate the Voigt line shape at all atmospheric levels. The thresholding technique described in Letchworth et al. (2007) is implemented in order to filter out lines with small intensities and save computation time. We have checked that the absorption coefficients obtained in this way were in good agreement with those obtained with the well-validated code LBLRTM (Clough et al., 1992,1995).

### 3.1.4. Absorption cross-section determination

We rely on the correlated-k approximation (Fu et al., 1992) to increase the speed of radiative transfer calculations. The correlated k-distributions (CKD) are defined by a Gauss-Legendre quadrature with 16 nodes on spectral intervals whose extents correspond to a fraction of the ILSF's FWHM (namely 0.9 and 0.4 cm⁻¹ respectively for the SPICAM and PFS). The CKD are tabulated as a function of the temperature (between 70 and 350 K by step of 10 K), pressure (5 points per decade between 10⁻² and 10 hPa) and (for NIR and TIR2 only) water vapor volume mixing ratio (58 points, regularly spaced on a log-scale between 10⁻⁷ and 4 × 10⁻²). The « kdistribution_light » package (Eymet, 2014) has been used to compute the CKD lookup tables (LUT).

## 3.2. Scattering by aerosols

In this study, scattering by water ice clouds has not been considered due to the small impact on the retrieved water vapor parameters as deduced from the tests performed with, a selected set of samples (see S6 for additional details). To describe the extinction by dust, we rely on an inherent optical properties (IOP) database provided to us by Wolff (Clancy et al., 2003; Wolff et al., 2003,2009). Starting from a set of experimental constraints on microphysical properties (complex refraction index, particle size) the absorption and angular scattering cross-section have been calculated, between 260 nm and 50 µm, using the Mie theory. In the database, the angular dependence of the scattering phase function is represented by its asymmetry parameter g. Seven models are available. They correspond to lognormal particle size distribution (dispersion: 0.3 µm) with an effective radius ranging from 0.1 and 3 µm for lognormal particle size distribution. The scattering is negligible in the TIR channels as a result of the low single scattering albedo of dust (SSA) in TIR2, TIR3 and TIR2 (this channel being dominated by the strong CO₂ v₂ band). This is illustrated in Section 5.1.

In the NIR domain, the IOP are rather independent of the wavelength. For small particle size (radius < 0.1 µm) the asymmetry parameter is close to 0.2 while for a coarser particle it is equal to roughly 0.7. Considering these values, the Henley-Greenstein function can be used within a good approximation to represent the scattering phase function. Our scattering model relies on the discrete ordinate method to solve the radiative transfer equation (RTE) (see Section 3.5). In this approach, the scattering phase function is described by Legendre polynomials whose truncation is chosen in accordance to the value of g. The Delta-M approach is used to reduce the dependence of the precision

of the calculations to the number of polynomial coefficients in case of a g value close to 1.

The vertical distribution of dust is described by the volume mixing ratio (vmr) profile predicted by the Mars climate database (MCD). Microphysical properties of dust are assumed independent of altitude. The dust absorption band in the TIR3 domain is used to define both the dust optical depth (DOD) at 9 µm and the best set of IOP (that is the best effective radius). We have assumed that the representation of dust in the MCD is sufficiently realistic to establish a correlation between the DOD predicted by the MCD and the DOD retrieved from the PFS data (see Sections 4.1 and 5.1 for details).

## 3.3. The solar source function

Two TOA solar irradiance model have been used for the modeling of the NIR spectra: (i) a high spectral resolution (with a resolving power of 500 000) synthetic spectrum calculated in the wavelength range from 150 nm to 0.3 cm with the ATLAS12 model by Kurucz (http://kurucz.harvard.edu/stars/sun), (ii) an experimental calibrated spectrum ranging from 1.2 and 5 µm with a spectral resolution of ~1 nm prepared for the exploitation of the PFS data by Fiorenza et al. (2005). This spectrum is derived from several sources: between 2 and 5 µm, the space based measurements from the ATMOS experiment were used while in the range 1.2–2 µm, which includes the NIR domain, ground based measurements from the Kitt Peak observatory have been used. Both parts of the spectrum have been cross-calibrated using the solar continuum calculated by Kurucz. These spectra have been averaged over the spectral intervals used to define the CKD and have been normalized according to the Sun-Mars distance.

In our study, the « PFS » spectrum is taken as the reference model. It gives slightly better fit and the impact on the retrieved water vapor column abundances is negligible compared to other sources of uncertainty (see S5 for details).

## 3.4. Surface properties

Surface properties are considered both at data selection (through FOV uniformity conditions, see Section 1) and during retrieval stages. The surface visible albedo and thermal emissivity have been established from TES observations (Bandfield et al., 2003; Christensen et al., 2001). Owing to its multi angular capabilities, spectral signatures associated to scattering in the atmosphere and to the surface reflection/emissions could be disentangled. Fig. 3 shows the TES surface emissivity for bright and dark surfaces as explained in Bandfield et al. (2003). A smoothing procedure has been applied in order to remove noise and fill the gap corresponding to the 15 µm CO₂ band that is totally absorbed. This defines the first guess emissivity used in the subsequent TIR modeling assuming isotropic emission from the surface.

In the NIR, the value of the continuum depends on the way the atmospheric scattering is handled. Our first guess albedo is defined through a first-degree polynomial fit of the low-pass filtered spectra. The forward model assumes a lambertian behavior of the reflection distribution function of the surface. In both cases (NIR and TIR) the non-parametric Bayesian approach explained in Section 4 allows for a flexible low frequency shape modulation of the continuum assuming prior auto-correlation functions.

## 3.5. Calculation of upwelling TOA radiances

The TOA radiance for a plane paralleled non-diffusing atmosphere is:
<figure><img src="image_5.png" /><figcaption>Fig. 3. Left: dark/bright region mask built from the TES albedo map. High albedo (>0.18) regions appear in white, low albedo in black. Right: the symbols display the spectral emissivity for dark and bright regions obtained from TES data. The lines show the interpolations used in the forward model.</figcaption></figure>

$$I(\lambda) = \frac{A(\lambda)}{\pi} \mu_S F(\lambda) \exp\left(-\frac{\tau(0, \lambda)}{\mu}\right) + \varepsilon(\lambda) B(T_S, \lambda) \exp\left(-\frac{\tau(0, \lambda)}{\mu_O}\right) + \int_0^{Z_{TOA}} B(T(z), \lambda) k(z, \lambda) \exp\left(-\frac{\tau(z, \lambda)}{\mu_O}\right) dz \quad (1)$$

with $\lambda$ the considered wavelength. A, $\varepsilon$, $T_S$ are the surface albedo, emissivity and temperature respectively, B the black-body emission, F the solar radiance. The integrated opacity of the layer located between altitude $z$ and TOA is given by:

$$\tau(z, \lambda) = \int_{z}^{Z_{TOA}} k(u, \lambda) du \quad (2)$$

$k(z, \sigma)$ is the total extinction coefficient corresponding to absorption by molecules and aerosols and scattering by aerosols. The molecular contributions are obtained by a tri-linear interpolation of the temperature, log-pressure and water vapor vmr space of the CKD nodes (see Section 3.1). The intra-layer vertical integration of Eq. (2) is done assuming a linear variation of the absorption coefficient with the altitude. The air density profile is computed for hydrostatic equilibrium from the high spatially resolved MOLA (Mars Orbiter Laser Altimeter) surface pressure map (accessed through the MCD) and the temperature vertical profile (taken from the MCD or estimated from the PFS-LW data).

To take into account the scattering, the single scattering albedo and scattering phase function of the atmospheric layers have to be defined as well. This is done using the Wolff database presented in Section 3.2 assuming a Legendre polynomial development of the Henley-Greenstein phase function. We use 30 terms with a number of half-space streams of 8.

The calculations of intensity taking into account the scattering rely on the LIDORT (Linearized pseudo-spherical scalar Discrete Ordinate Radiative Transfer model) RTE solver (Spurr, 2008). This model is based on a linearized formulation of the DISORT discrete ordinate approach (Stamnes et al., 1988) allowing for fast analytical determination of the radiance derivatives with the atmospheric layer's IOP and surface properties (i.e. lambertian albedo and emissivity). The derivatives of the layers IOP with temperature (respectively water vapor vmr) are obtained by symmetric finite difference with a 1% perturbation (respectively 50%). The derivatives of the radiance with respect to the other parameters are calculated by symmetric finite difference as well. The perturbation depends on the parameters: 1% for the stray light contribution to the SPICAM-IR's ILSF, 5% for the wavelength assignment shift of the PSF-LW spectra, 110% for the ILSF's FWHM of this instrument, 5% for the H₂O vmr at the surface in the parametric formulation discussed below and 100 m for the saturation level altitude.

## 4. Inversion technique

In this study, we have to consider the retrieval of a variety of quantities, some regarding instrument properties (stray light contribution in SPICAM-IR, the FWHM of the ILSF and the spectral assignment for PFS) while other relate to geophysical parameters (surface albedo/ emissivity, atmospheric temperature, water vapor content and dust properties). To achieve this, three steps are followed in the estimation process: first, the instrumental parameters are determined, second the parameters of the dust model, and finally the parameters related to water vapor, atmospheric temperature and continuum. The first stage has been discussed in Section 1. The two others are discussed below.

### 4.1. Retrieval of surface temperature and dust properties

Scattering by aerosols may have a significant impact on the H₂O retrieval from NIR data (Trokhimovskiy et al., 2015 and in S6). However, the distinction between spectral signatures associated to the scattering and to the surface reflexion is difficult to achieve and may result in large uncertainties on the DOD and column abundance. For this reason, the aerosol model has to be defined prior to the H₂O retrieval. SYN spectra in the TIR3 channel are used to retrieve i) the surface temperature, ii) the DOD at 9 µm assuming a wavelength-dependent IOP model taken from (Wolff et al., see Section 3.2) and iii) the “best” dust effective radius which controls the spectral dependence of aerosols.

Between 20 and 25 µm, spectral samples contaminated by water vapor absorption are rejected: a 1-σ threshold is applied to the difference between the considered spectrum and the same spectrum with a low-pass filter. In that spectral region, the dust SSA is lower than 0.2 and the extinction by water ice cloud is negligible. The second part of the TIR bands used in this analysis extends from 8.5 to 10.5 µm. It is characterized by a strong silicate absorption band of dust (SSA < 0.5). The region between 10.5 and 12.5 µm is excluded because it might be contaminated by water ice cloud extinction (see Fig. 1 of Smith, 2004). Finally, the small absorption by CO₂ between 9 and 9.5 µm is neglected.

The retrieval relies on a non-linear least square fitting Levenberg-Marquardt algorithm (as implemented in the IDL package “MPFIT”, Markwardt, 2009) using the atmospheric temperature predicted by the MCD and the TES model emissivity (Section 3.4). A lognormal variable substitution is used to force the retrieved DOD positivity. A low (0.01) value of the DOD is assumed as first guess. The surface temperature is defined with a first guess corresponding to the average brightness temperature between 12.1 and 12.5 µm.

The information content of the TIR domain is not sufficient to
constrain the microphysical properties of the aerosols. So, we use the MCD to select the best model within the Wolff database. Note that our fast-radiative transfer model neglects scattering. For each radius, a full retrieval is performed that includes dust opacity as a free parameter. The retained radius is the one providing the opacity value closest to that proposed by the MCD. MCD is just used as a criterion for selection, opacity being otherwise free to vary. Results are presented in Section 5.1.

## 4.2. Retrieval of temperature and water vapor content

The retrieval of the temperature and water vapor relies on the Bayesian formalism described in Rodgers (2000). We assume that data is linked to the forward model through the following relation:

$$d = g(\bar{m}) + \epsilon \quad (3)$$

where $\epsilon$ is an unbiased gaussian random variable representing the noise on the measurement vector $d$ and $\bar{m}$ correspond to the true value of the state vectors. The estimation $\tilde{m}$ is defined by the value of $m$ which minimizes the quadratic cost function:

$$C(m) = [d - g(m)]^T C_d^{-1} [d - g(m)] + [m - m_0]^T C_0^{-1} [m - m_0] \quad (4)$$

The first term corresponds to the maximum likelihood probability distribution. The second term defines the uncertainty associated to the prior $m_0$. Since $g$ is a non-linear function, the estimator can be computed via an iterative Raphson–Newton descent method:

$$\tilde{m}_{i+1} = m_0 + C_0 G_i^T (G_i C_0 G_i^T + C_d)^{-1} [d - g(\tilde{m}_i) + G_i(\tilde{m}_i - m_0)] \quad (5)$$

with $G_i$ the Jacobian matrix of $g(\tilde{m}_i)$ evaluated at iteration $i$. In this expression, the inversion is done (by Cholesky decomposition) in the data space, which has in our case, a smaller size than the parameter space (due to continuum parameters). The convergence criterion is based on successive variations of the normalized $\chi^2$:

$$\chi_i^2 = \frac{[d - g(\tilde{m}_i)]^T C_d^{-1} [d - g(\tilde{m}_i)]}{n_d} \quad (6)$$

where $n_d$ is the number of data. The covariance error on $\tilde{m}$ is given by:

$$C_m = (1 - \text{MRM}) C_0 \quad (7)$$

with MRM the averaging kernel (matrix resolution model):

$$\text{MRM} = C_0 G^T (G C_0 G^T + C_d)^{-1} G \quad (8)$$

computed at convergence. The trace of MRM defines the DOFFS which corresponds to the number of independent parameters which can be constrained by data considering an a priori knowledge ($m_0, C_0$).

In Eq. (4) the matrix inversion is regularized by the $C_0$ matrix. For continuum (that is surface albedo and/or emissivity) parameters, which depend on the wavelength, $C_0$ is defined by an exponential kernel with a sufficiently high correlation length so that only low frequency modulation of the prior is retained by the estimation. Gaussian kernels are used for atmospheric profile parameters (that is water vapor vmr fluctuations relative to MCD vmr and/or temperature).

Prior experiments performed on synthetic datasets indicate that the DOFFS of H$_2$O are roughly equal to 1 when retrievals are performed on separate spectral domains (slightly higher for NIR, slightly lower for TIR depending on geophysical conditions). This constitutes a theoretical confirmation that one spectral domain alone is not able to supply more than one independent information regarding water vapor. This information is usually assigned to the column-integrated abundance.

With synergy, one is capable to constrain an additional information for H$_2$O and therefore establish a two-parameter model of the water vapor vertical distribution. In doing so, a large prior uncertainty on the column abundance can be assumed, resulting in a column estimator close to the unbiased maximum likelihood estimator (in other words, an estimation not biased by the prior information). The two-parameter model of the water vapor vmr is represented by:

$$x(z; x_0, z_s) = Y(z - z_s)r(z) + [1 - Y(z - z_s)] \left[ (r(z_s) - x_0) \frac{z}{z_s} + x_0 \right] \quad (9)$$

which depends on the vmr at the surface, $x_0$, and the saturation level altitude $z_s$. In this expression $Y$ is the Heavyside function:

$$Y(u) = \begin{cases} 1, & u \ge 0 \\ 0, & u < 0 \end{cases} \quad (10)$$

and $r$ is the saturation ratio. The saturation vapor pressure depends only on temperature. It is calculated using the Goff-Gratch formula.

The two-parameter model assumes 100% saturation above $z_s$ and a linear variation with altitude down to the surface. Since it cannot be sampled on a discrete vertical grid, we substitute the Heaviside function by the error function:

$$x(z; x_0, z_s) = \text{erf}(z - z_s)r(z) + [1 - \text{erf}(z - z_s)] \left[ (r(z_s) - x_0) \frac{z}{z_s} + x_0 \right] \quad (11)$$

with

$$\text{erf}(u) = \frac{1}{2} + \frac{1}{\sqrt{2\pi}\Delta} \int_{0}^{u} \exp\left(-\frac{t^2}{2\Delta^2}\right) dt \quad (12)$$

$\Delta$ is chosen to be equal to twice the vertical sampling step (2 km). The first guess is initialized using the values of $z_s$ and of the column abundance $W$ predicted by the MCD. The former is obtained with the vmr profile of the MCD and

$$x_0 = W \frac{gM}{m} \frac{\int_0^{PS} [\text{erf}(z - z_s)r(z) + [1 - \text{erf}(z - z_s)]r(z_s) \frac{z}{z_s}] dP}{\int_0^{PS} [1 - \text{erf}(z - z_s)] (1 - \frac{z}{z_s}) dP} \quad (13)$$

Fig. 4 shows a comparison between an original MCD profile and the parametric model in a case where $W \sim 8$ pr.µm. In this approach, uncertainty on the column abundances is calculated from the covariance associated to $x_0$ and $z_s$ estimators.

A second, non-parametric, approach, has been investigated which uses the MCD predictions as prior for our Bayesian algorithm. In that case, the goal is to identify large vertical scale fluctuations compatible with the NIR + TIR DOFFS, which might be different from single-domain DOFFS due to the complementarity of the averaging kernels between NIR and TIR in most configurations (see Section 5.3.3). We infer two estimators for H$_2$O: (i) the column abundance from the surface up to the TOA, and (ii) the column abundance from the surface up to 5 km,

<figure><img src="image_6.png" /><figcaption>Fig. 4. An example of H$_2$O vmr profile from the MCD (red) compared with the corresponding parametric model (black). Both profiles yield the same column abundance.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** A typical fit of PFS data in the TIR3 band where surface temperature and dust optical depth are retrieved, blue: without scattering, red: full transfer with scattering. The values of the parameters are a DOD of 0.18 at 9 µm, an effective dust radius of 1.5 µm and a surface temperature of 253.2 K. The bottom plots zoom over the two edges of the TIR3 domain and display the corresponding histogram of the residuals given in %.</figcaption></figure>

both are obtained by vertical integration of the retrieved concentration.

# 5. Results

## 5.1. Step1: Retrieval of surface temperature and dust properties

The methodology chosen to retrieve the surface temperature and dust properties from the TIR3 SYN spectra has been applied to the MY27 SYN products as previously explained. An example of forward modeling is displayed in Fig. 5. Distributions of residuals are displayed separately for the two parts of the TIR3 band located on both sides of the CO₂ absorption band. In the shortest wavenumber region, the adjusted spectra have been smoothed in order to remove the main absorption features while preserving the continuum shape. As illustrated in this figure, data are equally fitted regardless of whether scattering has been accounted for or not. The differences between the corresponding modeling results appear small when compared to the radio-metric noise.

An additional set of empirical criteria was established to retain only the best SYN products that will receive the complete retrieval procedure. The criteria are: (i) the Levenberg-Marquardt has to converge to a solution comprising a plausible DOD value and a normalized $\chi^2$ lower than 3 (see Fig. 6-left), (ii) the retrieved effective dust radius is smaller than 3 µm. The aim of the second condition is to prevent side effects resulting from the dust model selection procedure, which would

translate into a suspicious increase at the upper edge of the effective dust radius distribution (Fig. 6-middle). We see that results behave approximately like a log-normal distribution.

Finally, the retained dataset comprises 232 SYN spectra. Fig. 6 shows the correlation between the DOD obtained while taking into account or not the scattering. The exact relation between the two depends on the dust properties (i.e. the dust radius) and the spectral domain used for the retrieval which controls the average value of the SSA. However, to some good approximation a linear law might be considered. We note that the DOD with scattering is on average 15%–20% higher than the DOD without scattering. In the sections that follow, the DOD without scattering will be used for the TIR channels and the DOD with scattering for the NIR channel.

The estimated DOD values as well as their geographical distribution are shown on Fig. 7. Zonal averaged values are presented on Fig. 8. When compared with the THEMIS (THermal Emission Imaging System) climatology (Smith et al., 2009, their Fig. 6), we find a slightly systematic excess of DOD values however with the correct trends. It is to be noted that our sparse longitude sampling does not allow for an unbiased determination of zonal averages.

The correlations between our results with the predictions of the MCD are shown on Fig. 9. The good correlation between the DOD is due to the retrieval approach itself since the best fit is always leaning towards the MCD predictions (see Section 4.1 on the IOP model selection procedure).
<figure><img src="image_8.png" /><figcaption>Fig. 6. (From left to right) reduced $\chi^2$, number of occurrences of the dust effective radius, correlation between the DOD obtained with scattering (y-axis) and the DOD obtained without scattering (x-axis). Color is coded according to the effective radius value. DOD values are normalized to a reference pressure of 610 Pa.</figcaption></figure>

The DOD values obtained without scattering are found to be higher by <10% than those of the MCD, a result consistent with the outcome of our comparison with the THEMIS climatology. The cause for this small discrepancy was not searched for and it was decided instead to assign a 10% systematic uncertainty to the retrieved DOD.

## 5.2. Step 2: Retrieval of atmospheric temperature profiles

The retrieval of the atmospheric temperature from the TIR1 channel is based on the Bayesian formalism described in Section 4.2. The aerosol model estimated in the previous section is used and the DOD is kept constant. In that step, the retrieved parameters are the surface temperature, whose first guess value and its associated uncertainty are extracted from the dust retrieval (Section 5.1). The spectral emissivity of the surface is also inferred with a first guess built from the TES emissivity model (see Section 3.4). A wavelength independent prior uncertainty of 1% is applied with an exponential correlation kernel (correlation “length”: $10^3$ cm$^{-1}$). The emissivity parameter space is bounded between 0 and 1 by substituting $\epsilon$ by $p \equiv \sin(2\epsilon - 1)$. A wavelength shift parameter is also adjusted, taking advantage of the highly structured spectral shape of the CO$_2$ 15 µm band. The derivation of the atmospheric temperature profile starts with a first guess given by the MCD. The prior covariance is defined arbitrarily, being equivalent to constant uncertainty of 30 K throughout the profile, whereas a Gaussian correlation kernel with a standard deviation of 10 km is assumed. The

FWHM of the ILSF is kept constant and is determined as explained in S1. A typical data modeling of the TIR 1 spectral interval is shown on Fig. 10.

Although a good $\chi^2$ is achieved in this case, the fit is imperfect as revealed by the residual distribution displayed on Fig. 11 where some absorption features appear in the wings of the strong CO$_2$ line at 16.3 µm (615 cm$^{-1}$). The Fig. 12 displays the retrieval of the continuum parameters. The large prior correlation “length” ensures that no steep absorption feature is assigned to emissivity fluctuation. However, the “length” is chosen small enough to tolerate moderate departure from the prior assumption. The retrieval of the atmospheric temperature is shown in Fig. 13. The temperature DOFFS is equal to 6.4 whereas the DOFFS for the entire retrieval is 10.8. This means that 6 to 7 independent parameters can be extracted from the temperature profile, while the other 4–5 parameters pertain to surface emissivity / temperature as well as to the wavelength shift.

## 5.3. Step 3: Retrieval of water vapor

In order to gauge the potential benefit of the synergistic retrieval compared to traditional single-domain inversion techniques, the same SYN dataset was used to perform retrieval in three different manners: first, using the NIR domain only, then solely the TIR domain and finally the combination of both which establishes the synergy between the two domains.

<figure><img src="image_9.png" /><figcaption>Fig. 7. Map of DOD scaled at 610 Pa and derived without accounting for scattering from the PFS data.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Zonal averages of normalized DOD. The error bars correspond to the standard deviation within the latitude/solar longitude bin.</figcaption></figure>

### 5.3.1. NIR-only inversion

In this subsection, we analyze the results of the retrieval performed while assuming the H₂O parametric vmr model of Eq. (11) and while restricting ourselves to the NIR part of the SYN spectra. Still, we employ the aerosol model derived from the TIR3 band (see Section 4.1) while the DOD is kept constant. In order to limit computing time, the fast-radiative transfer model presented in Eq. (1) is used. Scattering was neglected, an assumption supported by the study presented in S6.

In the “NIR only” case, the retrieved parameters are the spectral albedo of the surface whose first guess is built as explained in 3.4. A wavelength independent prior uncertainty of 10% is applied with an exponential correlation kernel with a “length” of 10⁷ cm⁻¹. The positivity of these parameters is ensured by a lognormal variable substitution. The other parameters are the H₂O vmr at the surface, x₀ and the saturation level altitude zₛ. The first guess on these parameters are defined based on the MCD predictions. The positivity of x₀ is again ensured by a lognormal variable substitution. A 100% prior uncertainty on x₀ and 200 m prior uncertainty on zₛ were chosen. The ILSF model presented in Section 2.1 has been used. A typical data modeling is shown on Fig. 14.

An examination of the residual distributions (Fig. 15) whose fluctuations show distinct wavelength dependent structures suggests that a better fit might possibly be achieved by adding a spectral shift in the state vector (as done for the PFS spectra).

The retrieval of the continuum parameters for the same SYN product is shown on Fig. 16. Since atmospheric scattering is neglected, the continuum incorporates both the surface and the atmospheric contributions. In this particular example, atmospheric scattering accounts for 15% of the continuum as deduced from an estimation performed with scattering. The large correlation length of the prior covariance kernel ensures that a good separation with the H₂O absorption spectral features is achieved while some continuum shape distortion is allowed.

The water vapor column abundances retrieved for the entire dataset are shown in Fig. 17. The Fig. 18-left shows a comparison with Trokhimovskiy et al. (2015). Trends are found to be in good agreement between the two despite the fact that scattering was neglected in our retrievals.

It is worth comparing retrievals based on the parametric approach to those based on the non-parametric approach (see Section 4.2 for a definition of the two models). This comparison is shown on Fig. 18-

<figure><img src="image_11.png" /><figcaption>Fig. 9. Retrievals of DOD without scattering (left) and of the surface temperature (right) from PFS data and compared to the predictions of the MCD. Correlation coefficients (R) between PFS and MCD datasets are given in the plots.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10. Forward modeling of the TIR1 spectra for orbit 205, SYN spectrum #8.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Residual (differences between data and model divided by the radiometric noise at 1-σ) distribution corresponding to the fit displayed in Fig. 10.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Retrieval of the spectral emissivity in the TIR1 band corresponding to the fit displayed in Fig. 10.</figcaption></figure>

right (see Section 5.3.3 for details on the non-parametric inversion settings). Column-integrated abundances obtained with these two approaches are strongly correlated yet the parametric inversion gives roughly 10% more water vapor than the non-parametric one. Due to the small prior uncertainty assigned to the saturation level altitude, only one H₂O parameter can be constrained on the vertical profile. However, the non-parametric approach gives a DOFFS greater than 1 (that is 1.4 on average, see Fig. 25 in Section 5.3.3), suggesting that the parametric model is too rigid, resulting in a systematic positive bias on the column

abundances. However, no firm conclusion should be made out of this remark since the $\chi^2$ improvement provided by the non-parametric approach is low.

Here, we present the results of the H₂O parametric inversion (11) applied to the TIR spectra only. Dust scattering is treated as described in Section 5.3.1. In addition to the parameters listed in Section 5.2 the retrieval also concerns the H₂O vmr at the surface, $x_0$ and the saturation level altitude $z_s$, with the same inversion parameters setting as for the NIR only case. The spectral emissivity of the surface in the TIR2 band (see Section 5.2 for the TIR1 emissivity) is also a derived product of the inversion.

Fig. 19 shows an example of TIR2 spectra modeling. In general, the residuals are not evenly distributed (Fig. 20). The features that we see in the continuum parameters solution are present in the prior (built from the TES data, see Section 3.4). As in the case of the TIR1 band (Section 5.2), the inversion parameters for emissivity ensures a good separation between absorption and surface spectral signatures.

The retrieved H₂O column abundances are found to be in good agreement with the SPICAM measurements presented in Trokhimovskiy et al. (2015) (Fig. 21, left). Owing to the lack of statistics in the latitude/solar longitude bins, it is not possible to claim a possible improvement resulting from the simultaneous temperature retrieval.

Unlike the NIR case, the inverse model formulation has little impact on the estimated H₂O column abundances (given the errors bars which are approximatively equal to 15% on average for both the parametric and the non-parametric approaches).
<figure><img src="image_15.png" /><figcaption>Fig. 13. (Left) Retrieval of the temperature profile. The prior corresponds to the MCD predictions for the spectra of Fig. 10. (Right): Error reduction factor on the atmospheric temperature parameters.</figcaption></figure>

### 5.3.3. Synergistic NIR + TIR inversion

Here, the retrievals of the H₂O vmr profile from the SYN spectra are presented. Inversion conditions are the same as the other inversion cases except for the parameters related to water vapor. We use the vmr profile predicted by the MCD as the reference profile. The estimated parameters correspond to an altitude profile of a scaling factor, whose first guess is set to one, and for which we consider a 100% prior uncertainty. Its positivity is forced by a lognormal variable substitution. A Gaussian kernel with a correlation length of 2 km ensures the vertical regularization. A comparison with the latest published SPICAM column abundances is presented on Fig. 22.

In the following, results obtained while using only NIR spectra are labeled by NIR, those obtained from TIR spectra only are labeled TIR and those obtained with both are labeled NIR + TIR. The averaging kernels obtained for a dry and extended ( ~ 6 pr-µm at $L_s = 20°$, latitude ~ 11°) or a wet and confined ( ~ 40 pr-µm at $L_s = 120°$, latitude ~ 60°) atmosphere are shown in Fig. 23. With respect to the NIR case, the synergistic approach produces an increase of the H₂O DOFFS by 30% in the dry configuration (1.77 for NIR + TIR vs. 1.35 for NIR and 1.15 for TIR) and by 40% in the wet configuration (1.63 for NIR + TIR vs. 1.18 for NIR and 1.08 for TIR and) thereby establishing the significant complementarity between the two spectral domains. This increase of “vertical resolution” is also reflected in the fact that lower altitude kernels (below 5 km) exhibit stronger peaks in the NIR + TIR case (like for the 4 and 6 km kernels) regardless of the wet or dry configuration.

The H₂O profile inversion for the dry and wet configurations are shown on Fig. 24-left. We note that the NIR and TIR retrievals give the same column abundances in the dry case while the NIR gives a significantly higher result in the wet case. This difference may be due to the fact that scattering is neglected in the NIR (see S6). The vmr estimators are marginally biased by the prior at the atmospheric levels where the posterior uncertainty is small compared to the prior uncertainty. This is illustrated by Fig. 24-right, which indicates that the uncertainty reduction is minimal near 2–3 km for the NIR case and around 3–4 km for the TIR case. Also, the vmr sensitivity region extends to higher altitudes in the wet case compared to the dry case.

A global assessment of the spectral synergy is performed (Fig. 25) by comparing the H₂O DOFFS in the NIR, TIR and NIR + TIR cases. As expected, the DOFFS differs between NIR and TIR cases as a result of the different information brought by the solar scattered and the thermally emitted radiances.

Combining the Bayesian and the maximum likelihood approaches one is able to define two unbiased estimators for the partial and the total column. This statement contrasts with the TIR inversion, where the DOFFS below 5 km is significantly lower than 1. This implies that the TIR estimator of the partial column is biased and thus depends on its initial first guess. Yet, the TIR inversion should provide an unbiased estimator of the total column. As a consequence, it is possible to make a relevant comparison between the NIR and TIR Bayesian estimators of the total column only when water vapor is confined in a 5-km thick layer near the surface (a case akin to the “wet” case investigated). The

<figure><img src="image_16.png" /><figcaption>Fig. 14. Modeling of a NIR spectrum for orbit 1485, SYN spectrum #1. The estimated column abundance is 25 ± 2 pr-µm.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 15. Residual distribution corresponding to the fit displayed in Fig. 14.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Continuum retrieval corresponding to the fit displayed in Fig. 14. The solution found while accounting for scattering is shown as well (DOD = 0.08, mean radius = 1 µm).</figcaption></figure>

DOFFS of the NIR + TIR inversion gives two main facts. First, the partial column estimator of the synergy is robust against the first guess, regardless of the vertical distribution. Second, the total column

abundance is better constrained than in the NIR- or TIR-only inversions (the difference between the total and partial DOFFS are greater in the NIR + TIR case).

However, if significant amount of water vapor is located above 5 km, NIR-only will provide a biased estimate. To illustrate these remarks, we compare the total and partial column estimates using two different first guesses for the water vapor profile (Fig. 27). Both correspond to the same total prior column. One is taken from the MCD predictions, while the other corresponds to the parametric model. From this comparison, it can be concluded that the NIR-only estimate of the total column is biased for wet cases (unless the right prior is used but in that case data does not bring additional information content), whereas it is the partial column estimate that is biased for TIR-only. The synergistic retrieval is robust in any configuration, “wet” or “dry.”

Fig. 26 shows zonal averages of the DOFFS for H₂O parameters up to TOA or below 5 km (referred as “partial column” hereinafter). We see a small latitudinal modulation, particularly in the NIR case. We check that it is uncorrelated with the local time or with the column abundance. NIR-only can yield an unbiased estimator of the H₂O partial column provided that its DOFFS is ≥1, which is the case for latitudes between 20°S and 40–60°N.

<figure><img src="image_19.png" /><figcaption>Fig. 17. Map of retrieved H₂O column abundances scaled at 6.1 hPa.TIR-only inversion.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 18. (Left) zonal average of scaled H₂O column abundances obtained for a parametric inversion. The inner error bars represent the statistical uncertainty, the outer error bar show the quadratic sum of statistical uncertainty and standard deviation in the latitudinal bin. The lines display the scaled column abundances extracted from Trokhimovskiy et al. (2015). (Right) scatter plot of the column abundances corresponding to parametric and non-parametric inversions. Statistic error bars on column abundances are ~6 and 10% on average in the parametric and non-parametric cases respectively.</figcaption></figure>

## 6. Discussion

In this study, we have developed and implemented a Bayesian approach to estimate water vapor abundances in the martian atmosphere from NIR and TIR spectra collected in a nadir viewing geometry. For this purpose, we have set up two processors that are operated

sequentially. The first processor generates co-localized spectra which are averaged over similar horizontal scales and which are characterized by high SNRs. The second processor retrieves water vapor from these spectra using a variety of configurations (i.e. NIR, TIR, NIR + TIR, parametric or non-parametric model, etc.). This set of tools have permitted an analysis of the different sources of instrumental uncertainties

<figure><img src="image_21.png" /><figcaption>Fig. 19. (Top) simulated TIR2 spectra for orbit 1040 (SYN spectra #2). (Bottom): Result obtained after the inversion of the continuum.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 20. Residual (expressed in %) distribution corresponding to the fit displayed in Fig. 19.</figcaption></figure>

(i.e. transfer function, spectral calibration) and a verification of the validity of some of the assumptions made in the forward model (that is those related to spectroscopy, scattering and solar source function).

Several improvements could be added, in particular the scattering by dust in the NIR domain even though results obtained in that case exhibit deviations from a no scattering case that generally remain within the radiometric noise. Nevertheless, when analysed at the annual scale, water vapor results obtained while neglecting scattering may produce some seasonal biases that may slightly distort the scientific interpretation. Another improvement could be designed for the parametric model which, in its current configuration, lacks flexibility and probably forces results too strongly around its intrinsic assumptions (for instance, the constant H₂O vmr set a saturation vapor pressure above the saturation level). Another limit of the current method is the neglect of the contribution of water ice clouds. For the time being, retrievals should avoid locations and seasons where thick water ice louds are susceptible to form. Some parameters indicating their presence could be easily derived from several spectral bands covered by the MEX instruments where water ice possesses diagnostic signatures. This supplementary selection stage should however be replaced at some point by a proper modelling of the cloud impact on the radiative transfer and thus on the collected spectra. This would create additional parameters for the retrieval and some of them might be difficult to constrain with sufficient confidence, like the cloud altitude and its vertical extension.

Regarding the performances of the synergistic method, an analysis

<figure><img src="image_23.png" /><figcaption>Fig. 22. Zonal averages of scaled H₂O column abundances obtained from a non-parametric NIR + TIR retrieval (see Fig. 18).</figcaption></figure>

of the measurements information content shows that, on average, the spectral synergy induces an increase of the H₂O DOFFS by 20 or >50% compared to a NIR- or TIR-only case. In that context, we also demonstrate that the NIR + TIR synergy makes it possible to establish an unbiased estimate of the partial column lying in the first 5 km above the surface, yielding at the same time an information on the total column

<figure><img src="image_24.png" /><figcaption>Fig. 21. (Left) zonal average of scaled H₂O column abundances obtained for parametric TIR estimation including or (Right) not including the simultaneous retrieval of the atmospheric temperature. As in Fig. 18, the inner error bars represent the statistical uncertainty, the outer error bar show the quadratic sum of statistical uncertainty and standard deviation in the latitudinal bin. The lines display the scaled column abundances extracted from Trokhimovskiy et al. (2015).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 23. Averaging kernels (see (8)) for NIR, TIR and NIR + TIR setups obtained in a dry (upper line of plots) and a wet (lower line of plots) configuration. Only the kernels computed for the altitude located at and below 10 km (see legend) are shown for the sake of clarity.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 24. (Left) H₂O vmr estimations vs. altitude. The retrieved column abundances are given in the legend box for all the types of inversion, compared to the prior estimate. (Right) post. to prior uncertainty ratio vs. altitude. The upper figures are for the “dry” atmosphere case, while the lower figures concern the “wet” atmosphere case.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 25. DOFFS for water vapor profile obtained in the NIR, TIR and NIR + TIR cases.</figcaption></figure>

abundance which is less tied to the partial column estimate than in the NIR or TIR only cases.

The results obtained while exploring two end-member situations (“wet” and “dry”) tend to suggest that the synergy performs optimally in a configuration with significant water confinement near the surface, like in the “wet” case. This characteristic of the synergy addresses situations where water vapor is expected to have maximum interactions with the surface. Therefore, an in-depth analysis of the synergistically generated H₂O dataset has the potential to shed new light on the processes potentially at work in the boundary layer for water vapor (for instance, it may potentially provide clues regarding the suspected role of adsorption/desorption to/from the martian regolith). Fig. 28 shows how the processed dataset permits a deduction of the vertical

partitioning of water as a function of latitude and for different periods of Ls. While the dataset seems too restricted to yield firm conclusion regarding the vertical behaviour of water, it is noteworthy that away from the equator, water vapor exhibits less confinement near the surface of Mars. Around the equator, the fraction kept in the first 5 km above the surface varies between 50 and 70%, whereas this fraction drops to ~40% at 30°S and 70°N. It could be indicative of a mechanism trapping water close to the surface around the equator and poleward of it, such as exchanges with the regolith. To be confirmed, such hypothesis would need to be tested against a more extended and more representative dataset.

Finally, one of the crucial advantages of the synergy lies in its robustness against changing atmospheric scenes and configurations. As

<figure><img src="image_28.png" /><figcaption>Fig. 26. (Upper) Zonal averages of H₂O DOFFS in the NIR. (Middle) TIR and (Lower) NIR + TIR cases for the total column (filled circle) and up to 5 km above the surface (filled square). The error bars give the standard deviation in the considered latitude/Ls bin.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Fig. 27.** Comparison between total (left) and partial (right) column abundances obtained using a parametric or a predicted (MCD) 1st guess of the water vapor profiles. From top to bottom: NIR, TIR, NIR + TIR cases. On each plot the value of the linear correlation factors, R, are indicated. The color code is the same as in Fig. 26.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. 28.** Plot showing the relative vertical partitioning of water as deduced from the NIR + TIR synergetic method. The quantity displayed, expressed in %, is the 5-km partial column ratioed with the total column abundance of water. Higher ratio values indicate stronger confinement of water near the surface.</figcaption></figure>
demonstrated in the previous subsection, the synergy produces re-
trievals that are essentially free from biases, contrarily to the single-
domain retrievals which cannot provide two nearly-independent in-
formation, one of them being necessarily biased. In addition, it appears
that the synergy not only supplies an additional information on the H₂O
profile, but it also provides column abundances with smaller error bars.

It is therefore reasonable to expect that the synergistic approach
could definitely settle the debate regarding the differences of water
vapor retrievals between MEX instruments which has created some
confusion as to which instrument lies closer to the truth. Saying that,
we believe that a complete reprocessing of the entire MEX mission using
a synergistic approach could establish a new reference climatology for
not only the horizontal and temporal distribution of water column
abundances but also for the partitioning between the near-surface layer
and the rest of the column above it.

7. Conclusion

A new method for retrieving information on water vapor in the
atmosphere of Mars has been developed and tested in the context of the
MEX mission. It is found that combining various spectral intervals from
the near- to the thermal infrared where water vapor possesses diag-
nostic absorption bands, brings additional information as well as ro-
bustness to the retrieval of the water vapor concentration when com-
pared to a single-interval retrieval approach which has been the
standard retrieval mode for water vapor so far. Whereas this so-called
"synergy" methodology cannot produce a detailed characterization of
the vertical profile of water vapor, it nonetheless sets fundamental
constraints on the partitioning of the water column between a near-
surface layer of 5 km depth and the rest of the column located above the
layer. This characteristic offers an optimized reassessment of the ex-
perimental data on Mars' atmospheric water collected from the orbit in
a nadir-looking mode. Because the latter is the prime mode of ob-
servations of the instruments used to explore Mars, the corresponding
data volume and its space-time coverage outmatch those of data col-
lected in dedicated viewing geometry like solar occultations (Fedorova
et al., 2009; Maltagiati et al., 2011) or limb staring (Todd Clancy et al.,
2017), and thus offers a rich opportunity to study a yet poorly explored
"territory" of martian observations.

We stress here that our guiding philosophy in accomplishing this
work was first and foremost the elaboration of a “proof of concept” to
be carefully validated against actual MEX data. The fact that only a
small sample (~0.1%) of the MEX dataset was tested to evaluate the
scientific benefit of this method imposes us to remain cautious and to
not draw too firm conclusions apart from the demonstration that a
significant potential for synergy exists between PFS and SPICAM ob-
servations and that this first exploration motivates a deeper scrutiny of
the entire MEX dataset.

Among the perspectives identified in this work, two emerge as the
main priorities. The first one is to modify the number of atmospheric
layers in the forward model to simplify the inclusion of aerosol scat-
tering. The 2-km vertical discretization employed in this analysis was
intended to make it possible to extract an information content of the
measurements by reducing the effect related to the discretization itself.
Since accounting for scattering is highly computer-time consuming,
only a reduced number of tests including this effect has been performed.
Inspired by our conclusions, it will be possible to reduce considerably
the number of layers in the model by considering, for example, a layer
between the surface and 5 km and typically two other ones above, one
of which could be used to simulate the presence of a water ice cloud and
its effect on radiative transfer, as discussed in the previous section.

As indicated above, the second priority is to extend our study to the
entire dataset of the Mars Express mission (more than 7 martian years
and ~200 000 SYN spectra to be processed) and to see to what extent
the data quality criteria currently defined could be relaxed to increase
the seasonal and geographical coverage of the synergistic analysis. This

work constitutes the final objective of this effort for H₂O since it will
subsequently permit an analysis of the spatial, seasonal and multi-an-
nual variations of the water vapor vertical layering which is known to
reflect the effects of transport and other physical processes prevailing in
the boundary layer.

This extension of the work could eventually encompass the sy-
nergistic retrieval of all the observed atmospheric species that possess
spectral signatures in multiple and distant wavelength ranges in the UV
to IR domain. The requirement to extract information on the vertical
from nadir is to have, for the same species, weighting functions (al-
ternatively averaging kernels) peaking at different altitudes depending
on the wavelength range. This is how a temperature profile can be
deduced from the observations of the CO₂ v₂ band at 15 µm. This is the
case for CO and a theoretical work has been conducted already by
Robert et al. (2017) who also addressed CH₄. Both species have indeed
separate signatures in the NIR and TIR ranges and thus fulfil the basic
requirements for applying a synergistic approach upon them. However,
the tests were performed on synthetic datasets and tests on experi-
mental data remain to be done.

A more complex synergistic retrieval is anticipated for dust and water ice clouds. In the case of dust; the combination of spectral intervals would rely on the 9 µm band of silicate for the TIR domain, while continuum values in the visible and in the near-infrared could form the complementing NIR domain. However, the difficulty encountered in NIR is the absence of distinct absorption bands of dust, whose effect on spectra is only a contribution to the continuum radiances and/or reflectances at a level commensurate with the contribution of surface albedo. Even more important, an ambiguity remains for particulate component such as mineral dust since variable physical or optical properties (e.g. *r*<sub>eff</sub>, *ν*<sub>eff</sub>) can have a similar signature on spectroscopy as a changing vertical distribution. Disentangling all these aspects might prove too challenging in a retrieval context, yet this option would certainly deserve some investigation since like for water, constraining the vertical distribution of concentration is key to a better understanding of the dust cycle.

Regarding water ice clouds, which have been neglected in this
study, the analysis of PFS data could include the 12 µm absorption band
of water ice to infer a cloud opacity that could be validated against
existing climatologies (that of THEMIS for example). Additional ab-
sorption features of water ice in the NIR domain (like the absorption
band at 1.5 µm studied by Madeleine et al., 2012, as well as the broader
one at 3 µm) could complete the synergy and provide a more reliable
estimate of the r<sub>eff</sub> of the particles composing the water ice clouds.
However, constraining cloud altitude can probably not leverage on
synergy considering the limited vertical extent of the cloud and the
variability of its altitude formation; unless some additional considera-
tion is included, such as assuming some coincidence between the cloud
altitude and that of water vapor saturation. Enlarging the synergistic
dataset to encompass cloudy scenes is nonetheless motivated both by
the importance clouds have for the global cycle of water on Mars and by
the fact their occurrence if sufficient to affect a significant fraction of
the synergistic information contained in the MEX dataset.

For the reasons described above, it appears that the synergistic
methodology opens new paths for the exploration of the atmosphere of
Mars, as it does for the Earth's atmosphere. With the ExoMars Trace Gas
Orbiter (TGO) of ESA, options for synergy are reduced compared to
MEX since no instrument onboard covers the thermal range of the band
of H₂O longward of 20 µm that is used in our study (the TIRVIM in-
strument onboard TGO which is similar to PFS has a cut-off wavelength
at 17 µm). However, species like CO, O₃ (the latter established from a
combination of UV and TIR bands), could be attempted in addition to
particulate components like dust or water ice whose potentials for sy-
nergy are yet to be demonstrated. Furthermore, synergy might be taken
into consideration in the design of future remote sensing instruments
and in the future selection of orbital payloads as it enhances the sci-
entific return and overcomes the intrinsic limitations imposed by a
constrained viewing geometry like nadir.

## Acknowledgments

This work was funded by the European Union's Horizon 2020 Programme under grant agreement UPWARDS-633127. No new data were used in producing this manuscript.

## References

Bandfield, J.L., Smith, M.D., 2003. Multiple emission angle surface-atmosphere separations of thermal emission spectrometer data. Icarus 161, 47–65. https://doi.org/10.1016/S0019-1035(02)00025-8.

Bertaux, J.-L., Korablev, O., Perrier, S., Quémerais, E., Montmessin, F., Leblanc, F., Lebonnois, S., Rannou, P., Lefèvre, F., Forget, F., Fedorova, A., Dimarellis, E., Reberac, A., Fonteyn, D., Chaufray, J.Y., Guibert, S., 2006. SPICAM on Mars Express: Observing modes and overview of UV spectrometer data and scientific results. J. Geophys. Res. (Planets) 111, 10. https://doi.org/10.1029/2006JE002690.

Brown, L.R., Humphrey, C.M., Gamache, R.R., 2007. CO₂-broadened water in the pure rotation and O₂ fundamental regions. J. Mol. Spectrosc. 246. https://doi.org/10.1016/j.jms.2007.07.010.

Christensen, P.R., et al., 2001. Mars Global Surveyor Thermal Emission Spectrometer experiment: Investigation description and surface science results. J. Geophys. Res. 106, 23823–23872. https://doi.org/10.1029/2000JE001370.

Christi, M., Stephens, G., 2004. Retrieving profiles of atmospheric CO₂ in clear sky and in the presence of thin cloud using spectroscopy from the near and thermal infrared: A preliminary case study. J. Geophys. Res. 109, D04316. https://doi.org/10.1029/2003JD004058.

Clancy, R.T., Wolff, M.J., Christensen, P.R., 2003. Mars aerosol studies with the MGS TES emission phase function observations: Optical depths, particle sizes, and ice cloud types versus latitude and solar longitude. J. Geophys. Res. 108 (E9), 5098. https://doi.org/10.1029/2003JE002058.

Clough, S.A., Iacono, M.J., Moncet, J.-L., 1992. Line-by-line calculation of atmospheric fluxes and cooling rates: application to water vapor. J. Geophys. Res. 97, 15761–15785.

Clough, S.A., Iacono, M.J., 1995. Line-by-line calculations of atmospheric fluxes and cooling rates II: application to carbon dioxide, ozone, methane, nitrous oxide, and the halocarbons. J. Geophys. Res. 100 (16).

Costantino, L., Cuesta, J., Emili, E., Coman, A., Foret, G., Dufour, G., Eremenko, M., Chailleux, Y., Beekmann, M., Flaud, J.-M., 2017. Potential of multispectral synergism for observing ozone pollution by combining IASI-NG and UVNS measurements from the EPS-SG satellite. Atmos. Meas. Tech. 10, 1281–1298. https://doi.org/10.5194/amt-10-1281-2017.

Encrenaz, T., Melchiorri, R., Fouchet, T., Drossart, P., Lellouch, E., Gondet, B., Bibring, J.-P., Langevin, Y., Titov, D., Ignatiev, N., Forget, F., 2005. A mapping of martian water sublimation during early northern summer using OMEGA/Mars Express. Astron. Astrophys. 441, L9–L12.

Eymet V. 2014. KDISTRIBUTION version 1.1.1 User Manual, available from http://www.lmd.jussieu.fr/~eymet.

Fedorova, A., Korablev, O., Bertaux, J.-L., Rodin, A., Kiselev, A., Perrier, S., 2006. Mars water vapor abundance from SPICAM IR spectrometer: seasonal and geographic distributions. J. Geophys. Res. (Planets) 111 (E10), 9–26. https://doi.org/10.1029/20 06JE002694. E09S07.

Fedorova, A., Korablev, O.I., Bertaux, J.-L., Rodin, A.V., Montmessin, F., Belyaev, D.A., Reberac, A., 2009. Solar infrared occultation observations by SPICAM experiment on Mars-Express: simultaneous measurements of the vertical distributions of H₂O, CO₂ and aerosol. Icarus 200, 96–117. https://doi.org/10.1016/j.icarus.2008.11.006.

Fedorova, A., Bertaux, J.-L., Betsis, D., Montmessin, F., Korablev, O., Maltagliati, L., Clarke, J., 2018. Water vapor in the middle atmosphere of Mars during the 2007 global dust storm. Icarus 300, 440–457. https://doi.org/10.1016/j.icarus.2017.09.025.

Fiorenza, C., Formisano, V., 2005. A solar spectrum for PFS data analysis. Planet. Space Sci. 53 (10), 1009–1016.

Forget, F., Haberle, R.M., Montmessin, F., Levrard, B., Head, J.W., 2006. Formation of glaciers on Mars by atmospheric precipitation at high obliquity. Science 311, 368–371.

Formisano, V., et al., 2005. The Planetary Fourier Spectrometer (PFS) onboard the European Mars Express mission. Planet. Space Sci. 53, 963–974.

Fouchet, T., Lellouch, E., Ignatiev, N.I., Forget, F., Titov, D.V., Tschimmel, M., Montmessin, F., Formisano, V., Giuranna, M., Maturilli, A., Encrenaz, T., 2007. Martian water vapor: Mars Express PFS/LW observations. Icarus 190 (1), 32–49.

Fu, Q., Liou, K.N., 1992. On the correlated k-distribution method for radiative transfer in nonhomogeneous atmospheres. J. Atmos. Sci. 49 (22).

Haberle, R.M., Jakosky, B.M., 1990. Sublimation and transport of water from the north residual polar cap on Mars. J. Geophys. Res. (Planets) 95, 1423–1437.

Humlíček, J., 1982. Optimized computation of the voigt and complex probability functions. J. Quant. Spectrosc. Radiat. Transf. 27 (12), 437–444.

Jakosky, B.M., 1983. The role of seasonal reservoir in Mars water cycle: II. Coupled model of the regolith, the polar caps and atmospheric transport. Icarus 55, 19–39.

Jakosky, B.M., Farmer, C.B., 1982. The seasonal and global behavior of water vapor in the Mars atmosphere: complete global results of the Viking atmospheric water detector experiment. J. Geophys. Res. 87, 2999–3019.

Korablev, O., et al., 2006. SPICAM IR acousto-optic spectrometer experiment on Mars

Express. J. Geophys. Res. 111, E09S03. https://doi.org/10.1029/2006JE002696.

Korablev, O., et al., 2013. Characterization of the stray light in a space borne atmospheric AOTF spectrometer. Opt. Express 21 (15), 18354–18360. https://doi.org/10.1364/ OE.21.018354.

Levrard, B., Forget, F., Montmessin, F., Laskar, J., 2004. Recent ice-rich deposits formed at high latitudes on Mars by sublimation of unstable equatorial ice during low ob-liquity. Nature 431, 1072–1075.

Laraia, A.L., Gamache, RobertR., Lamouroux, Julien, Gordon, JouliE., Rothman, LaurenceS., 2011. Total internal partition sums to support planetary remote sensing. Icarus 215, 391–400.

Letchworth, K.L., Benner, D.C., 2007. Rapid and accurate calculation of the Voigt func- tion. J. Quant. Spectrosc. Radiat. Transf. 107, 173–192.

Madeleine, J.-B., Forget, F., Head, J.W., Levrard, B., Montmessin, F., Millour, E., 2009. Amazonian northern mid-latitude glaciation on Mars: a proposed climate scenario. Icarus 203, 390–405. https://doi.org/10.1016/j.icarus.2009.04.037.

Madeleine, J.-B., et al., 2012. Aphelion water-ice cloud mapping and property retrieval using the OMEGA imaging spectrometer onboard Mars Express. J. Geophys. Res. 117, E00J07. https://doi.org/10.1029/2011JE003940.

Maltagliati, L., Montmessin, F., Fedorova, A., Korablev, O., Forget, F., Bertaux, J.-L., 2011. Evidence of water vapor in excess of saturation in the atmosphere of Mars. Science 333 (1868). https://doi.org/10.1126/science.1207957.

Maltagliati, L., Montmessin, F., Korablev, O., Fedorova, A., Forget, F., Määttänen, A., Lefèvre, F., Bertaux, J.-L., 2013. Annual survey of water vapor vertical distribution and water-aerosol coupling in the martian atmosphere observed by SPICAM/M XEs solar occultations. Icarus 223, 942–962. https://doi.org/10.1016/j.icarus.2012.12.

Markwardt, C.B. 2009, "Non-linear least squares fitting in IDL with MPFIT" in proc. Astronomical Data Analysis Software and Systems XVIII, Quebec, Canada, ASP Conference Series, Vol. 411, eds. D. Bohlender, P. Dowler & D. Durand (Astronomical Society of the Pacific: San Francisco), p. 251–254.

Montmessin, F., Forget, F., Rannou, P., Cabane, M., Haberle, R.M., 2004. Origin and role of water ice clouds in the martian water cycle as inferred from a general circulation model. J. Geophy. Res. (Planets) 109, 10004. https://doi.org/10.1029/2004JE02284.

Navarro, T., Madeleine, J.-B., Forget, F., Spiga, A., Millour, E., Montmessin, F., Määttänen, A., 2014. Global climate modeling of the martian water cycle with improved mi-crophysics and radiatively active water ice clouds. J. Geophys. Res. (Planets) 119. doi:10.10 02/2013JE04550.

Niro, F., Boulet, C., Hartmann, J.-M., 2004. Spectra calculations in central and wing re- gions of CO₂ IR bands between 10 and 20µm. I: model and laboratory measurements. J. Quant. Spectrosc. Radiat. Transf. 88, 483–497.

Pan, L.W., Edwards, D.P., Gille, J.C., Smith, M.W., Drummond, J.R., 1995. Satellite re- mote sensing of tropospheric CO and CH4 – forward model studies of the MOPITT instrument. Appl. Tech. 34, 6976–6988.

Pan, L.W., Gille, J.C., Edwards, D.P., Bailey, P.L., Rodgers, C.D., 1998. Retrieval of tro- pospheric carbon monoxide for the MOPITT experiment. J. Geophys. Res. 103, 32277–32290.

Perrier, S., Bertaux, J.L., Lefèvre, F., Lebonnois, S., Korablev, O., Fedorova, A., Montmessin, F., 2006. Global distribution of total ozone on Mars from SPICAM/MEX UV measurements. J. Geophys. Res. 111, E09S06. https://doi.org/10.1029/ 2006JE002681.

Razavi, A., Clerbaux, C., Wespes, C., Clarisse, L., Hurtmans, D., Payan, S., Camy-Peyret, C., Coheur, P.F., 2009. Characterization of methane retrievals from the IASI space-borne sounder. Atmos. Chem. Phys. 9, 7889–7899.

Richardson, M.I., Wilson, R.J., 2002. Investigation of the nature and stability of the martian seasonal water cycle with a general circulation model. J. Geophys. Res. 107 (E5), 5031. https://doi.org/10.1029/2001JE001536.

Robert, S., Camy-Peyret, C., Daerden, F., De Mazière, M., De Wachter, E., Neary, L., Vandenbussche, S., Vandaele, A.C., 2017. Two test-cases for synergistic detections in the Martian atmosphere: Carbon monoxide and methane. J. Quant. Spectrosc. Radiat. Transf. 189, 86–104. https://doi.org/10.1016/j.jqsrt.2016.11.003.

Rodgers, C., 2000. Inverse Methods for Atmospheric Sounding: Theory and Practice. World Sci., London.

Rothman, L.S., et al., 2013. The HITRAN2012 molecular spectroscopic database. J. Quant. Spectrosc. Radiat. Transf. 130, 4–50. https://doi.org/10.1016/j.jqsrt. 2013.07.002.

Smith, M.D., 2002. The annual cycle of water vapor on Mars as observed by the Thermal Emission Spectrometer. J. Geophys. Res. 107 (E11), 5115. https://doi.org/10.1029/ 2001JE001522.

Smith, M.D., 2004. Interannual variability in TES atmospheric observations of Mars during 1999-2003. Icarus 167, 148–165. https://doi.org/10.1016/j.jqsrt.2016.11.003.

Smith, M.D., Wolff, M.J., Clancy, R.T., Murchie, S.L., 2009. CRISM observations of water vapor and carbon monoxide. J. Geophys. Res. (Planets). https://doi.org/10.1029/ 2008JE003288.

Spurr, R., 2008. LIDORT and VLIDORT: linearized pseudo-spherical scalar and vector discrete ordinate radiative transfer models for use in remote sensing retrieval pro- blems. In: Kokhanovsky, A. (Ed.), Light Scattering Reviews 3 Sensing. Springer.

Stamnes, K., Tsay, SC., Wiscombe, W., Jayawaera, K., 1988. Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Appl. Opt. 27 (12), 2502–2509.

Todd Clancy, R., Smith, M.D., Lefèvre, F., McConnochie, T.H., Sandor, B.J., Wolff, M.J., Lee, S.W., Murchie, S.L., Toigo, A.D., Nair, H., Navarro, T., 2017. Vertical profiles of Mars 1.27 µm O₂ dayglow from MRO CRISM limb spectra: Seasonal/global behaviors, comparisons to LMDGCM simulations, and a global definition for Mars water vapor profiles. Icarus 293, 132–156. https://doi.org/10.1016/j.icarus.2017.04.011.

Trokhimovskiy, A.V., Fedorova, A., Korablev, O., Montmessin, F., Bertaux, J.-L., Rodin,
A., Smith, M.D., 2015. Mars' water vapor mapping by the SPICAM IR spectrometer: five martian years of observations. Icarus 251, 50–64.

Tschimmel, M., Ignatiev, N.I., Titov, D.V., Lellouch, E., Fouchet, T., Giuranna, M., Formisano, V., 2008. Investigation of water vapor on Mars with PFS/SW of Mars Express. Icarus 195 (2), 557–575.

Wells, R.J., 1999. Rapid Approximation to the Voigt/Faddeeva function and its derivatives. J. Quant. Spectrosc. Radiat. Transf. 62, 29–48.

Wolff, M.J., Todd Clancy, R., 2003. Constraints on the size of martian aerosols from thermal emission spectrometer observations. J. Geophys. Res 108 (E9), 5097. https://doi.org/10.1029/2003JE002057.

Wolff, M.J., et al., 2009. Wavelength dependence of dust aerosol single scattering albedo as observed by the Compact Reconnaissance Imaging Spectrometer. J. Geophys. Res. 114, E00D04. https://doi.org/10.1029/2009JE003350.

Worden, H.M., Deeter, M.N., Edwards, D.P., Gille, J.C., Drummond, J.R., Nédélec, P., 2010. Observations of near-surface carbon monoxide from space using MOPITT multispectral retrievals. J. Geophys. Res. 115, D18314. https://doi.org/10.1029/2010JD014242.