# Constraining the global composition of D/H and 18O/16O in Martian water using SOFIA/EXES

J. Alday, S. Aoki, C. Dewitt, Franck Montmessin, J. A. Holmes, Manish R.

Patel, J. P. Mason, Thérèse Encrenaz, M. J. Richter, P. G. J. Irwin, et al.

► To cite this version:

J. Alday, S. Aoki, C. Dewitt, Franck Montmessin, J. A. Holmes, et al. Constraining the global composition of D/H and 18O/16O in Martian water using SOFIA/EXES. Monthly Notices of the Royal Astronomical Society, In press. insu-04561301

HAL Id: insu-04561301

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

Submitted on 26 Apr 2024

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.
# Constraining the global composition of D/H and 18O/16O in Martian water using SOFIA/EXES

Juan Alday, Shohei Aoki, Curtis Dewitt, Franck Montmessin, James Holmes, Manish R. Patel, J. P. Mason, Thérèse Encrenaz, Matthew Richter, Patrick Gerard Joseph Irwin, et al.

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

## To cite this version:

Juan Alday, Shohei Aoki, Curtis Dewitt, Franck Montmessin, James Holmes, et al. Constraining the global composition of D/H and 18O/16O in Martian water using SOFIA/EXES. Monthly Notices of the Royal Astronomical Society, 2024, 530 (3), pp.2919-2932. 10.1093/mnras/stae1067. insu-04561301v2

HAL Id: insu-04561301

https://insu.hal.science/insu-04561301v2

Submitted on 8 Sep 2024

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

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

<figure><img src="image_2.png" /></figure>
# Constraining the global composition of D/H and ¹⁸O/¹⁶O in Martian water using SOFIA/EXES

J. Alday<fnref n="1" />,<fnref n="2" />*, S. Aoki<fnref n="1" />,<fnref n="3" />, C. DeWitt<fnref n="4" />,<fnref n="5" />, F. Montmessin<fnref n="6" />, J. A. Holmes<fnref n="2" />, M. R. Patel<fnref n="2" />, J. P. Mason<fnref n="2" />, T. Encrenaz<fnref n="7" />, M. J. Richter<fnref n="8" />, P. G. J. Irwin<fnref n="9" />, F. Daerden<fnref n="3" />, N. Terada<fnref n="10" />, and H. Nakagawa<fnref n="10" />

Accepted XXX. Received YYY; in original form ZZZ

## ABSTRACT

Isotopic ratios in water vapour carry important information about the water reservoir on Mars. Localised variations in these ratios can inform us about the water cycle and surface-atmosphere exchanges. On the other hand, the global isotopic composition of the atmosphere carries the imprints of the long-term fractionation, providing crucial information about the early water reservoir and its evolution throughout history. Here, we report the analysis of measurements of the D/H and ¹⁸O/¹⁶O isotopic ratios in water vapour in different seasons (Ls = 15°, 127°, 272°, 305°) made with SOFIA/EXES. These measurements, free of telluric absorption, provide a unique tool for constraining the global isotopic composition of Martian water vapour. We find the maximum planetary D/H ratio in our observations during the northern summer (D/H = 5.2 ± 0.2 with respect to the Vienna Standard Mean Ocean Water, VSMOW) and to exhibit relatively small variations throughout the year (D/H = 5.0 ± 0.2 and 4.3 ± 0.4 VSMOW during the northern winter and spring, respectively), which are to first order consistent though noticeably larger than the expectations from condensation-induced fractionation. Our measurements reveal the annually-averaged isotopic composition of water vapour to be consistent with D/H = 5.0 ± 0.2 and ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW. In addition, based on a comparison between the SOFIA/EXES measurements and the predictions from a Global Climate Model, we estimate the D/H in the northern polar ice cap to be ~5% larger than that in the atmospheric reservoir (D/Hice = 5.3 ± 0.3 VSMOW).

**Key words:** planets and satellites: atmospheres – planets and satellites: terrestrial planets – radiative transfer

## 1 INTRODUCTION

The water reservoir on Mars is thought to have evolved substantially in the past four billion years. Geomorphological and mineralogical features on the Martian surface suggest that the size of the water reservoir was much larger in the past than it is today. In addition, these surface features suggest that water could once have existed in its liquid form, in contrast to what the present-day climatic conditions enable (Carr & Clow 1981; Baker 2001). Therefore, reconciling the past and present evidence requires that at some point in history the water inventory available on the surface/atmosphere of Mars migrated to the subsurface, was incorporated to the crust, or escaped to space (e.g., Alsaeed & Jakosky (2019); Scheller et al. (2021)).

The isotopic composition of water vapour in Mars' atmosphere

carries the imprints of processes that have fractionated the water reservoir throughout history and may therefore help elucidate how Mars lost most of its water inventory. Among all the scenarios that can explain the loss of water, escape to space is the only process that produces isotopic differentiation. Therefore, the overall enrichment of D/H in Mars' atmosphere with respect to Earth of approximately a factor of six suggests that a substantial fraction of the atmosphere was lost to space (Owen et al. 1988; Bjoraker et al. 1989; Krasnopolsky et al. 1997). The enrichment in atmospheric D/H not only tells us that atmospheric escape has occurred, but also, if coupled with photochemical and evolution models, enables quantitative estimations of the amount of water lost to space throughout history (Yung et al. 1988; Krasnopolsky 2000; Cangi et al. 2023).

The ¹⁸O/¹⁶O isotopic ratio in water vapour is also enriched in the heavy isotopes with respect to Earth, although to a much lesser extent than D/H (¹⁸O/¹⁶O = 1.084 ± 0.010 with respect to Earth, Webster

<footnote n="1">¹Graduate School of Frontier Sciences, University of Tokyo, Kashiwa, Chiba, Japan</footnote>

<footnote n="2">²School of Physical Sciences, The Open University, Milton Keynes, United Kingdom</footnote>

<footnote n="3">³Royal Belgian Institute for Space Aeronomy, BIRA-IASB, Brussels, Belgium</footnote>

<footnote n="4">⁴SOFIA-USRA, NASA Ames Research Center, Moffett Field, CA, USA</footnote>

<footnote n="5">⁵Space Science Institute, Boulder, CO, USA</footnote>

<footnote n="6">⁶LATMOS/CNRS, Guyancourt, France</footnote>

<footnote n="7">⁷LESIA, Observatoire de Paris, PSL University, CNRS, UPMC, Univ. Denis Diderot, 92195 Meudon, France</footnote>

<footnote n="8">⁸Physics Department, University of California, Davis, CA, USA</footnote>

<footnote n="9">⁹AOPP, Department of Physics, University of Oxford, Oxford, United Kingdom</footnote>

<footnote n="10">¹⁰Graduate School of Science, Tohoku University, Sendai, Miyagi, Japan</footnote>

<footnote>* E-mail: juan.alday@open.ac.uk</footnote>
# Constraining the global composition of D/H and ¹⁸O/¹⁶O in Martian water using SOFIA/EXES

J. Alday<sup>1,2*</sup> S. Aoki,<sup>1,3</sup> C. DeWitt,<sup>4,5</sup> F. Montmessin,<sup>6</sup> J. A. Holmes,<sup>2</sup> M. R. Patel,<sup>2</sup> J. P. Mason,<sup>2</sup> T. Encrenaz,<sup>7</sup> M. J. Richter,<sup>8</sup> P. G. J. Irwin<sup>1,9</sup> F. Daerden,<sup>3</sup> N. Terada<sup>10</sup> and H. Nakagawa<sup>10</sup>

Accepted 2024 April 17. Received 2024 March 21; in original form 2023 October 12

## ABSTRACT

Isotopic ratios in water vapour carry important information about the water reservoir on Mars. Localized variations in these ratios can inform us about the water cycle and surface-atmosphere exchanges. On the other hand, the global isotopic composition of the atmosphere carries the imprints of the long-term fractionation, providing crucial information about the early water reservoir and its evolution throughout history. Here, we report the analysis of measurements of the D/H and ¹⁸O/¹⁶O isotopic ratios in water vapour in different seasons (L<sub>S</sub> = 15°, 127°, 272°, and 305°) made with the Echelon-Cross-Echelle Spectrograph (EXES) aboard the Stratospheric Observatory for Infrared Astronomy (SOFIA). These measurements, free of telluric absorption, provide a unique tool for constraining the global isotopic composition of Martian water vapour. We find the maximum planetary D/H ratio in our observations during the northern summer (D/H = 5.2 ± 0.2 with respect to the Vienna Standard Mean Ocean Water, VSMOW) and to exhibit relatively small variations throughout the year (D/H = 5.0 ± 0.2 and 4.3 ± 0.4 VSMOW during the northern winter and spring, respectively), which are to first order consistent though noticeably larger than the expectations from condensation-induced fractionation. Our measurements reveal the annually averaged isotopic composition of water vapour to be consistent with D/H = 5.0 ± 0.2 and ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW. In addition, based on a comparison between the SOFIA/EXES measurements and the predictions from a Global Climate Model, we estimate the D/H in the northern polar ice cap to be ~5 per cent larger than that in the atmospheric reservoir (D/H<sub>ice</sub> = 5.3 ± 0.3 VSMOW).

**Key words:** radiative transfer – planets and satellites: atmospheres – planets and satellites: terrestrial planets.

## 1. INTRODUCTION

The water reservoir on Mars is thought to have evolved substantially in the past four billion years. Geomorphological and mineralogical features on the Martian surface suggest that the size of the water reservoir was much larger in the past than it is today. In addition, these surface features suggest that water could once have existed in its liquid form, in contrast to what the present-day climatic conditions enable (Carr & Clow 1981; Baker 2001). Therefore, reconciling the past and present evidence requires that at some point in history the water inventory available on the surface/atmosphere of Mars migrated to the subsurface, was incorporated to the crust, or escaped to space (e.g. Alsaeed & Jakosky 2019; Scheller et al. 2021).

The isotopic composition of water vapour in Mars' atmosphere carries the imprints of processes that have fractionated the water reservoir throughout history and may therefore help elucidate how Mars lost most of its water inventory. Among all the scenarios that can explain the loss of water, escape to space is the only process that produces isotopic differentiation. Therefore, the overall enrichment of D/H in Mars' atmosphere with respect to Earth of approximately a factor of six suggests that a substantial fraction of the atmosphere was lost to space (Owen et al. 1988; Bjoraker, Mumma & Larson 1989; Krasnopolsky et al. 1997). The enrichment in atmospheric D/H not only tells us that atmospheric escape has occurred, but also, if coupled with photochemical and evolution models, enables quantitative estimations of the amount of water lost to space throughout history (Yung et al. 1988; Krasnopolsky 2000; Cangi et al. 2023).

The ¹⁸O/¹⁶O isotopic ratio in water vapour is also enriched in the heavy isotopes with respect to Earth, although to a much lesser extent

<footnote n="1"><sup>1</sup>Graduate School of Frontier Sciences, University of Tokyo, Kashiwa, Chiba 277-8561, Japan</footnote>

<footnote n="2"><sup>2</sup>School of Physical Sciences, The Open University, Milton Keynes MK7 6AA, UK</footnote>

<footnote n="3"><sup>3</sup>Royal Belgian Institute for Space Aeronomy, BIRA-IASB, Brussels B-1180, Belgium</footnote>

<footnote n="4"><sup>4</sup>SOFIA-USRA, NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

<footnote n="5"><sup>5</sup>Space Science Institute, Boulder, CO 80301, USA</footnote>

<footnote n="6"><sup>6</sup>LATMOS/CNRS, Guyancourt 78280, France</footnote>

<footnote n="7"><sup>7</sup>LESIA, Observatoire de Paris, PSL University, CNRS, UPMC, Université Denis Diderot, Meudon 92195, France</footnote>

<footnote n="8"><sup>8</sup>Physics Department, University of California, Davis, CA 95616., USA</footnote>

<footnote n="9"><sup>9</sup>AOPP, Department of Physics, University of Oxford, Oxford OX1 3PU, UK</footnote>

<footnote n="10"><sup>10</sup>Graduate School of Science, Tohoku University, Sendai, Miyagi 980-0845, Japan</footnote>

<footnote>*E-mail: juan.alday@open.ac.uk</footnote>
et al. (2013)). While atmospheric escape is also expected to enrich the atmosphere in $^{18}$O/$^{16}$O interactions between the oxygen reservoirs in H$_2$O, CO$_2$ and minerals on the crust make this connection much more complex (Jakosky & Phillips 2001). Indeed, water vapour in the atmosphere of Mars is more enriched in $^{18}$O/$^{16}$O than CO$_2$, which might indicate a photochemical transfer of depleted $^{18}$O/$^{16}$O from H$_2$O to CO$_2$ (Alday et al. 2021a,b).

Apart from the long-term isotopic fractionation that affects the overall isotopic composition of the atmospheric bulk, atmospheric processes at present can also produce localised variations of the D/H and $^{18}$O/$^{16}$O ratios (Cheng et al. 1999; Montmessin et al. 2005; Hu 2019). The variations of D/H in the atmosphere of Mars have been the subject of several observing campaigns using ground-based telescopes (Novak et al. 2011; Villanueva et al. 2015; Aoki et al. 2015; Krasnopolsky 2015, 2021; Encrenaz et al. 2016, 2018; Khayat et al. 2019), as well as recent observations of its vertical structure with the ExoMars Trace Gas Orbiter (Villanueva et al. 2021, 2022; Alday et al. 2021a). In general, the measurements reveal variability of D/H ranging from values as low as 1-2 with respect to the Earth's standard (Vienna Standard Mean Ocean Water; D/H = 155.76 × 10⁻⁶ and $^{18}$O/$^{16}$O = 2005.20 × 10⁻⁶) to values as high as 10 VSMOW, with averaged results typically ranging from 4-7 VSMOW. The localised variations of D/H revealed by the measurements are to first order consistent with the predictions from condensation-induced fractionation, although the measurements tend to show greater variability than predicted by the models (Montmessin et al. 2005; Rossi et al. 2022; Vals et al. 2022; Daerden et al. 2022).

Encrenaz et al. (2016, 2018) analysed observations made using the Echelon-Cross-Echelle Spectrograph (EXES) aboard the Stratospheric Observatory for Infrared Astronomy (SOFIA) in April 2014, March 2016 and January 2017 to map D/H in the atmosphere of Mars. Thanks to the high altitude of SOFIA, these measurements are mostly free of terrestrial contamination and therefore provide a unique means to characterise the D/H ratio on Mars. The analysis of the SOFIA/EXES measurements revealed relatively low disk-integrated values ranging from 4.0-4.5 VSMOW, and showed less variability than other ground-based observatories, possibly due to the relatively low spatial resolution of the instrument (~3") in comparison with the size of the Mars' disk (~5-15").

Here, we continue the investigations of the isotopic composition of water vapour using SOFIA/EXES. In particular, we analyse measurements made at different seasons on Mars to better understand the annual changes in the global values of D/H and $^{18}$O/$^{16}$O, aiming to disentangle the effect of condensation-induced fractionation from the measurements and infer the overall isotopic composition of the atmospheric bulk. In this study, we reprocess two of the previous observations reported by Encrenaz et al. (2018) using a full radiative transfer model, and include the analysis of two new observations made in October 2018 and April 2019. Additionally, we include the characterisation of the $^{18}$O/$^{16}$O ratio into our analysis of the data.

In the following sections, we describe in detail the characteristics of the measurements and the analysis of the data. Section 2 presents the observational setup of the SOFIA/EXES measurements and their main characteristics. The radiative transfer analysis and retrieval scheme are presented in Section 3. The retrieved data products are then presented in Section 4, as well as an in-depth analysis of the implications of these measurements. Finally, Section 5 summarises the main conclusions of this study.

## 2 SUMMARY OF THE SOFIA/EXES OBSERVATIONS

In this study, we analyse spectroscopic observations of Mars made with the EXES instrument aboard SOFIA during four different observing campaigns between 2016 and 2019 (see Table 1). EXES is an infrared spectrograph that operates between 4.5 and 28.3 µm in high (R = λ/Δλ = 50,000-100,000), medium (R = 5,000-20,000) and low (R = 1,000-3,000) spectral resolution modes (Richter et al. 2018). In the high-resolution mode used for the observations of Mars, 11-12 different diffraction orders are dispersed onto the detector frame covering an instantaneous spectral range between approximately 1383-1392 cm⁻¹. Additionally, the cross dispersion design of the instrument enables the characterisation of spatial variations along the height of the instrument's slit. In particular, the diffraction orders are projected along the detector's y-axis, and the 72-pixel width of each of these corresponds to the height of the slit (11.95"). On the other hand, the spectral information of the measurement is projected along the x-axis. The measured spectra in each diffraction order comes together with a quantification of the uncertainty provided by the SOFIA/EXES team derived from the Poisson noise of the detector counts propagated through all the nod-subtraction and flat fielding steps. Additionally, the EXES team report an uncertainty of ~10% in the flux calibration of the spectra (DeWitt et al. 2023). However, given that the retrievals in this study use normalised spectra, the derived isotopic ratios will not be affected by the uncertainty in the absolute flux calibration.

The observational setup of all four campaigns analysed in this study was similar, with the main difference being whether the observations were performed in mapping or non-mapping modes:

* In the observations made in March 2016 and October 2018, the apparent size of Mars in the sky (~11-13") was relatively big compared to the spatial resolution of the instrument (~2.35"). In this case, several scans were performed while sequentially moving the slit in constant steps (0.72") in a direction perpendicular to the slit height (i.e., along the direction of the slit width). The combination of the spatial information from the slit height together with the several scans along the slit width allows the extraction of radiance maps at various points across the Mars' disk.

* In the observations made in January 2017 and April 2019, since the apparent size of Mars was small compared to the spatial resolution, the measurements were not performed aiming to map the whole disk. Instead, several scans were performed, but only centering the slit at one or two locations on the planet. Therefore, the spatial information in these observations is only provided along the slit length, which is aligned approximately perpendicular to the equator, allowing the analysis of the latitudinal variations.

It must be noted that there is some uncertainty in the pointing, and the geometry of the observations must be accurately determined in the processing of the data. In particular, we calibrate the geometry of the observations by comparing the measured continuum radiance (ν = 1390.96 cm⁻¹) with that simulated with the forward model (see Supplementary Material). Figure 1 shows an example of the final geometry determined for the observations made during 2018 and 2019. For the observations made in 2017 and 2019 (i.e., non-mapping mode), there is an uncertainty in the geometry of the observations along the direction of the slit width. In order to estimate the impact of this uncertainty in our results, we performed the retrievals perturbing the geometry of the observations. The results from these tests suggest that the impact of this in the calculation of the isotopic ratios is small because of the relatively small size of the planet compared to the spatial resolution of the instrument (see Supplementary Material).
<table><caption>**Table 1.** Observational parameters of the SOFIA/EXES measurements analysed in this study.</caption><thead><tr><th>Date (YYYY-MM-DD)</th><th>Time (UTC)</th><th>L<sub>S</sub> (°)</th><th>Subobserver point (°, °)</th><th>Subsolar point (°, °)</th><th>Doppler velocity (km s<sup>-1</sup>)</th><th>Mars diameter (arcsec)</th><th>SOFIA altitude (km)</th></tr></thead><tbody><tr><td>2016-03-24</td><td>12:06:05</td><td>127</td><td>6.9°N –168.6°E</td><td>20.0°N 159.7°E</td><td>-15.41</td><td>10.9</td><td>13.7</td></tr><tr><td>2017-01-25</td><td>01:10:56</td><td>305</td><td>-26.4°N 20.1°E</td><td>-20.5°N 55.6°E</td><td>11.71</td><td>5.19</td><td>11.9</td></tr><tr><td>2018-10-20</td><td>02:58:17</td><td>272</td><td>-18.1°N –48.6°E</td><td>-24.4°N –2.9°E</td><td>11.19</td><td>13.19</td><td>11.7</td></tr><tr><td>2019-04-24</td><td>04:18:19</td><td>15</td><td>-4.30°N –46.4°E</td><td>-6.6°N –23.9°E</td><td>12.17</td><td>4.26</td><td>11.6</td></tr></tbody></table>

than D/H ($^{18}$O/$^{16}$O = 1.084 ± 0.010 with respect to Earth, Webster et al. 2013). While atmospheric escape is also expected to enrich the atmosphere in $^{18}$O/$^{16}$O, interactions between the oxygen reservoirs in H$_2$O, CO$_2$, and minerals on the crust make this connection much more complex (Jakosky & Phillips 2001). Indeed, water vapour in the atmosphere of Mars is more enriched in $^{18}$O/$^{16}$O than CO$_2$, which might indicate a photochemical transfer of depleted $^{18}$O/$^{16}$O from H$_2$O to CO$_2$ (Alday et al. 2021a, b).

Apart from the long-term isotopic fractionation that affects the overall isotopic composition of the atmospheric bulk, atmospheric processes at present can also produce localized variations of the D/H and $^{18}$O/$^{16}$O ratios (Cheng et al. 1999; Montmessin, Fouchet & Forget 2005; Hu 2019). The variations of D/H in the atmosphere of Mars have been the subject of several observing campaigns using ground-based telescopes (Novak, Mumma & Villanueva 2011; Aoki et al. 2015; Krasnopolsky 2015, 2021; Villanueva et al. 2015; Encrenaz et al. 2016, 2018; Khayat et al. 2019), as well as recent observations of its vertical structure with the ExoMars Trace Gas Orbiter (TGO, Alday et al. 2021a; Villanueva et al. 2021, 2022). In general, the measurements reveal variability of D/H ranging from values as low as 1–2 with respect to the Earth’s standard [Vienna Standard Mean Ocean Water (VSMOW), D/H = 155.76 × 10$^{-6}$ and $^{18}$O/$^{16}$O = 2005.20 × 10$^{-6}$] to values as high as 10 VSMOW, with averaged results typically ranging 4–7 VSMOW. The localized variations of D/H revealed by the measurements are to first order consistent with the predictions from condensation-induced fractionation, although the measurements tend to show greater variability than predicted by the models (Montmessin et al. 2005; Daerden et al. 2022; Rossi et al. 2022; Vals et al. 2022).

Encrenaz et al. (2016, 2018) analysed observations made using the Echelon-Cross-Echelle Spectrograph (EXES) aboard the Stratospheric Observatory for Infrared Astronomy (SOFIA) in 2014 April, 2016 March, and 2017 January to map D/H in the atmosphere of Mars. Thanks to the high altitude of SOFIA, these measurements are mostly free of terrestrial contamination and therefore provide a unique means to characterize the D/H ratio on Mars. The analysis of the SOFIA/EXES measurements revealed relatively low disc-integrated values ranging 4.0–4.5 VSMOW, and showed less variability than other ground-based observatories, possibly due to the relatively low spatial resolution of the instrument (~3 arcsec) in comparison with the size of the Mars’ disc (~5–15 arcsec).

Here, we continue the investigations of the isotopic composition of water vapour using SOFIA/EXES. In particular, we analyse measurements made at different seasons on Mars to better understand the annual changes in the global values of D/H and $^{18}$O/$^{16}$O, aiming to disentangle the effect of condensation-induced fractionation from the measurements and infer the overall isotopic composition of the atmospheric bulk. In this study, we reprocess two of the previous observations reported by Encrenaz et al. (2018) using a full radiative transfer model, and include the analysis of two new observations made in 2018 October and 2019 April. Additionally, we include the characterization of the $^{18}$O/$^{16}$O ratio into our analysis of the data.

In the following sections, we describe in detail the characteristics of the measurements and the analysis of the data. Section 2 presents the observational set-up of the SOFIA/EXES measurements and their main characteristics. The radiative transfer analysis and retrieval scheme are presented in Section 3. The retrieved data products are then presented in Section 4, as well as an in-depth analysis of the implications of these measurements. Finally, Section 5 summarizes the main conclusions of this study.

## 2. SUMMARY OF THE SOFIA/EXES OBSERVATIONS

In this study, we analyse spectroscopic observations of Mars made with the EXES instrument aboard SOFIA during four different observing campaigns between 2016 and 2019 (see Table 1). EXES is an infrared spectrograph that operates between 4.5 and 28.3 µm in high- ($R = \lambda/\Delta\lambda = 50000–100000$), medium- ($R = 5000–20000$), and low-spectral resolution ($R = 1000–3000$) modes (Richter et al. 2018). In the high-resolution mode used for the observations of Mars, 11–12 different diffraction orders are dispersed on to the detector frame covering an instantaneous spectral range between approximately 1383 and 1392 cm$^{-1}$. Additionally, the cross-dispersion design of the instrument enables the characterization of spatial variations along the height of the instrument’s slit. In particular, the diffraction orders are projected along the detector’s y-axis, and the 72-pixel width of each of these corresponds to the height of the slit (11.95 arcsec). On the other hand, the spectral information of the measurement is projected along the x-axis. The measured spectra in each diffraction order come together with a quantification of the uncertainty provided by the SOFIA/EXES team derived from the Poisson noise of the detector counts propagated through all the nod-subtraction and flat-fielding steps. Additionally, the EXES team report an uncertainty of ~10 per cent in the flux calibration of the spectra (DeWitt, Montiel & Rashman 2023). However, given that the retrievals in this study use normalized spectra, the derived isotopic ratios will not be affected by the uncertainty in the absolute flux calibration.

The observational set-up of all four campaigns analysed in this study was similar, with the main difference being whether the observations were performed in mapping or non-mapping modes:

(i) In the observations made in 2016 March and 2018 October, the apparent size of Mars in the sky (~11–13 arcsec) was relatively big compared with the spatial resolution of the instrument (~2.35 arcsec). In this case, several scans were performed while sequentially moving the slit in constant steps (0.72 arcsec) in a direction perpendicular to the slit height (i.e. along the direction of the slit width). The combination of the spatial information from the slit height together with the several scans along the slit width allows the extraction of radiance maps at various points across the Mars’ disc.

(ii) In the observations made in 2017 January and 2019 April, since the apparent size of Mars was small compared with the spatial resolution, the measurements were not performed aiming to map the
<table><caption>**Table 1.** Observational parameters of the SOFIA/EXES measurements analysed in this study.</caption><thead><tr><th>Date (YYYY-MM-DD)</th><th>Time (UTC)</th><th>L<sub>S</sub> (°)</th><th>Sub-observer point (°, °)</th><th>Sub-solar point (°, °)</th><th>Doppler velocity (km/s)</th><th>Mars diameter (arcsec)</th><th>SOFIA altitude (km)</th></tr></thead><tbody><tr><td>2016-03-24</td><td>12:06:05</td><td>127</td><td>6.9°N -168.6°E</td><td>20.0°N 159.7°E</td><td>-15.41</td><td>10.9</td><td>13.7</td></tr><tr><td>2017-01-25</td><td>01:10:56</td><td>305</td><td>-26.4°N 20.1°E</td><td>-20.5°N 55.6°E</td><td>11.71</td><td>5.19</td><td>11.9</td></tr><tr><td>2018-10-20</td><td>02:58:17</td><td>272</td><td>-18.1°N -48.6°E</td><td>-24.4°N -2.9°E</td><td>11.19</td><td>13.19</td><td>11.7</td></tr><tr><td>2019-04-24</td><td>04:18:19</td><td>15</td><td>-4.30°N -46.4°E</td><td>-6.6°N -23.9°E</td><td>12.17</td><td>4.26</td><td>11.6</td></tr></tbody></table>

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Example of forward model and comparison with the SOFIA/EXES data at $v = 1390.96 \text{ cm}^{-1}$ for the observations made in October 2018 (top) and April 2019 (bottom). The left panels show the results of the forward model described in the manuscript using the extracted parameters from the OpenMARS dataset. The middle panels show the measurements from SOFIA/EXES in the two campaigns, which were made using the mapping and non-mapping modes, respectively. The right panels show the comparison between the normalised radiances between the model and the measurement, where each of the lines represents the radiance along the slit height for each of the scans along the slit width direction. The difference between the measured and modelled radiance is shown as thin lines with an offset of -0.2 for the clarity of the figure.</figcaption></figure>

# 3 RADIATIVE TRANSFER ANALYSIS

In the following sections we describe in detail the radiative transfer analysis of the SOFIA/EXES observations, starting from the description of the forward model, the retrieval pipeline, and the calculation of the isotopic ratios in the Martian atmosphere.

## 3.1 Description of the forward model

In this section, we describe the forward model developed for the analysis of the SOFIA/EXES observations of Mars. First of all, the geometry of the observations is calculated using the SPICE toolkit (Acton et al. 2018), which we use to determine the observation angles at any location on the planet, as well as some important geometric parameters summarised in Table 1. Once the geometry has been determined, the surface temperature and vertical profiles of pressure, temperature, water vapour mixing ratio and dust and water ice optical depths are extracted from the OpenMARS database (Holmes et al.

2020) in a grid of locations throughout the Mars disk at the exact time of the observations, which serve as the input parameters in our reference atmosphere.

The OpenMARS database assimilates different atmospheric parameters derived from several spacecraft into a Mars Global Circulation Model (GCM) to obtain an accurate and global representation of the atmosphere of Mars. The characteristics of the assimilated parameters into the GCM vary depending on the availability of measurements. In this study, we use the OpenMARS database v2 (Holmes et al. 2023), which includes the assimilation of temperature profiles and dust column optical depths measured by the Mars Climate Sounder (MCS) onboard the Mars Reconnaissance Orbiter (MRO) (Kleinböhl et al. 2017). In addition, the dataset also assimilates observations from the ExoMars Trace Gas Orbiter (TGO) after March 2018, including measurements of the water vapour abundances (Holmes et al. 2022).

The reference surface and atmospheric parameters extracted from the OpenMARS dataset are then used to generate high-resolution
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of forward model and comparison with the SOFIA/EXES data at $v = 1390.96 \text{ cm}^{-1}$ for the observations made in 2018 October (top panel) and 2019 April (bottom panel). The left panels show the results of the forward model described in the manuscript using the extracted parameters from the OpenMARS data set. The middle panels show the measurements from SOFIA/EXES in the two campaigns, which were made using the mapping and non-mapping modes, respectively. The right panels show the comparison between the normalized radiances between the model and the measurement, where each of the lines represents the radiance along the slit height for each of the scans along the slit width direction. The difference between the measured and modelled radiance is shown as thin lines with an offset of -0.2 for the clarity of the figure.</figcaption></figure>

whole disc. Instead, several scans were performed, but only centring the slit at one or two locations on the planet. Therefore, the spatial information in these observations is only provided along the slit length, which is aligned approximately perpendicular to the equator, allowing the analysis of the latitudinal variations.

It must be noted that there is some uncertainty in the pointing, and the geometry of the observations must be accurately determined in the processing of the data. In particular, we calibrate the geometry of the observations by comparing the measured continuum radiance ($v = 1390.96 \text{ cm}^{-1}$) with that simulated with the forward model (see Supplementary Material). Fig. 1 shows an example of the final geometry determined for the observations made during 2018 and 2019. For the observations made in 2017 and 2019 (i.e. non-mapping mode), there is an uncertainty in the geometry of the observations along the direction of the slit width. In order to estimate the impact of this uncertainty in our results, we performed the retrievals perturbing the geometry of the observations. The results from these tests suggest that the impact of this in the calculation of the isotopic ratios is small because of the relatively small size of the planet compared with the spatial resolution of the instrument (see Supplementary Material).

### 3. RADIATIVE TRANSFER ANALYSIS

In the following sections, we describe in detail the radiative transfer analysis of the SOFIA/EXES observations, starting from the description of the forward model, the retrieval pipeline, and the calculation of the isotopic ratios in the Martian atmosphere.

#### 3.1 Description of the forward model

In this section, we describe the forward model developed for the analysis of the SOFIA/EXES observations of Mars. First of all, the geometry of the observations is calculated using the SPICE toolkit (Acton et al. 2018), which we use to determine the observation angles at any location on the planet, as well as some important geometric parameters summarized in Table 1. Once the geometry has been determined, the surface temperature and vertical profiles of pressure, temperature, water vapour mixing ratio and dust and water ice optical depths are extracted from the OpenMARS data base (Holmes, Lewis & Patel 2020) in a grid of locations throughout the Mars disc at the exact time of the observations, which serve as the input parameters in our reference atmosphere.

The OpenMARS data base assimilates different atmospheric parameters derived from several spacecraft into a Mars Global Circulation Model (GCM) to obtain an accurate and global representation of the atmosphere of Mars. The characteristics of the assimilated parameters into the GCM vary depending on the availability of measurements. In this study, we use the OpenMARS data base v2 (Holmes, Lewis & Patel 2023), which includes the assimilation of temperature profiles and dust column optical depths measured by the Mars Climate Sounder onboard the Mars Reconnaissance Orbiter (Kleinböhl, Friedson & Schofield 2017). In addition, the data set also assimilates observations from the ExoMars TGO after 2018 March, including measurements of the water vapour abundances (Holmes et al. 2022).
spectra in each grid location on the Mars disk using the NEMESIS algorithm (Irwin et al. 2008). The radiative transfer calculations are performed using line-by-line modelling and account for the thermal emission of the atmosphere and surface including multiple scattering effects, and for the Doppler shift between the Earth and Mars. The spectroscopic parameters of the water vapour isotopologues are taken from the 2020 version of the HITRAN database (Gordon et al. 2022), including the pressure-broadening coefficients suitable for a CO₂-dominated atmosphere reported in Devi et al. (2017) and Régalia et al. (2019). The spectroscopic parameters of CO₂ are also taken from the 2020 version of the HITRAN database, using in this case complete self-broadening for the calculation of the pressure-broadening coefficients. The optical properties of Martian dust are taken from Connour et al. (2022), assuming an effective radius of r<sub>eff</sub> = 1.5 µm. In the case of water ice, the optical properties are calculated using Mie Theory following the refractive indices reported in Warren & Brandt (2008), and assuming a log-normal distribution of with r<sub>eff</sub> = 3.5 µm and ν<sub>eff</sub> = 0.1 (Clancy 2003).

The modelled spectra in each of the locations along the Mars disk are then multiplied by the Earth's atmospheric transmission. Earth's vertical profiles of pressure, temperature and volume mixing ratios above the altitude of SOFIA are extracted from the NASA/MERRA-2 model using the online extraction tool from the Planetary Spectrum Generator (Villanueva et al. 2018). These profiles are then used to calculate the line-of-sight density for each gas, their optical depth, and finally the Earth transmittance. It must be noted that SOFIA, as an airborne observatory making its observations in the stratosphere, is above the atmospheric bulk and allows observing other celestial objects with low telluric contamination. This, together with the high-spectral resolution and Doppler shift, generally allows the selection of spectral windows with no telluric contamination (see section 3.2).

Once the radiative transfer calculations have been performed (i.e., Mars spectra and telluric transmittance), the modelled spectra are convolved with the Point Spread Function (PSF) to account for the instrumental response of SOFIA/EXES in terms of spatial resolution. In this study, we determined empirically the PSF of the instrument using observations of a red giant star made during the same night as the Mars' measurements in October 2018. Our calculations suggest that the PSF shows no dependence with wavelength and can be modelled using a Voigt function with parameters σ = 0.57 and γ = 0.65 (see Supplementary Material). Finally, the spatially-convolved spectra are convolved with the Instrument Lineshape (ILS) to account for the spectral response of the instrument. In this study, the ILS is determined empirically by using measurements of spectral emissions in a C₂H₂ gas-cell (see Supplementary Material).

Figure 1 shows a comparison of the forward model described here and the Mars measurements from SOFIA/EXES made in October 2018 and April 2019. While the absolute calibration of the measurements is unknown, the spatial distribution of the normalised modelled and measured radiance throughout the images may be compared. While some offsets between the measured and modelled radiance might depict small departures in the surface temperature in the model, we find the overall agreement in the radiance distribution to be satisfactory, suggesting that the instrument's spatial response is accurately accounted for in our forward model. Similarly, Figure 2 shows examples of the comparison between measured and modelled strong CO₂ absorption lines for each of the analysed observations in this study. This comparison suggests that the described forward model satisfactorily reproduces the CO₂ absorption features, which supports the choice of parameters in our reference atmosphere and the suitability of our forward model.

## 3.2 Selection of spectral windows

The spectroscopic observations made by SOFIA/EXES achieve high-spectral resolution in a relatively wide spectral range encompassing several absorption lines of the different water isotopes and carbon dioxide. Figure 3 shows an example of a spectrum of Mars at the resolution of the SOFIA/EXES measurements, highlighting the spectral features of different gases and isotopes. Since the main goal of this study is to measure the isotopic composition of water vapour, the absorption lines for targeted isotopes ($H_2^{16}O$, HD$^{16}O$ and $H_2^{18}O$) should be ideally comparable in optical depth and free of telluric and CO₂ contamination. Based on these requirements, we select two spectral windows to derive the column densities of each of the water isotopologues (see Figures 3 and 4).

The lines selected to retrieve the column abundances of $H_2^{16}O$ and HD$^{16}O$ are the same ones used by Encrenaz et al. (2018). These lines are ideal for deriving the D/H ratio because they are of similar optical depth, do not overlap with any CO₂ spectral features, and are free of telluric contamination for relatively large Doppler shifts (i.e., free of telluric contamination for all observations analysed in this study). On the other hand, the different sensitivities of the line strengths of each absorption line to the temperature field can potentially induce systematic errors in the derivation of the D/H ratio. In particular, the strength of the HDO absorption lines is weakly dependent on temperature, but a departure of $\Delta T = 10$ K in the temperature field can induce changes of approximately 50% and 20% in the strengths of the H₂O lines. Therefore, the D/H ratio derived from the measurements may be systematically under or overestimated if the temperature field is not accurately represented in the reference atmosphere. This dependence of the absorption line on the temperature field likely represents the major source of uncertainty in the analysis of these measurements.

In the case of the $H_2^{18}O$ isotopologue, the selection of suitable spectral windows within the measured spectral range is more challenging. The strongest absorption line of $H_2^{18}O$ occurs at 1389.91 cm⁻¹ and its optical depth is similar to the transitions of HDO and H₂O. However, this absorption line completely overlaps with absorption by carbon dioxide within the spectral resolution of SOFIA/EXES (see Figures 3 and 4). Since the overlapping absorption of $H_2^{18}O$ and CO₂ can lead to degeneracies in the spectral fitting that can impact the accuracy and reliability of the retrieval, we use instead another absorption line of $H_2^{18}O$, free of telluric and CO₂ contamination, located close to 1383.33 cm⁻¹. Unfortunately, the absorption line at 1383.33 cm⁻¹, used for the derivation of the $^{18}O/^{16}O$ ratio in this study, is located at the edge of the detector frame and was not measured in the observations made in October 2018 and April 2019. Therefore, we only report measurements of $^{18}O/^{16}O$ for the observations made in March 2016 and January 2017.

## 3.3 Description of the retrieval scheme

One of the main challenges to model the spectra measured in each pixel of the SOFIA/EXES images is the reconstruction of the field-of-view (FOV): due to the relatively poor spatial resolution of the measurements, spectra at different locations all over the Mars disk have to be calculated and convolved with the PSF in each iteration to accurately model the measured spectrum in each pixel of the detector.

The reconstruction of the FOV thus makes the analysis of the spectra computationally expensive. In order to overcome this problem, and given that most of the parameters in our atmospheric model are fixed, we modify the forward model described in section 3.1 to make use of tabulated look-up tables of high-resolution spectra that we
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Example of forward model and comparison with the SOFIA/EXES observations analysed in this study. The points represent the SOFIA/EXES measurements in a spectral range encompassing two strong absorption features of CO₂, while the solid lines correspond to the forward models calculated using the atmospheric parameters extracted from the OpenMARS data base. The cross on the radiance maps highlights the location of the spectra shown in the different panels.</figcaption></figure>

The reference surface and atmospheric parameters extracted from the OpenMARS data set are then used to generate high-resolution spectra in each grid location on the Mars disc using the NEMESIS algorithm (Irwin et al. 2008). The radiative transfer calculations are performed using line-by-line modelling and account for the thermal emission of the atmosphere and surface including multiple scattering effects, and for the Doppler shift between the Earth and Mars. The spectroscopic parameters of the water vapour isotopologues are taken from the 2020 version of the HITRAN data base (Gordon et al. 2022), including the pressure-broadening coefficients suitable for a CO₂-dominated atmosphere reported in Devi et al. (2017) and Régalia et al. (2019). The spectroscopic parameters of CO₂ are also taken from the 2020 version of the HITRAN data base, using in this case complete self-broadening for the calculation of the pressure-broadening coefficients. The optical properties of Martian dust are taken from Connour et al. (2022), assuming an effective radius of $r_{\text{eff}} = 1.5$ µm. In the case of water ice, the optical properties are calculated using Mie Theory following the refractive indices reported in Warren & Brandt (2008), and assuming a log-normal distribution of with $r_{\text{eff}} = 3.5$ µm and $v_{\text{eff}} = 0.1$ (Clancy 2003).

The modelled spectra in each of the locations along the Mars disc are then multiplied by the Earth's atmospheric transmission. Earth's vertical profiles of pressure, temperature, and volume mixing ratios above the altitude of SOFIA are extracted from the NASA/MERRA-2 model using the online extraction tool from the Planetary Spectrum Generator (Villanueva et al. 2018). These profiles are then used to calculate the line-of-sight density for each gas, their optical depth, and finally the Earth transmittance. It must be noted that SOFIA, as an airborne observatory making its observations in the stratosphere, is above the atmospheric bulk and allows observing other celestial objects with low telluric contamination. This, together with the high-spectral resolution and Doppler shift, generally allows the selection of spectral windows with no telluric contamination (see Section 3.2).

Once the radiative transfer calculations have been performed (i.e. Mars spectra and telluric transmittance), the modelled spectra are convolved with the point spread function (PSF) to account for the

instrumental response of SOFIA/EXES in terms of spatial resolution. In this study, we determined empirically the PSF of the instrument using observations of a red giant star made during the same night as the Mars' measurements in 2018 October. Our calculations suggest that the PSF shows no dependence with wavelength and can be modelled using a Voigt function with parameters $\sigma = 0.57$ and $\gamma = 0.65$ (see Supplementary Material). Finally, the spatially convolved spectra are convolved with the instrument lineshape (ILS) to account for the spectral response of the instrument. In this study, the ILS is determined empirically by using measurements of spectral emissions in a C₂H₂ gas-cell (see Supplementary Material).

Fig. 1 shows a comparison of the forward model described here and the Mars measurements from SOFIA/EXES made in 2018 October and 2019 April. While the absolute calibration of the measurements is unknown, the spatial distribution of the normalized modelled and measured radiance throughout the images may be compared. While some offsets between the measured and modelled radiance might depict small departures in the surface temperature in the model, we find the overall agreement in the radiance distribution to be satisfactory, suggesting that the instrument's spatial response is accurately accounted for in our forward model. Similarly, Fig. 2 shows examples of the comparison between measured and modelled strong CO₂ absorption lines for each of the analysed observations in this study. This comparison suggests that the described forward model satisfactorily reproduces the CO₂ absorption features, which supports the choice of parameters in our reference atmosphere and the suitability of our forward model.

## 3.2 Selection of spectral windows

The spectroscopic observations made by SOFIA/EXES achieve high-spectral resolution in a relatively wide spectral range encompassing several absorption lines of the different water isotopes and carbon dioxide. Fig. 3 shows an example of a spectrum of Mars at the resolution of the SOFIA/EXES measurements, highlighting the spectral features of different gases and isotopes. Since the main goal
<figure><img src="image_2.png" /></figure>

**Figure 2.** Example of forward model and comparison with the SOFIA/EXES observations analysed in this study. The black points represent the SOFIA/EXES measurements in a spectral range encompassing two strong absorption features of CO₂, while the red lines correspond to the forward models calculated using the atmospheric parameters extracted from the OpenMARS database. The white cross on the radiance maps highlights the location of the spectra shown in the different panels.

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

**Figure 3.** Modelled spectrum of the atmosphere of Mars at the spectral resolution of SOFIA/EXES. The coloured shadowed regions indicate the position and strength of the absorption lines of the water isotopes, while the black solid line indicates the absorption by CO₂. The vertical dashed lines show the position of the absorption lines of each isotopologue considered in this study.

interpolate over instead of performing the radiative transfer calculations in each iteration. In particular, we generate high-resolution spectra on a grid of locations all over the Mars disk by fixing all parameters in our reference atmosphere (i.e. pressure, temperature and aerosol densities) but the total abundance of the water vapour isotopes. Specifically, the vertical profiles of water vapour from OpenMARS in each location are multiplied by scaling factors ranging from 0 to 3 with a step of 0.1, which is equivalent to generating spectra for different water column densities. Similarly, spectra for the minor isotopes are also generated using the water vapour profiles from OpenMARS and scaling them following an isotopic composition given by $D/H = 5 \text{ VSMOW}$ and $¹⁸\text{O}/¹⁶\text{O} = 1 \text{ VSMOW}$. Once the look-up tables for the different isotopes have been generated, the spectra at an arbitrary column density and location can be calculated by interpolation.

Once the Mars high-resolution spectra are calculated at the different locations, the forward model follows the steps described in section 3.1: calculation of the telluric transmittance, spatial convolution with the PSF, and spectral convolution with the ILS. Finally, we multiply the modelled spectrum with a second order polynomial to accurately fit the continuum radiance. In summary, for each measurement, we fit the following model parameters:

* Scaling factor to the reference column abundance for any of the water isotopes. Note that one single scaling factor is applied to all locations on the planet.

* Scaling factor for the line-of-sight density of H₂O, CO₂ and CH₄ in the Earth's atmosphere. Note that while this capability is implemented in our forward model, the selection of spectral windows
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Modelled spectrum of the atmosphere of Mars at the spectral resolution of SOFIA/EXES. The coloured shadowed regions indicate the position and strength of the absorption lines of the water isotopes, while the black solid line indicates the absorption by CO₂. The vertical dashed lines show the position of the absorption lines of each isotopologue considered in this study.</figcaption></figure>

of this study is to measure the isotopic composition of water vapour, the absorption lines for targeted isotopes ($H_2^{16}O$, $HD^{16}O$, and $H_3^{18}O$) should be ideally comparable in optical depth and free of telluric and CO₂ contamination. Based on these requirements, we select two spectral windows to derive the column densities of each of the water isotopologues (see Figs 3 and 4).

The lines selected to retrieve the column abundances of $H_2^{16}O$ and $HD^{16}O$ are the same ones used by Encrenaz et al. (2018). These lines are ideal for deriving the D/H ratio because they are of similar optical depth, do not overlap with any CO₂ spectral features, and are free of telluric contamination for relatively large Doppler shifts (i.e. free of telluric contamination for all observations analysed in this study). On the other hand, the different sensitivities of the line strengths of each absorption line to the temperature field can potentially induce systematic errors in the derivation of the D/H ratio. In particular, the strength of the HDO absorption lines is weakly dependent on temperature, but a departure of $\Delta T = 10$ K in the temperature field can induce changes of approximately 50 per cent and 20 per cent in the strengths of the H₂O lines. Therefore, the D/H ratio derived from the measurements may be systematically under or overestimated if the temperature field is not accurately represented in the reference atmosphere. This dependence of the absorption line on the temperature field likely represents the major source of uncertainty in the analysis of these measurements.

In the case of the $H_2^{18}O$ isotopologue, the selection of suitable spectral windows within the measured spectral range is more challenging. The strongest absorption line of $H_2^{18}O$ occurs at 1389.91 $cm^{-1}$ and its optical depth is similar to the transitions of HDO and H₂O. However, this absorption line completely overlaps with absorption by carbon dioxide within the spectral resolution of SOFIA/EXES (see Figs 3 and 4). Since the overlapping absorption of $H_2^{18}O$ and CO₂ can lead to degeneracies in the spectral fitting that can impact the accuracy and reliability of the retrieval, we use instead another absorption line of $H_2^{18}O$, free of telluric and CO₂ contamination, located close to 1383.33 $cm^{-1}$. Unfortunately, the absorption line at 1383.33 $cm^{-1}$, used for the derivation of the $^{18}O/^{16}O$ ratio in this study, is located at the edge of the detector frame and was not measured in the observations made in 2018 October and 2019 April. Therefore, we only report measurements of $^{18}O/^{16}O$ for the observations made in 2016 March and 2017 January.

### 3.3 Description of the retrieval scheme

One of the main challenges to model the spectra measured in each pixel of the SOFIA/EXES images is the reconstruction of the field of view (FOV): Due to the relatively poor spatial resolution of the measurements, spectra at different locations all over the Mars disc have to be calculated and convolved with the PSF in each iteration to accurately model the measured spectrum in each pixel of the detector.

The reconstruction of the FOV thus makes the analysis of the spectra computationally expensive. In order to overcome this problem, and given that most of the parameters in our atmospheric model are fixed, we modify the forward model described in Section 3.1 to make use of tabulated look-up tables of high-resolution spectra that we interpolate over instead of performing the radiative transfer calculations in each iteration. In particular, we generate high-resolution spectra on a grid of locations all over the Mars disc by fixing all parameters in our reference atmosphere (i.e. pressure, temperature, and aerosol densities) but the total abundance of the water vapour isotopes. Specifically, the vertical profiles of water vapour from OpenMARS in each location are multiplied by scaling factors ranging from 0 to 3 with a step of 0.1, which is equivalent to generating spectra for different water column densities. Similarly, spectra for the minor isotopes are also generated using the water vapour profiles from OpenMARS and scaling them following an isotopic composition given by $D/H = 5$ VSMOW and $^{18}O/^{16}O = 1$ VSMOW. Once the look-up tables for the different isotopes have been generated, the spectra at an arbitrary column density and location can be calculated by interpolation.

Once the Mars high-resolution spectra are calculated at the different locations, the forward model follows the steps described in Section 3.1: calculation of the telluric transmittance, spatial convolution with the PSF, and spectral convolution with the ILS. Finally, we multiply the modelled spectrum with a second-order polynomial to accurately fit the continuum radiance. In summary, for each measurement, we fit the following model parameters:

(i) Scaling factor to the reference column abundance for any of the water isotopes. Note that one single scaling factor is applied to all locations on the planet.

(ii) Scaling factor for the line-of-sight density of H₂O, CO₂, and CH₄ in the Earth's atmosphere. Note that while this capability is
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Spectroscopy of the water vapour isotopes in the spectral windows used for the analysis of the SOFIA/EXES observations. The panels in the left and middle columns show the cross-sections of different absorption lines of H<sub>2</sub><sup>16</sup>O (top), HD<sup>16</sup>O (middle), H<sub>2</sub><sup>18</sup>O (bottom). The black dashed line in the mid-bottom panel shows the overlapping cross-section of CO<sub>2</sub>. The panels on the right show the integrated cross-sections as a function of temperature, which are a proxy for the dependence of the absorption line strengths with temperature.</figcaption></figure>

for this study avoids any telluric contamination and therefore this capability is not currently used.

* Coefficients of second order polynomial to fit the continuum. Note that this factor is implemented to correct for the absence of absolute calibration of the measurements, but the impact on the retrieved scaling factors for the water isotopes is minimal.

The model parameters are iteratively modified using an Optimal Estimation approach (Rodgers 2000; Irwin et al. 2008) to find the best fit between the modelled and measured spectra. Figure 5 shows several examples of the fits for the different spectral windows and observation campaigns. The developed retrieval scheme generally shows a good convergence and the residuals are well within the uncertainty of the measurements, which is crucial to properly characterise the random error of the retrieved parameters.

This retrieval scheme is applied to every spectral window and every pixel of the SOFIA/EXES images, which enables us to derive a scaling factor for each of those. It must be noted that the derived scaling factor is in reality a convolution of scaling factors from different locations of the Mars disk, or the convolution of different column densities from different locations. In order to derive the effective col-

umn density for each measurement we calculate the weighted average following

$$ \bar{m} = s \cdot \frac{\sum_i K_i m_i}{\sum_i K_i} \qquad (1) $$

where $\bar{m}$ represents the effective column density of a given isotope, $m_i$ is the column density of the isotope in each location of the planet, $s$ represents the retrieved scaling factor from a given measurement, and $K_i$ are the diagonal elements of the averaging kernel matrix, which can be used as a proxy for the relative sensitivity to the water abundance across the Martian disk.

For each pixel in the SOFIA/EXES images (i.e., every measurement), the averaging kernel matrix is calculated with the forward model following the Optimal Estimation formalism (Rodgers 2000), as explained in more detail in the Supplementary Material. The values of the diagonal elements of the averaging kernel matrix provide a useful diagnostic to characterise the relative sensitivity to the water abundance across the Martian disk, which can have important implications for our interpretation of the derived isotopic ratios from the SOFIA/EXES measurements. Figure 6 shows examples of the diagonal elements of the averaging kernel matrix calculated for three
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Spectroscopy of the water vapour isotopes in the spectral windows used for the analysis of the SOFIA/EXES observations. The panels in the left and middle columns show the cross-sections of different absorption lines of H<sub>2</sub><sup>16</sup>O (top panel), HD<sup>16</sup>O (middle panel), and H<sub>2</sub><sup>18</sup>O (bottom panel). The dashed line in the mid-bottom panel shows the overlapping cross-section of CO<sub>2</sub>. The panels on the right show the integrated cross-sections as a function of temperature, which are a proxy for the dependence of the absorption line strengths with temperature.</figcaption></figure>

implemented in our forward model, the selection of spectral windows for this study avoids any telluric contamination and therefore this capability is not currently used.

(iii) Coefficients of second-order polynomial to fit the continuum. Note that this factor is implemented to correct for the absence of absolute calibration of the measurements, but the impact on the retrieved scaling factors for the water isotopes is minimal.

The model parameters are iteratively modified using an Optimal Estimation approach (Rodgers 2000; Irwin et al. 2008) to find the best fit between the modelled and measured spectra. Fig. 5 shows several examples of the fits for the different spectral windows and observation campaigns. The developed retrieval scheme generally shows a good convergence and the residuals are well within the uncertainty of the measurements, which is crucial to properly characterize the random error of the retrieved parameters.

This retrieval scheme is applied to every spectral window and every pixel of the SOFIA/EXES images, which enables us to derive a scaling factor for each of those. It must be noted that the derived scaling factor is in reality a convolution of scaling factors from different locations of the Mars disc, or the convolution of different column densities from different locations. In order to derive the effective column density for each measurement we calculate the

weighted average following

$$ \bar{m} = s \cdot \frac{\sum_{i} K_{i} m_{i}}{\sum_{i} K_{i}} \qquad (1) $$

where $\bar{m}$ represents the effective column density of a given isotope, $m_i$ is the column density of the isotope in each location of the planet, s represents the retrieved scaling factor from a given measurement, and $K_i$ are the diagonal elements of the averaging kernel matrix, which can be used as a proxy for the relative sensitivity to the water abundance across the Martian disc.

For each pixel in the SOFIA/EXES images (i.e. every measurement), the averaging kernel matrix is calculated with the forward model following the Optimal Estimation formalism (Rodgers 2000), as explained in more detail in the Supplementary Material. The values of the diagonal elements of the averaging kernel matrix provide a useful diagnostic to characterize the relative sensitivity to the water abundance across the Martian disc, which can have important implications for our interpretation of the derived isotopic ratios from the SOFIA/EXES measurements. Fig. 6 shows examples of the diagonal elements of the averaging kernel matrix calculated for three measurements during the observations made in 2018 October. On one hand, the calculated values for the measurements pointing close to the subsolar point (i.e. warmer temperatures; see panels d and e in Fig. 6) closely represent the projection of the PSF on the
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Example of spectral fits to the SOFIA/EXES data. The different rows of panels represent examples of the spectral fits in the observations made in March 2016 and January 2017. The different columns represent the spectral fits for each of the spectral windows selected to derive the abundance of the different water isotopes. In the upper panels of each row, the black dots are the SOFIA/EXES measurements, while the coloured lines represent the best fit to the data, with the colours representing the isotopologue the absorption lines correspond to. Finally, the lower panels show the residuals between the best fit and the data (coloured lines), while the black dashed lines represent the uncertainties of the measurement.</figcaption></figure>

measurements during the observations made in October 2018. On one hand, the calculated values for the measurements pointing close to the subsolar point (i.e., warmer temperatures; see panels d and e in Figure 6) closely represent the projection of the PSF on the Martian disk, with the peak values of the diagonal elements of the averaging kernel matrix coinciding with the centre of the FOV and quickly decreasing with increasing distance from this point. On the other hand, the values for the measurements pointing away from the subsolar point (i.e., colder temperatures; see panel f in Figure 6) show a more complicated behaviour: the sensitivity to the water abundance peaks at the centre of the FOV but the distribution is much broader, showing substantial sensitivity to the water abundance from the subsolar region. This behaviour occurs because while the values of the PSF away from the centre of the FOV are small, the signal contribution from the subsolar region is much greater, with the sensitivity to the water abundance being a trade-off between the two.

After obtaining the effective column densities from each measurement, we enhance the precision of the derived quantities by averaging the retrieved values with uncertainties lower than 25% into 0.15"-wide bins. Figure 7 shows a comparison of the derived H₂O and HDO effective column densities from each spectral window. In particular, while the HDO densities derived from the two spectral windows are on average consistent with each other, we find that those in H₂O are systematically different. We interpret this difference to be arising from a small departure in the assumed temperatures in our reference atmosphere, which creates a systematic bias on the retrieved column densities due to the different sensitivity of the H₂O absorption lines to the temperature field (see Figure 4). On the other hand, given the low sensitivity of the HDO lines to the temperature, a small bias in the assumed temperatures does not induce any large systematic errors in the derived HDO column densities.

In order to correct for such a temperature bias, we perform the retrievals several times by shifting the surface and atmospheric temperatures by constant offsets with respect to the values extracted from the OpenMARS dataset. Panels a and b in Figure 8 show that, as could be expected from the sensitivity of the absorption lines in Figure 4, the retrieved water column densities from spectral window #1 are more sensitive to the reference atmospheric temperature than those from spectral window #2. Theoretically, the correct atmospheric temperature will be such that the retrieved column densities of H₂O from both spectral windows are equal. Therefore, we can estimate the bias in our reference atmosphere $\Delta T_b$ by finding the temperature offset at which both spectral windows agree (see panels c and d in Figure 8). The calculated distribution of $\Delta T_b$ is then spatially smoothed by fitting a polynomial function along the slit height direction for each of the measurement scans. Finally, we calculate the corrected column densities of the water isotopologues and derive the D/H and ¹⁸O/¹⁶O isotopic ratios. A detailed description of this methodology, as well as the impact of this correction in the derived isotopic ratios, is provided in the Supplementary Material.

We consider this approach to provide a reasonable estimation of potential temperature biases in our reference atmosphere that can affect the derivation of the isotopic ratios. Indeed, we find that this correction removes most of the unexpected variations in the retrieved isotopic ratios that would otherwise be difficult to reconcile with our current understanding of isotopic fractionation in the Martian atmosphere. However, while we consider the results to be generally satisfactory, it must also be noted that there are limitations to the applicability of this correction. In particular, since each measurement is analysed independently, we need to consider a constant temperature bias for all altitude levels and across the whole disk for the correction of the measurement. While this assumption is reasonable when
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example of spectral fits to the SOFIA/EXES data. The different rows of panels represent examples of the spectral fits in the observations made in 2016 March and 2017 January. The different columns represent the spectral fits for each of the spectral windows selected to derive the abundance of the different water isotopes. In the upper panels of each row, the dots are the SOFIA/EXES measurements, while the coloured lines represent the best fit to the data, with the colours representing the isotopologue the absorption lines correspond to. Finally, the lower panels show the residuals between the best fit and the data (coloured lines), while the dashed lines represent the uncertainties of the measurement.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Sensitivity to the water abundance across the Martian disc for the observations performed in 2018 October. The top panels represent the brightness distribution measured by SOFIA/EXES (a), as well as the surface temperature (b), and water vapour column densities (c) extracted from the OpenMARS data base. The bottom panels show the diagonal elements of the averaging kernel matrix for the three measurements indicated with crosses in panel a, which can be used as a proxy for the sensitivity to the water vapour abundance across the Martian disc. The sensitivity maps in panels d–f correspond to the measurements indicated with crosses in panel a from south to north, respectively.</figcaption></figure>

Martian disc, with the peak values of the diagonal elements of the averaging kernel matrix coinciding with the centre of the FOV and quickly decreasing with increasing distance from this point. On the other hand, the values for the measurements pointing away from the subsolar point (i.e. colder temperatures; see panel f in Fig. 6) show a more complicated behaviour: The sensitivity to the water abundance peaks at the centre of the FOV but the distribution is much broader, showing substantial sensitivity to the water abundance from the subsolar region. This behaviour occurs because while the values of the PSF away from the centre of the FOV are small, the signal contribution from the subsolar region is much greater, with the sensitivity to the water abundance being a trade-off between the two.

After obtaining the effective column densities from each measurement, we enhance the precision of the derived quantities by averaging the retrieved values with uncertainties lower than 25 per cent into 0.15-arcsec-wide bins. Fig. 7 shows a comparison of the derived H₂O and HDO effective column densities from each spectral window. In particular, while the HDO densities derived from the two spectral windows are on average consistent with each other, we find that those in H₂O are systematically different. We interpret this difference to be arising from a small departure in the assumed temperatures in our reference atmosphere, which creates a systematic bias on the retrieved column densities due to the different sensitivity of the H₂O absorption lines to the temperature field (see Fig. 4). On the other hand, given the low sensitivity of the HDO lines to the temperature, a small bias in the assumed temperatures does not induce any large systematic errors in the derived HDO column densities.

In order to correct for such a temperature bias, we perform the retrievals several times by shifting the surface and atmospheric temperatures by constant offsets with respect to the values extracted from the OpenMARS data set. Panels a and b in Fig. 8 show that, as
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Sensitivity to the water abundance across the Martian disk for the observations performed in October 2018. The top panels represent the brightness distribution measured by SOFIA/EXES (a), as well as the surface temperature (b) and water vapour column densities (c) extracted from the OpenMARS database. The bottom panels show the diagonal elements of the averaging kernel matrix for the three measurements indicated with black crosses in panel a, which can be used as a proxy for the sensitivity to the water vapour abundance across the Martian disk. The sensitivity maps in panels d-f correspond to the measurements indicated with black crosses in panel a from south to north, respectively.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Example of the retrieved H₂O and HDO column densities derived from the measurements made in October 2018. The retrieval scheme allows us to derive the effective column density of H₂O and HDO (top panels) and their uncertainties (bottom panels) from each spectral window. While the column densities of HDO derived from spectral windows #3 and #4 are consistent with each other, we find those in H₂O derived from spectral window #1 to be systematically lower than those from #2.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Calculation of the temperature bias $\Delta T_b$ from the retrievals of the H₂O column density in the October 2018 observations. The top panels (a,b) show the values of the H₂O column densities retrieved from the two spectral windows when varying the temperatures in the reference atmosphere with respect to those extracted from the OpenMARS dataset. Panel c shows the difference between the retrieved column densities from both spectral windows at the different temperature offsets. Panel d shows the calculated temperature bias $\Delta T_b$ in our reference atmosphere.</figcaption></figure>

analysing the measurements close to the subsolar region, where the sensitivity to the water abundance is confined to a narrow region on the disk (and the temperature bias in our reference atmosphere is likely similar), this can be an oversimplification for the measurements further away from the subsolar region, where the sensitivity to the water abundance is spread over a wide region on the disk (and the temperature bias might be very different in different regions of the Mars disk) (see Figure 6).

## 4 RESULTS

Using the retrieval approach explained in the previous section, we report the abundances of H₂O, HDO, H₂¹⁸O and the D/H and ¹⁸O/¹⁶O isotopic ratios in the atmosphere of Mars for the four SOFIA/EXES observations selected in our study (see Figure 9). In the following sections, we discuss the results of the measurements in more detail, focusing on describing the observed latitudinal distribution and the seasonal variability. In addition, in order to better evaluate the implications of our measurements to the current understanding of isotopic fractionation in the atmosphere of Mars, we compare the results of the measurements with the simulations of the Martian D/H cycle modelled using a Global Climate Model and reported in Vals et al. (2022) and Rossi et al. (2022).

### 4.1 Latitudinal distribution of water vapour and its isotopic ratios

The seasonal condensation/sublimation cycles of the polar ice caps on Mars produce substantial localised variations in the water column abundance and its isotopic composition. In particular, the lower saturation vapour pressure of HDO favours its condensation over H₂O, which induces a decrease of the D/H ratio in water vapour when condensation occurs (Merlivat & Nief 1967). Seasonally, this effect can produce approximately a two-fold decrease of D/H in the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Example of the retrieved H₂O and HDO column densities derived from the measurements made in 2018 October. The retrieval scheme allows us to derive the effective column density of H₂O and HDO (top panels) and their uncertainties (bottom panels) from each spectral window. While the column densities of HDO derived from spectral windows #3 and #4 are consistent with each other, we find those in H₂O derived from spectral window #1 to be systematically lower than those from #2.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Calculation of the temperature bias $\Delta T_b$ from the retrievals of the H₂O column density in the 2018 October observations. The top panels (a and b) show the values of the H₂O column densities retrieved from the two spectral windows when varying the temperatures in the reference atmosphere with respect to those extracted from the OpenMARS data set. Panel c shows the difference between the retrieved column densities from both spectral windows at the different temperature offsets. Panel d shows the calculated temperature bias $\Delta T_b$ in our reference atmosphere.</figcaption></figure>

could be expected from the sensitivity of the absorption lines in Fig. 4, the retrieved water column densities from spectral window #1 are more sensitive to the reference atmospheric temperature than those from spectral window #2. Theoretically, the correct atmospheric temperature will be such that the retrieved column densities of H₂O from both spectral windows are equal. Therefore, we can estimate the bias in our reference atmosphere $\Delta T_b$ by finding the temperature offset at which both spectral windows agree (see panels c and d in Fig. 8). The calculated distribution of $\Delta T_b$ is then spatially smoothed by fitting a polynomial function along the slit height direction for each of the measurement scans. Finally, we calculate the corrected column densities of the water isotopologues and derive the D/H and ¹⁸O/¹⁶O isotopic ratios. A detailed description of this methodology,

as well as the impact of this correction in the derived isotopic ratios, is provided in the Supplementary Material.

We consider this approach to provide a reasonable estimation of potential temperature biases in our reference atmosphere that can affect the derivation of the isotopic ratios. Indeed, we find that this correction removes most of the unexpected variations in the retrieved isotopic ratios that would otherwise be difficult to reconcile with our current understanding of isotopic fractionation in the Martian atmosphere. However, while we consider the results to be generally satisfactory, it must also be noted that there are limitations to the applicability of this correction. In particular, since each measurement is analysed independently, we need to consider a constant temperature bias for all altitude levels and across the whole disc for the correction of the measurement. While this assumption is reasonable when analysing the measurements close to the subsolar region, where the sensitivity to the water abundance is confined to a narrow region on the disc (and the temperature bias in our reference atmosphere is likely similar), this can be an oversimplification for the measurements further away from the subsolar region, where the sensitivity to the water abundance is spread over a wide region on the disc (and the temperature bias might be very different in different regions of the Mars disc) (see Fig. 6).

# 4. RESULTS

Using the retrieval approach explained in the previous section, we report the abundances of H₂O, HDO, H₂¹⁸O, and the D/H and ¹⁸O/¹⁶O isotopic ratios in the atmosphere of Mars for the four SOFIA/EXES observations selected in our study (see Fig. 9). In the following sections, we discuss the results of the measurements in more detail, focusing on describing the observed latitudinal distribution and the seasonal variability. In addition, in order to better evaluate the implications of our measurements to the current understanding of isotopic fractionation in the atmosphere of Mars, we compare the results of the measurements with the simulations of the Martian D/H cycle modelled using a Global Climate Model and reported in Vals et al. (2022) and Rossi et al. (2022).

## 4.1 Latitudinal distribution of water vapour and its isotopic ratios

The seasonal condensation/sublimation cycles of the polar ice caps on Mars produce substantial localized variations in the water column abundance and its isotopic composition. In particular, the lower saturation vapour pressure of HDO favours its condensation over H₂O, which induces a decrease of the D/H ratio in water vapour when condensation occurs (Merlivat & Nief 1967). Seasonally, this effect can produce approximately a two-fold decrease of D/H in the winter hemispheres, where most of the water vapour condenses on to the seasonal water ice frosts (Montmessin et al. 2005). In the case of ¹⁸O/¹⁶O, there are no quantitative estimations for the amount of fractionation on Mars, but based on the relation between the variations of D/H (i.e. $\delta$D) and ¹⁸O/¹⁶O (i.e. $\delta$¹⁸O) on Earth, we may expect $\delta$D ~ 8 × $\delta$¹⁸O (Craig 1961).

Fig. 10 shows the latitudinal distribution of the water vapour column density and the associated D/H and ¹⁸O/¹⁶O isotopic ratios derived from the SOFIA/EXES measurements. The variations of the water column abundance follow the expectations from the current knowledge of the water cycle (Montmessin et al. 2004; Navarro et al. 2014). The annual peak of the water column abundance is found during the summer season in the Northern hemisphere (i.e. 2016 March, $L_s$ = 127°), coinciding with the sublimation of the
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Summary of the data products derived from the SOFIA/EXES measurements, plotted in order of increasing $L_s$. The panels in this figure show the retrieved column density of H₂O (top) and the associated D/H (middle) and ¹⁸O/¹⁶O (bottom) isotopic ratios in the atmosphere of Mars derived from the SOFIA/EXES measurements made in different campaigns, as indicated by the title in each column. The black cross in each of the panels represents the location of the subsolar point in each of the observations.</figcaption></figure>

winter hemispheres, where most of the water vapour condenses onto the seasonal water ice frosts (Montmessin et al. 2005). In the case of ¹⁸O/¹⁶O, there are no quantitative estimations for the amount of fractionation on Mars, but based on the relation between the variations of D/H (i.e., δD) and ¹⁸O/¹⁶O (i.e., δ¹⁸O) on Earth, we may expect δD ~ 8 × δ¹⁸O (Craig 1961).

Figure 10 shows the latitudinal distribution of the water vapour column density and the associated D/H and ¹⁸O/¹⁶O isotopic ratios derived from the SOFIA/EXES measurements. The variations of the water column abundance follow the expectations from the current knowledge of the water cycle (Montmessin et al. 2004; Navarro et al. 2014). The annual peak of the water column abundance is found during the summer season in the northern hemisphere (i.e. March 2016, $L_s = 127°$), coinciding with the sublimation of the northern polar cap. On the other hand, the water column density peaks in the southern hemisphere during its summer season (i.e., October 2018 and January 2017, $L_s = 272°$ and $L_s = 305°$), although the peak in column abundance is lower than that found in the northern summer, in agreement with the observations from previous orbiters (Trokhimovskiy et al. 2015; Smith et al. 2018; Crismani et al. 2021; Knutsen et al. 2022). Close to the spring equinox period (i.e., April 2019, $L_s = 15°$), the latitude distribution is relatively constant in the equatorial region, showing a slight increase towards the northern hemisphere, with the values of the column abundance being overall lower than during the solstice periods. It must be noted that the SOFIA/EXES measurements are smoothed with the relatively large PSF but we find overall a good agreement with the values of the water column abundance retrieved from the Mars orbiters.

As opposed to the latitudinal variations observed in the water vapour column abundances, we find the latitudinal distribution of the D/H ratio to be almost constant in all observations. The only exception to this is found in the observations in January 2017, where the retrievals suggest a decrease of D/H from approximately 5 to 4 VSMOW from the equator to 50°N. However, this decrease in D/H coincides with a strong and rapidly changing increase of the calculated temperature error in the observations (see Supplementary Figure 5) and we suspect this decrease might be produced by a systematic bias caused by the assumption of a constant temperature offset for the whole disk. The measurements made in March 2016 and October 2018, which show essentially constant latitude trends, were made using the mapping mode and are likely the most accurate observations, given that we can average many points to reduce the random error. Finally, the observations made in April 2019 appear to be noisier than the rest, mainly because fewer measurements were performed, which directly impacts the effective SNR. Within the observed variability in each latitude bin in this campaign, we do not observe evidence for latitudinal variations above the level of the uncertainty.

In order to investigate whether the absence of latitudinal variability is in agreement with the expectations from condensation-induced fractionation, we compare the observed latitudinal trends with the GCM simulations (Vals et al. 2022; Rossi et al. 2022) (see Figure 11). However, in order to perform a one-to-one comparison between both datasets, we need to smooth the GCM predictions in a similar way as SOFIA/EXES smooths the retrieved column densities. Therefore, we calculate the effective D/H ratio from the GCM following
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Summary of the data products derived from the SOFIA/EXES measurements, plotted in order of increasing $L_s$. The panels in this figure show the retrieved column density of H₂O (top panel) and the associated D/H (middle panel) and ¹⁸O/¹⁶O (bottom panel) isotopic ratios in the atmosphere of Mars derived from the SOFIA/EXES measurements made in different campaigns, as indicated by the title in each column. The cross in each of the panels represents the location of the subsolar point in each of the observations.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Latitudinal distribution of the H₂O column density and the D/H and ¹⁸O/¹⁶O isotopic ratios from the analysis of the SOFIA/EXES observations. The dots in these panels represent the values derived from the measurements, as well as their associated uncertainties. The solid lines represent the latitudinal trends of the D/H (middle) and ¹⁸O/¹⁶O (bottom) isotope ratios, where the error bars represent the standard deviation of the measurements in each latitude band.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Latitudinal distribution of the H₂O column density and the D/H and ¹⁸O/¹⁶O isotopic ratios from the analysis of the SOFIA/EXES observations. The dots in these panels represent the values derived from the measurements, as well as their associated uncertainties. The red and green lines represent the latitudinal trends of the D/H and ¹⁸O/¹⁶O isotope ratios, where the error bars represent the standard deviation of the measurements in each latitude band.</figcaption></figure>

equation 1 (see red dots in Figure 11). When performing this convection, we observe that most of the variations predicted by the GCM essentially disappear. Therefore, we conclude that the spatial resolution of SOFIA/EXES smooths most variations predicted from the GCM and therefore impedes us to perform a quantitative evaluation on the effect of condensation-induced fractionation.

The analysis of the SOFIA/EXES measurements in March 2016 and January 2017 by Encrenaz et al. (2018) also revealed smooth variations of the D/H ratio along the Mars disk. Particularly, they found an absence of strong latitudinal variations for the observations in March 2016 ($L_S = 127°$), which is consistent with our analysis of the data. In their analysis of the observations in January 2017 ($L_S = 305°$) they also reported smooth variations in D/H, but noted the presence of a slight increase in D/H with latitude, with values of D/H~4.5 near 90°S and values of D/H~5.0 near 45°N. Based on our analysis of the data, we believe the source of this trend might be caused by inaccuracies of the assumed temperature in the reference atmosphere: we observe a similar increase of D/H with latitude that disappears when correcting for the temperature bias (see Supplementary Figure 5). The difference between the assumed temperature fields in both analysis of the data might also explain the difference in the disk-averaged results (D/H ~4 and ~4.5 (Encrenaz et al. 2018) compared to D/H=5.16 and D/H=5.04 in this study for the observations in March 2016 and January 2017, respectively). In that aspect, we consider that the method used in this study, combining the varying temperature field from the OpenMARS database and the convolution with the PSF of SOFIA/EXES, likely provides a more accurate estimation of the atmospheric structure and minimises the systematic uncertainties in the derived D/H.

Finally, the latitudinal distribution of ¹⁸O/¹⁶O was also derived from the observations made in March 2016 and January 2017. Over-

all, we do not find any significant latitudinal variations above the level of uncertainty in the observations. There are no quantitative calculations for the amount of expected fractionation in the ¹⁸O/¹⁶O isotopic ratio in water vapour in the atmosphere of Mars. Nevertheless, considering that the spatial resolution of the instrument smooths most of the variations in D/H, and that the variations in ¹⁸O/¹⁶O are expected to be smaller than those in D/H, it is probably safe to assume that no latitudinal variations of the ¹⁸O/¹⁶O ratio are expected within the level of uncertainties of the data.

Interestingly, the observations made in January 2017 do not show a decrease of the ¹⁸O/¹⁶O isotopic ratio for latitudes higher than 0°N, as it is observed with D/H. Nevertheless, we still suspect that this feature in the D/H ratio might be caused by a bias in the temperature field: the temperature dependence of the H₂O and H₂¹⁸O absorption lines is similar, and a bias in the temperature field can average out when computing the isotopic ratio (see Figure 4). In the next sections, we exclude the D/H measurements made in January 2017 in the northern hemisphere from the rest of the analysis.

## 4.2 Seasonal variations of the D/H and ¹⁸O/¹⁶O ratios

The analysis of the latitudinal distributions in the previous section shows that the SOFIA/EXES measurements are, to a large extent, insensitive to localised variations in the isotopic ratios, but instead provide a good measure of the global D/H and ¹⁸O/¹⁶O ratios in the atmospheric reservoir. Using the simulations of D/H from a Global Climate Model accounting for condensation-induced fractionation, Montmessin et al. (2005) suggested that while local values of the D/H ratio might follow strong variations, the global inventory of D/H should be much more constant throughout a Martian year. In particular, their simulations suggest that the planetary D/H ratio varies
northern polar cap. On the other hand, the water column density peaks in the Southern hemisphere during its summer season (i.e. 2018 October and 2017 January, $L_s = 272^\circ$ and $L_s = 305^\circ$), although the peak in column abundance is lower than that found in the northern summer, in agreement with the observations from previous orbiters (Trokhimovskiy et al. 2015; Smith et al. 2018; Crismani et al. 2021; Knutsen et al. 2022). Close to the spring equinox period (i.e. 2019 April, $L_s = 15^\circ$), the latitude distribution is relatively constant in the equatorial region, showing a slight increase towards the Northern hemisphere, with the values of the column abundance being overall lower than during the solstice periods. It must be noted that the SOFIA/EXES measurements are smoothed with the relatively large PSF but we find overall a good agreement with the values of the water column abundance retrieved from the Mars orbiters.

As opposed to the latitudinal variations observed in the water vapour column abundances, we find the latitudinal distribution of the D/H ratio to be almost constant in all observations. The only exception to this is found in the observations in 2017 January, where the retrievals suggest a decrease of D/H from approximately 5 to 4 VSMOW from the equator to $50^\circ$N. However, this decrease in D/H coincides with a strong and rapidly changing increase of the calculated temperature error in the observations (see Supplementary Fig. 5) and we suspect this decrease might be produced by a systematic bias caused by the assumption of a constant temperature offset for the whole disc. The measurements made in 2016 March and 2018 October, which show essentially constant latitude trends, were made using the mapping mode and are likely the most accurate observations, given that we can average many points to reduce the random error. Finally, the observations made in 2019 April appear to be noisier than the rest, mainly because fewer measurements were performed, which directly impacts the effective signal-to-noise ratio (SNR). Within the observed variability in each latitude bin in this campaign, we do not observe evidence for latitudinal variations above the level of the uncertainty.

In order to investigate whether the absence of latitudinal variability is in agreement with the expectations from condensation-induced fractionation, we compare the observed latitudinal trends with the GCM simulations (Rossi et al. 2022; Vals et al. 2022, see Fig. 11). However, in order to perform a one-to-one comparison between both data sets, we need to smooth the GCM predictions in a similar way as SOFIA/EXES smooths the retrieved column densities. Therefore, we calculate the effective D/H ratio from the GCM following equation (1) (see red dots in Fig. 11). When performing this convolution, we observe that most of the variations predicted by the GCM essentially disappear. Therefore, we conclude that the spatial resolution of SOFIA/EXES smooths most variations predicted from the GCM and therefore impedes us to perform a quantitative evaluation on the effect of condensation-induced fractionation.

The analysis of the SOFIA/EXES measurements in 2016 March and 2017 January by Encrenaz et al. (2018) also revealed smooth variations of the D/H ratio along the Mars disc. Particularly, they found an absence of strong latitudinal variations for the observations in March 2016 ($L_s = 127^\circ$), which is consistent with our analysis of the data. In their analysis of the observations in 2017 January ($L_s = 305^\circ$) they also reported smooth variations in D/H, but noted the presence of a slight increase in D/H with latitude, with values of D/H ~ 4.5 near $90^\circ$S and values of D/H ~ 5.0 near $45^\circ$N. Based on our analysis of the data, we believe the source of this trend might be caused by inaccuracies of the assumed temperature in the reference atmosphere: We observe a similar increase of D/H with latitude that disappears when correcting for the temperature bias (see Supplementary Fig. 5). The difference between the assumed

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Latitudinal distribution of the D/H ratio modelled using a Global Climate Model accounting for condensation-induced fractionation and comparison with the retrieved D/H distribution measured with SOFIA/EXES. The grey dots show the values of the D/H ratio predicted by the GCM simulations, evidencing substantial variability of D/H in the atmosphere. When the model predictions are convolved with the calculated averaging kernels for SOFIA/EXES (red dots), this variability disappears. The solid lines represent the latitudinal trends derived from the SOFIA/EXES measurements, where the error bars represent the standard deviation of the measurements within each latitude band.</figcaption></figure>

temperature fields in both analysis of the data might also explain the difference in the disc-averaged results [D/H ~ 4 and ~4.5 (Encrenaz et al. 2018) compared with D/H = 5.16 and 5.04 in this study for the observations in 2016 March and 2017 January, respectively). In that aspect, we consider that the method used in this study, combining the varying temperature field from the OpenMARS data base and the convolution with the PSF of SOFIA/EXES, likely provides a more accurate estimation of the atmospheric structure and minimizes the systematic uncertainties in the derived D/H.

Finally, the latitudinal distribution of $^{18}$O/$^{16}$O was also derived from the observations made in 2016 March and 2017 January. Overall, we do not find any significant latitudinal variations above the level of uncertainty in the observations. There are no quantitative calculations for the amount of expected fractionation in the $^{18}$O/$^{16}$O isotopic ratio in water vapour in the atmosphere of Mars. Nevertheless, considering that the spatial resolution of the instrument smooths most of the variations in D/H, and that the variations in $^{18}$O/$^{16}$O are expected to be smaller than those in D/H, it is probably safe to assume that no latitudinal variations of the $^{18}$O/$^{16}$O ratio are expected within the level of uncertainties of the data.

Interestingly, the observations made in 2017 January do not show a decrease of the $^{18}$O/$^{16}$O isotopic ratio for latitudes higher than $0^\circ$N, as it is observed with D/H. Nevertheless, we still suspect that this feature in the D/H ratio might be caused by a bias in the temperature field: The temperature dependence of the H$_2$O and H$_2^{18}$O absorption lines is similar, and a bias in the temperature field can average out when computing the isotopic ratio (see Fig. 4). In the next sections, we exclude the D/H measurements made in 2017 January in the Northern hemisphere from the rest of the analysis.
<figure><img src="image_11.png" /></figure>

**Figure 11.** Latitudinal distribution of the D/H ratio modelled using a Global Climate Model accounting for condensation-induced fractionation and comparison with the retrieved D/H distribution measured with SOFIA/EXES. The grey dots show the values of the D/H ratio predicted by the GCM simulations, evidencing substantial variability of D/H in the atmosphere. When the model predictions are convolved with the calculated averaging kernels for SOFIA/EXES (red dots), this variability disappears. The black lines represent the latitudinal trends derived from the SOFIA/EXES measurements, where the error bars represent the standard deviation of the measurements within each latitude band.

within 2%, finding maximum peaks during the summer/winter seasons and minima during spring/autumn. In this section, we update the calculations of Montmessin et al. (2005) with the more recent version of the GCM presented by Vals et al. (2022) and Rossi et al. (2022), and compare them with the temporal variability found in the SOFIA/EXES observations.

Table 2 and Figure 12 show the weighted averaged values of the D/H and ¹⁸O/¹⁶O isotopic ratios in each of the SOFIA/EXES observations, as well as the values predicted from the GCM simulations. The GCM predicts the planetary D/H ratio to peak at $L_S = 135^\circ$ with a value of D/H = 4.84, coinciding with the period when the water content in the atmospheric reservoir peaks. The planetary D/H ratio has another local maximum (D/H = 4.81) around $L_S = 300^\circ$ during the summer season in the southern hemisphere. In agreement with Montmessin et al. (2005), the minima are found during the spring/autumn seasons, where the planetary D/H ratio decreases to D/H = 4.60 around $L_S = 40^\circ$ and D/H = 4.62 around $L_S = 190^\circ$. Therefore, the GCM simulations suggest seasonal variations in the D/H ratio of approximately 5%.

The temporal variations in the global D/H ratios derived from the measurements are qualitatively similar to the predictions from the GCM. We find the maximum D/H ratio in our observations during the northern summer (D/H = 5.16 ± 0.23 VSMOW at $L_S = 127^\circ$). The values derived for the southern summer are approximately 3% smaller, although this is similar to the derived uncertainties. Finally, we find the minimum value in the measurements made in April 2019 (D/H = 4.31 ± 0.38) during the spring season in the northern

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

**Figure 12.** Annual cycle of the D/H ratio in the atmosphere of Mars. Simulations of the D/H ratio in water vapour in a GCM (and scaled by a factor of 1.05) suggest that it is subject to strong localised variations induced by the condensation of water vapour on the polar ice caps and ice clouds (top panel). On the other hand, the simulations suggest the planetary D/H ratio only varies up to 5% throughout a full Martian year (red line in bottom panel). The measurements of D/H from SOFIA/EXES (black dots in bottom panel) evidence the absence of strong variability in the global inventory of D/H, with the exception of the measurements made in $L_S = 15^\circ$, which depict a substantially lower value than the rest (~15%).

hemisphere. While the GCM predicts an annual minimum of the global D/H ratio around this season, the variation derived from the measurements is substantially larger than predicted by the GCM (see Table 2). While the uncertainties derived for these observations are somewhat larger, as discussed in the previous section, we do not have any reasons to disprove these values. Interestingly, Krasnopolsky (2015) also reported significantly lower values than the GCM during the spring season ($L_S = 42^\circ$), which might depict inaccuracies of the model during this season.

Krasnopolsky (2015, 2021) measured the seasonal and latitudinal variations of the D/H ratio during several observing campaigns with CSHELL at NASA Infrared Telescope Facility (IRTF), covering approximately the first half of a Martian year ($L_S = 20-160^\circ$). They report a global-mean D/H ratio of 4.6 ± 0.7 VSMOW during this period, consistent with our analysis of the SOFIA/EXES measurements. One of their observations was performed on 28 March 2016, three days after the observations made by SOFIA/EXES during northern summer ($L_S = 127^\circ$). Their observations reveal an essentially constant D/H ratio in the northern hemisphere with D/H~5, similar to our measured value during this campaign (D/H = 5.16 ± 0.23 VSMOW). In contrast to the GCM, Krasnopolsky (2015) find the D/H ratio to peak during aphelion ($L_S=71^\circ$) and to be generally more variable than predicted by the GCM, with a minimum during northern spring at $L_S=42^\circ$ with values ranging between ~2-4 VSMOW.

The measurements from both SOFIA/EXES and IRTF/CSHELL appear to reveal more temporal variability in D/H than predicted by the GCM, which might indicate the presence of missing fractionating processes apart from condensation-induced fractionation that are not currently captured in the model. Krasnopolsky (2021) estimated that photolysis-induced fractionation should only have an effect above
<table><caption>**Table 2.** Global isotopic ratios derived from the SOFIA/EXES observations. The D/H and ¹⁸O/¹⁶O ratios represent the weighted average for each observation, where the uncertainty is calculated as the standard deviation from the measurements. The values of δD correspond to the variations in the global isotopic ratios with respect to that at Ls = 127° (2016 March).</caption><thead><tr><th>Date (YYYY-MM-DD)</th><th>Time (UTC)</th><th>L<sub>S</sub> (°)</th><th>D/H (VSMOW)</th><th>δD (per cent)</th><th>D/H (GCM) (VSMOW)</th><th>δD (GCM) (per cent)</th><th>¹⁸O/¹⁶O (VSMOW)</th></tr></thead><tbody><tr><td>2016-03-24</td><td>12:06:05</td><td>127</td><td>5.16 ± 0.23</td><td>0.00</td><td>4.84</td><td>0.00</td><td>1.12 ± 0.07</td></tr><tr><td>2017-01-25</td><td>01:10:56</td><td>305</td><td>5.04 ± 0.18</td><td>-2.3</td><td>4.81</td><td>-0.62</td><td>1.05 ± 0.06</td></tr><tr><td>2018-10-20</td><td>02:58:17</td><td>272</td><td>4.95 ± 0.21</td><td>-4.1</td><td>4.76</td><td>-1.7</td><td>-</td></tr><tr><td>2019-04-24</td><td>04:18:19</td><td>15</td><td>4.31 ± 0.38</td><td>-16.5</td><td>4.66</td><td>-3.7</td><td>-</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Annual cycle of the D/H ratio in the atmosphere of Mars. Simulations of the D/H ratio in water vapour in a GCM (and scaled by a factor of 1.05) suggest that it is subject to strong localized variations induced by the condensation of water vapour on the polar ice caps and ice clouds (top panel). On the other hand, the simulations suggest the planetary D/H ratio only varies up to 5 per cent throughout a full Martian year (solid line in bottom panel). The measurements of D/H from SOFIA/EXES (dots in bottom panel) evidence the absence of strong variability in the global inventory of D/H, with the exception of the measurements made in Ls = 15°, which depict a substantially lower value than the rest (~15 per cent).</figcaption></figure>

## 4.2 Seasonal variations of the D/H and ¹⁸O/¹⁶O ratios

The analysis of the latitudinal distributions in the previous section shows that the SOFIA/EXES measurements are, to a large extent, insensitive to localized variations in the isotopic ratios, but instead provide a good measure of the global D/H and ¹⁸O/¹⁶O ratios in the atmospheric reservoir. Using the simulations of D/H from a Global Climate Model accounting for condensation-induced fractionation, Montmessin et al. (2005) suggested that while local values of the D/H ratio might follow strong variations, the global inventory of D/H should be much more constant throughout a Martian year. In particular, their simulations suggest that the planetary D/H ratio varies within 2 per cent, finding maximum peaks during the summer/winter seasons and minima during spring/autumn. In this section, we update the calculations of Montmessin et al. (2005) with the more recent version of the GCM presented by Vals et al. (2022) and Rossi et al. (2022), and compare them with the temporal variability found in the SOFIA/EXES observations.

Table 2 and Fig. 12 show the weighted averaged values of the D/H and ¹⁸O/¹⁶O isotopic ratios in each of the SOFIA/EXES observations,

as well as the values predicted from the GCM simulations. The GCM predicts the planetary D/H ratio to peak at Ls = 135° with a value of D/H = 4.84, coinciding with the period when the water content in the atmospheric reservoir peaks. The planetary D/H ratio has another local maximum (D/H = 4.81) around Ls = 300° during the summer season in the Southern hemisphere. In agreement with Montmessin et al. (2005), the minima are found during the spring/autumn seasons, where the planetary D/H ratio decreases to D/H = 4.60 around Ls = 40° and D/H = 4.62 around Ls = 190°. Therefore, the GCM simulations suggest seasonal variations in the D/H ratio of approximately 5 per cent.

The temporal variations in the global D/H ratios derived from the measurements are qualitatively similar to the predictions from the GCM. We find the maximum D/H ratio in our observations during the northern summer (D/H = 5.16 ± 0.23 VSMOW at Ls = 127°). The values derived for the southern summer are approximately 3 per cent smaller, although this is similar to the derived uncertainties. Finally, we find the minimum value in the measurements made in 2019 April (D/H = 4.31 ± 0.38) during the spring season in the Northern hemisphere. While the GCM predicts an annual minimum of the global D/H ratio around this season, the variation derived from the measurements is substantially larger than predicted by the GCM (see Table 2). While the uncertainties derived for these observations are somewhat larger, as discussed in the previous section, we do not have any reasons to disprove these values. Interestingly, Krasnopolsky (2015) also reported significantly lower values than the GCM during the spring season (Ls = 42°), which might depict inaccuracies of the model during this season.

Krasnopolsky (2015, 2021) measured the seasonal and latitudinal variations of the D/H ratio during several observing campaigns with CSHELL at NASA Infrared Telescope Facility (IRTF), covering approximately the first half of a Martian year (Ls = 20–160°). They report a global-mean D/H ratio of 4.6 ± 0.7 VSMOW during this period, consistent with our analysis of the SOFIA/EXES measurements. One of their observations was performed on 2016 March 28, 3 d after the observations made by SOFIA/EXES during northern summer (Ls = 127°). Their observations reveal an essentially constant D/H ratio in the Northern hemisphere with D/H ~ 5, similar to our measured value during this campaign (D/H = 5.16 ± 0.23 VSMOW). In contrast to the GCM, Krasnopolsky (2015) find the D/H ratio to peak during aphelion (Ls = 71°) and to be generally more variable than predicted by the GCM, with a minimum during northern spring at Ls = 42° with values ranging between ~2 and 4 VSMOW.

The measurements from both SOFIA/EXES and IRTF/CSHELL appear to reveal more temporal variability in D/H than predicted by the GCM, which might indicate the presence of missing fractionating processes apart from condensation-induced fractionation that are not currently captured in the model. Krasnopolsky (2021) estimated that photolysis-induced fractionation should only have an effect above ~50 km. This was later confirmed by Vals et al. (2022) and
<table><caption>**Table 2.** Global isotopic ratios derived from the SOFIA/EXES observations. The D/H and ¹⁸O/¹⁶O ratios represent the weighted average for each observation, where the uncertainty is calculated as the standard deviation from the measurements. The values of δD correspond to the variations in the global isotopic ratios with respect to that at Ls = 127° (March 2016).</caption><thead><tr><th>Date (YYYY-MM-DD)</th><th>Time (UTC)</th><th>L<sub>S</sub> (°)</th><th>D/H (VSMOW)</th><th>δD (%)</th><th>D/H (GCM) (VSMOW)</th><th>δD (GCM) (%)</th><th>¹⁸O/¹⁶O (VSMOW)</th></tr></thead><tbody><tr><td>2016-03-24</td><td>12:06:05</td><td>127</td><td>5.16 ± 0.23</td><td>0.00</td><td>4.84</td><td>0.00</td><td>1.12 ± 0.07</td></tr><tr><td>2017-01-25</td><td>01:10:56</td><td>305</td><td>5.04 ± 0.18</td><td>-2.3</td><td>4.81</td><td>-0.62</td><td>1.05 ± 0.06</td></tr><tr><td>2018-10-20</td><td>02:58:17</td><td>272</td><td>4.95 ± 0.21</td><td>-4.1</td><td>4.76</td><td>-1.7</td><td>-</td></tr><tr><td>2019-04-24</td><td>04:18:19</td><td>15</td><td>4.31 ± 0.38</td><td>-16.5</td><td>4.66</td><td>-3.7</td><td>-</td></tr></tbody></table>

~50 km. This was later confirmed by Vals et al. (2022) and Rossi et al. (2022), who introduced this effect in the most updated version of the GCM and showed that it has little impact on the horizontal distribution of D/H. One potential candidate for explaining the enhanced variability in the observations is the fractionation during the adsorption of water vapour on the Martian regolith. Hu (2019) showed that this process can produce variations in D/H of about ~15% near the surface and might account for the extra fractionation required to reconcile the variability in the measurements and models.

In the case of the ¹⁸O/¹⁶O, if we assume that the variations induced by condensation are proportional to those in D/H following the meteoric line (δD ~ 8 × δ¹⁸O), we might expect variations in the global ¹⁸O/¹⁶O ratio lower than 1%, which is well below the uncertainty of the measurements. The calculations of the averaged ¹⁸O/¹⁶O ratios from the observations in March 2016 and January 2017 indicate an enrichment of ¹⁸O/¹⁶O = 1.12 ± 0.07 VSMOW and ¹⁸O/¹⁶O = 1.05 ± 0.06 VSMOW, respectively (see Table 2). In order to investigate whether these differences are significant, we compute some statistical parameters on the observed distribution of measurements. Figure 13 shows a histogram with all the measured values of ¹⁸O/¹⁶O, as well as the relation of these with the calculated uncertainties of each measurement. It can be observed in Figure 13 that the vast majority of the points show variability within the uncertainties (1σ). The mean value of the observed distribution is μ = 1.09 VSMOW and the standard deviation is σ = 0.08 VSMOW. The standard deviation of the measurements is indeed very similar to the level of the calculated uncertainty, which suggests that the width of the distribution is a consequence of the measurement noise. Therefore, we conclude that the observed variations are consistent with a constant value of ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW.

## 4.3 The global inventory of D/H and ¹⁸O/¹⁶O in Mars water

In the previous sections, we explored the spatial and temporal variations of the D/H and ¹⁸O/¹⁶O isotopic ratios with respect to the model expectations from condensation-induced fractionation. In this section, we aim at further discussing the effects of condensation in our measurements and infer what the overall isotopic composition of the present-day exchangeable water reservoir on Mars is. These values carry the information about the long-term fractionation on Mars and are therefore crucial to better understand the evolution of the water reservoir throughout history.

The annual cycle of water vapour is currently dominated by the condensation and sublimation of the north polar cap, which hosts a vast amount of water ice that can exchange with the atmospheric reservoir (Montmessin et al. 2017). When exposed to sunlight during the spring/summer seasons, the amount of ice sublimed from the north polar cap controls the density of water in the atmosphere. Similarly, the D/H ratio released from the north polar cap controls

<figure><img src="image_13.png" /><figcaption>**Figure 13.** Statistical analysis of the ¹⁸O/¹⁶O measurements. The distribution of measured points (top panel) resembles a Gaussian distribution with mean and standard deviations of μ = 1.09 and σ = 0.08. When comparing the variations in ¹⁸O/¹⁶O with the uncertainties of the measurement (bottom panel), we observe that the observed variability is lower than the 1-σ uncertainties (black dashed line), with only a few exceptions that exceed this line.</figcaption></figure>

the D/H ratio of water in the atmospheric reservoir. Therefore, the isotopic composition characteristic of the surface/atmosphere water reservoir at present will be determined by the isotopic ratios in the north polar cap (Montmessin et al. 2005; Daerden et al. 2022).

Modelling the annual D/H cycle in a GCM, Montmessin et al. (2005) calculated that the D/H ratio of the atmospheric reservoir (~5 VSMOW) was about 15% less than that prescribed in the simulations in the north polar cap (5.6 VSMOW). Using the most updated version of the GCM (Rossi et al. 2022; Vals et al. 2022), we calculate that this dichotomy between the two reservoirs decreases to 5.5% (D/H prescribed in the north polar cap is 5 VSMOW; the annually averaged D/H in the atmospheric reservoir is 4.72 VSMOW). This value is also close to that reported by Daerden et al. (2022), where there is a difference of 5% between the value prescribed on the north polar cap and that in the atmosphere.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Statistical analysis of the ¹⁸O/¹⁶O measurements. The distribution of measured points (top panel) resembles a Gaussian distribution with mean and standard deviations of μ = 1.09 and σ = 0.08. When comparing the variations in ¹⁸O/¹⁶O with the uncertainties of the measurement (bottom panel), we observe that the observed variability is lower than the 1σ uncertainties (1σ dashed line), with only a few exceptions that exceed this line.</figcaption></figure>

Rossi et al. (2022), who introduced this effect in the most updated version of the GCM and showed that it has little impact on the horizontal distribution of D/H. One potential candidate for explaining the enhanced variability in the observations is the fractionation during the adsorption of water vapour on the Martian regolith. Hu (2019) showed that this process can produce variations in D/H of about ~15 per cent near the surface and might account for the extra fractionation required to reconcile the variability in the measurements and models.

In the case of the ¹⁸O/¹⁶O, if we assume that the variations induced by condensation are proportional to those in D/H following the meteoric line (δD ~ 8 × δ¹⁸O), we might expect variations in the global ¹⁸O/¹⁶O ratio lower than 1 per cent, which is well below the uncertainty of the measurements. The calculations of the averaged ¹⁸O/¹⁶O ratios from the observations in 2016 March and 2017 January indicate an enrichment of ¹⁸O/¹⁶O = 1.12 ± 0.07 VSMOW and ¹⁸O/¹⁶O = 1.05 ± 0.06 VSMOW, respectively (see Table 2). In order to investigate whether these differences are significant, we compute some statistical parameters on the observed distribution of measurements. Fig. 13 shows a histogram with all the measured values of ¹⁸O/¹⁶O, as well as the relation of these with the calculated uncertainties of each measurement. It can be observed in Fig. 13 that the vast majority of the points show variability within the uncertainties (1σ). The mean value of the observed distribution is μ = 1.09 VSMOW and the standard deviation is σ = 0.08 VSMOW. The standard deviation of the measurements is indeed very similar to the level of the calculated uncertainty, which suggests that the width of the distribution is a consequence of the measurement noise.

Therefore, we conclude that the observed variations are consistent with a constant value of ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW.

## 4.3 The global inventory of D/H and ¹⁸O/¹⁶O in Mars water

In the previous sections, we explored the spatial and temporal variations of the D/H and ¹⁸O/¹⁶O isotopic ratios with respect to the model expectations from condensation-induced fractionation. In this section, we aim at further discussing the effects of condensation in our measurements and infer what the overall isotopic composition of the present-day exchangeable water reservoir on Mars is. These values carry the information about the long-term fractionation on Mars and are therefore crucial to better understand the evolution of the water reservoir throughout history.

The annual cycle of water vapour is currently dominated by the condensation and sublimation of the north polar cap, which hosts a vast amount of water ice that can exchange with the atmospheric reservoir (Montmessin et al. 2017). When exposed to sunlight during the spring/summer seasons, the amount of ice sublimed from the north polar cap controls the density of water in the atmosphere. Similarly, the D/H ratio released from the north polar cap controls the D/H ratio of water in the atmospheric reservoir. Therefore, the isotopic composition characteristic of the surface/atmosphere water reservoir at present will be determined by the isotopic ratios in the north polar cap (Montmessin et al. 2005; Daerden et al. 2022).

Modelling the annual D/H cycle in a GCM, Montmessin et al. (2005) calculated that the D/H ratio of the atmospheric reservoir (~5 VSMOW) was about 15 per cent less than that prescribed in the simulations in the north polar cap (5.6 VSMOW). Using the most updated version of the GCM (Rossi et al. 2022; Vals et al. 2022), we calculate that this dichotomy between the two reservoirs decreases to 5.5 per cent (D/H prescribed in the north polar cap is 5 VSMOW; the annually averaged D/H in the atmospheric reservoir is 4.72 VSMOW). This value is also close to that reported by Daerden et al. (2022), where there is a difference of 5 per cent between the value prescribed on the north polar cap and that in the atmosphere.

Using this relation and the information about the global isotopic ratios derived from the SOFIA/EXES measurements, we can estimate the D/H of the north polar cap. We use two approaches to make this calculation. In the first approach, we assume that the measurements made in the northern summer (i.e. 2016 March, $L_s = 127°$) are the least affected by condensation-induced fractionation, since a large portion of the water mass is located in the atmospheric reservoir. Using this value, we estimate that the GCM simulations need to be corrected by 6.7 ± 4.7 per cent, which implies that the annually averaged D/H ratio in the atmospheric reservoir is $D/H = 5.03 ± 0.22$ VSMOW, and that in the north polar cap is $D/H = 5.34 ± 0.24$ VSMOW. In the second approach, we calculate the ratio between the model and measurements using all measurements, and calculate that the model should be scaled by 4.1 ± 2.3 per cent instead. In this case, the D/H of the atmospheric reservoir and the north polar cap would be given by 4.91 ± 0.11 and 5.20 ± 0.12 VSMOW, respectively. The results are indeed very similar when making both assumptions. The average between both cases would be $D/H = 5.0 ± 0.2$ VSMOW and $D/H = 5.3 ± 0.3$ VSMOW, where we kept the largest uncertainty.

In the case of ¹⁸O/¹⁶O, we expect variations in the global isotopic ratio to be smaller than the uncertainties. Therefore, we estimate the averaged isotopic composition of the atmospheric reservoir to be consistent with ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW (see Fig. 13). This value is consistent with previous measurements of the ¹⁸O/¹⁶O from the Curiosity Rover and the ExoMars TGO (Webster et al.
Using this relation and the information about the global isotopic ratios derived from the SOFIA/EXES measurements, we can estimate the D/H of the north polar cap. We use two approaches to make this calculation. In the first approach, we assume that the measurements made in the northern summer (i.e., March 2016, $L_S = 127°$) are the least affected by condensation-induced fractionation, since a large portion of the water mass is located in the atmospheric reservoir. Using this value, we estimate that the GCM simulations need to be corrected by $6.7 ± 4.7\%$, which implies that the annually averaged D/H ratio in the atmospheric reservoir is $D/H = 5.03 ± 0.22$ VSMOW, and that in the north polar cap is $D/H = 5.34 ± 0.24$ VSMOW. In the second approach, we calculate the ratio between the model and measurements using all measurements, and calculate that the model should be scaled by $4.1 ± 2.3\%$ instead. In this case, the D/H of the atmospheric reservoir and the north polar cap would be given by $4.91 ± 0.11$ and $5.20 ± 0.12$ VSMOW, respectively. The results are indeed very similar when making both assumptions. The average between both cases would be $D/H = 5.0 ± 0.2$ VSMOW and $D/H = 5.3 ± 0.3$ VSMOW, where we kept the largest uncertainty.

In the case of ¹⁸O/¹⁶O, we expect variations in the global isotopic ratio to be smaller than the uncertainties. Therefore, we estimate the averaged isotopic composition of the atmospheric reservoir to be consistent with ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW (see Figure 13). This value is consistent with previous measurements of the ¹⁸O/¹⁶O from the Curiosity Rover and the ExoMars Trace Gas Orbiter (Webster et al. 2013; Alday et al. 2019, 2021a). Similar to the case of D/H, the ¹⁸O/¹⁶O composition of the atmospheric water reservoir might also be controlled by that in the north polar cap. However, it must also be noted that the oxygen reservoir in water vapour might exchange with that in carbon dioxide through photochemical reactions, which makes the connection much more complex. In particular, Alday et al. (2021a) suggested that oxygen depleted in ¹⁸O/¹⁶O could be transferred from H₂O to CO₂ because of the preferential photolysis of H₂¹⁶O over H₂¹⁸O, which might enrich the bulk ¹⁸O/¹⁶O composition in H₂O with respect to that in CO₂ with time. This is similar to the differences in the C and O isotopic composition between CO₂ and CO, which are fractionated through photochemical reactions (Alday et al. 2023; Aoki et al. 2023; Yoshida et al. 2023). Because of this, without further photochemical modelling, it is difficult to say whether the global ¹⁸O/¹⁶O composition of water vapour is controlled by the value in the north polar cap, or is also buffered with the isotopic composition of CO₂.

## 5 CONCLUSIONS

Isotopic ratios in water vapour can provide important information about the characteristics of the water reservoir on early Mars and its evolution throughout history. Nevertheless, reconstructing the history of water from its isotopic ratios first requires disentangling the isotopic signatures induced by long-term and short-term fractionation.

Here, we analyse measurements of H₂O, HDO and H₂¹⁸O made with SOFIA/EXES throughout different periods of a Martian year ($L_S = 15°, 127°, 272°, 305°$) to constrain the global inventories of D/H and ¹⁸O/¹⁶O in the Martian water reservoir. In particular, we retrieve the abundances of the water isotopes from several absorption lines in the spectral range $v = 1382-1392$ cm⁻¹, using a full radiative transfer code that accounts for the instrumental characteristics of SOFIA/EXES. In addition, we develop a retrieval scheme to refine the accuracy of the retrieved column abundances when the assumed temperatures in the reference atmosphere depart from the real ones.

Our analysis of the SOFIA/EXES measurements reveals a mostly uniform spatial distribution of D/H and ¹⁸O/¹⁶O in the atmosphere of Mars, similar to previous studies using these data (Encrenaz et al. 2016, 2018). Nevertheless, we show that the absence of latitudinal variations is most likely caused by the poor spatial resolution of the instrument, which makes these measurements largely insensitive to the spatial variations. On the other hand, the high altitude of SOFIA during the observations and the high spectral resolution of the measurements removes the contamination from the terrestrial atmosphere, which enables accurate measurements of the global isotopic ratios in the Martian atmosphere.

We find the measured temporal variations of the D/H ratio to be consistent to first order with the predictions from a Global Climate Model accounting for condensation-induced fractionation (Vals et al. 2022; Rossi et al. 2022): the maximum global D/H ratio in water vapour in our observations is found during the northern summer ($D/H = 5.16 ± 0.23$ VSMOW), coinciding with the sublimation of the northern polar cap and the period with the highest atmospheric water content; the global D/H ratio during the southern summer is smaller than the annual maximum by just a few percent (~3% in the SOFIA/EXES observations, ~1% in the model); we find the lowest D/H ratio during the spring season in the northern hemisphere ($D/H = 4.31 ± 0.38$ VSMOW), which is substantially lower than the model predicts (~15% decrease in SOFIA/EXES with respect to northern summer, while the model predicts a decrease of only ~4%).

We combine the information about the seasonal trends in the planetary D/H ratio from the GCM simulations together with the measured values with SOFIA/EXES to estimate the global D/H representative of the water reservoir. In particular, we estimate the annually averaged D/H ratio in the atmospheric reservoir to be $D/H = 5.0 ± 0.2$ VSMOW, while the water sublimating from the north polar cap is about ~5% larger than in the atmosphere ($D/H = 5.3 ± 0.3$ VSMOW).

In the case of ¹⁸O/¹⁶O, we estimate the global isotopic ratio in water vapour to be ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW, which is consistent with previous investigations. While the value of ¹⁸O/¹⁶O in water vapour may be controlled by that of water ice in the northern polar cap, the oxygen isotopes in water vapour may also interact with those in carbon dioxide, which complicates the relation between the isotopic composition in different reservoirs.

## ACKNOWLEDGEMENTS

This article is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 2002 to the University of Stuttgart. Financial support for this work was provided by the Japan Society for the Promotion of Science (PE22045), the UK Space Agency and Science and Technology Facilities Council (ST/Y000234/1, ST/Y005414/1, ST/X001180/1, ST/X006549/1). SA is supported by JSPS KAKENHI Grant Number 22K03709, 22H05151, 22H00164, and 19H00707.

## DATA AVAILABILITY

The SOFIA/EXES observations analysed in this study are publicly available and can be accessed through the SOFIA Science Archive (https://irsa.ipac.caltech.edu/Missions/)
2013; Alday et al. 2019, 2021a). Similar to the case of D/H, the ¹⁸O/¹⁶O composition of the atmospheric water reservoir might also be controlled by that in the north polar cap. However, it must also be noted that the oxygen reservoir in water vapour might exchange with that in carbon dioxide through photochemical reactions, which makes the connection much more complex. In particular, Alday et al. (2021a) suggested that oxygen depleted in ¹⁸O/¹⁶O could be transferred from H₂O to CO₂ because of the preferential photolysis of H₂¹⁶O over H₂¹⁸O, which might enrich the bulk ¹⁸O/¹⁶O composition in H₂O with respect to that in CO₂ with time. This is similar to the differences in the C and O isotopic composition between CO₂ and CO, which are fractionated through photochemical reactions (Alday et al. 2023; Aoki et al. 2023; Yoshida et al. 2023). Because of this, without further photochemical modelling, it is difficult to say whether the global ¹⁸O/¹⁶O composition of water vapour is controlled by the value in the north polar cap, or is also buffered with the isotopic composition of CO₂.

## 5. CONCLUSIONS

Isotopic ratios in water vapour can provide important information about the characteristics of the water reservoir on early Mars and its evolution throughout history. Nevertheless, reconstructing the history of water from its isotopic ratios first requires disentangling the isotopic signatures induced by long- and short-term fractionation.

Here, we analyse measurements of H₂O, HDO, and H₂¹⁸O made with SOFIA/EXES throughout different periods of a Martian year (Ls = 15°, 127°, 272°, and 305°) to constrain the global inventories of D/H and ¹⁸O/¹⁶O in the Martian water reservoir. In particular, we retrieve the abundances of the water isotopes from several absorption lines in the spectral range v = 1382–1392 cm⁻¹, using a full radiative transfer code that accounts for the instrumental characteristics of SOFIA/EXES. In addition, we develop a retrieval scheme to refine the accuracy of the retrieved column abundances when the assumed temperatures in the reference atmosphere depart from the real ones.

Our analysis of the SOFIA/EXES measurements reveals a mostly uniform spatial distribution of D/H and ¹⁸O/¹⁶O in the atmosphere of Mars, similar to previous studies using these data (Encrenaz et al. 2016, 2018). Nevertheless, we show that the absence of latitudinal variations is most likely caused by the poor spatial resolution of the instrument, which makes these measurements largely insensitive to the spatial variations. On the other hand, the high altitude of SOFIA during the observations and the high spectral resolution of the measurements remove the contamination from the terrestrial atmosphere, which enables accurate measurements of the global isotopic ratios in the Martian atmosphere.

We find the measured temporal variations of the D/H ratio to be consistent to first order with the predictions from a Global Climate Model accounting for condensation-induced fractionation (Rossi et al. 2022; Vals et al. 2022): The maximum global D/H ratio in water vapour in our observations is found during the northern summer (D/H = 5.16 ± 0.23 VSMOW), coinciding with the sublimation of the northern polar cap and the period with the highest atmospheric water content; the global D/H ratio during the southern summer is smaller than the annual maximum by just a few per cent (~3 per cent in the SOFIA/EXES observations, ~1 per cent in the model); we find the lowest D/H ratio during the spring season in the Northern hemisphere (D/H = 4.31 ± 0.38 VSMOW), which is substantially lower than the model predicts (~15 per cent decrease in SOFIA/EXES with respect to northern summer, while the model predicts a decrease of only ~4 per cent).

We combine the information about the seasonal trends in the planetary D/H ratio from the GCM simulations together with the measured values with SOFIA/EXES to estimate the global D/H representative of the water reservoir. In particular, we estimate the annually averaged D/H ratio in the atmospheric reservoir to be D/H = 5.0 ± 0.2 VSMOW, while the water sublimating from the north polar cap is about ~5 per cent larger than in the atmosphere (D/H = 5.3 ± 0.3 VSMOW).

In the case of ¹⁸O/¹⁶O, we estimate the global isotopic ratio in water vapour to be ¹⁸O/¹⁶O = 1.09 ± 0.08 VSMOW, which is consistent with previous investigations. While the value of ¹⁸O/¹⁶O in water vapour may be controlled by that of water ice in the northern polar cap, the oxygen isotopes in water vapour may also interact with those in carbon dioxide, which complicates the relation between the isotopic composition in different reservoirs.

## ACKNOWLEDGEMENTS

This article is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 2002 to the University of Stuttgart. Financial support for this work was provided by the Japan Society for the Promotion of Science (PE22045), the UK Space Agency, and Science and Technology Facilities Council (ST/Y000234/1, ST/Y005414/1, ST/X001180/1, and ST/X006549/1). SA is supported by JSPS KAKENHI Grant Number 22K03709, 22H05151, 22H00164, and 19H00707.

## DATA AVAILABILITY

The SOFIA/EXES observations analysed in this study are publicly available and can be accessed through the SOFIA Science Archive (https://irsa.ipac.caltech.edu/Missions/sofia.html). The NEMESIS algorithm, used to perform the radiative transfer calculations, can be downloaded from Irwin (2022). The D/H data from the Global Climate Model used in this study can be accessed from Vals & Rossi (2022). The parameters from the OpenMARS data base used for defining the reference atmosphere in the analysis of the measurements can be downloaded from Holmes et al. (2023). The data products derived from the analysis of the measurements in this study are available in Alday (2023).

## REFERENCES

* Acton C., Bachman N., Semenov B., Wright E., 2018, Planet. Space Sci., 150, 9

* Alday J., 2023, Isotopic composition of Martian water vapour from SOFIA/EXES, The Open University, Version 1, https://doi.org/10.219 54/ou.rd.24580579.v1

* Alday J. et al., 2019, A&A, 630, A91

* Alday J. et al., 2021a, Nature Astron., 5, 943

* Alday J. et al., 2021b, J. Geophys. Res. Planets, 126, e06992

* Alday J. et al., 2023, Nature Astron., 7, 867

* Alsaed N. R., Jakosky B. M., 2019, J. Geophys. Res. Planets, 124, 3344

* Aoki S., Nakagawa H., Sagawa H., Giuranna M., Sindoni G., Aronica A., Kasaba Y., 2015, Icarus, 260, 7

* Aoki S. et al., 2023, Planet. Sci. J., 4, 97

* Baker V. R., 2001, Nature, 412, 228

* Bjoraker G. L., Mumma M. J., Larson H. P., 1989, BAAS, 21, 991

* Cangi E., Chaffin M., Yelle R., Gregory B., Deighan J., 2023, J. Geophys. Res. Planets, 128, e2022JE007713
sofia.html). The NEMESIS algorithm, used to perform the radiative transfer calculations, can be downloaded from Irwin (2022). The D/H data from the Global Climate Model used in this study can be accessed from Vals & Rossi (2022). The parameters from the OpenMARS database used for defining the reference atmosphere in the analysis of the measurements can be downloaded from Holmes et al. (2023). The data products derived from the analysis of the measurements in this study are available in Alday (2023).

## REFERENCES

Acton C., Bachman N., Semenov B., Wright E., 2018, *Planetary and Space Science*, 150, 9

Alday J., 2023, Isotopic composition of Martian water vapour from SOFIA/EXES, doi:10.21954/OU.RD.24580579.V1, https://ordo.open.ac.uk/articles/dataset/Isotopic_composition_of_Martian_water_vapour_from_SOFIA_EXES/24580579/1

Alday J., et al., 2019, Astronomy & Astrophysics, 630, A91

Alday J., et al., 2021a, Nature Astronomy, 5, 943

Alday J., et al., 2021b, Journal of Geophysical Research: Planets, 126

Alday J., et al., 2023, Nature Astronomy, 7, 867

Alsaeed N. R., Jakosky B. M., 2019, *Journal of Geophysical Research: Planets*, 124, 3344

Aoki S., Nakagawa H., Sagawa H., Giuranna M., Sindoni G., Aronica A., Kasaba Y., 2015, Icarus, 260, 7

Aoki S., et al., 2023, *The Planetary Science Journal*, 4, 97

Baker V. R., 2001, Nature, 412, 228

Bjoraker G. L., Mumma M. J., Larson H. P., 1989, in *Bulletin of the American Astronomical Society*. p. 991, http://adsabs.harvard.edu/abs/1989BAAS...21Q.991B

Cangi E., Chaffin M., Yelle R., Gregory B., Deighan J., 2023, *Journal of Geophysical Research: Planets*, 128, e2022JE007713

Carr M. H., Clow G. D., 1981, Icarus, 48, 91

Cheng B.-M., Chew E. P., Liu C.-P., Bahou M., Lee Y.-P., Yung Y. L., Gerstell M. F., 1999, *Geophysical Research Letters*, 26, 3657

Clancy R. T., 2003, *Journal of Geophysical Research*, 108, 5098

Connour K., Wolff M. J., Schneider N. M., Deighan J., Lefèvre F., Jain S. K., 2022, Icarus, 387, 115177

Craig H., 1961, Science, 133, 1702

Crismani M. M. J., et al., 2021, *Journal of Geophysical Research: Planets*, 126

Daerden F., et al., 2022, *Journal of Geophysical Research: Planets*, 127

DeWitt C., Montiel E., Rashman M., 2023, *EXES Handbook for Archive Users* v1.0, https://irsa.ipac.caltech.edu/data/SOFIA/docs/instruments/handbooks/EXES_Handbook_for_Archive_Users_Ver1.0.pdf

Devi V. M., et al., 2017, *Journal of Quantitative Spectroscopy and Radiative Transfer*, 203, 158

Encrenaz T., et al., 2016, *Astronomy & Astrophysics*, 586, A62

Encrenaz T., et al., 2018, *Astronomy & Astrophysics*, 612, A112

Gordon I. E., et al., 2022, *Journal of Quantitative Spectroscopy and Radiative Transfer*, 277, 107949

Holmes J. A., Lewis S. R., Patel M. R., 2020, *Planetary and Space Science*, 188, 104962

Holmes J. A., et al., 2022, *Journal of Geophysical Research: Planets*, 127

Holmes J., Lewis S., Patel M., 2023, *OpenMARS* database, doi:10.21954/OU.RD.C.4278950.V2, https://ordo.open.ac.uk/collections/OpenMARS_database/4278950/2

Hu R., 2019, *Earth and Planetary Science Letters*, 519, 192

Irwin P., 2022, *nemesiscode/radtrancode: NEMESIS*, doi:10.5281/ZENODO.5816714, https://zenodo.org/record/5816714

Irwin P., et al., 2008, *Journal of Quantitative Spectroscopy and Radiative Transfer*, 109, 1136

Jakosky B. M., Phillips R. J., 2001, *Nature*, 412, 237

Khayat A. S., Villanueva G. L., Smith M. D., Guzewich S. D., 2019, *Icarus*, 330, 204

Kleinböhl A., Friedson A. J., Schofield J. T., 2017, *Journal of Quantitative Spectroscopy and Radiative Transfer*, 187, 511

Knutsen E. W., et al., 2022, *Journal of Geophysical Research: Planets*, 127, e2022JE007252

Krasnopolsky V., 2000, *Icarus*, 148, 597

Krasnopolsky V. A., 2015, *Icarus*, 257, 377

Krasnopolsky V. A., 2021, *Planetary and Space Science*, 208, 105345

Krasnopolsky V. A., Bjoraker G. L., Mumma M. J., Jennings D. E., 1997, *Journal of Geophysical Research: Planets*, 102, 6525

Merlivat L., Nief G., 1967, *Tellus*, 19, 122

Montmessin F., Forget F., Rannou P., Cabane M., Haberle R. M., 2004, *Journal of Geophysical Research: Planets*, 109, E10004

Montmessin F., Fouchet T., Forget F., 2005, *Journal of Geophysical Research: Planets*, 110, E03006

Montmessin F., Smith M. D., Langevin Y., Mellon M. T., Fedorova A., 2017, in Haberle R. M., Clancy R. T., Forget F., Smith M. D., Zurek R. W., eds., *The Atmosphere and Climate of Mars*, 1 edn, Cambridge University Press, pp 338–373, doi:10.1017/9781139060172.011, https://www.cambridge.org/core/product/identifier/9781139060172%23CN-bp-11/type/book_part

Navarro T., Madeleine J.-B., Forget F., Spiga A., Millour E., Montmessin F., Määttänen A., 2014, *Journal of Geophysical Research: Planets*, 119, 1479

Novak R., Mumma M., Villanueva G., 2011, *Planetary and Space Science*, 59, 163

Owen T., Maillard J. P., de Bergh C., Lutz B. L., 1988, *Science*, 240, 1767

Richter M. J., DeWitt C. N., McKelvey M., Montiel E., McMurray R., Case M. E., 2018, *Journal of Astronomical Instrumentation*, 07, 1840013

Rodgers C. D., 2000, *Inverse Methods for Atmospheric Sounding: Theory and Practice. Series on Atmospheric, Oceanic and Planetary Physics Vol. 2, World Scientific*, doi:10.1142/3171, https://www.worldscientific.com/worldscibooks/10.1142/3171

Rossi L., et al., 2022, *Journal of Geophysical Research: Planets*, 127

Régalia L., Cousin E., Gamache R., Vispoel B., Robert S., Thomas X., 2019, *Journal of Quantitative Spectroscopy and Radiative Transfer*, 231, 126

Scheller E. L., Ehlmann B. L., Hu R., Adams D. J., Yung Y. L., 2021, *Science*, p. eabc7717

Smith M. D., Daerden F., Neary L., Khayat A., 2018, *Icarus*, 301, 117

Trokhimovskiy A., Fedorova A., Korablev O., Montmessin F., Bertaux J.-L., Rodin A., Smith M. D., 2015, *Icarus*, 251, 50

Vals M., Rossi L., 2022, *Simulations of the HDO cycle on Mars with the upgraded LMD Mars GCM*, doi:10.14768/9BA77F68-29A3-4BD6-8966-0900E2316A54, https://data.ipsl.fr/catalog/metadata/9ba77f68-29a3-4bd6-8966-0900e2316a54

Vals M., et al., 2022, preprint, *Improved modeling of Mars’ HDO cycle using a Mars’ Global Climate Model*, http://www.essoar.org/doi/10.1002/essoar.10510380.1. Planetology, doi:10.1002/essoar.10510380.1, http://www.essoar.org/doi/10.1002/essoar.10510380.1

Villanueva G. L., et al., 2015, *Science*, 348, 218

Villanueva G., Smith M., Protopapa S., Faggi S., Mandell A., 2018, *Journal of Quantitative Spectroscopy and Radiative Transfer*, 217, 86

Villanueva G. L., et al., 2021, *Science Advances*, 7, eabc8843

Villanueva G. L., et al., 2022, *Geophysical Research Letters*, 49

Warren S. G., Brandt R. E., 2008, *Journal of Geophysical Research*, 113, D14220

Webster C. R., et al., 2013, *Science*, 341, 260

Yoshida T., et al., 2023, *The Planetary Science Journal*, 4, 53

Yung Y. L., Wen J.-S., Pinto J. P., Allen M., Pierce K. K., Paulson S., 1988, *Icarus*, 76, 146

This paper has been typeset from a TeX/LA TeX file prepared by the author.
Carr M. H., Clow G. D., 1981, *Icarus*, **48**, 91

Cheng B.-M., Chew E. P., Liu C.-P., Bahou M., Lee Y.-P., Yung Y. L., Gerstell M. F., 1999, *Geophys. Res. Lett.*, **26**, 3657

Clancy R. T., 2003, *J. Geophys. Res.*, **108**, 5098

Connour K., Wolff M. J., Schneider N. M., Deighan J., Lefèvre F., Jain S. K., 2022, *Icarus*, **387**, 115177

Craig H., 1961, *Science*, **133**, 1702

Crismani M. M. J. et al., 2021, *J. Geophys. Res. Planets*, **126**, e06878

Daerden F. et al., 2022, *J. Geophys. Res. Planets*, **127**, e07079

DeWitt C., Montiel E., Rashman M., 2023, *EXES Handbook for Archive Users, SOFIA - Stratospheric Observatory for Infrared Astronomy, Rev1.0*, https://irsa.ipac.caltech.edu/data/SOFIA/docs/instruments/handbooks/EXES_Handbook_for_Archive_Users_Ver1.0.pdf

Devi V. M. et al., 2017, *J. Quant. Spectrosc. Radiat. Transfer*, **203**, 158

Encrenaz T. et al., 2016, *A&A*, **586**, A62

Encrenaz T. et al., 2018, *A&A*, **612**, A112

Gordon I. E. et al., 2022, *J. Quant. Spectrosc. Radiat. Transfer*, **277**, 107949

Holmes J. A., Lewis S. R., Patel M. R., 2020, *Planet. Space Sci.*, **188**, 104962

Holmes J. A. et al., 2022, *J. Geophys. Res. Planets*, **127**, e2022JE007203

Holmes J., Lewis S., Patel M., 2023, *OpenMARS* database, The Open University, https://doi.org/10.21954/ou.rd.c.4278950

Hu R., 2019, *Earth Planet. Sci. Lett.*, **519**, 192

Irwin P. et al., 2008, *J. Quant. Spectrosc. Radiat. Transfer*, **109**, 1136

Irwin P. et al., 2022, *nemesiscode/radtrancode: NEMESIS*. Zenodo, https://doi.org/10.5281/zenodo.581671

Jakošky B. M., Phillips R. J., 2001, *Nature*, **412**, 237

Khayat A. S., Villanueva G. L., Smith M. D., Guzewich S. D., 2019, *Icarus*, **330**, 204

Kleinböhl A., Friedson A. J., Schofield J. T., 2017, *J. Quant. Spectrosc. Radiat. Transfer*, **187**, 511

Knutsen E. W. et al., 2022, *J. Geophys. Res. Planets*, **127**, e2022JE007252

Krasnopolsky V., 2000, *Icarus*, **148**, 597

Krasnopolsky V. A., 2015, *Icarus*, **257**, 377

Krasnopolsky V. A., 2021, *Planet. Space Sci.*, **208**, 105345

Krasnopolsky V. A., Bjoraker G. L., Mumma M. J., Jennings D. E., 1997, *J. Geophys. Res. Planets*, **102**, 6525

Merlivat L., Nief G., 1967, *Tellus*, **19**, 122

Montmessin F., Forget F., Rannou P., Cabane M., Haberle R. M., 2004, *J. Geophys. Res. Planets*, **109**, E10004

Montmessin F., Fouchet T., Forget F., 2005, *J. Geophys. Res. Planets*, **110**, E03006

Montmessin F., Smith M. D., Langevin Y., Mellon M. T., Fedorova A., 2017, in *Haberle R. M., Clancy R. T., Forget F., Smith M. D., Zurek R. W., eds, The Atmosphere and Climate of Mars*, 1st edn. Cambridge Univ. Press, Cambridge, p. 338

Navarro T., Madeleine J.-B., Forget F., Spiga A., Millour E., Montmessin F., Määttänen A., 2014, *J. Geophys. Res. Planets*, **119**, 1479

Novak R., Mumma M., Villanueva G., 2011, *Planet. Space Sci.*, **59**, 163

Owen T., Maillard J. P., de Bergh C., Lutz B. L., 1988, *Science*, **240**, 1767

Régalia L., Cousin E., Gamache R., Vispoel B., Robert S., Thomas X., 2019, *J. Quant. Spectrosc. Radiat. Transfer*, **231**, 126

Richter M. J., DeWitt C. N., McKelvey M., Montiel E., McMurray R., Case M. E., 2018, *J. Astron. Instrum.*, **07**, 1840013

Rodgers C. D., 2000, *Inverse Methods for Atmospheric Sounding: Theory and Practice, Series on Atmospheric, Oceanic and Planetary Physics*, Vol. 2. World Scientific Press, Singapore

Rossi L. et al., 2022, *J. Geophys. Res. Planets*, **127**, e07201

Scheller E. L., Ehlmann B. L., Hu R., Adams D. J., Yung Y. L., 2021, *Science*, **372**, 56

Smith M. D., Daerden F., Neary L., Khayat A., 2018, *Icarus*, **301**, 117

Trokhimovskiy A., Fedorova A., Korablev O., Montmessin F., Bertaux J.-L., Rodin A., Smith M. D., 2015, *Icarus*, **251**, 50

Vals M., Rossi L., 2022, *Simulations of the HDO cycle on Mars with the upgraded LMD Mars GCM, ESPRI/IPSL*, https://doi.org/10.14768/9BA77F68-29A3-4BD6-8966-0900E2316A54

Vals M. et al., 2022, *ESS Open Arch. eprints*, 105, essoar.10510380

Villanueva G. L. et al., 2015, *Science*, **348**, 218

Villanueva G. L., Smith M., Protopapa S., Faggi S., Mandell A., 2018, *J. Quant. Spectrosc. Radiat. Transfer*, **217**, 86

Villanueva G. L. et al., 2021, *Sci. Adv.*, **7**, eabc8843

Villanueva G. L. et al., 2022, *Geophys. Res. Lett.*, **49**, e98161

Warren S. G., Brandt R. E., 2008, *J. Geophys. Res.*, **113**, D14220

Webster C. R. et al., 2013, *Science*, **341**, 260

Yoshida T. et al., 2023, *Planet. Sci. J.*, **4**, 53

Yung Y. L., Wen J.-S., Pinto J. P., Allen M., Pierce K. K., Paulson S., 1988, *Icarus*, **76**, 146

## SUPPORTING INFORMATION

Supplementary data are available at *MNRAS* online.

### Alday2023_MNRAS_SOFIA_SM_v4.pdf

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.