# Simulations of the hydrogen and deuterium thermal and non-thermal escape at Mars at Spring Equinox

Jean-Yves Chaufray, F. Gonzalez-Galindo, François Leblanc, Ronan Modolo, Margaux Vals, Franck Montmessin, Franck Lefèvre, F. Forget, M.

Lopez-Valverde, G. Gilli

► To cite this version:

Jean-Yves Chaufray, F. Gonzalez-Galindo, François Leblanc, Ronan Modolo, Margaux Vals, et al. Simulations of the hydrogen and deuterium thermal and non-thermal escape at Mars at Spring Equinox. Icarus, 2024, 418, pp.116152. 10.1016/j.icarus.2024.116152. insu-04591622

HAL Id: insu-04591622

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

Submitted on 31 May 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_1.png" /></figure>

# Simulations of the hydrogen and deuterium thermal and non-thermal escape at Mars at Spring Equinox

J.-Y. Chaufray<sup>a, *</sup>, F. Gonzalez-Galindo<sup>b</sup>, F. Leblanc<sup>a</sup>, R. Modolo<sup>a</sup>, M. Vals<sup>a</sup>, F. Montmessin<sup>a</sup>, F. Lefèvre<sup>a</sup>, F. Forget<sup>c</sup>, M. Lopez-Valverde<sup>b</sup>, G. Gilli<sup>b</sup>

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

## ARTICLE INFO

Keywords:

Mars (1007)

Planetary ionospheres (2185)

Atmospheric evolution (2301)

## ABSTRACT

We present simulations of the thermal and nonthermal escape processes for H and D, under atomic, molecular and ion forms at Mars during spring equinox. These processes include Jeans escape, several photochemical reactions and the escape associated to the solar wind interaction with Mars. While the hydrogen escape is dominated by the atomic Jeans escape, we find that the deuterium escape is dominated by the photochemical atomic escape. Ions escape represent only 10% of the total escape for both species and is mostly due to charge exchange between neutral and solar wind protons. Including all the processes, we find a D/H fractionation factor (D/H escape ratio divided by the D/H atmospheric ratio) f = 0.04, with a main uncertainty associated to the elastic collisional cross sections needed to accurately derive the photochemical escape rate. Using this fractionation factor and considering a 30 m exchangeable reservoir of water, the average hydrogen escape rate needed to fractionate the Martian water from its primordial value to its current D/H value during the last 4.5 Gyr is ~1.0 × 10<sup>28</sup> s<sup>-1</sup> which is larger than the current average escape rate (~ 2 × 10<sup>26</sup> s<sup>-1</sup>).

## 1. Introduction

Numerous evidences suggest that liquid water was present in the Martian past not only during the Noachian, eroding the surface and forming the valley networks (Carr, 1996) and the phyllosilicate minerals (Bibring et al., 2006), but also during the Hesperian, forming the outflow channels and a possible ocean during the late Hesperian (~3 Gyr ago) (Dickeson and Davis, 2020; Schmidt et al., 2022). The fate of this water is still unknown, but a substantial part of it could have escaped into the interplanetary medium as suggested by the D/H ratio measured in the Martian water vapor (Owen et al., 1988) or in the Hesperian clay minerals (Mahaffy et al., 2015). The D/H ratio measured in the current Martian atmosphere is 5 to 6 times the Standard Mean Ocean Water (SMOW) value (Owen et al., 1988; Krasnopolsky et al., 1997; Webster et al., 2013) and has also been observed to vary spatially (Novak et al., 2011; Villanueva et al., 2015; Aoki et al., 2015; Krasnopolsky, 2015). The D/H ratio has been estimated to 1.9 after the end of the hydrodynamic phase (Leshin, 2000) and Curiosity/SAM measured a D/H ratio of 3 ± 0.2 SMOW in strongly bound water or hydroxyl in

ancient (hesperian) martian clays (Mahaffy et al., 2015), indicating an increase of the D/H fractionation from the end of the hydrodynamics phase (~ 4 Gyr ago) (Jakosky and Phillips, 2001) to the Hesperian and another increase from Hesperian to now, most likely driven by the preferential atmospheric escape of H compared to D due to its lower atomic mass (Jakosky et al., 2018).

Let us assume a model with two water reservoirs, consisting of the atmosphere and the seasonally exchanging polar ices (reservoir 1) and the permanent polar ices (reservoir 2). If we assume that the D/H ratio is the same for the two reservoirs, and, noting N<sub>H</sub>(t) and N<sub>D</sub>(t) the total number of H and D atoms (under atomic or molecular forms) at a given time t, the evolution of the ratio [D/H](t) = N<sub>H</sub>(t)/N<sub>D</sub>(t) during an infinitesimal time variation dt is given by

$$ \left[ \frac{D}{H} \right] (t+dt) = \frac{N_D(t) - \Phi_D(t)dt}{N_H(t) - \Phi_H(t)dt} \qquad (1) $$

where Φ<sub>H</sub>(t) and Φ<sub>D</sub>(t) are the hydrogen and deuterium atmospheric loss rate (in atoms per time units) at time t.

A first order development of this equation leads to the differential

<footnote><sup>a</sup> LATMOS-IPSL, CNRS, UVSQ, Sorbonne Université, Guyancourt, France</footnote>

<footnote><sup>b</sup> IAA, CSIC, Granada, Spain</footnote>

<footnote><sup>c</sup> LMD-IPSL, CNRS, Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: chaufray@latmos.ipsl.fr (J.-Y. Chaufray).</footnote>
<figure><img src="image_1.png" /></figure>

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

# Simulations of the hydrogen and deuterium thermal and non-thermal escape at Mars at Spring Equinox

J.-Y. Chaufray<sup>a, *</sup>, F. Gonzalez-Galindo<sup>b</sup>, F. Leblanc<sup>a</sup>, R. Modolo<sup>a</sup>, M. Vals<sup>a</sup>, F. Montmessin<sup>a</sup>, F. Lefèvre<sup>a</sup>, F. Forget<sup>c</sup>, M. Lopez-Valverde<sup>b</sup>, G. Gilli<sup>b</sup>

## ARTICLE INFO

Keywords:

Mars (1007)

Planetary ionospheres (2185)

Atmospheric evolution (2301)

## ABSTRACT

We present simulations of the thermal and nonthermal escape processes for H and D, under atomic, molecular and ion forms at Mars during spring equinox. These processes include Jeans escape, several photochemical reactions and the escape associated to the solar wind interaction with Mars. While the hydrogen escape is dominated by the atomic Jeans escape, we find that the deuterium escape is dominated by the photochemical atomic escape. Ions escape represent only 10% of the total escape for both species and is mostly due to charge exchange between neutral and solar wind protons. Including all the processes, we find a D/H fractionation factor (D/H escape ratio divided by the D/H atmospheric ratio) f = 0.04, with a main uncertainty associated to the elastic collisional cross sections needed to accurately derive the photochemical escape rate. Using this fractionation factor and considering a 30 m exchangeable reservoir of water, the average hydrogen escape rate needed to fractionate the Martian water from its primordial value to its current D/H value during the last 4.5 Gyr is ~1.0 × 10<sup>28</sup> s<sup>-1</sup> which is larger than the current average escape rate (~ 2 × 10<sup>26</sup> s<sup>-1</sup>).

## 1. Introduction

Numerous evidences suggest that liquid water was present in the Martian past not only during the Noachian, eroding the surface and forming the valley networks (Carr, 1996) and the phyllosilicate minerals (Bibring et al., 2006), but also during the Hesperian, forming the outflow channels and a possible ocean during the late Hesperian (~3 Gyr ago) (Dickeson and Davis, 2020; Schmidt et al., 2022). The fate of this water is still unknown, but a substantial part of it could have escaped into the interplanetary medium as suggested by the D/H ratio measured in the Martian water vapor (Owen et al., 1988) or in the Hesperian clay minerals (Mahaffy et al., 2015). The D/H ratio measured in the current Martian atmosphere is 5 to 6 times the Standard Mean Ocean Water (SMOW) value (Owen et al., 1988; Krasnopolsky et al., 1997; Webster et al., 2013) and has also been observed to vary spatially (Novak et al., 2011; Villanueva et al., 2015; Aoki et al., 2015; Krasnopolsky, 2015). The D/H ratio has been estimated to 1.9 after the end of the hydrodynamic phase (Leshin, 2000) and Curiosity/SAM measured a D/H ratio of 3 ± 0.2 SMOW in strongly bound water or hydroxyl in

ancient (hesperian) martian clays (Mahaffy et al., 2015), indicating an increase of the D/H fractionation from the end of the hydrodynamics phase (~ 4 Gyr ago) (Jakosky and Phillips, 2001) to the Hesperian and another increase from Hesperian to now, most likely driven by the preferential atmospheric escape of H compared to D due to its lower atomic mass (Jakosky et al., 2018).

Let us assume a model with two water reservoirs, consisting of the atmosphere and the seasonally exchanging polar ices (reservoir 1) and the permanent polar ices (reservoir 2). If we assume that the D/H ratio is the same for the two reservoirs, and, noting N<sub>H</sub>(t) and N<sub>D</sub>(t) the total number of H and D atoms (under atomic or molecular forms) at a given time t, the evolution of the ratio [D/H](t) = N<sub>H</sub>(t)/N<sub>D</sub>(t) during an infinitesimal time variation dt is given by

$$ \left[ \frac{D}{H} \right] (t+dt) = \frac{N_D(t) - \Phi_D(t)dt}{N_H(t) - \Phi_H(t)dt} \quad (1) $$

where Φ<sub>H</sub>(t) and Φ<sub>D</sub>(t) are the hydrogen and deuterium atmospheric loss rate (in atoms per time units) at time t.

A first order development of this equation leads to the differential

<footnote><sup>a</sup> LATMOS-IPSL, CNRS, UVSQ, Sorbonne Université, Guyancourt, France</footnote>

<footnote><sup>b</sup> IAA, CSIC, Granada, Spain</footnote>

<footnote><sup>c</sup> LMD-IPSL, CNRS, Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: chaufray@latmos.ipsl.fr (J.-Y. Chaufray).</footnote>
equation of the evolution of D/H:

$$ \frac{d[D/H]}{dt} = \left[ \frac{D/H}{H} \right] (t) \frac{1 - f(t)}{\tau_H(t)} \qquad (2) $$

where $\tau_H(t) = N_H(t)/\Phi_H(t)$ is the time scale of the hydrogen loss and $f(t)$ is the fractionation factor given by

$$ f(t) = \frac{\Phi_D(t)/\Phi_H(t)}{N_D(t)/N_H(t)} \qquad (3) $$

Two obvious conditions must be fulfilled to produce fractionation during escape (Eq. (2)). The time scale $\tau_H$ of the hydrogen escape rate should not be much larger than the planet's age, the fractionation factor should be different from one. If $f(t) = 1$ all along the planet history, then H and D escape in proportion of their relative abundance leading to no fractionation. Because the mass of H is lower than D, the atmospheric escape is more efficient for H than D and then $f(t) < 1$ most of the time leading to an increase of D/H with time. The time scale of the D/H fractionation is $\tau_H(t)/(1-f(t)) \le \tau_H(t)$ and would be equal to the hydrogen loss timescale if the escape of D was negligible ($f(t) = 0$).

The two reservoirs model can be expanded to a 3 reservoirs model (Krasnopolsky, 2000) where a third reservoir (reservoir 3) with the primordial D/H ratio of 1 SMOW is added to the previous system. The D/H ratio of this reservoir is denoted (D/H)3. This reservoir does not exchange water directly to the reservoir 1, but can supply the second reservoir (e.g. permanent ice caps) in water to maintain a constant hydrogen abundance in reservoir 2. The reservoir 2 can exchange (on long time scale) with the atmosphere and is replenished from reservoir 3. The evolution of the D/H ratio in reservoir 2 relative to the primordial D/H ratio in reservoir 3 ($r = (D/H)_2/(D/H)_3$) is solution of the differential equation (Krasnopolsky, 2000):

$$ \frac{dr}{dt} = \frac{1 - r(t)f(t)}{\tau_H(t)} \qquad (4) $$

where the evolution of r is still determined by the two parameters $\tau_H(t)$ and $f(t)$ defined previously. At $t = 0$, $r(t) = 1$ so if $f(t) = 1$ all the time, we still find no fractionation of the water as expected. But once the reservoir 2 has been fractionated ($r(t) > 1$), $f(t)$ should be smaller than $1/r(t)$ in order to continue having an increase of $r(t)$. If at a time t, $f(t) = 1/r(t)$, the D/H increase of reservoir 2 due to the preferential loss of H is balanced by its decrease due to the lower primordial D/H flux coming from reservoir 3.

Currently, hydrogen escape, which is dominated by Jeans escape, has been estimated from a diverse set of observations (Anderson and Hord, 1971; Chaufray et al., 2008; Chaffin et al., 2014, 2018; Halekas, 2017; Rahmati et al., 2018; Qin, 2021; Mayyasi et al., 2023) and modeling (Krasnopolsky, 2002, 2010, 2019; Chaufray et al., 2015, 2021a), showing large seasonal variations. Atomic deuterium is more difficult to measure but past measurements and more recent MAVEN measurements have shown also large seasonal variations of its concentration in the upper atmosphere (Clarke et al., 2017; Mayyasi et al., 2017; Chaufray et al., 2021b), and thus of its escape. But contrary to H, D escape is expected to be mostly non-thermal (Krasnopolsky, 2002, 2010) due to its larger mass. The main photochemical escape processes for D have been studied in details only recently (Cangi et al., 2023) using photochemical modeling of deuterium ion chemistry and Monte Carlo modeling of hot escape probabilities (Gregory et al., 2023a). These simulations have shown that DCO$^+$ recombination and the reaction between CO$_2^+$ and HD are the two most important sources of deuterium escape. In this study, we simulate the most important escape processes for the major form of H and D (H, H$_2$, D, HD, H$^+$, D$^+$) in the upper atmosphere of Mars including: thermal, photochemical and ions escape, the latter driven by the solar wind interaction with Mars' atmosphere. The photochemical escape includes not only the ion-neutral reactions

able to produce hot hydrogen and deuterium (Gregory et al., 2023b; Cangi et al., 2023) but also the collisions between hot oxygen produced from the O$_2^+$ dissociative recombination (e.g. Deighan et al., 2015, Leblanc et al., 2017, Lillis et al., 2017, Jakosky et al., 2018) with atomic and molecular hydrogen and deuterium (Gacesa et al., 2012; Shematovich, 2013). The different models used to estimate the different escape processes are described in Section 2 and the simulated values are presented in Section 3. These values of f(t) and $\tau_H(t)$ as well as the current main limits of our approach are discussed in Section 4 before a conclusion and perspectives of the study (Section 5).

## 2. Models

### 2.1. Model of the Martian atmosphere and the thermal escape

To compute the thermal escape of H, D, H$_2$ and HD, we use the Planetary Climate Model (PCM) in its Martian version (Mars-PCM). Mars-PCM is a 3D hydrodynamic model of the Martian atmosphere and ionosphere from the surface to the exobase near ~200 km, the altitude at which thermal escape is calculated (Forget et al., 1999; Gonzalez-Galindo et al., 2009, 2013, 2015). The vertical coordinate is the atmospheric pressure, and the upper altitude is not constant but can vary between ~200 to 250 km.

This model has been used in the past to compute the H and H$_2$ thermal escape (Chaufray et al. 2015, 2018, 2021). The spatial resolution is 5.625° in longitude, 3.75° latitude and a pressure vertical coordinate is used leading to a variable vertical resolution, which is typically of ~10 km in the thermosphere.

The main updates used in this study that affect the H and D species in the thermosphere, are the inclusion of detailed cloud microphysics that considers the details of the formation of clouds and their microphysics (nucleation, ice growth, scavenging) and their radiative effect, needed to simulate the supersaturation (Navarro et al., 2014; Vals et al., 2022; Naar, 2023), HDO as a tracer in its vapor and icy phases (Vals et al., 2022; Rossi et al., 2022) and the incorporation of the photochemistry of deuterated species, including ions (Table 1). The inclusion of the supersaturation makes the hygropause more porous, allowing H$_2$O and HDO to reach the upper atmosphere (Vals et al., 2022).

The coefficients of the reactions between HD and an ion i, $k_{HD,i}$, are approximated from the Langevin frequency (e.g. Krasnopolsky, 2002) assuming the same polarizability between HD and H$_2$ (Ishiguro et al., 1952), yet considering the difference in the mass, leading to (Krasnopolsky, 2002)

$$ k_{HD,i} = \sqrt{2/3}k_{H2,i} $$

These reactions can produce atomic H or D and we assume each channel to be equiprobable. For example, the important reaction H$_2$ + CO$_2^+$ has a reaction rate $9.5 \times 10^{10}$ cm$^3$.s$^{-1}$ in Mars-PCM (Table 1). The equivalent reaction with HD: HD + CO$_2^+$ has a total reaction rate of $(2/3)^{1/2} \times 9.5 \times 10^{10}$ cm$^3$.s$^{-1} = 7.8 \times 10^{10}$ cm$^3$.s$^{-1}$ and the two channels: HD + CO$_2^+$ → HCO$_2^+$ + D or HD + CO$_2^+$ → DCO$_2^+$ + H are equiprobable and then both channels have a coefficient equal to $3.9 \times 10^{10}$ cm$^3$.s$^{-1}$.

A similar approximation is used for the reaction between D and ions, while because the root square of mass ratio between H$_2$O and HDO is close to 1 (~0.97), we use the same coefficient for the two species.

Coefficients of the recombination of deuterated ions DCO$_2^+$, DCO$^+$ and HDO$^+$ are equal to the coefficients of the HCO$_2^+$, HCO$^+$ and H$_2$O$^+$ ions) and taken from the UMIST database as for the reaction between deuterated ions and neutrals species. The diverse measured values of the coefficient rate of the HCO$^+$ dissociative recombination and its temperature variations present important disagreement (e.g. Fox, 2015). Experimental measurements of the reaction rate of HCO$^+$ and DCO$^+$ between 150 and 270 K done by Korolov et al., 2009 indicate a possible 15% lower coefficient rate at 300 K for the DCO$^+$ recombination compared to the HCO$^+$ recombination which is inside the error bars of
equation of the evolution of D/H:

$$ \frac{d[D/H]}{dt} = \left[ \frac{D/H}{} \right] (t) \frac{1 - f(t)}{\tau_H(t)} \qquad (2) $$

where $\tau_H(t) = N_H(t)/\Phi_H(t)$ is the time scale of the hydrogen loss and $f(t)$ is the fractionation factor given by

$$ f(t) = \frac{\Phi_D(t)/\Phi_H(t)}{N_D(t)/N_H(t)} \qquad (3) $$

Two obvious conditions must be fulfilled to produce fractionation during escape (Eq. (2)). The time scale $\tau_H$ of the hydrogen escape rate should not be much larger than the planet's age, the fractionation factor should be different from one. If $f(t) = 1$ all along the planet history, then H and D escape in proportion of their relative abundance leading to no fractionation. Because the mass of H is lower than D, the atmospheric escape is more efficient for H than D and then $f(t) < 1$ most of the time leading to an increase of D/H with time. The time scale of the D/H fractionation is $\tau_H(t)/(1-f(t)) \le \tau_H(t)$ and would be equal to the hydrogen loss timescale if the escape of D was negligible ($f(t) = 0$).

The two reservoirs model can be expanded to a 3 reservoirs model (Krasnopolsky, 2000) where a third reservoir (reservoir 3) with the primordial D/H ratio of 1 SMOW is added to the previous system. The D/H ratio of this reservoir is denoted (D/H)3. This reservoir does not exchange water directly to the reservoir 1, but can supply the second reservoir (e.g. permanent ice caps) in water to maintain a constant hydrogen abundance in reservoir 2. The reservoir 2 can exchange (on long time scale) with the atmosphere and is replenished from reservoir 3. The evolution of the D/H ratio in reservoir 2 relative to the primordial D/H ratio in reservoir 3 ($r = (D/H)_2/(D/H)_3$) is solution of the differential equation (Krasnopolsky, 2000):

$$ \frac{dr}{dt} = \frac{1 - r(t)f(t)}{\tau_H(t)} \qquad (4) $$

where the evolution of r is still determined by the two parameters $\tau_H(t)$ and $f(t)$ defined previously. At $t = 0$, $r(t) = 1$ so if $f(t) = 1$ all the time, we still find no fractionation of the water as expected. But once the reservoir 2 has been fractionated ($r(t) > 1$), $f(t)$ should be smaller than $1/r(t)$ in order to continue having an increase of $r(t)$. If at a time t, $f(t) = 1/r(t)$, the D/H increase of reservoir 2 due to the preferential loss of H is balanced by its decrease due to the lower primordial D/H flux coming from reservoir 3.

Currently, hydrogen escape, which is dominated by Jeans escape, has been estimated from a diverse set of observations (Anderson and Hord, 1971; Chaufray et al., 2008; Chaffin et al., 2014, 2018; Halekas, 2017; Rahmati et al., 2018; Qin, 2021; Mayyasi et al., 2023) and modeling (Krasnopolsky, 2002, 2010, 2019; Chaufray et al., 2015, 2021a), showing large seasonal variations. Atomic deuterium is more difficult to measure but past measurements and more recent MAVEN measurements have shown also large seasonal variations of its concentration in the upper atmosphere (Clarke et al., 2017; Mayyasi et al., 2017; Chaufray et al., 2021b), and thus of its escape. But contrary to H, D escape is expected to be mostly non-thermal (Krasnopolsky, 2002, 2010) due to its larger mass. The main photochemical escape processes for D have been studied in details only recently (Cangi et al., 2023) using photochemical modeling of deuterium ion chemistry and Monte Carlo modeling of hot escape probabilities (Gregory et al., 2023a). These simulations have shown that DCO$^+$ recombination and the reaction between CO$_2^+$ and HD are the two most important sources of deuterium escape. In this study, we simulate the most important escape processes for the major form of H and D (H, H$_2$, D, HD, H$^+$, D$^+$) in the upper atmosphere of Mars including: thermal, photochemical and ions escape, the latter driven by the solar wind interaction with Mars' atmosphere. The photochemical escape includes not only the ion-neutral reactions

able to produce hot hydrogen and deuterium (Gregory et al., 2023b; Cangi et al., 2023) but also the collisions between hot oxygen produced from the O$_2^+$ dissociative recombination (e.g. Deighan et al., 2015, Leblanc et al., 2017, Lillis et al., 2017, Jakosky et al., 2018) with atomic and molecular hydrogen and deuterium (Gacesa et al., 2012; Shematovich, 2013). The different models used to estimate the different escape processes are described in Section 2 and the simulated values are presented in Section 3. These values of f(t) and $\tau_H(t)$ as well as the current main limits of our approach are discussed in Section 4 before a conclusion and perspectives of the study (Section 5).

## 2. Models

### 2.1. Model of the Martian atmosphere and the thermal escape

To compute the thermal escape of H, D, H$_2$ and HD, we use the Planetary Climate Model (PCM) in its Martian version (Mars-PCM). Mars-PCM is a 3D hydrodynamic model of the Martian atmosphere and ionosphere from the surface to the exobase near ~200 km, the altitude at which thermal escape is calculated (Forget et al., 1999; Gonzalez-Galindo et al., 2009, 2013, 2015). The vertical coordinate is the atmospheric pressure, and the upper altitude is not constant but can vary between ~200 to 250 km.

This model has been used in the past to compute the H and H$_2$ thermal escape (Chaufray et al., 2015, 2018, 2021). The spatial resolution is 5.625° in longitude, 3.75° latitude and a pressure vertical coordinate is used leading to a variable vertical resolution, which is typically of ~10 km in the thermosphere.

The main updates used in this study that affect the H and D species in the thermosphere, are the inclusion of detailed cloud microphysics that considers the details of the formation of clouds and their microphysics (nucleation, ice growth, scavenging) and their radiative effect, needed to simulate the supersaturation (Navarro et al., 2014; Vals et al., 2022; Naar, 2023), HDO as a tracer in its vapor and icy phases (Vals et al., 2022; Rossi et al., 2022) and the incorporation of the photochemistry of deuterated species, including ions (Table 1). The inclusion of the supersaturation makes the hygropause more porous, allowing H$_2$O and HDO to reach the upper atmosphere (Vals et al., 2022).

The coefficients of the reactions between HD and an ion i, $k_{HD,i}$, are approximated from the Langevin frequency (e.g. Krasnopolsky, 2002) assuming the same polarizability between HD and H$_2$ (Ishiguro et al., 1952), yet considering the difference in the mass, leading to (Krasnopolsky, 2002)

$$ k_{HD,i} = \sqrt{2/3k_{H2,i}} $$

These reactions can produce atomic H or D and we assume each channel to be equiprobable. For example, the important reaction H$_2$ + CO$_2^+$ has a reaction rate $9.5 \times 10^{10}$ cm$^3$.s$^{-1}$ in Mars-PCM (Table 1). The equivalent reaction with HD: HD + CO$_2^+$ has a total reaction rate of (2/3)$^{1/2} \times 9.5 \times 10^{10}$ cm$^3$.s$^{-1} = 7.8 \times 10^{10}$ cm$^3$.s$^{-1}$ and the two channels: HD + CO$_2^+$ → HCO$_2^+$ + D or HD + CO$_2^+$ → DCO$_2^+$ + H are equiprobable and then both channels have a coefficient equal to $3.9 \times 10^{10}$ cm$^3$.s$^{-1}$.

A similar approximation is used for the reaction between D and ions, while because the root square of mass ratio between H$_2$O and HDO is close to 1 (~0.97), we use the same coefficient for the two species.

Coefficients of the recombination of deuterated ions DCO$_2^+$, DCO$^+$ and HDO$^+$ are equal to the coefficients of the HCO$_2^+$, HCO$^+$ and H$_2$O$^+$ ions) and taken from the UMIST database as for the reaction between deuterated ions and neutrals species. The diverse measured values of the coefficient rate of the HCO$^+$ dissociative recombination and its temperature variations present important disagreement (e.g. Fox, 2015). Experimental measurements of the reaction rate of HCO$^+$ and DCO$^+$ between 150 and 270 K done by Korolov et al., 2009 indicate a possible 15% lower coefficient rate at 300 K for the DCO$^+$ recombination compared to the HCO$^+$ recombination which is inside the error bars of
<caption>**Table 1**</caption>

hydrogen ionospheric chemical reactions included in PCM-Mars and their reaction rates (in cm³·s⁻¹). The third column indicate if the reaction is included to describe deuterium or not, or if only one of the two possible products of the reaction is included.

<table><thead><tr><th>Reaction</th><th>Rate (cm³·s⁻¹)</th><th>Used for deuterium</th></tr></thead><tbody><tr><td>HCO<sub>2</sub><sup>+</sup> + e<sup>-</sup> → H + O + CO</td><td>8.1 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.64</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + e<sup>-</sup> → OH + CO</td><td>3.2 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.64</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + e<sup>-</sup> → H + CO<sub>2</sub></td><td>6.0 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.64</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + O → HCO<sup>+</sup> + O<sub>2</sub></td><td>1.0 × 10<sup>-9</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + CO → HCO<sup>+</sup> + CO<sub>2</sub></td><td>7.8 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>CO<sub>2</sub><sup>+</sup> + H → HCO<sup>+</sup> + O</td><td>4.5 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>H<sup>+</sup> + CO<sub>2</sub> → HCO<sup>+</sup> + O</td><td>3.5 × 10<sup>-9</sup></td><td>No</td></tr><tr><td>CO<sup>+</sup> + H<sub>2</sub> → HCO<sub>2</sub><sup>+</sup> + H</td><td>7.5 × 10<sup>-10</sup></td><td>Yes (2 channels)</td></tr><tr><td>HCO<sup>+</sup> + e<sup>-</sup> → CO + H</td><td>2.4 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.69</sup></td><td>Yes</td></tr><tr><td>CO<sub>2</sub><sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + CO<sub>2</sub></td><td>2.04 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>CO<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + CO</td><td>1.72 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>O<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + O</td><td>3.2 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>N<sub>2</sub><sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + N<sub>2</sub></td><td>2.3 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>N<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + N</td><td>2.8 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>H<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + H</td><td>6.9 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>HCO<sup>+</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + CO</td><td>2.5 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + O<sub>2</sub> → O<sub>2</sub><sup>+</sup> + H<sub>2</sub>O</td><td>4.6 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + CO → HCO<sup>+</sup> + OH</td><td>5.0 × 10<sup>-10</sup></td><td>Yes (2 channels)</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + O → O<sub>2</sub><sup>+</sup> + H<sub>2</sub></td><td>4.0 × 10<sup>-11</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + NO → NO<sup>+</sup> + H<sub>2</sub>O</td><td>2.7 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + OH</td><td>2.1 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + H<sub>2</sub> → H<sub>3</sub>O<sup>+</sup> + H</td><td>6.4 × 10<sup>-10</sup></td><td>Yes (1 channel only)</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → H + H + O</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → H + OH</td><td>8.6 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes (2 channels)</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub></td><td>3.9 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H + H</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → H<sub>2</sub>O + H</td><td>7.09 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H<sub>2</sub></td><td>5.37 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub> + H</td><td>5.6 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>O<sup>+</sup> + H<sub>2</sub> → OH<sup>+</sup> + H</td><td>1.7 × 10<sup>-9</sup></td><td>Yes (1 channel only)</td></tr><tr><td>OH<sup>+</sup> + O → O<sub>2</sub><sup>+</sup> + H</td><td>7.1 × 10<sup>-10</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + CO<sub>2</sub> → HCO<sub>2</sub><sup>+</sup> + O</td><td>1.44 × 10<sup>-9</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + CO → HCO<sup>+</sup> + O</td><td>1.05 × 10<sup>-9</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + NO → NO<sup>+</sup> + OH</td><td>3.59 × 10<sup>-10</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + H<sub>2</sub> → H<sub>2</sub>O<sup>+</sup> + H</td><td>1.01 × 10<sup>-9</sup></td><td>Yes (1 channel)</td></tr><tr><td>OH<sup>+</sup> + O<sub>2</sub> → O<sub>2</sub><sup>+</sup> + OH</td><td>5.9 × 10<sup>-10</sup></td><td>No</td></tr><tr><td>O<sup>+</sup> + H → H<sup>+</sup> + O</td><td>5.66 × 10<sup>-10</sup>(T<sub>e</sub>/300)<sup>0.36</sup>exp(8.6/T<sub>e</sub>)</td><td>No</td></tr><tr><td>H<sup>+</sup> + O → O<sup>+</sup> + H</td><td>6.86 × 10<sup>-10</sup>(T<sub>e</sub>/300)<sup>0.26</sup>exp(-224.3/T<sub>e</sub>)</td><td>No</td></tr><tr><td>CO<sub>2</sub><sup>+</sup> + H<sub>2</sub> → HCO<sub>2</sub><sup>+</sup> + H</td><td>9.5 × 10<sup>-10</sup></td><td>Yes (2 channels)</td></tr><tr><td>CO<sup>+</sup> + H → H<sup>+</sup> + CO</td><td>4.0 × 10<sup>-10</sup></td><td>No</td></tr></tbody></table>

the measured coefficients. They also find a slightly different temperature dependence (power $-1.3 \pm 0.2$ for HCO<sup>+</sup> dissociative recombination and $-1.1 \pm 0.2$ for DCO<sup>+</sup> dissociative recombination) which is different from the value of $-0.69$ given in the UMIST database. At Te = 1000 K (representative of the electron temperature near the exobase), the reaction rate of the HCO<sup>+</sup> recombination in our study is $\sim 1.0 \times 10^{-7}$ cm³·s⁻¹ while the expressions from Korolov et al. (2009) (used by Cangi et al., 2023), leads to $4.2 \times 10^{-8}$ cm³·s⁻¹ and the expressions used by Fox (2015) gives an intermediate value of $6 \times 10^{-8}$ cm³·s⁻¹. Using the coefficient rate from Korolov et al. (2009) should therefore reduce our

simulated production and escape rate by a factor $\sim 2$.

In the current version, only the two major deuterated ions DCO<sup>+</sup> and HCO<sub>2</sub><sup>+</sup>, which are needed to describe the hot atomic deuterium (Cangi et al., 2023), and one minor ion HDO<sup>+</sup> are included. The inclusion of other minor deuterated ions as D<sup>+</sup>, H<sub>2</sub>DO<sup>+</sup> and OD<sup>+</sup> should not change the results presented in this paper.

The simulation presented here is done for $L_s = 0°$, and the solar average scenario, supposed to be representative of the current average conditions (Gonzalez-Galindo et al., 2015). The ion dynamics presented in Chaufray et al. (2014) is included and extended to describe HCO<sup>+</sup>, DCO<sup>+</sup>, HCO<sub>2</sub><sup>+</sup> and DCO<sub>2</sub><sup>+</sup>. In the model, the computation of the molecular diffusion of light escaping species requires a boundary condition on the effusion velocity (Chaufray et al., 2015). We choose the Jeans effusion velocity for H, H<sub>2</sub> and D, and no escape for HD. This assumption is discussed in Section 4. The output of the model is a 3D field of atmospheric (composition, temperature) and ionospheric (composition) variables as well as the Jeans escape rate of the light species (even if the boundary condition for HD is a zero-effusion velocity, its Jeans escape rate can be calculated from the simulated temperature and density at the exobase). The ion and electron temperatures from Viking measurements are imposed (Chaufray et al., 2014). Electron and Ion temperatures measured by MAVEN/LPW (Ergun et al., 2021; Hanley et al., 2021) at low SZA and MAVEN/STATIC are not expected to affect our calculations appreciably, but should be included in the future to confirm this statement.

## 2.2. Model of the exosphere and the photochemical escape

The exosphere content is composed by thermal and non-thermal populations. The thermal population is the extension of the atmosphere resulting from the ballistic motions of the atoms and molecules, above the exobase, assuming a Maxwell-Boltzmann velocity distribution function at the exobase temperature. The non-thermal population is the population created by the photochemical reactions near the exobase with a large initial energy and not fully thermalized due to the small number of collisions near the exobase. This separation between thermal and non-thermal populations is a usual assumption for atomic oxygen on Mars but is more questionable for hydrogen and discussed in Section 4.2.

The 3D thermospheric neutral densities simulated with PCM-Mars are reinterpolated into an altitude grid and extended above the exobase (fixed at 200 km) using the formalism presented by Vidal-Madjar and Bertaux (1972) (see also appendix) to compute the density of the thermal population of the different species (H, O, CO<sub>2</sub>, D, H<sub>2</sub> and HD). This approach can lead to small discontinuities at the exobase (especially for light species) because the exospheric transport is not included in the calculation of the density at the exobase in the PCM-Mars simulations presented here. This assumption is discussed in Section 4.2.

The non-thermal components of the exosphere and the photochemical escape rates are computed using a 3D Monte Carlo model. In Monte Carlo simulations, contrary to hydrodynamic code, we use a particular approach, where the hot species are described by several particles (called test particles) each one representative of a large number of real particles. The background atmosphere is considered as a fluid characterized by its density, composition and temperature, not modified by the presence of these hot particles.

This model works in Mars-centered Solar Orbital (MSO) frame and tracks the hot particles, produced in the Martian upper atmosphere between 120 and 200 km. The hot particles are followed between 120 and 20,000 km.

The model has been used to simulate the hot oxygen produced by the dissociative recombination of O<sub>2</sub><sup>+</sup> ions and to derive the oxygen escape rate. The collisions of this hot population with the background atmosphere can produce other hotter species (CO<sub>2</sub>, N<sub>2</sub>, Ar) and were included to interpret the observations of a non-thermal argon population detected in the exosphere of Mars (Leblanc et al., 2019). In this paper, we have added the collisions of the hot oxygen with four new species: H, H<sub>2</sub>, D
<caption>**Table 1**</caption>

hydrogen ionospheric chemical reactions included in PCM-Mars and their reaction rates (in cm³·s⁻¹). The third column indicate if the reaction is included to describe deuterium or not, or if only one of the two possible products of the reaction is included.

<table><thead><tr><th>Reaction</th><th>Rate (cm³·s⁻¹)</th><th>Used for deuterium</th></tr></thead><tbody><tr><td>HCO<sub>2</sub><sup>+</sup> + e<sup>-</sup> → H + O + CO</td><td>8.1 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.64</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + e<sup>-</sup> → OH + CO</td><td>3.2 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.64</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + e<sup>-</sup> → H + CO<sub>2</sub></td><td>6.0 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.64</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + O → HCO<sup>+</sup> + O<sub>2</sub></td><td>1.0 × 10<sup>-9</sup></td><td>Yes</td></tr><tr><td>HCO<sub>2</sub><sup>+</sup> + CO → HCO<sup>+</sup> + CO<sub>2</sub></td><td>7.8 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>CO<sub>2</sub><sup>+</sup> + H → HCO<sup>+</sup> + O</td><td>4.5 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>H<sup>+</sup> + CO<sub>2</sub> → HCO<sup>+</sup> + O</td><td>3.5 × 10<sup>-9</sup></td><td>No</td></tr><tr><td>CO<sup>+</sup> + H<sub>2</sub> → HCO<sub>2</sub><sup>+</sup> + H</td><td>7.5 × 10<sup>-10</sup></td><td>Yes (2 channels)</td></tr><tr><td>HCO<sup>+</sup> + e<sup>-</sup> → CO + H</td><td>2.4 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.69</sup></td><td>Yes</td></tr><tr><td>CO<sub>2</sub><sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + CO<sub>2</sub></td><td>2.04 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>CO<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + CO</td><td>1.72 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>O<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + O</td><td>3.2 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>N<sub>2</sub><sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + N<sub>2</sub></td><td>2.3 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>N<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + N</td><td>2.8 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>H<sup>+</sup> + H<sub>2</sub>O → H<sub>2</sub>O<sup>+</sup> + H</td><td>6.9 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>HCO<sup>+</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + CO</td><td>2.5 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + O<sub>2</sub> → O<sub>2</sub><sup>+</sup> + H<sub>2</sub>O</td><td>4.6 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + CO → HCO<sup>+</sup> + OH</td><td>5.0 × 10<sup>-10</sup></td><td>Yes (2 channels)</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + O → O<sub>2</sub><sup>+</sup> + H<sub>2</sub></td><td>4.0 × 10<sup>-11</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + NO → NO<sup>+</sup> + H<sub>2</sub>O</td><td>2.7 × 10<sup>-10</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + H<sub>2</sub>O → H<sub>3</sub>O<sup>+</sup> + OH</td><td>2.1 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + H<sub>2</sub> → H<sub>3</sub>O<sup>+</sup> + H</td><td>6.4 × 10<sup>-10</sup></td><td>Yes (1 channel only)</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → H + H + O</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → H + OH</td><td>8.6 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes (2 channels)</td></tr><tr><td>H<sub>2</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub></td><td>3.9 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>Yes</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H + H</td><td>3.05 × 10<sup>-7</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → H<sub>2</sub>O + H</td><td>7.09 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → OH + H<sub>2</sub></td><td>5.37 × 10<sup>-8</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup> → O + H<sub>2</sub> + H</td><td>5.6 × 10<sup>-9</sup>(300/T<sub>e</sub>)<sup>0.5</sup></td><td>No</td></tr><tr><td>O<sup>+</sup> + H<sub>2</sub> → OH<sup>+</sup> + H</td><td>1.7 × 10<sup>-9</sup></td><td>Yes (1 channel only)</td></tr><tr><td>OH<sup>+</sup> + O → O<sub>2</sub><sup>+</sup> + H</td><td>7.1 × 10<sup>-10</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + CO<sub>2</sub> → HCO<sub>2</sub><sup>+</sup> + O</td><td>1.44 × 10<sup>-9</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + CO → HCO<sup>+</sup> + O</td><td>1.05 × 10<sup>-9</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + NO → NO<sup>+</sup> + OH</td><td>3.59 × 10<sup>-10</sup></td><td>No</td></tr><tr><td>OH<sup>+</sup> + H<sub>2</sub> → H<sub>2</sub>O<sup>+</sup> + H</td><td>1.01 × 10<sup>-9</sup></td><td>Yes (1 channel)</td></tr><tr><td>OH<sup>+</sup> + O<sub>2</sub> → O<sub>2</sub><sup>+</sup> + OH</td><td>5.9 × 10<sup>-10</sup></td><td>No</td></tr><tr><td>O<sup>+</sup> + H → H<sup>+</sup> + O</td><td>5.66 × 10<sup>-10</sup>(T<sub>e</sub>/300)<sup>0.36</sup>exp(8.6/T<sub>e</sub>)</td><td>No</td></tr><tr><td>H<sup>+</sup> + O → O<sup>+</sup> + H</td><td>6.86 × 10<sup>-10</sup>(T<sub>e</sub>/300)<sup>0.26</sup>exp(-224.3/T<sub>e</sub>)</td><td>No</td></tr><tr><td>CO<sub>2</sub><sup>+</sup> + H<sub>2</sub> → HCO<sub>2</sub><sup>+</sup> + H</td><td>9.5 × 10<sup>-10</sup></td><td>Yes (2 channels)</td></tr><tr><td>CO<sup>+</sup> + H → H<sup>+</sup> + CO</td><td>4.0 × 10<sup>-10</sup></td><td>No</td></tr></tbody></table>

the measured coefficients. They also find a slightly different temperature dependence (power $-1.3 \pm 0.2$ for HCO<sup>+</sup> dissociative recombination and $-1.1 \pm 0.2$ for DCO<sup>+</sup> dissociative recombination) which is different from the value of $-0.69$ given in the UMIST database. At Te = 1000 K (representative of the electron temperature near the exobase), the reaction rate of the HCO<sup>+</sup> recombination in our study is $\sim 1.0 \times 10^{-7}$ cm³·s⁻¹ while the expressions from Korolov et al. (2009) (used by Cangi et al., 2023), leads to $4.2 \times 10^{-8}$ cm³·s⁻¹ and the expressions used by Fox (2015) gives an intermediate value of $6 \times 10^{-8}$ cm³·s⁻¹. Using the coefficient rate from Korolov et al. (2009) should therefore reduce our

simulated production and escape rate by a factor $\sim 2$.

In the current version, only the two major deuterated ions DCO<sup>+</sup> and HCO<sub>2</sub><sup>+</sup>, which are needed to describe the hot atomic deuterium (Cangi et al., 2023), and one minor ion HDO<sup>+</sup> are included. The inclusion of other minor deuterated ions as D<sup>+</sup>, H<sub>2</sub>DO<sup>+</sup> and OD<sup>+</sup> should not change the results presented in this paper.

The simulation presented here is done for $L_s = 0°$, and the solar average scenario, supposed to be representative of the current average conditions (Gonzalez-Galindo et al., 2015). The ion dynamics presented in Chaufray et al. (2014) is included and extended to describe HCO<sup>+</sup>, DCO<sup>+</sup>, HCO<sub>2</sub><sup>+</sup> and DCO<sub>2</sub><sup>+</sup>. In the model, the computation of the molecular diffusion of light escaping species requires a boundary condition on the effusion velocity (Chaufray et al., 2015). We choose the Jeans effusion velocity for H, H<sub>2</sub> and D, and no escape for HD. This assumption is discussed in Section 4. The output of the model is a 3D field of atmospheric (composition, temperature) and ionospheric (composition) variables as well as the Jeans escape rate of the light species (even if the boundary condition for HD is a zero-effusion velocity, its Jeans escape rate can be calculated from the simulated temperature and density at the exobase). The ion and electron temperatures from Viking measurements are imposed (Chaufray et al., 2014). Electron and Ion temperatures measured by MAVEN/LPW (Ergun et al., 2021; Hanley et al., 2021) at low SZA and MAVEN/STATIC are not expected to affect our calculations appreciably, but should be included in the future to confirm this statement.

## 2.2. Model of the exosphere and the photochemical escape

The exosphere content is composed by thermal and non-thermal populations. The thermal population is the extension of the atmosphere resulting from the ballistic motions of the atoms and molecules, above the exobase, assuming a Maxwell-Boltzmann velocity distribution function at the exobase temperature. The non-thermal population is the population created by the photochemical reactions near the exobase with a large initial energy and not fully thermalized due to the small number of collisions near the exobase. This separation between thermal and non-thermal populations is a usual assumption for atomic oxygen on Mars but is more questionable for hydrogen and discussed in Section 4.2.

The 3D thermospheric neutral densities simulated with PCM-Mars are reinterpolated into an altitude grid and extended above the exobase (fixed at 200 km) using the formalism presented by Vidal-Madjar and Bertaux (1972) (see also appendix) to compute the density of the thermal population of the different species (H, O, CO<sub>2</sub>, D, H<sub>2</sub> and HD). This approach can lead to small discontinuities at the exobase (especially for light species) because the exospheric transport is not included in the calculation of the density at the exobase in the PCM-Mars simulations presented here. This assumption is discussed in Section 4.2.

The non-thermal components of the exosphere and the photochemical escape rates are computed using a 3D Monte Carlo model. In Monte Carlo simulations, contrary to hydrodynamic code, we use a particular approach, where the hot species are described by several particles (called test particles) each one representative of a large number of real particles. The background atmosphere is considered as a fluid characterized by its density, composition and temperature, not modified by the presence of these hot particles.

This model works in Mars-centered Solar Orbital (MSO) frame and tracks the hot particles, produced in the Martian upper atmosphere between 120 and 200 km. The hot particles are followed between 120 and 20,000 km.

The model has been used to simulate the hot oxygen produced by the dissociative recombination of O<sub>2</sub><sup>+</sup> ions and to derive the oxygen escape rate. The collisions of this hot population with the background atmosphere can produce other hotter species (CO<sub>2</sub>, N<sub>2</sub>, Ar) and were included to interpret the observations of a non-thermal argon population detected in the exosphere of Mars (Leblanc et al., 2019). In this paper, we have added the collisions of the hot oxygen with four new species: H, H<sub>2</sub>, D
and HD to estimate their escape rate induced by these collisions.

Other chemical reactions in the Martian upper atmosphere produce atomic hydrogen and deuterium at high energy. When these reactions occur deeply below the exobase, all the atoms produced are thermalized by collisions with the main species (CO₂, O). However, a small fraction of atoms produced slightly below or above the exobase will not be fully thermalized and may have an energy larger than the escape energy. This small fraction can be neglected in the simulation of the density below the exobase as assumed in the PCM-Mars, but for species like atomic D it can be a non-negligible fraction of the escape (e.g. Cangi et al., 2023). The photochemical reactions included in our models to simulate the non-thermal escape of H and D are the main chemical source of atomic hydrogen and deuterium in the ionosphere (Krasnopolsky, 2002, 2019; Fox, 2015; Gregory et al., 2023b):

$$H_2 + CO_2^+ \rightarrow HCO_2^+ + H + \Delta E = 1 \text{ eV} \quad r_1 = 9.5 \times 10^{-10} \text{ cm}^3 \cdot \text{s}^{-1} \quad (R1)$$

$$HD + CO_2^+ \rightarrow HCO_2^+ + D + \Delta E = 1 \text{ eV} \quad r_2 = 3.9 \times 10^{-10} \text{ cm}^3 \cdot \text{s}^{-1} \quad (R2)$$

$$\begin{aligned} HCO^+ + e^- \rightarrow CO(X^1\Sigma^+) + H + \Delta E &= 7.31 \text{ eV} \quad r_{3a} \\ &= p_{3a} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad (R3a)$$

$$\begin{aligned} HCO^+ + e^- \rightarrow CO(a^3\Pi_r) + H + \Delta E &= 1.30 \text{ eV} \quad r_{3b} \\ &= p_{3b} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad (R3b)$$

$$\begin{aligned} DCO^+ + e^- \rightarrow CO(X^1\Sigma^+) + D + \Delta E &= 7.31 \text{ eV} \quad r_{4a} \\ &= p_{3a} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad (R4a)$$

$$\begin{aligned} DCO^+ + e^- \rightarrow CO(a^3\Pi_r) + D + \Delta E &= 1.30 \text{ eV} \quad r_{4b} \\ &= p_{4b} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad (R4b)$$

The rates of the two reactions producing H atoms (R1 and R3) are from Krasnopolsky (2019). For R3, we consider the two channels 3a and 3b with a probability $p_{3a}$ and $p_{3b} = 0.77$ and 0.23 respectively (Gregory et al., 2023a). The reaction rate for R2 is the one of R1 multiplied by 0.82, based on the Langevin formula (Krasnopolsky, 2002) and divided by 2 because we assumed that the other channel leading to ($DCO_2^+ + H$), not included in our simulation is equally possible. For R4, we consider two channels as for R3 with $p_{4a} = p_{3a}$ and $p_{4b} = p_{3b}$.

The inputs for the exospheric model are the density of the ions needed to calculate the production rates of the hot oxygen ($O_2^+$) and the hot hydrogen and deuterium ($CO_2^+$, $HCO^+$ and $DCO^+$), as well as the electron density and temperature. The 3D ion density is taken from the PCM-Mars (see Section 2.1) from 120 km to 200 km.

The main species needed to describe the sources and the energy degradation of the hot atoms are CO₂, N₂, Ar, O, H₂, H, D and HD. The neutral density of the main atmospheric species (CO₂, N₂, O, Ar) is extrapolated above the top boundary of Mars-PCM up to 20,000 km by a simple 1D hydrostatic law at the local (latitude, local time) exospheric temperature and considering the mass of each species. This approach is sufficient to describe accurately the density of these heavy species up to altitudes where their density become negligible (e.g. Kallio et al., 2011). This approach is less valid for light species (e.g. Kallio et al., 2011), and the density of the light species (H, H₂, D and HD) are calculated using the kinetic approach from Vidal-Madjar and Bertaux (1972) described partly in appendix. The collisions of the hot oxygen atoms with all the species are calculated at all altitudes between 120 and 20,000 km.

The electron temperature is the one used as input in Mars-PCM (see section 2.1). Because the ions and electron density and temperatures are not extrapolated, the production of hot species is limited below 200 km. The test particles produced between 120 and 200 km are followed in the full simulated domain between 120 and 20,000 km. For the simulations of the reactions R1, R2, R3 and R4, only the test particles with a kinetic energy larger than the potential energy divided by 10 are described (i.e with a velocity $v < v_{esc}/3.2$). The test particle is considered as

thermalized if its energy is lower than this limit. This choice is arbitrary but does not impact the simulated escape rate.

To simulate the hot hydrogen and deuterium produced by collisions of thermal H and D with the hot oxygen, because the number of simulated test particles is very large, we only considered test particles (whatever the species) with an energy larger than the local hydrogen escape energy multiplied by 0.8 (~ 0.1 eV at 200 km). Collisions of hot oxygen with CO₂, O (cold), N₂, Ar and the studied species (H, D, H₂ and HD) are considered. Given the large number of possible collisions, the collision cross sections of the different hot population and atmospheric species included in this study are taken from the simple universal model of Lewkow and Kharchenko (2014) and collisions between hot particles are neglected.

For all mechanisms, we use a temporal loop and eject regularly a constant number of test particles of different weights (the number of real particles represented by a test particle). This weight is proportional to the local production rate. The particles are removed from the simulation if they are thermalized, if they reach the lower boundary (120 km) or the upper boundary (20,000 km) of the model. Only the particles reaching the upper boundary with a velocity larger than the escape velocity are considered to compute the escape rate. The number of particles reaching the upper boundary with a velocity lower than the escape velocity is much smaller than the escaping particles. The simulated escape flux is given by the sum of the weights of the escaping test particles divided by the simulation time. We perform long simulations, such that at the end of the simulations, the simulated escape flux doesn't vary significantly and the total weight of the test particles remaining in the simulation divided by the simulation time is much lower (< 1%) than the derived escape rate. So, even if these remaining test particles were escaping, their contribution to the estimated escape rate would be negligible.

## 2.3. Model of the solar wind interaction with Mars

To simulate the solar wind interaction with the Martian upper atmosphere, we used the LatHys model (Modolo et al., 2016). LatHys is a parallel 3D, multispecies hybrid code, where ions are described as individual charged particles and electrons as a mass-less charge neutralizing fluid. In the hybrid approach, ions trajectory is computed from the equation of motion, including only the Lorenz force, while electric and magnetic (EM) fields are derived from the Maxwell equations (e.g. Kallio et al., 2011). The 3D Martian atmosphere simulated with Mars-PCM and extended into the exosphere (CO₂, O, H densities) was used as input of the model. At this stage the D corona is not included in the LatHyS simulations since the D⁺ production is not expected to play a role on the global interaction.

The hot oxygen corona is not considered in this simulation since its effect on the electromagnetic environment and on the H⁺ and D⁺ escape rate is expected to be of only few % (Dong et al., 2018). The electromagnetic environment is computed on a cartesian MSO grid with a spatial resolution of 80 km. The size of the simulation box is [-2.4; +2.4] Martian radius along X (Mars - Sun direction) and [-4.6; +4.6] Martian radius along Y and Z directions. The model is used to compute the electric and magnetic fields around Mars only. The solar wind (SW) and interplanetary magnetic field (IMF) parameters used in this simulation are summarized in Table 2, and are typical values at the Martian orbit. Crustal magnetic fields are not considered in the performed simulations to not complexify the global SW - Mars interaction.

In a second step, the production of H⁺ and D⁺ planetary ions by photoionization is computed as the products of the local neutral (H or D)

<table><caption>**Table 2**

Parameters of the solar wind and magnetic field used to simulate the interaction between Mars and the solar wind.</caption><thead><tr><td>SW Density</td><td>SW Velocity</td><td>IMF magnitude</td><td>IMF direction</td></tr></thead><tbody><tr><td>2.7 cm<sup>-3</sup></td><td>480 km</td><td>3 nT</td><td>Parker (57°)</td></tr></tbody></table>
and HD to estimate their escape rate induced by these collisions.

Other chemical reactions in the Martian upper atmosphere produce atomic hydrogen and deuterium at high energy. When these reactions occur deeply below the exobase, all the atoms produced are thermalized by collisions with the main species (CO₂, O). However, a small fraction of atoms produced slightly below or above the exobase will not be fully thermalized and may have an energy larger than the escape energy. This small fraction can be neglected in the simulation of the density below the exobase as assumed in the PCM-Mars, but for species like atomic D it can be a non-negligible fraction of the escape (e.g. Cangi et al., 2023). The photochemical reactions included in our models to simulate the non-thermal escape of H and D are the main chemical source of atomic hydrogen and deuterium in the ionosphere (Krasnopolsky, 2002, 2019; Fox, 2015; Gregory et al., 2023b):

$$H_2 + CO_2^+ \rightarrow HCO_2^+ + H + \Delta E = 1 \text{ eV} \quad r_1 = 9.5 \times 10^{-10} \text{ cm}^3 \cdot \text{s}^{-1} \quad \text{(R1)}$$

$$HD + CO_2^+ \rightarrow HCO_2^+ + D + \Delta E = 1 \text{ eV} \quad r_2 = 3.9 \times 10^{-10} \text{ cm}^3 \cdot \text{s}^{-1} \quad \text{(R2)}$$

$$\begin{aligned} HCO^+ + e^- \rightarrow CO(X^1\Sigma^+) + H + \Delta E &= 7.31 \text{ eV} \quad r_{3a} \\ &= p_{3a} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad \text{(R3a)}$$

$$\begin{aligned} HCO^+ + e^- \rightarrow CO(a^3\Pi_r) + H + \Delta E &= 1.30 \text{ eV} \quad r_{3b} \\ &= p_{3b} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad \text{(R3b)}$$

$$\begin{aligned} DCO^+ + e^- \rightarrow CO(X^1\Sigma^+) + D + \Delta E &= 7.31 \text{ eV} \quad r_{4a} \\ &= p_{3a} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad \text{(R4a)}$$

$$\begin{aligned} DCO^+ + e^- \rightarrow CO(a^3\Pi_r) + D + \Delta E &= 1.30 \text{ eV} \quad r_{4b} \\ &= p_{4b} \times 2.4 \times 10^{-7} \text{ 300/T}_e \end{aligned} \quad \text{(R4b)}$$

The rates of the two reactions producing H atoms (R1 and R3) are from Krasnopolsky (2019). For R3, we consider the two channels 3a and 3b with a probability $p_{3a}$ and $p_{3b} = 0.77$ and 0.23 respectively (Gregory et al., 2023a). The reaction rate for R2 is the one of R1 multiplied by 0.82, based on the Langevin formula (Krasnopolsky, 2002) and divided by 2 because we assumed that the other channel leading to ($DCO_2^+ + H$), not included in our simulation is equally possible. For R4, we consider two channels as for R3 with $p_{4a} = p_{3a}$ and $p_{4b} = p_{3b}$.

The inputs for the exospheric model are the density of the ions needed to calculate the production rates of the hot oxygen ($O_2^+$) and the hot hydrogen and deuterium ($CO_2^+$, $HCO^+$ and $DCO^+$), as well as the electron density and temperature. The 3D ion density is taken from the PCM-Mars (see Section 2.1) from 120 km to 200 km.

The main species needed to describe the sources and the energy degradation of the hot atoms are CO₂, N₂, Ar, O, H₂, H, D and HD. The neutral density of the main atmospheric species (CO₂, N₂, O, Ar) is extrapolated above the top boundary of Mars-PCM up to 20,000 km by a simple 1D hydrostatic law at the local (latitude, local time) exospheric temperature and considering the mass of each species. This approach is sufficient to describe accurately the density of these heavy species up to altitudes where their density become negligible (e.g. Kallio et al., 2011). This approach is less valid for light species (e.g. Kallio et al., 2011), and the density of the light species (H, H₂, D and HD) are calculated using the kinetic approach from Vidal-Madjar and Bertaux (1972) described partly in appendix. The collisions of the hot oxygen atoms with all the species are calculated at all altitudes between 120 and 20,000 km.

The electron temperature is the one used as input in Mars-PCM (see section 2.1). Because the ions and electron density and temperatures are not extrapolated, the production of hot species is limited below 200 km. The test particles produced between 120 and 200 km are followed in the full simulated domain between 120 and 20,000 km. For the simulations of the reactions R1, R2, R3 and R4, only the test particles with a kinetic energy larger than the potential energy divided by 10 are described (i.e with a velocity $v < v_{esc}/3.2$). The test particle is considered as

thermalized if its energy is lower than this limit. This choice is arbitrary but does not impact the simulated escape rate.

To simulate the hot hydrogen and deuterium produced by collisions of thermal H and D with the hot oxygen, because the number of simulated test particles is very large, we only considered test particles (whatever the species) with an energy larger than the local hydrogen escape energy multiplied by 0.8 (~0.1 eV at 200 km). Collisions of hot oxygen with CO₂, O (cold), N₂, Ar and the studied species (H, D, H₂ and HD) are considered. Given the large number of possible collisions, the collision cross sections of the different hot population and atmospheric species included in this study are taken from the simple universal model of Lewkow and Kharchenko (2014) and collisions between hot particles are neglected.

For all mechanisms, we use a temporal loop and eject regularly a constant number of test particles of different weights (the number of real particles represented by a test particle). This weight is proportional to the local production rate. The particles are removed from the simulation if they are thermalized, if they reach the lower boundary (120 km) or the upper boundary (20,000 km) of the model. Only the particles reaching the upper boundary with a velocity larger than the escape velocity are considered to compute the escape rate. The number of particles reaching the upper boundary with a velocity lower than the escape velocity is much smaller than the escaping particles. The simulated escape flux is given by the sum of the weights of the escaping test particles divided by the simulation time. We perform long simulations, such that at the end of the simulations, the simulated escape flux doesn't vary significantly and the total weight of the test particles remaining in the simulation divided by the simulation time is much lower (<1%) than the derived escape rate. So, even if these remaining test particles were escaping, their contribution to the estimated escape rate would be negligible.

## 2.3. Model of the solar wind interaction with Mars

To simulate the solar wind interaction with the Martian upper atmosphere, we used the LatHys model (Modolo et al., 2016). LatHys is a parallel 3D, multispecies hybrid code, where ions are described as individual charged particles and electrons as a mass-less charge neutralizing fluid. In the hybrid approach, ions trajectory is computed from the equation of motion, including only the Lorenz force, while electric and magnetic (EM) fields are derived from the Maxwell equations (e.g. Kallio et al., 2011). The 3D Martian atmosphere simulated with Mars-PCM and extended into the exosphere (CO₂, O, H densities) was used as input of the model. At this stage the D corona is not included in the LatHyS simulations since the D⁺ production is not expected to play a role on the global interaction.

The hot oxygen corona is not considered in this simulation since its effect on the electromagnetic environment and on the H⁺ and D⁺ escape rate is expected to be of only few % (Dong et al., 2018). The electromagnetic environment is computed on a cartesian MSO grid with a spatial resolution of 80 km. The size of the simulation box is [-2.4; +2.4] Martian radius along X (Mars - Sun direction) and [-4.6; +4.6] Martian radius along Y and Z directions. The model is used to compute the electric and magnetic fields around Mars only. The solar wind (SW) and interplanetary magnetic field (IMF) parameters used in this simulation are summarized in Table 2, and are typical values at the Martian orbit. Crustal magnetic fields are not considered in the performed simulations to not complexify the global SW - Mars interaction.

In a second step, the production of H⁺ and D⁺ planetary ions by photoionization is computed as the products of the local neutral (H or D)

<table><caption>**Table 2**
Parameters of the solar wind and magnetic field used to simulate the interaction between Mars and the solar wind.</caption><thead><tr><td>SW Density</td><td>SW Velocity</td><td>IMF magnitude</td><td>IMF direction</td></tr></thead><tbody><tr><td>2.7 cm<sup>-3</sup></td><td>480 km</td><td>3 nT</td><td>Parker (57°)</td></tr></tbody></table>
density and the ionization frequency $v = 3.1 \times 10^{-8} s^{-1}$ for both species. The production of $H^+$ and $D^+$ planetary ions by charge exchange is computed with a test particles code in which $12 \times 10^6$ test particles representative of the solar wind protons are followed in the LatHyS Martian electromagnetic fields. The simulated production of planetary ions ($H^+$ and $D^+$) is recorded on a spherical grid.

The escape rates of $H^+$ and $D^+$, produced either by photoionization or by charge exchange are also computed using the test particle code. For each process and each species, $12 \times 10^6$ test particles (planetary $H^+$ and $D^+$) are followed until they escape from the simulation box or reach the lower boundary at 200 km.

# 3. Results

## 3.1. Thermal population and Jeans escape rate

### 3.1.1. Horizontal distribution near the exobase

The simulated temperature near the exobase (200 km) shows the expected day/night asymmetry (Fig. 1). For this average solar condition scenario, the dayside exospheric temperature is ~210 K and ~150 K on the nightside.

The atomic hydrogen and deuterium density at the exobase (Fig. 2), as well as $H_2$ and HD (not shown) present a pattern similar to those simulated by Chaufray et al. (2015) at $Ls = 0^{\circ}-30^{\circ}$, with a dawn/dusk asymmetry and a density bulge at 4–5 am. This enhancement of light species results from the convergence and downwelling of the thermospheric winds in the morning nightside. Because these thermospheric winds are enriched in light species, they produce a local bulge in light species (e.g. Elrod et al., 2017). Such an enhancement has been observed for helium and molecular hydrogen by MAVEN (Elrod et al., 2017; Stone et al., 2022) and is expected for other light species.

For H, the exospheric ballistic transport can reduce this enhancement by redistributing the hydrogen at the exobase from the nightside to the dayside (Chaufray et al., 2018) and leading to a distribution closer to the theoretical solution $nT^{5/2} = \text{Constant}$ (Hodges and Johnson, 1968). Such an exospheric transport is numerically difficult to extend to heavier species like D or $H_2$ because the computation of the net flux (i.e. the difference between the upward and downward exospheric flux) must be more accurate when the species is heavier. Therefore, it is not included in these simulations for any species (see Section 4).

<figure><img src="image_3.png" /><figcaption>Fig. 1. PCM generated map of the simulated temperature at 200 km (near the exobase) in Mars body-fixed rotating frame for average solar condition. Noon is at longitude = 0°, dawn at longitude = 270°, dusk at longitude = 90° and midnight at longitude = 180°. The north pole is at co-latitude = 0° and south pole at co-latitude = 180°.</figcaption></figure>

### 3.1.2. Thermospheric and ionospheric density vertical profiles

Mars-PCM was used to simulate the neutral and ion density in the Martian thermosphere and ionosphere at $Ls = 0^{\circ}$ for an average solar activity. The density of H, D, $H_2$, HD, O and $CO_2$ (the main neutral species useful for this work) are shown in Fig. 3 (left panel), and the density of electrons, $O_2^+$, $CO_2^+$, $HCO^+$, and $DCO^+$ ions (the main ion species useful for this work) are shown in Fig. 3 (right panel).

The vertical variations of the neutral species are controlled by the local temperature (~200 K at the subsolar point in this simulation). The simulated density of $H_2$ near 200 km is ~2 × $10^5$ cm$^{-3}$, while the transition between the $CO_2$ dominant species to the O dominant species is at ~200 km (when the $CO_2$ density is ~10$^7$ cm$^{-3}$), in agreement with the observations by MAVEN/NGIMS (Stone et al., 2022).

The simulated $HCO^+$ density is in the range of the density simulated by Fox (2015), Krasnopolsky (2019) and Cangi et al. (2023) and observed by MAVEN/NGIMS (Benna et al., 2015), but the altitude of the peak is slightly below (170 km) and less pronounced than the peak simulated by Fox (2015) and Krasnopolsky (2019) while in better agreement with Cangi et al. (2023) who find also a secondary peak around 170 km with a density ~$10^3$ cm$^{-3}$. This difference may be due to the limited altitude range of Mars-PCM and Cangi et al. (2023) compared to the models by Fox (2015) and Krasnopolsky (2019) but should not modify our conclusions on the relative importance of the different escape processes.

### 3.1.3. Thermal escape

The thermal escape of H, $H_2$, D and HD is computed by integration of the Jeans equation over the exobase as done by Chaufray et al. (2015). The simulated escape rates of the four species due to the different simulated processes are given in Table 3 as well as the effusion velocities. The effusion velocity $u_{eff}$ is a 1D concept (e.g. Chassefiere and Leblanc, 2004; Chaufray, 2021) and in our 3D simulation, we will define it, for all processes, as the simulated escape rate $\Phi$ divided by the average density <n>$_{200}$ of the species at 200 km (calculated with PCM-Mars) and by the surface of the exobase $S_{200}$ in order to compare our values to those reported by Krasnopolsky (2010):

$$u_{eff} = \frac{\Phi}{S_{200}\langle n \rangle_{200}}$$

The thermal effusion velocity is controlled by the temperature at the exobase and our values are close to the derived values by Krasnopolsky (2010) for T = 200 K ($3.9 \times 10^2$, $4.8 \times 10^{-1}$, $4.8 \times 10^{-1}$, and $5 \times 10^{-4}$ cm/s for H, $H_2$, D and HD). Note that for a 1D model, the effusion velocity of $H_2$ and D is equal due to their equal mass while in our simulation, slight differences in the horizontal distribution of the density of the two species lead to a small (<5%) difference in their effusion velocity.

## 3.2. Photochemical escape

### 3.2.1. Production by dissociative recombination of $HCO^+$ and $DCO^+$

The simulated production rate of H and D by dissociative recombination of $HCO^+$ and $DCO^+$ (noted $XCO^+$ with X = H or D) above 80 km are $6.1 \times 10^{25}$ s$^{-1}$ and $2.8 \times 10^{22}$ s$^{-1}$. Our simulated production rate of hot H is in good agreement with the production calculated by Gregory et al. (2023a) above 80 km. Because almost all the atoms produced between 80 and 120 km are thermalized by collisions and their probability to escape is null (Gregory et al., 2023a), we only simulate the production above 120 km. The simulated production of hot H and D above 120 km is $3.0 \times 10^{25}$ and $1.2 \times 10^{22}$ s$^{-1}$ respectively. The simulated escape rates and the associated effusion velocity are given in Table 4.

Half of the hydrogen atoms produced by this reaction above 120 km escapes, while one third of the deuterium atoms escape. The difference in the efficiency is due to the mass of the atom and the larger energy needed for D to escape. The simulated effusion velocity is lower than the
density and the ionization frequency $\nu = 3.1 \times 10^{-8} s^{-1}$ for both species. The production of $H^+$ and $D^+$ planetary ions by charge exchange is computed with a test particles code in which $12 \times 10^6$ test particles representative of the solar wind protons are followed in the LatHyS Martian electromagnetic fields. The simulated production of planetary ions ($H^+$ and $D^+$) is recorded on a spherical grid.

The escape rates of $H^+$ and $D^+$, produced either by photoionization or by charge exchange are also computed using the test particle code. For each process and each species, $12 \times 10^6$ test particles (planetary $H^+$ and $D^+$) are followed until they escape from the simulation box or reach the lower boundary at 200 km.

# 3. Results

## 3.1. Thermal population and Jeans escape rate

### 3.1.1. Horizontal distribution near the exobase

The simulated temperature near the exobase (200 km) shows the expected day/night asymmetry (Fig. 1). For this average solar condition scenario, the dayside exospheric temperature is ~210 K and ~150 K on the nightside.

The atomic hydrogen and deuterium density at the exobase (Fig. 2), as well as $H_2$ and HD (not shown) present a pattern similar to those simulated by Chaufray et al. (2015) at $Ls = 0^{\circ}-30^{\circ}$, with a dawn/dusk asymmetry and a density bulge at 4–5 am. This enhancement of light species results from the convergence and downwelling of the thermospheric winds in the morning nightside. Because these thermospheric winds are enriched in light species, they produce a local bulge in light species (e.g. Elrod et al., 2017). Such an enhancement has been observed for helium and molecular hydrogen by MAVEN (Elrod et al., 2017; Stone et al., 2022) and is expected for other light species.

For H, the exospheric ballistic transport can reduce this enhancement by redistributing the hydrogen at the exobase from the nightside to the dayside (Chaufray et al., 2018) and leading to a distribution closer to the theoretical solution $nT^{5/2} = \text{Constant}$ (Hodges and Johnson, 1968). Such an exospheric transport is numerically difficult to extend to heavier species like D or $H_2$ because the computation of the net flux (i.e. the difference between the upward and downward exospheric flux) must be more accurate when the species is heavier. Therefore, it is not included in these simulations for any species (see Section 4).

<figure><img src="image_3.png" /><figcaption>Fig. 1. PCM generated map of the simulated temperature at 200 km (near the exobase) in Mars body-fixed rotating frame for average solar condition. Noon is at longitude = 0°, dawn at longitude = 270°, dusk at longitude = 90° and midnight at longitude = 180°. The north pole is at co-latitude = 0° and south pole at co-latitude = 180°.</figcaption></figure>

### 3.1.2. Thermospheric and ionospheric density vertical profiles

Mars-PCM was used to simulate the neutral and ion density in the Martian thermosphere and ionosphere at $Ls = 0^{\circ}$ for an average solar activity. The density of H, D, $H_2$, HD, O and $CO_2$ (the main neutral species useful for this work) are shown in Fig. 3 (left panel), and the density of electrons, $O_2^+$, $CO_2^+$, $HCO^+$, and $DCO^+$ ions (the main ion species useful for this work) are shown in Fig. 3 (right panel).

The vertical variations of the neutral species are controlled by the local temperature (~200 K at the subsolar point in this simulation). The simulated density of $H_2$ near 200 km is ~2 $\times 10^5$ cm$^{-3}$, while the transition between the $CO_2$ dominant species to the O dominant species is at ~200 km (when the $CO_2$ density is ~10$^7$ cm$^{-3}$), in agreement with the observations by MAVEN/NGIMS (Stone et al., 2022).

The simulated $HCO^+$ density is in the range of the density simulated by Fox (2015), Krasnopolsky (2019) and Cangi et al. (2023) and observed by MAVEN/NGIMS (Benna et al., 2015), but the altitude of the peak is slightly below (170 km) and less pronounced than the peak simulated by Fox (2015) and Krasnopolsky (2019) while in better agreement with Cangi et al. (2023) who find also a secondary peak around 170 km with a density ~$10^3$ cm$^{-3}$. This difference may be due to the limited altitude range of Mars-PCM and Cangi et al. (2023) compared to the models by Fox (2015) and Krasnopolsky (2019) but should not modify our conclusions on the relative importance of the different escape processes.

### 3.1.3. Thermal escape

The thermal escape of H, $H_2$, D and HD is computed by integration of the Jeans equation over the exobase as done by Chaufray et al. (2015). The simulated escape rates of the four species due to the different simulated processes are given in Table 3 as well as the effusion velocities. The effusion velocity $u_{eff}$ is a 1D concept (e.g. Chassefiere and Leblanc, 2004; Chaufray, 2021) and in our 3D simulation, we will define it, for all processes, as the simulated escape rate $\Phi$ divided by the average density $\langle n \rangle_{200}$ of the species at 200 km (calculated with PCM-Mars) and by the surface of the exobase $S_{200}$ in order to compare our values to those reported by Krasnopolsky (2010):

$$u_{eff} = \frac{\Phi}{S_{200}\langle n \rangle_{200}}$$

The thermal effusion velocity is controlled by the temperature at the exobase and our values are close to the derived values by Krasnopolsky (2010) for $T = 200$ K ($3.9 \times 10^2$, $4.8 \times 10^{-1}$, $4.8 \times 10^{-1}$, and $5 \times 10^{-4}$ cm/s for H, $H_2$, D and HD). Note that for a 1D model, the effusion velocity of $H_2$ and D is equal due to their equal mass while in our simulation, slight differences in the horizontal distribution of the density of the two species lead to a small (<5%) difference in their effusion velocity.

## 3.2. Photochemical escape

### 3.2.1. Production by dissociative recombination of $HCO^+$ and $DCO^+$

The simulated production rate of H and D by dissociative recombination of $HCO^+$ and $DCO^+$ (noted $XCO^+$ with X = H or D) above 80 km are $6.1 \times 10^{25}$ s$^{-1}$ and $2.8 \times 10^{22}$ s$^{-1}$. Our simulated production rate of hot H is in good agreement with the production calculated by Gregory et al. (2023a) above 80 km. Because almost all the atoms produced between 80 and 120 km are thermalized by collisions and their probability to escape is null (Gregory et al., 2023a), we only simulate the production above 120 km. The simulated production of hot H and D above 120 km is $3.0 \times 10^{25}$ and $1.2 \times 10^{22}$ s$^{-1}$ respectively. The simulated escape rates and the associated effusion velocity are given in Table 4.

Half of the hydrogen atoms produced by this reaction above 120 km escapes, while one third of the deuterium atoms escape. The difference in the efficiency is due to the mass of the atom and the larger energy needed for D to escape. The simulated effusion velocity is lower than the
<figure><img src="image_4.png" /><figcaption>Fig. 2. PCM generated maps of the simulated hydrogen density (left panel) and deuterium density (right panel) in Mars body-fixed rotating frame at 200 km.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Left: PCM simulated density of several neutral species at subsolar point from 120 km to the top of the atmosphere (~ 225 km). Right: PCM simulated density of several ions at subsolar point from 120 km to the top of the atmosphere (~ 225 km).</figcaption></figure>

<table><caption>**Table 3**
PCM Simulated thermal escape rates and effusion velocity for H, D, H₂ and HD.</caption><thead><tr><th></th><th>H</th><th>D</th><th>H<sub>2</sub></th><th>HD</th></tr></thead><tbody><tr><td>Thermal escape</td><td></td><td></td><td></td><td></td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>2.1 &times; 10<sup>26</sup></td><td>1.3 &times; 10<sup>20</sup></td><td>1.9 &times; 10<sup>23</sup></td><td>1.8 &times; 10<sup>17</sup></td></tr><tr><td>Effusion velocities (cm. s<sup>-1</sup>)</td><td>3.5 &times; 10<sup>2</sup></td><td>4.1 &times; 10<sup>-1</sup></td><td>4.0 &times; 10<sup>-1</sup></td><td>7.1 &times; 10<sup>-4</sup></td></tr></tbody></table>

Jeans effusion velocity by more than one order of magnitude for H, but larger by more than one order of magnitude for D. The simulated H escape rate is in good agreement with the H escape rate simulated by Gregory et al. (2023a) for their low solar activity case. Considering we simulate a similar production and used comparable elastic cross sections between H and CO₂ (in our case σ ~ 4.6 × 10⁻¹⁵ cm² while they use a value slightly lower ~3.9 × 10⁻¹⁵ cm²) it suggests a global agreement between our two models.

### 3.2.2. Production by reaction of H₂ and HD with CO₂⁺

The reactions R1 and R2 are another major source of hot H and D in the ionosphere of Mars. The simulated production rates and escape rates are given in Table 4. The reaction between H₂ and CO₂⁺ leads to an escape flux larger than the escape flux computed by Gregory et al. (2023b) by a factor 2, which may be due to the difference in the H₂ and CO₂⁺ density profiles. However, similarly to these authors, we find a lower escape rate for the reaction R1 compared to R3. The efficiency (ratio between escape and production) is slightly lower for the former reaction for both H and D than the dissociative recombination of XCO⁺ due to the lower energy transferred to the atoms. The effusion velocity of H produced by R1 and R3 is lower by one order of magnitude than the effusion velocity of H produced by thermal escape. Our value is two times larger than the total photochemical escape by photoproduction estimated by Krasnopolsky (2010) at solar minimum which could result from his approximation that only the particles produced above the exobase escape, while part of the particles produced below are not fully thermalized and also contribute to the escape. For D, our effusion velocity by R2 and R4 is four times the value of Krasnopolsky (2010) but it
<figure><img src="image_4.png" /><figcaption>Fig. 2. PCM generated maps of the simulated hydrogen density (left panel) and deuterium density (right panel) in Mars body-fixed rotating frame at 200 km.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Left: PCM simulated density of several neutral species at subsolar point from 120 km to the top of the atmosphere (~ 225 km). Right: PCM simulated density of several ions at subsolar point from 120 km to the top of the atmosphere (~ 225 km).</figcaption></figure>

<table><caption>**Table 3**
PCM Simulated thermal escape rates and effusion velocity for H, D, H₂ and HD.</caption><thead><tr><th></th><th>H</th><th>D</th><th>H<sub>2</sub></th><th>HD</th></tr></thead><tbody><tr><td>Thermal escape</td><td></td><td></td><td></td><td></td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>2.1 &times; 10<sup>26</sup></td><td>1.3 &times; 10<sup>20</sup></td><td>1.9 &times; 10<sup>23</sup></td><td>1.8 &times; 10<sup>17</sup></td></tr><tr><td>Effusion velocities (cm. s<sup>-1</sup>)</td><td>3.5 &times; 10<sup>2</sup></td><td>4.1 &times; 10<sup>-1</sup></td><td>4.0 &times; 10<sup>-1</sup></td><td>7.1 &times; 10<sup>-4</sup></td></tr></tbody></table>

Jeans effusion velocity by more than one order of magnitude for H, but larger by more than one order of magnitude for D. The simulated H escape rate is in good agreement with the H escape rate simulated by Gregory et al. (2023a) for their low solar activity case. Considering we simulate a similar production and used comparable elastic cross sections between H and CO₂ (in our case σ ~ 4.6 × 10⁻¹⁵ cm² while they use a value slightly lower ~3.9 × 10⁻¹⁵ cm²) it suggests a global agreement between our two models.

### 3.2.2. Production by reaction of H₂ and HD with CO₂⁺

The reactions R1 and R2 are another major source of hot H and D in the ionosphere of Mars. The simulated production rates and escape rates are given in Table 4. The reaction between H₂ and CO₂⁺ leads to an escape flux larger than the escape flux computed by Gregory et al. (2023b) by a factor 2, which may be due to the difference in the H₂ and CO₂⁺ density profiles. However, similarly to these authors, we find a lower escape rate for the reaction R1 compared to R3. The efficiency (ratio between escape and production) is slightly lower for the former reaction for both H and D than the dissociative recombination of XCO⁺ due to the lower energy transferred to the atoms. The effusion velocity of H produced by R1 and R3 is lower by one order of magnitude than the effusion velocity of H produced by thermal escape. Our value is two times larger than the total photochemical escape by photoproduction estimated by Krasnopolsky (2010) at solar minimum which could result from his approximation that only the particles produced above the exobase escape, while part of the particles produced below are not fully thermalized and also contribute to the escape. For D, our effusion velocity by R2 and R4 is four times the value of Krasnopolsky (2010) but it
<caption>**Table 4**</caption>

MC simulated photochemical production, escape rates and effusion velocities for H, D, $H_2$ and HD.

<table><thead><tr><th></th><th>H</th><th>D</th><th>H<sub>2</sub></th><th>HD</th></tr></thead><tbody><tr><td colspan="5">Dissociative Recombination of XCO<sup>+</sup></td></tr><tr><td>Production above 120 km<br>(s<sup>-1</sup>)</td><td>3.0 ×<br>10<sup>25</sup></td><td>1.2 ×<br>10<sup>22</sup></td><td>X</td><td>X</td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>1.5 ×<br>10<sup>25</sup></td><td>4.0 ×<br>10<sup>21</sup></td><td>X</td><td>X</td></tr><tr><td>Effusion velocities (cm.s<sup>-1</sup>)</td><td>23.3</td><td>12.4</td><td>X</td><td>X</td></tr><tr><td colspan="5">XH + CO<sub>2</sub><sup>+</sup></td></tr><tr><td>Production above 120 km</td><td>2.3 ×<br>10<sup>25</sup></td><td>6.2 ×<br>10<sup>21</sup></td><td>X</td><td>X</td></tr><tr><td>Escape rates (s<sup>-1</sup>)</td><td>9.8 ×<br>10<sup>24</sup></td><td>1.6 ×<br>10<sup>21</sup></td><td>X</td><td>X</td></tr><tr><td>Effusion velocity (cm.s<sup>-1</sup>)</td><td>16.3</td><td>5.0</td><td>X</td><td>X</td></tr><tr><td colspan="5">Collision with hot O</td></tr><tr><td>Production above 120 km<br>(s<sup>-1</sup>)</td><td>4.9 ×<br>10<sup>25</sup></td><td>1.6 ×<br>10<sup>22</sup></td><td>3.2 ×<br>10<sup>25</sup></td><td>1.4 ×<br>10<sup>22</sup></td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>1.9 ×<br>10<sup>25</sup></td><td>1.8 ×<br>10<sup>21</sup></td><td>3.4 ×<br>10<sup>24</sup></td><td>1.0 ×<br>10<sup>21</sup></td></tr><tr><td>Effusion velocity (cm.s<sup>-1</sup>)</td><td>31.7</td><td>5.5</td><td>7.3</td><td>3.9</td></tr></tbody></table>

is not clear if the dissociative recombination of $DCO^+$ was included in the estimate of Krasnopolsky (2010). If this reaction was not included then the effusion velocity computed from R2 is only 1.4 larger than the effusion velocity calculated by Krasnopolsky (2010).

### 3.2.3. Production by collisions with hot oxygen

Hot oxygen produced by dissociative recombination of the $O_2^+$, the main ion in the Martian ionosphere, can transfer part of its energy to the background species, producing hot neutrals (e.g. Leblanc et al., 2019) and contributing to the escape of light species (Gacesa et al., 2012; Shematovich, 2013). The simulated oxygen escape rate is $9 \times 10^{25} \text{ s}^{-1}$ slightly larger than the escape rates reported by Lillis et al. (2017) (between 1.2 and $5.5 \times 10^{25} \text{ s}^{-1}$), who used a model based on the thermospheric/ionospheric MAVEN observations, but in very good agreement with the average value derived from pickup ions measurements by Rahmati et al. (2018) (mean value $\sim 9 \times 10^{25} \text{ s}^{-1}$). However, as discussed in several studies, this value is sensitive to the collisional cross section (Leblanc et al., 2017), and should vary almost inversely with the $O-CO_2$ cross section (Cravens et al., 2017). The $O-CO_2$ cross section used in our simulation $\sim 4.6 \times 10^{-15} \text{ cm}^2$ (Lewkow and Kharchenko, 2014) is $\sim$ four times lower than the cross section used in Lillis et al. (2017), which can explain the difference. Our simulated production and escape rate of hot H and D by collisions with hot oxygen (or reprints produced by hot oxygen) are also indicated in Table 4. These values should also be very sensitive to the collisional cross section. Our derived effusion velocity for H is larger than the value estimated by Krasnopolsky (2010) by one order of magnitude but comparable for $H_2$, D and HD, although in our simulations, the simulated effusion velocity (i.e. efficiency of the process) decreases with the mass of the species. The simulated effusion velocity of $H_2$ is larger by 40% than that of D because the “universal” cross sections of Lewkow and Kharchenko (2014) are 40% larger for molecule-atom collisions that atom-atom collisions (parameter $\gamma$ in their Eq. (2)). Therefore, the energy transfer from hot O to $H_2$ is more efficient than the energy transfer from hot O to D. The ratio between the escape rate and the production rate is $\sim 10\%$ for both species, confirming that the difference is due to a more efficient production of hot $H_2$ than hot D by collisions. The difference in the effusion velocity of $H_2$ and D due to collisions with hot oxygen should be directly scalable with their relative collisional cross sections with O. Our effusion velocities for $H_2$ and HD are one order of magnitude larger than the effusion velocities computed by Gacesa et al. (2012) for an exospheric temperature of 240 K (derived from their simulated escape flux given in

their Table 1 to be 0.59 cm/s for $H_2$ and 0.27 cm/s for HD at 200 km). The exospheric temperature is not expected to have a large effect on this process, and the elastic cross section between O and $H_2$ for the energy range considered is not very different. The main difference is due to a lower collisional cross section between $H_2$ and $CO_2$, since in our simulations we used the value given by Lewkow and Kharchenko (2014) leading to $\sigma \sim 4.6 \times 10^{-15} \text{ cm}^2$, while Gacesa et al. (2012) used the $O-N_2$ cross section from Balakrishnan et al. (1998) $\sim 2 \times 10^{-14} \text{ cm}^2$ so their thermalization of hot $H_2$ with $CO_2$ is four times larger. Assuming a simple inverse relation between this cross section and the escape rate, similar to the approximation considered by Cravens et al. (2017) for oxygen escape, a difference by a factor of $\sim 2$ still remains that could be due to the fact that we neglect the inelastic collisions (i.e the vibrational excitation of these two molecules induced by collisions) which should reduce the kinetic energy transferred to $H_2$ and HD. The ratio between the effusion velocity of $H_2$ and HD ($\sim 1.9$) is in good agreement with the ratio of $\sim 2.1$ computed by Gacesa et al. (2012). The ratio between the thermal escape of H and the escape of H by collisions with hot O is equal to $\sim 11$. This value is slightly lower than the ratio obtained by Shematovich (2013) which is equal to $\sim 18$. Here again, this difference by a factor $\sim 1.6$ is most likely due to different collisional cross sections.

### 3.3. Solar wind loss

To simulate the ion hydrogen and deuterium loss, we use the LatHyS model and the associated particle test. The electromagnetic environment of Mars is simulated with LatHyS using the neutral atmosphere – exosphere obtained with the thermal population only. The effect of the non-thermal component on the electromagnetic environment is neglected in this study. The neutral D and H can be ionized by different processes: photoionization, charge exchange with solar wind protons and electron impact. The last process is not considered in this study since it was found negligible for atomic oxygen (Chaufray et al., 2007).

#### 3.3.1. Photoionization

The photoionization rate of H is computed from 200 km to 10 Martian radii (henceforth $R_m$), while it is computed from 200 km to 3 $R_m$ for D. Due to its low scale height, the ionization of D above 3 Martian radii can be neglected but the ionization of H between 3 $R_m$ and 10 $R_m$ represent $\sim 20\%$ of the total production. We do not consider the attenuation of the solar flux since the medium is optically thin above 200 km. Then, the photoionization production is simply given by the ionization frequency (solar average scenario) multiplied by the local density except in the shadow of Mars where the production is null. The total production of $H^+$ by photoionization is $3.9 \times 10^{24} \text{ s}^{-1}$ and the total production of $D^+$ is $4.3 \times 10^{20} \text{ s}^{-1}$. We use the test particle code to describe the produced ions and compute their trajectory. The ions can be partly lost by charge exchange with atomic H or O, or reimpact the planet (at 200 km) or escape. The escape flux is then estimated from the escaping fraction (Table 5). Since the spatial domain of LatHyS is limited (see Section 2.3), some ions are produced outside the LatHyS domain, where the electromagnetic fields are not computed. We assume that all ions produced outside the domain of LatHyS are escaping. This is a good approximation since these pick-up ions, produced at very high altitudes, are convected by the solar wind after a few gyromotions.

<caption>**Table 5**</caption>

LatHyS simulated production rates, escape rates and effusion velocity for the ion $H^+$ and $D^+$ escape. hv: photoionization, CE: charge exchange.

<table><thead><tr><th>Solar wind interaction</th><th>H (hv)</th><th>D (hv)</th><th>H (CE)</th><th>D (CE)</th></tr></thead><tbody><tr><td>Production rates (s<sup>-1</sup>)</td><td>3.9 × 10<sup>24</sup></td><td>4.3 × 10<sup>20</sup></td><td>2.3 × 10<sup>25</sup></td><td>1.2 × 10<sup>21</sup></td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>1.9 × 10<sup>24</sup></td><td>1.0 × 10<sup>20</sup></td><td>1.9 × 10<sup>25</sup></td><td>0.8 × 10<sup>21</sup></td></tr><tr><td>Effusion velocities (cm. s<sup>-1</sup>)</td><td>3.2</td><td>3.1 × 10<sup>-1</sup></td><td>31.7</td><td>2.5</td></tr></tbody></table>
<caption>**Table 4**</caption>

MC simulated photochemical production, escape rates and effusion velocities for H, D, H₂ and HD.

<table><thead><tr><th></th><th>H</th><th>D</th><th>H₂</th><th>HD</th></tr></thead><tbody><tr><td colspan="5">Dissociative Recombination of XCO<sup>+</sup></td></tr><tr><td>Production above 120 km<br>(s<sup>-1</sup>)</td><td>3.0 ×<br>10<sup>25</sup></td><td>1.2 ×<br>10<sup>22</sup></td><td>X</td><td>X</td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>1.5 ×<br>10<sup>25</sup></td><td>4.0 ×<br>10<sup>21</sup></td><td>X</td><td>X</td></tr><tr><td>Effusion velocities (cm.s<sup>-1</sup>)</td><td>23.3</td><td>12.4</td><td>X</td><td>X</td></tr><tr><td colspan="5">XH + CO<sub>2</sub><sup>+</sup></td></tr><tr><td>Production above 120 km</td><td>2.3 ×<br>10<sup>25</sup></td><td>6.2 ×<br>10<sup>21</sup></td><td>X</td><td>X</td></tr><tr><td>Escape rates (s<sup>-1</sup>)</td><td>9.8 ×<br>10<sup>24</sup></td><td>1.6 ×<br>10<sup>21</sup></td><td>X</td><td>X</td></tr><tr><td>Effusion velocity (cm.s<sup>-1</sup>)</td><td>16.3</td><td>5.0</td><td>X</td><td>X</td></tr><tr><td colspan="5">Collision with hot O</td></tr><tr><td>Production above 120 km<br>(s<sup>-1</sup>)</td><td>4.9 ×<br>10<sup>25</sup></td><td>1.6 ×<br>10<sup>22</sup></td><td>3.2 ×<br>10<sup>25</sup></td><td>1.4 ×<br>10<sup>22</sup></td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>1.9 ×<br>10<sup>25</sup></td><td>1.8 ×<br>10<sup>21</sup></td><td>3.4 ×<br>10<sup>24</sup></td><td>1.0 ×<br>10<sup>21</sup></td></tr><tr><td>Effusion velocity (cm.s<sup>-1</sup>)</td><td>31.7</td><td>5.5</td><td>7.3</td><td>3.9</td></tr></tbody></table>

is not clear if the dissociative recombination of DCO⁺ was included in the estimate of Krasnopolsky (2010). If this reaction was not included then the effusion velocity computed from R2 is only 1.4 larger than the effusion velocity calculated by Krasnopolsky (2010).

### 3.2.3. Production by collisions with hot oxygen

Hot oxygen produced by dissociative recombination of the O₂⁻, the main ion in the Martian ionosphere, can transfer part of its energy to the background species, producing hot neutrals (e.g. Leblanc et al., 2019) and contributing to the escape of light species (Gacesa et al., 2012; Shematovich, 2013). The simulated oxygen escape rate is 9 × 10²⁵ s⁻¹ slightly larger than the escape rates reported by Lillis et al. (2017) (between 1.2 and 5.5 × 10²⁵ s⁻¹), who used a model based on the thermospheric/ionospheric MAVEN observations, but in very good agreement with the average value derived from pickup ions measurements by Rahmati et al. (2018) (mean value ~9 × 10²⁵ s⁻¹). However, as discussed in several studies, this value is sensitive to the collisional cross section (Leblanc et al., 2017), and should vary almost inversely with the O-CO₂ cross section (Cravens et al., 2017). The O-CO₂ cross section used in our simulation ~4.6 × 10⁻¹⁵ cm² (Lewkow and Kharchenko, 2014) is ~ four times lower than the cross section used in Lillis et al. (2017), which can explain the difference. Our simulated production and escape rate of hot H and D by collisions with hot oxygen (or reprints produced by hot oxygen) are also indicated in Table 4. These values should also be very sensitive to the collisional cross section. Our derived effusion velocity for H is larger than the value estimated by Krasnopolsky (2010) by one order of magnitude but comparable for H₂, D and HD, although in our simulations, the simulated effusion velocity (i.e. efficiency of the process) decreases with the mass of the species. The simulated effusion velocity of H₂ is larger by 40% than that of D because the “universal” cross sections of Lewkow and Kharchenko (2014) are 40% larger for molecule-atom collisions that atom-atom collisions (parameter γ in their Eq. (2)). Therefore, the energy transfer from hot O to H₂ is more efficient than the energy transfer from hot O to D. The ratio between the escape rate and the production rate is ~10% for both species, confirming that the difference is due to a more efficient production of hot H₂ than hot D by collisions. The difference in the effusion velocity of H₂ and D due to collisions with hot oxygen should be directly scalable with their relative collisional cross sections with O. Our effusion velocities for H₂ and HD are one order of magnitude larger than the effusion velocities computed by Gacesa et al. (2012) for an exospheric temperature of 240 K (derived from their simulated escape flux given in

their Table 1 to be 0.59 cm/s for H₂ and 0.27 cm/s for HD at 200 km). The exospheric temperature is not expected to have a large effect on this process, and the elastic cross section between O and H₂ for the energy range considered is not very different. The main difference is due to a lower collisional cross section between H₂ and CO₂, since in our simulations we used the value given by Lewkow and Kharchenko (2014) leading to σ ~ 4.6 × 10⁻¹⁵ cm², while Gacesa et al. (2012) used the O-N₂ cross section from Balakrishnan et al. (1998) ~ 2 × 10⁻¹⁴ cm² so their thermalization of hot H₂ with CO₂ is four times larger. Assuming a simple inverse relation between this cross section and the escape rate, similar to the approximation considered by Cravens et al. (2017) for oxygen escape, a difference by a factor of ~2 still remains that could be due to the fact that we neglect the inelastic collisions (i.e the vibrational excitation of these two molecules induced by collisions) which should reduce the kinetic energy transferred to H₂ and HD. The ratio between the effusion velocity of H₂ and HD (~1.9) is in good agreement with the ratio of ~2.1 computed by Gacesa et al. (2012). The ratio between the thermal escape of H and the escape of H by collisions with hot O is equal to ~11. This value is slightly lower than the ratio obtained by Shematovich (2013) which is equal to ~18. Here again, this difference by a factor ~ 1.6 is most likely due to different collisional cross sections.

### 3.3. Solar wind loss

To simulate the ion hydrogen and deuterium loss, we use the LatHyS model and the associated particle test. The electromagnetic environment of Mars is simulated with LatHyS using the neutral atmosphere – exosphere obtained with the thermal population only. The effect of the non-thermal component on the electromagnetic environment is neglected in this study. The neutral D and H can be ionized by different processes: photoionization, charge exchange with solar wind protons and electron impact. The last process is not considered in this study since it was found negligible for atomic oxygen (Chaufray et al., 2007).

#### 3.3.1. Photoionization

The photoionization rate of H is computed from 200 km to 10 Martian radii (henceforth R<sub>m</sub>), while it is computed from 200 km to 3 R<sub>m</sub> for D. Due to its low scale height, the ionization of D above 3 Martian radii can be neglected but the ionization of H between 3 R<sub>m</sub> and 10 R<sub>m</sub> represent ~20% of the total production. We do not consider the attenuation of the solar flux since the medium is optically thin above 200 km. Then, the photoionization production is simply given by the ionization frequency (solar average scenario) multiplied by the local density except in the shadow of Mars where the production is null. The total production of H⁺ by photoionization is 3.9 × 10²⁴ s⁻¹ and the total production of D⁺ is 4.3 × 10²⁰ s⁻¹. We use the test particle code to describe the produced ions and compute their trajectory. The ions can be partly lost by charge exchange with atomic H or O, or reimpact the planet (at 200 km) or escape. The escape flux is then estimated from the escaping fraction (Table 5). Since the spatial domain of LatHyS is limited (see Section 2.3), some ions are produced outside the LatHyS domain, where the electromagnetic fields are not computed. We assume that all ions produced outside the domain of LatHyS are escaping. This is a good approximate since these pick-up ions, produced at very high altitudes, are convected by the solar wind after a few gyromotions.

<caption>**Table 5**</caption>

LatHyS simulated production rates, escape rates and effusion velocity for the ion H⁺ and D⁺ escape. hv: photoionization, CE: charge exchange.

<table><thead><tr><th>Solar wind interaction</th><th>H (hv)</th><th>D (hv)</th><th>H (CE)</th><th>D (CE)</th></tr></thead><tbody><tr><td>Production rates (s<sup>-1</sup>)</td><td>3.9 × 10<sup>24</sup></td><td>4.3 × 10<sup>20</sup></td><td>2.3 × 10<sup>25</sup></td><td>1.2 × 10<sup>21</sup></td></tr><tr><td>Escape rate (s<sup>-1</sup>)</td><td>1.9 × 10<sup>24</sup></td><td>1.0 × 10<sup>20</sup></td><td>1.9 × 10<sup>25</sup></td><td>0.8 × 10<sup>21</sup></td></tr><tr><td>Effusion velocities (cm. s<sup>-1</sup>)</td><td>3.2</td><td>3.1 × 10<sup>-1</sup></td><td>31.7</td><td>2.5</td></tr></tbody></table>
### 3.3.2. Charge exchange with solar wind protons

To simulate the ions loss by charge exchange with solar wind protons, we use a test particle code and proceed in two steps. First, we inject solar wind protons on a square area of length $\Delta Y = \Delta Z = 9R_m$, perpendicular to the x-axis placed at $x = +2.4R_m$, upstream of the planet. We solve the equation of particle's motions assuming the electric and magnetic environment determined by the LatHyS simulation. At each time step, we compute the charge exchange reactions between the solar wind protons and the neutral H and D of the Martian exosphere, producing planetary $H^+$ and $D^+$ ions around Mars (Fig. 4).

The total production rate for both ions is larger than the production rate by photoionization: $2.3 \times 10^{25} \text{ s}^{-1}$ for $H^+$ and $1.2 \times 10^{21} \text{ s}^{-1}$ for $D^+$. The ratio is slightly larger for H than D because of the larger scale height of H, and the fact that charge exchange is not important below the magnetic pile-up boundary, while the photoionization production is maximum at 200 km. The production of $D^+$ and $H^+$ are largest at the magnetic pile-up boundary and at the bow shock (Fig. 4). Far from the planet, the solar wind flux is uniform and then production rates variations are due to neutral hydrogen/deuterium variations resulting in a lower scale height of production for $D^+$ compared to $H^+$. Escape rates of these ions are computed with a particle test code for the second step of the analysis. We inject $H^+$ (and $D^+$ respectively), with a zero velocity,

according to the charge exchange or photoionization productions derived from the first step. The numerical weight of the planetary ion reflects the local production (by charge exchange or photoionization). The dynamic of these planetary $H^+$ and $D^+$ is thus determined by using the test-particle code with the same electric and magnetic field, allowing to characterize the escaping particles. Charge exchange between $H^+$ and $D^+$ ions and the neutral atoms are also considered, re-neutralizing such species, and forming Energetic Neutral Atoms (ENAs), during their trajectory. Most of the produced ions escape in our simulations (Table 5) and therefore the ENA escape can be neglected.

The escape rate of $H^+$ and $D^+$ produced by charge exchange is larger than the escape rate of ions produced by photoionization, in agreement with the simple estimate from Krasnopolsky (2010). As found by this author, we also find that the photoionization channel is comparable to the thermal escape for $D^+$ at $T_{\text{exo}} \sim 200 \text{ K}$ (Krasnopolsky, 2010). We derived an effusion velocity by charge exchange similar to Krasnopolsky (2010) for H ($23.9 \text{ cm.s}^{-1}$ for Texo $\sim 200 \text{ K}$) but lower by a factor of 4 for D. This difference is certainly due to the differences in the simulated solar wind proton flux around Mars, showing structures (Fig. 4) in our simulations, while Krasnopolsky (2010) consider a uniform solar wind flux above the ionopause. The difference is less important for H, because a larger fraction of H is produced at high altitude above the bow-shock

<figure><img src="image_6.png" /><figcaption>Fig. 4. production rate of planetary H+ (left panel) and D+ (right panel) ions around Mars by charge exchange between H and D with the solar wind protons (X-Y MSO plane).</figcaption></figure>
### 3.3.2. Charge exchange with solar wind protons

To simulate the ions loss by charge exchange with solar wind protons, we use a test particle code and proceed in two steps. First, we inject solar wind protons on a square area of length $\Delta Y = \Delta Z = 9R_m$, perpendicular to the x-axis placed at $x = +2.4R_m$, upstream of the planet. We solve the equation of particle's motions assuming the electric and magnetic environment determined by the LatHyS simulation. At each time step, we compute the charge exchange reactions between the solar wind protons and the neutral H and D of the Martian exosphere, producing planetary $H^+$ and $D^+$ ions around Mars (Fig. 4).

The total production rate for both ions is larger than the production rate by photoionization: $2.3 \times 10^{25} \text{ s}^{-1}$ for $H^+$ and $1.2 \times 10^{21} \text{ s}^{-1}$ for $D^+$. The ratio is slightly larger for H than D because of the larger scale height of H, and the fact that charge exchange is not important below the magnetic pile-up boundary, while the photoionization production is maximum at 200 km. The production of $D^+$ and $H^+$ are largest at the magnetic pile-up boundary and at the bow shock (Fig. 4). Far from the planet, the solar wind flux is uniform and then production rates variations are due to neutral hydrogen/deuterium variations resulting in a lower scale height of production for $D^+$ compared to $H^+$. Escape rates of these ions are computed with a particle test code for the second step of the analysis. We inject $H^+$ (and $D^+$ respectively), with a zero velocity,

according to the charge exchange or photoionization productions derived from the first step. The numerical weight of the planetary ion reflects the local production (by charge exchange or photoionization). The dynamic of these planetary $H^+$ and $D^+$ is thus determined by using the test-particle code with the same electric and magnetic field, allowing to characterize the escaping particles. Charge exchange between $H^+$ and $D^+$ ions and the neutral atoms are also considered, re-neutralizing such species, and forming Energetic Neutral Atoms (ENAs), during their trajectory. Most of the produced ions escape in our simulations (Table 5) and therefore the ENA escape can be neglected.

The escape rate of $H^+$ and $D^+$ produced by charge exchange is larger than the escape rate of ions produced by photoionization, in agreement with the simple estimate from Krasnopolsky (2010). As found by this author, we also find that the photoionization channel is comparable to the thermal escape for $D^+$ at $T_{\text{exo}} \sim 200 \text{ K}$ (Krasnopolsky, 2010). We derived an effusion velocity by charge exchange similar to Krasnopolsky (2010) for H ($23.9 \text{ cm.s}^{-1}$ for Texo $\sim 200 \text{ K}$) but lower by a factor of 4 for D. This difference is certainly due to the differences in the simulated solar wind proton flux around Mars, showing structures (Fig. 4) in our simulations, while Krasnopolsky (2010) consider a uniform solar wind flux above the ionopause. The difference is less important for H, because a larger fraction of H is produced at high altitude above the bow-shock

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** production rate of planetary H+ (left panel) and D+ (right panel) ions around Mars by charge exchange between H and D with the solar wind protons (X-Y MSO plane).</figcaption></figure>
(due to its larger scale height) where the solar wind is uniform. Ionization by charge exchange of H₂ and HD are much lower due to a cross section lower by ~2 orders of magnitude compared to H and D (Nakai et al., 1987) and therefore is not considered in this study.

# 4. Discussion

## 4.1. Fractionation and timescale

The dominant atmospheric loss process for hydrogen is the thermal loss of atomic hydrogen due to its low mass. The total simulated escape rate (including H and H₂) is 2.8 × 10²⁶ s⁻¹, the Jeans escape rate representing ~75% of this escape rate. This value agrees with the observed value at equinox reported by Chaufray et al., 2008, Chaffin et al., 2014) from Mars Express but is larger than the value derived from MAVEN by a factor of ~2–3 (Mayyasi et al., 2023) and the Jeans escape rate is larger than our previous estimate at Ls = 0° by a factor of ~4 (Chaufray et al., 2015, 2021). This difference is due to the new microphysics description of the clouds by Navarro et al. (2014) allowing water vapor to reach altitudes >100 km where it can be photolyzed or decomposed by ion reactions, thereby producing more atomic hydrogen, in agreement with observations. This access of water vapor to such high altitudes is supported by observations (Maltagliati et al., 2013; Belyaev et al., 2021).

For deuterium, the total simulated escape rate is 9.4 × 10²¹ s⁻¹ (including D and HD), the Jeans escape being much lower than the photochemical escape. The DCO⁺ dissociative recombination, and the collisions between hot oxygen and deuterium (atomic and molecular) are the two dominant processes and represent ~70% of the total escape. This result differs from Krasnopolsky (2010) who found that charge exchange between solar wind ions and atomic deuterium was the main loss process. In our simulation, the charge exchange is the main escape process in the ion form but represents only 10% of the total escape for both H and D. The simulated escape rates of H and D can be used to derive the two parameters controlling the D/H fractionation on Mars presented in the introduction. The fractionation factor f deduced from our simulations, assuming a D/H ratio ~ 5.5 Standard Mean Ocean Water (SMOW) in the atmosphere, is f = 0.04. This value is << 1 and therefore the timescale of the D/H fractionation (τH/(1-f)) is close to the hydrogen loss timescale. This value is similar to the value obtained by Cangi et al. (2023) who also find that thermal H escape is dominant while D escape is dominated by photochemical loss. These authors did not include the HD and H₂ escape produced by collisions with hot oxygen, explaining their very low escape value compared to our study or to other studies (Krasnopolsky, 2010; Gacesa et al., 2012).

If we only consider the water content of the current Martian atmosphere, that is ~10 pr.µm, then the loss timescale of hydrogen is τH = 11,000 years. This requires that water ice in the polar caps balances the escape rate to maintain a stable abundance of water vapor in the Martian atmosphere. Therefore, the calculation of τH(t) should account for the full size of the exchangeable reservoir of water which is very uncertain (Krasnopolsky and Feldman, 2001) and time dependent due to the variations of the orbital and inclination of Mars along its history (Touma and Wisdom, 1993; Laskar et al., 2004; Forget et al., 2006; Naar, 2023).

Considering the three reservoirs model of Krasnopolsky (2000) presented in Section 1, the solution of Eq. (4) is

$$r(t) = e^{-\int_{0}^{t} \frac{f(t')}{\tau(t')} dt'} \left( 1 + \int_{0}^{t} \frac{dt'}{\tau(t')} e^{\int_{0}^{t} \frac{f(t')}{\tau(t')} dt'} \right)$$

If we suppose f(t) to be constant, this solution can be simplified to

$$r(t) = \frac{1}{f} \left( 1 - (1 - f) e^{-\int_{0}^{t} \frac{dt'}{\tau(t')}} \right)$$

and then, following Krasnopolsky (2000, 2015), we can constrain the integral:

$$\alpha = \int_{0}^{t} \frac{dt'}{\tau(t')} = \frac{1}{f} \ln \left( \frac{1 - f}{1 - f'} \right)$$

Using r = 5.5 and f = 0.04, we found α ~ 5.2. This can be used to deduce an average escape rate during a period of T = 4.5 Gyr.

$$\langle \Phi \rangle = \alpha \frac{N_H}{T} \approx 5.2 \frac{N_H}{T}$$

This average flux will be equal to our simulated hydrogen escape rate for a water exchangeable reservoir of 0.8 m global equivalent layer (GEL) which is much smaller than the water contained in the ice caps (e. g. Lasue et al., 2013). In the frame of this model, even a time dependent f (t) will lead to a very small reservoir since α is minimum for f < < 1 and equal to r-1 = 4.5, which would lead to a water exchangeable reservoir of 0.9 m GEL. A more likely explanation is that the average escape rate during the last 4.5 Gyr is larger than our hydrogen escape rate simulated for the current conditions.

If we consider an exchangeable reservoir of 30 m GEL and the simulated value of f(t), the computed value for <Φ> is equal to 1.0 × 10²⁸ s⁻¹ which is larger than our simulated value and the hydrogen escape rate derived from Lyman-α observation near northern solstice (e. g. Chaffin et al., 2014; Mayyasi et al., 2023). However, this value is more consistent with the H Jeans escape rate simulated at larger Martian obliquity during dusty scenarios, even though still larger by a factor ~ 3 (Gilli et al., 2022).

We can also consider that the current escape is only representative of the last 3.5 Gyr and that at t₁ = 3.5 Gyr ago, the D/H ratio of the exchangeable reservoir was already r₁ = 3 SMOW as measured by Cúriosity/SAM in strongly bound water or hydroxyl in ancient (hesperian) martian clays (Mahaffy et al., 2015). In that case

$$\alpha(t) - \alpha(t_1) = \int_{t_1}^{t} \frac{dt'}{\tau(t')} = \frac{1}{f} \ln \left( \frac{1 - f r_1}{1 - f r} \right) \approx 3.0$$

and the average escape rate <Φ> needed to fractionate the exchangeable reservoir during the last 3.5 Gyr is equal to 7.4 × 10²⁷ s⁻¹, only reduced by 25%.

## 4.2. Current limits of this study

Several assumptions done in this study will need to be refined in the future to better describe the hydrogen fractionation in the upper atmosphere of Mars.

As mentioned in previous sections, the elastic cross sections between H, D, H₂ and HD and other atmospheric species (CO₂, O, Ar, N₂) are not well known whereas the simulated photochemical escape is strongly sensitive to these values. Such uncertainty has a limited impact on the simulated H loss because the thermal escape is dominant. However, the photochemical escape is the major loss process for D and therefore the D escape, and fractionation factor f are very sensitive to the elastic cross sections. Based on the different estimate of the collision cross sections (Lewkow and Kharchenko, 2014, Gacesa et al., 2012, see also Cravens et al., 2017), we estimate the uncertainty on the cross section to be equivalent to a factor ~ 10 at most for H₂ (HD) and CO₂, probably less for H(D) and CO₂. Consequently, a similar uncertainty on the simulated photochemical escape rates, assuming an inverse relation between the escape flux of a species and the collisional cross section of this species with CO₂ (Cravens et al., 2017) is expected. The differences of the simulated hot H₂ and HD escape with Gacesa et al. (2012) (see Section 3.1.3) seem consistent with this simple relation.

Another limit of the current study is the upper boundary condition (at P ~ 10⁻⁸ Pa, ~ 200 km) used in the Mars-PCM for D, H₂ and HD. We only consider the thermal escape while we show in this paper that non-thermal escape is more important for these species. Including the non-
(due to its larger scale height) where the solar wind is uniform. Ionization by charge exchange of H₂ and HD are much lower due to a cross section lower by ~2 orders of magnitude compared to H and D (Nakai et al., 1987) and therefore is not considered in this study.

# 4. Discussion

## 4.1. Fractionation and timescale

The dominant atmospheric loss process for hydrogen is the thermal loss of atomic hydrogen due to its low mass. The total simulated escape rate (including H and H₂) is 2.8 × 10²⁶ s⁻¹, the Jeans escape rate representing ~75% of this escape rate. This value agrees with the observed value at equinox reported by Chaufray et al., 2008, Chaffin et al., 2014) from Mars Express but is larger than the value derived from MAVEN by a factor of ~2–3 (Mayyasi et al., 2023) and the Jeans escape rate is larger than our previous estimate at Ls = 0° by a factor of ~4 (Chaufray et al. 2015, 2021). This difference is due to the new microphysics description of the clouds by Navarro et al. (2014) allowing water vapor to reach altitudes >100 km where it can be photolyzed or decomposed by ion reactions, thereby producing more atomic hydrogen, in agreement with observations. This access of water vapor to such high altitudes is supported by observations (Maltagliati et al., 2013; Belyaev et al., 2021).

For deuterium, the total simulated escape rate is 9.4 × 10²¹ s⁻¹ (including D and HD), the Jeans escape being much lower than the photochemical escape. The DCO⁺ dissociative recombination, and the collisions between hot oxygen and deuterium (atomic and molecular) are the two dominant processes and represent ~70% of the total escape. This result differs from Krasnopolsky (2010) who found that charge exchange between solar wind ions and atomic deuterium was the main loss process. In our simulation, the charge exchange is the main escape process in the ion form but represents only 10% of the total escape for both H and D. The simulated escape rates of H and D can be used to derive the two parameters controlling the D/H fractionation on Mars presented in the introduction. The fractionation factor f deduced from our simulations, assuming a D/H ratio ~ 5.5 Standard Mean Ocean Water (SMOW) in the atmosphere, is f = 0.04. This value is << 1 and therefore the timescale of the D/H fractionation (τH/(1-f)) is close to the hydrogen loss timescale. This value is similar to the value obtained by Cangi et al. (2023) who also find that thermal H escape is dominant while D escape is dominated by photochemical loss. These authors did not include the HD and H₂ escape produced by collisions with hot oxygen, explaining their very low escape value compared to our study or to other studies (Krasnopolsky, 2010; Gacesa et al., 2012).

If we only consider the water content of the current Martian atmosphere, that is ~10 pr.µm, then the loss timescale of hydrogen is τH = 11,000 years. This requires that water ice in the polar caps balances the escape rate to maintain a stable abundance of water vapor in the Martian atmosphere. Therefore, the calculation of τH(t) should account for the full size of the exchangeable reservoir of water which is very uncertain (Krasnopolsky and Feldman, 2001) and time dependent due to the variations of the orbital and inclination of Mars along its history (Touma and Wisdom, 1993; Laskar et al., 2004; Forget et al., 2006; Naar, 2023).

Considering the three reservoirs model of Krasnopolsky (2000) presented in Section 1, the solution of Eq. (4) is

$$r(t) = e^{-\int_{0}^{t} \frac{f(t')}{\tau(t')} dt'} \left( 1 + \int_{0}^{t} \frac{dt'}{\tau(t')} e^{\int_{0}^{t} \frac{f(t')}{\tau(t')} dt'} \right)$$

If we suppose f(t) to be constant, this solution can be simplified to

$$r(t) = \frac{1}{f} \left( 1 - (1 - f) e^{-\int_{0}^{t} \frac{dt'}{\tau(t')}} \right)$$

and then, following Krasnopolsky (2000, 2015), we can constrain the integral:

$$\alpha = \int_{0}^{t} \frac{dt'}{\tau(t')} = \frac{1}{f} \ln \left( \frac{1 - f}{1 - f'} \right)$$

Using r = 5.5 and f = 0.04, we found α ~ 5.2. This can be used to deduce an average escape rate during a period of T = 4.5 Gyr.

$$\langle \Phi \rangle = \alpha \frac{N_H}{T} \approx 5.2 \frac{N_H}{T}$$

This average flux will be equal to our simulated hydrogen escape rate for a water exchangeable reservoir of 0.8 m global equivalent layer (GEL) which is much smaller than the water contained in the ice caps (e. g. Lasue et al., 2013). In the frame of this model, even a time dependent f (t) will lead to a very small reservoir since α is minimum for f < < 1 and equal to r-1 = 4.5, which would lead to a water exchangeable reservoir of 0.9 m GEL. A more likely explanation is that the average escape rate during the last 4.5 Gyr is larger than our hydrogen escape rate simulated for the current conditions.

If we consider an exchangeable reservoir of 30 m GEL and the simulated value of f(t), the computed value for <Φ> is equal to 1.0 × 10²⁸ s⁻¹ which is larger than our simulated value and the hydrogen escape rate derived from Lyman-α observation near northern solstice (e. g. Chaffin et al., 2014; Mayyasi et al., 2023). However, this value is more consistent with the H Jeans escape rate simulated at larger Martian obliquity during dusty scenarios, even though still larger by a factor ~ 3 (Gilli et al., 2022).

We can also consider that the current escape is only representative of the last 3.5 Gyr and that at t₁ = 3.5 Gyr ago, the D/H ratio of the exchangeable reservoir was already r₁ = 3 SMOW as measured by Cúriosity/SAM in strongly bound water or hydroxyl in ancient (hesperian) martian clays (Mahaffy et al., 2015). In that case

$$\alpha(t) - \alpha(t_1) = \int_{t_1}^{t} \frac{dt'}{\tau(t')} = \frac{1}{f} \ln \left( \frac{1 - f r_1}{1 - f r} \right) \approx 3.0$$

and the average escape rate <Φ> needed to fractionate the exchangeable reservoir during the last 3.5 Gyr is equal to 7.4 × 10²⁷ s⁻¹, only reduced by 25%.

## 4.2. Current limits of this study

Several assumptions done in this study will need to be refined in the future to better describe the hydrogen fractionation in the upper atmosphere of Mars.

As mentioned in previous sections, the elastic cross sections between H, D, H₂ and HD and other atmospheric species (CO₂, O, Ar, N₂) are not well known whereas the simulated photochemical escape is strongly sensitive to these values. Such uncertainty has a limited impact on the simulated H loss because the thermal escape is dominant. However, the photochemical escape is the major loss process for D and therefore the D escape, and fractionation factor f are very sensitive to the elastic cross sections. Based on the different estimate of the collision cross sections (Lewkow and Kharchenko, 2014, Gacesa et al., 2012, see also Cravens et al., 2017), we estimate the uncertainty on the cross section to be equivalent to a factor ~ 10 at most for H₂ (HD) and CO₂, probably less for H(D) and CO₂. Consequently, a similar uncertainty on the simulated photochemical escape rates, assuming an inverse relation between the escape flux of a species and the collisional cross section of this species with CO₂ (Cravens et al., 2017) is expected. The differences of the simulated hot H₂ and HD escape with Gacesa et al. (2012) (see Section 3.1.3) seem consistent with this simple relation.

Another limit of the current study is the upper boundary condition (at P ~ 10⁻⁸ Pa, ~ 200 km) used in the Mars-PCM for D, H₂ and HD. We only consider the thermal escape while we show in this paper that non-thermal escape is more important for these species. Including the non-
thermal escape should reduce the simulated density at the exobase (e.g. Krasnopolsky, 2010) and therefore reduce their escape rate, requiring an iterative procedure between PCM-Mars and the simulation of the non-thermal process to be consistent. Such iteration will be studied in the future, but we do not expect it will change the simulated escape rate by more than a factor of 2 (Krasnopolsky, 2010), much less than the uncertainty due to the elastic collisional cross sections.

We also neglect the ballistic transport in the exosphere in the upper boundary conditions in Mars-PCM that should smooth the horizontal variations of the density of the light species at the exobase. Such process was included for H by Chaufray et al. (2018), but the numerical procedure used for this study is difficult to extend to heavier species. However, the effect of the exospheric ballistic transport is expected to be less important for heavier species because the horizontal motion is shorter and therefore, the horizontal redistribution of the density at the exobase limited.

Due to the short timescale variations of the solar wind flux, the production of planetary ions by charge exchange should be also highly variable. Simple scalings are difficult to apply here, because for a much larger dynamic pressure, the magnitude of the magnetic field in the magnetic pileup region would be larger, reducing the efficiency of the escape process by charge exchange (deflecting the solar wind protons further from Mars) while increasing the planetary ion production. A study based on a larger range of solar wind parameter should be considered to better estimate its average contribution to the atmospheric escape.

We only consider the spring equinox season and mean solar conditions. This is assumed to be representative of the average escape rate, but seasonal variations, particularly the southern summer season ($L_s = 270°$) where large amount of water vapor can reach high altitudes, increasing the source of atomic hydrogen (Fedorova et al., 2021; Belyaev et al., 2021; Chaffin et al., 2021; Shaposhnikov et al., 2019) would be interesting to consider because the fractionation of D and H should be more sensitive to the difference in the $H_2O/HDO$ photolysis rate.

Finally, we summed up the escape rates obtained for each process to derive the total escape, assuming implicitly that the thermal and non-thermal population can be fully separated. For ions, the assumption of an independence between thermal and non-thermal escape is valid because the planetary ions are produced below the bow-shock, at an altitude where non-escaping atoms are dominant.

The photochemical reactions used to simulate the photochemical escape also produce thermal hydrogen and deuterium. Thus, part of the hydrogen and deuterium produced by photochemical reactions contribute to both thermal escape e and photochemical escape in our simulations. However, treating escape processes separately is probably a valid assumption because most of the hot (non-thermalized) escaping particles are produced at high altitudes (above ~150 km) where collisions become rare but the production is low, and therefore do not contribute much to the density at the exobase. The density at the exobase is more sensitive to the largest production at lower altitudes (<~150 km), where all the atoms are quickly thermalized.

For the collisions between hot oxygen and atomic hydrogen, this assumption is even less obvious and should be investigated in the future to be validated. However, since the thermal escape is dominant, an overestimate of the production of hot H by collision with hot oxygen atoms will not change the simulated total escape rate. For heavier species it is valid because the thermal escape is much lower and then negligible.

Considering only the largest uncertainty due to the elastic cross sections on the deuterium photochemical escape, we propagate the uncertainty on $f, \alpha$ and $<\Phi>$ (Table 6). To estimate the largest value of $f$, we consider only thermal escape for H, and rescale our simulated D escape flux considering a minimum $D-CO_2$ cross section and a maximum D—O cross section.

The maximum values for $D-CO_2$ and $HD-CO_2$ correspond to the value

<caption>**Table 6**</caption>

Estimated range of different elastic collisional cross sections and their propagation to f, $\alpha$ and $<\Phi>$.

<table><thead><tr><th></th><th>This simulation</th><th>Maximum value</th><th>Minimum value</th></tr></thead><tbody><tr><td>σD-CO<sub>2</sub> (cm<sup>2</sup>)</td><td>4.6 × 10<sup>-15</sup></td><td>2 × 10<sup>-14</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>σD<sub>O</sub> (cm<sup>2</sup>)</td><td>3.3 × 10<sup>-15</sup></td><td>5 × 10<sup>-15</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>σHD-O (cm<sup>2</sup>)</td><td>4.6 × 10<sup>-15</sup></td><td>5 × 10<sup>-15</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>σHD-CO<sub>2</sub> (cm<sup>2</sup>)</td><td>4.6 × 10<sup>-15</sup></td><td>2 × 10<sup>-14</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>Φ(D) (s<sup>-1</sup>)</td><td>8.4 × 10<sup>21</sup></td><td>20.2 × 10<sup>21</sup></td><td>2.6 × 10<sup>21</sup></td></tr><tr><td>Φ(HD) (s<sup>-1</sup>)</td><td>1.0 × 10<sup>21</sup></td><td>2.5 × 10<sup>21</sup></td><td>1.0 × 10<sup>20</sup></td></tr><tr><td>f</td><td>0.04</td><td>0.12</td><td>0.011</td></tr><tr><td>α</td><td>5.2</td><td>8.3</td><td>4.7</td></tr><tr><td>&lt;Φ&gt; for 30 m GEL (s<sup>-1</sup>)</td><td>1.0 × 10<sup>28</sup></td><td>1.6 × 10<sup>28</sup></td><td>8.9 × 10<sup>27</sup></td></tr></tbody></table>

between O and $N_2$ at 1 eV, used for example by Gacesa et al. (2012) to describe the collisions between $CO_2$ and $H_2$. The maximum value for D—O and HD-O correspond to the value for D—O from Zhang et al. (2009) near 1 eV. The minimum values for $D-CO_2$, $HD-CO_2$ correspond to the value for $H_2$—O from Gacesa et al. (2012) near 5 eV.

To calculate the minimum and maximum escape rate of D, we rescaled our simulated escape flux obtained for R2 and R4 assuming an escape rate inversely proportional to the $D-CO_2$ cross section, while for collisions between hot O and D we assume an escape flux proportional to the O—D cross section and inversely proportional to $CO_2$-D cross section. We also add the thermal and ion escape ($1.0 \times 10^{21} s^{-1}$).

For HD a similar rescale is done for collisions between hot O and HD. The range for f is large, between 0.011 and 0.12, while for the three reservoirs model, the effect on the average escape rate needed to fractionate a 30 m GEL exchangeable reservoir from the primordial value to the current value is limited (60% uncertainty).

## 5. Conclusion

In this paper, we estimate numerically the escape rates of H and D (as atoms and molecules) to derive the current fractionation factor of the hydrogen at $L_s = 0°$ and average solar conditions. The hydrogen escape is dominated by the thermal escape, while the deuterium escape is dominated by photochemical processes in the thermosphere/ionosphere. Several limitations, especially the inaccurate knowledge of the elastic and inelastic cross sections between (H, $H_2$, D, HD) and the atmospheric species (O, $CO_2$, $N_2$, Ar) might impact our estimates. Such uncertainty could be propagated to f, assuming a simple linear relation between the escape of the species X produced by collisions with hot oxygen ~ $\sigma_{O,X}$ and an inverse relation with $\sigma_{CO_2,X}$. The derived value of f agrees with previous studies but is smaller by a factor ~5–10 compared to the values deduced from MAVEN/IUVS observations near $L_s = 270°$. Increasing the elastic cross section of O—D by a factor of 5 or reducing the elastic-cross section between $D-CO_2$ should increase the D escape rate, and hence, the value of f but the value of the $D-CO_2$ elastic cross section in our study is already in the lower range of values considered in the literature (Shematovich, 2013). MAVEN/IUVS observations were performed during the dusty season, near $L_s = 270°$ while our simulations are done at $L_s = 0°$, illustrating that f could be seasonally dependent, driven by fractionation of the D/H in the Martian thermosphere due to different processes not investigated in this paper: HDO/$H_2O$ condensation, photodissociation, reaction rates of $H_2O$ and HDO (Montmessin et al., 2005), particularly important at $L_s = 270°$ where large amount of $H_2O/HDO$ could reach the thermosphere. A seasonally variable f is also supported by Cangi et al. (2024).

Other seasons and also other solar activity as well as specific events like flares or CMEs that could change the relative importance of the ion escape resulting from the solar wind interaction (Jakosky et al., 2015; Romanelli et al., 2018) will be considered in the future.
thermal escape should reduce the simulated density at the exobase (e.g. Krasnopolsky, 2010) and therefore reduce their escape rate, requiring an iterative procedure between PCM-Mars and the simulation of the non-thermal process to be consistent. Such iteration will be studied in the future, but we do not expect it will change the simulated escape rate by more than a factor of 2 (Krasnopolsky, 2010), much less than the uncertainty due to the elastic collisional cross sections.

We also neglect the ballistic transport in the exosphere in the upper boundary conditions in Mars-PCM that should smooth the horizontal variations of the density of the light species at the exobase. Such process was included for H by Chaufray et al. (2018), but the numerical procedure used for this study is difficult to extend to heavier species. However, the effect of the exospheric ballistic transport is expected to be less important for heavier species because the horizontal motion is shorter and therefore, the horizontal redistribution of the density at the exobase limited.

Due to the short timescale variations of the solar wind flux, the production of planetary ions by charge exchange should be also highly variable. Simple scalings are difficult to apply here, because for a much larger dynamic pressure, the magnitude of the magnetic field in the magnetic pileup region would be larger, reducing the efficiency of the escape process by charge exchange (deflecting the solar wind protons further from Mars) while increasing the planetary ion production. A study based on a larger range of solar wind parameter should be considered to better estimate its average contribution to the atmospheric escape.

We only consider the spring equinox season and mean solar conditions. This is assumed to be representative of the average escape rate, but seasonal variations, particularly the southern summer season (Ls = 270°) where large amount of water vapor can reach high altitudes, increasing the source of atomic hydrogen (Fedorova et al., 2021; Belyaev et al., 2021; Chaffin et al., 2021; Shaposhnikov et al., 2019) would be interesting to consider because the fractionation of D and H should be more sensitive to the difference in the H₂O/HDO photolysis rate.

Finally, we summed up the escape rates obtained for each process to derive the total escape, assuming implicitly that the thermal and non-thermal population can be fully separated. For ions, the assumption of an independence between thermal and non-thermal escape is valid because the planetary ions are produced below the bow-shock, at an altitude where non-escaping atoms are dominant.

The photochemical reactions used to simulate the photochemical escape also produce thermal hydrogen and deuterium. Thus, part of the hydrogen and deuterium produced by photochemical reactions contribute to both thermal escape e and photochemical escape in our simulations. However, treating escape processes separately is probably a valid assumption because most of the hot (non-thermalized) escaping particles are produced at high altitudes (above ~150 km) where collisions become rare but the production is low, and therefore do not contribute much to the density at the exobase. The density at the exobase is more sensitive to the largest production at lower altitudes (<~150 km), where all the atoms are quickly thermalized.

For the collisions between hot oxygen and atomic hydrogen, this assumption is even less obvious and should be investigated in the future to be validated. However, since the thermal escape is dominant, an overestimate of the production of hot H by collision with hot oxygen atoms will not change the simulated total escape rate. For heavier species it is valid because the thermal escape is much lower and then negligible.

Considering only the largest uncertainty due to the elastic cross sections on the deuterium photochemical escape, we propagate the uncertainty on f, α and <Φ> (Table 6). To estimate the largest value of f, we consider only thermal escape for H, and rescale our simulated D escape flux considering a minimum D-CO₂ cross section and a maximum D—O cross section.

The maximum values for D-CO₂ and HD-CO₂ correspond to the value

<table><caption>**Table 6**
Estimated range of different elastic collisional cross sections and their propagation to f, α and <Φ>.</caption><thead><tr><th></th><th>This simulation</th><th>Maximum value</th><th>Minimum value</th></tr></thead><tbody><tr><td>σD-CO₂ (cm²)</td><td>4.6 × 10<sup>-15</sup></td><td>2 × 10<sup>-14</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>σD_O (cm²)</td><td>3.3 × 10<sup>-15</sup></td><td>5 × 10<sup>-15</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>σHD-O (cm²)</td><td>4.6 × 10<sup>-15</sup></td><td>5 × 10<sup>-15</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>σHD-CO₂ (cm²)</td><td>4.6 × 10<sup>-15</sup></td><td>2 × 10<sup>-14</sup></td><td>2.0 × 10<sup>-15</sup></td></tr><tr><td>Φ(D) (s<sup>-1</sup>)</td><td>8.4 × 10<sup>21</sup></td><td>20.2 × 10<sup>21</sup></td><td>2.6 × 10<sup>21</sup></td></tr><tr><td>Φ(HD) (s<sup>-1</sup>)</td><td>1.0 × 10<sup>21</sup></td><td>2.5 × 10<sup>21</sup></td><td>1.0 × 10<sup>20</sup></td></tr><tr><td>f</td><td>0.04</td><td>0.12</td><td>0.011</td></tr><tr><td>α</td><td>5.2</td><td>8.3</td><td>4.7</td></tr><tr><td>&lt;Φ&gt; for 30 m GEL (s<sup>-1</sup>)</td><td>1.0 × 10<sup>28</sup></td><td>1.6 × 10<sup>28</sup></td><td>8.9 × 10<sup>27</sup></td></tr></tbody></table>

between O and N₂ at 1 eV, used for example by Gacesa et al. (2012) to describe the collisions between CO₂ and H₂. The maximum value for D—O and HD-O correspond to the value for D—O from Zhang et al. (2009) near 1 eV. The minimum values for D-CO₂, HD-CO₂ correspond to the value for H₂—O from Gacesa et al. (2012) near 5 eV.

To calculate the minimum and maximum escape rate of D, we rescaled our simulated escape flux obtained for R2 and R4 assuming an escape rate inversely proportional to the D-CO₂ cross section, while for collisions between hot O and D we assume an escape flux proportional to the O—D cross section and inversely proportional to CO₂-D cross section. We also add the thermal and ion escape (1.0 × 10²¹ s⁻¹).

For HD a similar rescale is done for collisions between hot O and HD. The range for f is large, between 0.011 and 0.12, while for the three reservoirs model, the effect on the average escape rate needed to fractionate a 30 m GEL exchangeable reservoir from the primordial value to the current value is limited (60% uncertainty).

## 5. Conclusion

In this paper, we estimate numerically the escape rates of H and D (as atoms and molecules) to derive the current fractionation factor of the hydrogen at Ls = 0° and average solar conditions. The hydrogen escape is dominated by the thermal escape, while the deuterium escape is dominated by photochemical processes in the thermosphere/ionosphere. Several limitations, especially the inaccurate knowledge of the elastic and inelastic cross sections between (H, H₂, D, HD) and the atmospheric species (O, CO₂, N₂, Ar) might impact our estimates. Such uncertainty could be propagated to f, assuming a simple linear relation between the escape of the species X produced by collisions with hot oxygen ~σO·X and an inverse relation with σCO₂·X. The derived value of f agrees with previous studies but is smaller by a factor ~5–10 compared to the values deduced from MAVEN/IUVS observations near Ls = 270. Increasing the elastic cross section of O—D by a factor of 5 or reducing the elastic-cross section between D-CO₂ should increase the D escape rate, and hence, the value of f but the value of the D-CO₂ elastic cross section in our study is already in the lower range of values considered in the literature (Shematovich, 2013). MAVEN/IUVS observations were performed during the dusty season, near Ls = 270° while our simulations are done at Ls = 0°, illustrating that f could be seasonally dependent, driven by fractionation of the D/H in the Martian thermosphere due to different processes not investigated in this paper: HDO/H₂O condensation, photodissociation, reaction rates of H₂O and HDO (Montmessin et al., 2005), particularly important at Ls = 270° where large amount of H₂O/HDO could reach the thermosphere. A seasonally variable f is also supported by Cangi et al. (2024).

Other seasons and also other solar activity as well as specific events like flares or CMEs that could change the relative importance of the ion escape resulting from the solar wind interaction (Jakosky et al., 2015; Romanelli et al., 2018) will be considered in the future.
### CRediT authorship contribution statement

**J.-Y. Chaufray:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **F. Gonzalez-Galindo:** Writing – original draft, Software, Methodology, Investigation, Funding acquisition. **F. Leblanc:** Writing – original draft, Software, Methodology, Investigation. **R. Modolo:** Writing – original draft, Software, Methodology, Investigation. **M. Vals:** Software, Methodology, Investigation. **F. Montmessin:** Writing – original draft, Software, Methodology. **F. Lefèvre:** Writing – original draft, Software, Methodology. **F. Forget:** Software, Methodology. **M. Lopez-Valverde:** Software, Methodology. **G. Gilli:** Writing – original draft, Methodology.

### Declaration of competing interest

The authors declare that they have no known competing financial

interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

No data was used for the research described in the article.

### Acknowledgments

This work is supported by the Programme National de Planétologie (PNP, France) of CNRS-INSU co-funded by CNES and Programme National Soleil Terre (PNST, France) of CNRS-INSU co-funded by CNES and CEA. The IAA team is financed by grant PREP2022-137579NB-100 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, and acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/ 501100011033.

## Appendix A. Extension of the thermal component into the exosphere

The method to extend the density of the thermal population from the thermosphere to the exosphere is based on Vidal-Madjar and Bertaux (1972) and was also partly presented in Chaufray et al. (2015, appendix). In this appendix, we derive the three equations needed to compute the latitude $δ_C$ and the longitude $α_C$ of the initial position $A_C$ of the atom at the exobase, knowing its position $(r, δ, α)$ and its velocity $(V, θ, φ)$ at A on the trajectory.

We use the same notation as Vidal-Madjar and Bertaux (1972). The position A of the atom on its trajectory is noted $(r, δ, α)$ and its velocity $(V, θ, φ)$ where $r, δ, α$ are the radial, latitude and longitude of A. $V$ is the magnitude of its velocity, $θ$ the angle between the velocity and the radial direction, and $φ$ the last angle defining the velocity vector. In the spherical basis, the position and velocity vectors at A are given by

$$ \vec{r} = r \vec{e}_r $$

$$ \vec{V} = V \cos\theta \vec{e}_r + V \sin\theta \cos\varphi \vec{e}_\theta + V \sin\theta \sin\varphi \vec{e}_\varphi $$

We define a new cartesian basis $(x, y, z')$ where $z'$ is perpendicular to the trajectory (i.e. perpendicular to $r$ and $V$), $x'$ is directed from the planet center toward the pericenter of the trajectory and $y'$ complete the system to have a direct basis:

These three vectors, expressed in $(e_r, e_θ, e_φ)$ are given by

$$ \vec{e}' = \frac{\vec{r} \times \vec{V}}{\| \vec{r} \times \vec{V} \|} = - \sin\varphi \vec{e}_\theta + \cos\varphi \vec{e}_\varphi $$

$$ \vec{x}' = \frac{\vec{e}'}{e} = \frac{1}{e} \left[ \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) \vec{e}_r - \frac{2u^2}{y} \sin\theta\cos\theta(\cos\varphi\vec{e}_\theta + \sin\varphi\vec{e}_\varphi) \right] $$

$$ \vec{y}' = \vec{z}' \times \vec{x}' = \frac{2u^2}{y} \sin\theta\cos\theta \vec{e}_r + \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) (\cos\varphi\vec{e}_\theta + \sin\varphi\vec{e}_\varphi) $$

where we used the eccentric vector (its norm is the eccentricity of the trajectory and its direction the direction of the pericenter) defined by

$$ \vec{e} = \frac{\vec{V} \times \vec{L}}{GMm} - \vec{e}_r = \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) \vec{e}_r - \frac{2u^2}{y} \sin\theta\cos\theta(\cos\varphi\vec{e}_\theta + \sin\varphi\vec{e}_\varphi) $$

and defined the adimensionned velocity $u$ and distance $y$ as Vidal-Madjar and Bertaux (1972) by

$$ u = \frac{V}{V_{\text{esc,c}}}; u^2 = \frac{r_c V^2}{2GM} $$

$$ y = \frac{r_c}{r} $$

$L$ is the angular momentum defined by

$$ \vec{L} = m \vec{r} \times \vec{V} = mr \sin\theta \vec{z}' $$

$M$ is the mass of the planet, $G$ the gravitational constant, $m$ the mass of the atom, and $V_{\text{esc,c}}$ the escape velocity at the exobase distance $r_c$.

From the norm of the eccentric vector, we can derive the expression of the eccentricity of the trajectory leading to eq. 12 of Vidal-Madjar and Bertaux (1972):

$$ e = \sqrt{1 + \frac{4u^2}{y^2} \sin^2\theta(u^2 - y)} $$
### CRediT authorship contribution statement

**J.-Y. Chaufray:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **F. Gonzalez-Galindo:** Writing – original draft, Software, Methodology, Investigation, Funding acquisition. **F. Leblanc:** Writing – original draft, Software, Methodology, Investigation. **R. Modolo:** Writing – original draft, Software, Methodology, Investigation. **M. Vals:** Software, Methodology, Investigation. **F. Montmessin:** Writing – original draft, Software, Methodology. **F. Lefèvre:** Writing – original draft, Software, Methodology. **F. Forget:** Software, Methodology. **M. Lopez-Valverde:** Software, Methodology. **G. Gilli:** Writing – original draft, Methodology.

### Declaration of competing interest

The authors declare that they have no known competing financial

interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

No data was used for the research described in the article.

### Acknowledgments

This work is supported by the Programme National de Planétologie (PNP, France) of CNRS-INSU co-funded by CNES and Programme National Soleil Terre (PNST, France) of CNRS-INSU co-funded by CNES and CEA. The IAA team is financed by grant PREP2022–137579NB-100 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, and acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033.

## Appendix A. Extension of the thermal component into the exosphere

The method to extend the density of the thermal population from the thermosphere to the exosphere is based on Vidal-Madjar and Bertaux (1972) and was also partly presented in Chaufray et al. (2015, appendix). In this appendix, we derive the three equations needed to compute the latitude $\delta_C$ and the longitude $\alpha_C$ of the initial position $A_C$ of the atom at the exobase, knowing its position $(r, \delta, \alpha)$ and its velocity $(V, \theta, \varphi)$ at A on the trajectory.

We use the same notation as Vidal-Madjar and Bertaux (1972). The position A of the atom on its trajectory is noted $(r, \delta, \alpha)$ and its velocity $(V, \theta, \varphi)$ where $r, \delta, \alpha$ are the radial, latitude and longitude of A. V is the magnitude of its velocity, $\theta$ the angle between the velocity and the radial direction, and $\varphi$ the last angle defining the velocity vector. In the spherical basis, the position and velocity vectors at A are given by

$$ \vec{r} = r \vec{e}_r $$

$$ \vec{V} = V \cos\theta \vec{e}_r + V \sin\theta \cos\varphi \vec{e}_\theta + V \sin\theta \sin\varphi \vec{e}_\varphi $$

We define a new cartesian basis $(x, y, z')$ where $z'$ is perpendicular to the trajectory (i.e. perpendicular to $r$ and $V$), $x'$ is directed from the planet center toward the pericenter of the trajectory and $y'$ complete the system to have a direct basis:

These three vectors, expressed in $(e_r, e_\theta, e_\varphi)$ are given by

$$ \vec{z'} = \frac{\vec{r} \times \vec{V}}{\| \vec{r} \times \vec{V} \|} = - \sin\varphi \vec{e}_\theta + \cos\varphi \vec{e}_\varphi $$

$$ \vec{x'} = \frac{\vec{e}'}{e} = \frac{1}{e} \left[ \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) \vec{e}_r - \frac{2u^2}{y} \sin\theta\cos\theta(\cos\varphi\vec{e}_\theta + \sin\varphi\vec{e}_\varphi) \right] $$

$$ \vec{y'} = \vec{z'} \times \vec{x'} = \frac{2u^2}{y} \sin\theta\cos\theta \vec{e}_r + \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) (\cos\varphi\vec{e}_\theta + \sin\varphi\vec{e}_\varphi) $$

where we used the eccentric vector (its norm is the eccentricity of the trajectory and its direction the direction of the pericenter) defined by

$$ \vec{e} = \frac{\vec{V} \times \vec{L}}{GMm} - \vec{e}_r = \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) \vec{e}_r - \frac{2u^2}{y} \sin\theta\cos\theta(\cos\varphi\vec{e}_\theta + \sin\varphi\vec{e}_\varphi) $$

and defined the adimensionned velocity $u$ and distance $y$ as Vidal-Madjar and Bertaux (1972) by

$$ u = \frac{V}{V_{\text{esc,c}}}; u^2 = \frac{r_c V^2}{2GM} $$

$$ y = \frac{r_c}{r} $$

$L$ is the angular momentum defined by

$$ \vec{L} = m \vec{r} \times \vec{V} = mr \sin\theta \vec{z'} $$

$M$ is the mass of the planet, $G$ the gravitational constant, $m$ the mass of the atom, and $V_{\text{esc,c}}$ the escape velocity at the exobase distance $r_c$.

From the norm of the eccentric vector, we can derive the expression of the eccentricity of the trajectory leading to eq. 12 of Vidal-Madjar and Bertaux (1972):

$$ e = \sqrt{1 + \frac{4u^2}{y^2} \sin^2\theta(u^2 - y)} $$
The true anomaly angle $\beta$ of A is the angle between the pericenter direction $x'$ and the position of A ($e_r$) and then the cosine and sinus of $\beta$ are given by

$$ \cos\beta = \vec{x}' \cdot \vec{e}_r = \frac{1}{e} \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) $$

$$ \sin\beta = \vec{y}' \cdot \vec{e}_r = \frac{2u^2}{ey} \sin\theta\cos\theta $$

The semi-lactus rectum $\omega$ of the trajectory can be deduced from the polar equation (Vidal-Madjar and Bertaux, 1972)

$$ \omega = 2r_c \frac{u^2 \sin^2\theta}{y^2} $$

Note that $\sin\beta < 0$ when $\theta > \pi/2$ (when the particle at A moves downward). For the hyperbolic trajectories (escaping atoms), only $\theta < \pi/2$ is permitted since, by assumption, no particles come back to the exobase with a velocity larger than the escape velocity.

The expressions of $x$, $y'$ and $z'$ can be simplified using $\beta$

$$ \vec{z}' = -\sin\varphi\vec{e}_{\theta} + \cos\varphi\vec{e}_{\varphi} $$

$$ \vec{x}' = [\cos\beta\vec{e}_r - \sin\beta(\cos\varphi\vec{e}_{\theta} + \sin\varphi\vec{e}_{\varphi})] $$

$$ \vec{y}' = \sin\beta\vec{e}_r + \cos\beta(\cos\varphi\vec{e}_{\theta} + \sin\varphi\vec{e}_{\varphi}) $$

The particle at A is coming from the origin point $A_c$ at the exobase ($r_c$), defined by its true anomaly angle $\beta_c$. In the $(x', y', z')$ frame, the radial vector $r_c$ is given by

$$ \vec{r}_c = r_c(\cos\beta_c\vec{x}' + \sin\beta_c\vec{y}') $$

and then, in the spherical frame ($e_r$, $e_\theta$, $e_\varphi$) we have

$$ \vec{r}_c = r_c(\cos\beta_c\cos\beta + \sin\beta\sin\beta_c)\vec{e}_r - (\sin\beta\cos\beta_c - \sin\beta_c\cos\beta)(\cos\varphi\vec{e}_{\theta} + \sin\varphi\vec{e}_{\varphi}) $$

Noting $\psi = \beta - \beta_c$, we can simplify this expression by

$$ \vec{r}_c = r_c(\cos\psi\vec{e}_r - \sin\psi\cos\varphi\vec{e}_{\theta} - \sin\psi\sin\varphi\vec{e}_{\varphi}) $$

In the cartesian basis xyz, the coordinates of $r$ and $r_c$ are given by

$$ \vec{r} = r \begin{pmatrix} \cos\delta\cos\alpha \\ \cos\delta\sin\alpha \\ \sin\delta \end{pmatrix} $$

$$ \vec{r}_c = r_c \begin{pmatrix} \cos\delta_c\cos\alpha_c \\ \cos\delta_c\sin\alpha_c \\ \sin\delta_c \end{pmatrix} $$

In the xyz basis, the spherical vectors can be expressed using the well-known expressions:

$$ \vec{e}_r = \begin{pmatrix} \cos\delta\cos\alpha \\ \cos\delta\sin\alpha \\ \sin\delta \end{pmatrix} $$

$$ \vec{e}_{\theta} = \begin{pmatrix} \sin\delta\cos\alpha \\ \sin\delta\sin\alpha \\ -\cos\delta \end{pmatrix} $$

$$ \vec{e}_{\varphi} = \begin{pmatrix} -\sin\alpha \\ \cos\alpha \\ 0 \end{pmatrix} $$

and then, the projections of $r_c$ along $x$, $y$ and $z$ lead to the three equations needed to fully determine the latitude $\delta_c$ and longitude $\alpha_c$ of $A_c$:

$$ \begin{pmatrix} \cos\delta_c\cos\alpha_c \\ \cos\delta_c\sin\alpha_c \\ \sin\delta_c \end{pmatrix} = \begin{pmatrix} \cos\psi\cos\delta\cos\alpha - \sin\psi\cos\varphi\sin\delta\cos\alpha + \sin\psi\sin\varphi\sin\alpha \\ \cos\psi\cos\delta\sin\alpha - \sin\psi\cos\varphi\sin\delta\sin\alpha - \sin\psi\sin\varphi\cos\alpha \\ \cos\psi\sin\delta + \sin\psi\cos\varphi\cos\delta \end{pmatrix} $$

The first equation corresponds to eq. 14 of Vidal-Madjar and Bertaux (1972). The third equation to eq. 13 of Vidal-Madjar and Bertaux (1972). The second equation, not given in Vidal-Madjar and Bertaux (1972), is needed to fully constrain the longitude $\alpha_c$ in the general case when there is no symmetry of the density and the temperature at the exobase.
The true anomaly angle $\beta$ of A is the angle between the pericenter direction $x'$ and the position of A ($e_r$) and then the cosine and sinus of $\beta$ are given by

$$ \cos\beta = \vec{x}' \cdot \vec{e}_r = \frac{1}{e} \left( \frac{2u^2}{y} \sin^2\theta - 1 \right) $$

$$ \sin\beta = \vec{y}' \cdot \vec{e}_r = \frac{2u^2}{ey} \sin\theta\cos\theta $$

The semi-lactus rectum $\omega$ of the trajectory can be deduced from the polar equation (Vidal-Madjar and Bertaux, 1972)

$$ \omega = 2r_c \frac{u^2 \sin^2\theta}{y^2} $$

Note that $\sin\beta < 0$ when $\theta > \pi/2$ (when the particle at A moves downward). For the hyperbolic trajectories (escaping atoms), only $\theta < \pi/2$ is permitted since, by assumption, no particles come back to the exobase with a velocity larger than the escape velocity.

The expressions of $x$, $y'$ and $z'$ can be simplified using $\beta$

$$ \vec{z}' = -\sin\varphi\vec{e}_{\theta} + \cos\varphi\vec{e}_{\varphi} $$

$$ \vec{x}' = [\cos\beta\vec{e}_r - \sin\beta(\cos\varphi\vec{e}_{\theta} + \sin\varphi\vec{e}_{\varphi})] $$

$$ \vec{y}' = \sin\beta\vec{e}_r + \cos\beta(\cos\varphi\vec{e}_{\theta} + \sin\varphi\vec{e}_{\varphi}) $$

The particle at A is coming from the origin point $A_c$ at the exobase ($r_c$), defined by its true anomaly angle $\beta_c$. In the $(x',y',z')$ frame, the radial vector $r_c$ is given by

$$ \vec{r}_c = r_c(\cos\beta_c\vec{x}' + \sin\beta_c\vec{y}') $$

and then, in the spherical frame ($e_r, e_{\theta}, e_{\varphi}$) we have

$$ \vec{r}_c = r_c(\cos\beta_c\cos\beta + \sin\beta\sin\beta_c)\vec{e}_r - (\sin\beta\cos\beta_c - \sin\beta_c\cos\beta)(\cos\varphi\vec{e}_{\theta} + \sin\varphi\vec{e}_{\varphi}) $$

Noting $\psi = \beta - \beta_c$, we can simplify this expression by

$$ \vec{r}_c = r_c(\cos\psi\vec{e}_r - \sin\psi\cos\varphi\vec{e}_{\theta} - \sin\psi\sin\varphi\vec{e}_{\varphi}) $$

In the cartesian basis xyz, the coordinates of $r$ and $r_c$ are given by

$$ \vec{r} = r \begin{pmatrix} \cos\delta\cos\alpha \\ \cos\delta\sin\alpha \\ \sin\delta \end{pmatrix} $$

$$ \vec{r}_c = r_c \begin{pmatrix} \cos\delta_c\cos\alpha_c \\ \cos\delta_c\sin\alpha_c \\ \sin\delta_c \end{pmatrix} $$

In the xyz basis, the spherical vectors can be expressed using the well-known expressions:

$$ \vec{e}_r = \begin{pmatrix} \cos\delta\cos\alpha \\ \cos\delta\sin\alpha \\ \sin\delta \end{pmatrix} $$

$$ \vec{e}_{\theta} = \begin{pmatrix} \sin\delta\cos\alpha \\ \sin\delta\sin\alpha \\ -\cos\delta \end{pmatrix} $$

$$ \vec{e}_{\varphi} = \begin{pmatrix} -\sin\alpha \\ \cos\alpha \\ 0 \end{pmatrix} $$

and then, the projections of $r_c$ along $x$, $y$ and $z$ lead to the three equations needed to fully determine the latitude $\delta_c$ and longitude $\alpha_c$ of $A_c$:

$$ \begin{pmatrix} \cos\delta_c\cos\alpha_c \\ \cos\delta_c\sin\alpha_c \\ \sin\delta_c \end{pmatrix} = \begin{pmatrix} \cos\psi\cos\delta\cos\alpha - \sin\psi\cos\varphi\sin\delta\cos\alpha + \sin\psi\sin\varphi\sin\alpha \\ \cos\psi\cos\delta\sin\alpha - \sin\psi\cos\varphi\sin\delta\sin\alpha - \sin\psi\sin\varphi\cos\alpha \\ \cos\psi\sin\delta + \sin\psi\cos\varphi\cos\delta \end{pmatrix} $$

The first equation corresponds to eq. 14 of Vidal-Madjar and Bertaux (1972). The third equation to eq. 13 of Vidal-Madjar and Bertaux (1972). The second equation, not given in Vidal-Madjar and Bertaux (1972), is needed to fully constrain the longitude $\alpha_c$ in the general case when there is no symmetry of the density and the temperature at the exobase.
## References

Anderson, D.E., Hord, C.W., 1971. Mariner 6 and 7 ultraviolet spectrometer experiment: analysis of hydrogen Lyman-Alpha data. J. Geophys. Res. 76, 6666-6673.

Aoki, S., Nakagawa, H., Sagawa, H., Giuranna, M., Sindoni, G., Aronica, A., Kasaba, Y., 2015. Seasonal variation of the HDO/H₂O ratio in the atmosphere of Mars at the middle of northern spring and beginning of northern summer. Icarus 260, 7-22.

Balakrishnan, N., Kharchenko, V., Dalgarno, A., 1998. Slowing of energetic O(3P) atom in collisions with N₂. J. Geophys. Res. 103, 23393-23398.

Belyaev, D.A., Fedorova, A.A., Trokhimovskiy, A., Alday, J., Montmessin, F., Korablev, O.I., et al., 2021. Revealing a high water abundance in the upper mesosphere of Mars with ACS onboard TGO. Geophys. Res. Lett. 48 (10) https://doi.org/10.1029/2021GL093411.

Benna, M., Mahaffy, P.R., Grebowsky, J.M., Fox, J.L., Yelle, R.V., Jakosky, B.M., 2015. First measurements of composition and dynamics of the Martian ionosphere by MAVEN's neutral gas and ion mass spectrometer. Geophys. Res. Lett. 42, 8958-8965. https://doi.org/10.1002/2015GL066146.

Bibring, J.-P.Y., Langevin, J.F. Mustard, Poulet, F., d'Arvidson, R., Gendrin, A., Gondet, B., Mangold, N., Pinet, P., Forget, F., the OMEGA team, 2006. Global mineralogical and aqueous Mars history derived from OMEGA/Mars express data. Science 312, 400-404.

Cangi, E., Chaffin, M., Yelle, R., Gregory, B., Deighan, J., 2023. Fully coupled photochemistry of the deuterated ionosphere of Mars and its effects on escape of H and D. J. Geophys. Res: Planet 128. https://doi.org/10.1029/2022JE007713 e2022JE07713.

Cangi, E., Chaffin, M.S., Yelle, R.V., Gregory, B.S., Deighan, J., 2024. Seasonal Enhancement in Upper Atmospheric D/H at Mars Driven by both Thermospheric Temperature and Mesospheric Water. Geophys. Res. Lett. p. 51.

Carr, M.H., 1996. Water on Mars. Oxford University Press, New York.

Chaffin, M.S., Chaufray, J.-Y., Stewart, I., Montmessin, F., Schneider, N., Bertaux, J.-L., 2014. Unexpected variability of Martian hydrogen escape. Geophys. Res. Lett. 41, 314-320. https://doi.org/10.1002/2013GL058578.

Chaffin, M.S., Chaufray, J.-Y., Deighan, J., Schneider, N.M., Mayyasi, M., Clarke, J.T., Thiemann, E., Jain, S.K., Crismani, M.M.J., Stiepen, A., Eparvier, F.G., McClintock, W.E., Stewart, A.I.F., Holswak, G.M., Montmessin, F., Jakosky, B., 2018. Mars H escape rates derived from MAVEN/IUVS Lyman-α brightness measurements and their dependence on model assumptions. J. Geophys. Res: Planets 123, 2192-2210.

Chaffin, M.S., Kass, D.M., Aoki, S., Fedorova, A.A., Deighan, J., Connour, K., et al., 2021. Regional dust storms on Mars enhance water loss to space. Nat. Astron. 5, 1036-1042.

Chassefière, E., Leblanc, F., 2004. Mars atmospheric escape and evolution: interaction with the solar wind. Planet. Space Sci. 52, 1039-1058 doi:1016/j.pss.204.07.002.

Chaufray, J.-Y., 2021. Departure of the thermal escape rate from the Jeans escape rate for atomic hydrogen at Earth, Mars and Pluto. Planet. Space Sci. 198, 105178 doi: 10.1016/j.pss.2021.105178.

Chaufray, J.-Y., Bertaux, J.-L., Leblanc, F., Quémérais, E., 2008. Observation of the hydrogen corona with SPICAM on Mars express. Icarus 195, 598-613.

Chaufray, J.-Y., Gonzalez-Galindo, F., Forget, F., Lopez-Valverde, M., Leblanc, F., Modolo, R., Hess, S., Yagi, M., Blelie, P.-L., Witasse, O., 2014. Three dimensional Martian ionosphere model: II effect of transport processes due to pressure gradients. J. Geophys. Res: Planets 119, 1614-1636. https://doi.org/10.1002/2013JE004551.

Chaufray, J.-Y., Gonzalez-Galindo, F., Forget, M.A. Lopez-Valverde, Leblanc, F., Modolo, R., Hess, S., 2015. Variability of the hydrogen in the martian upper atmosphere as simulated by a 3D atmosphere-exosphere coupling. Icarus 245, 282-294.

Chaufray, J.-Y., Modolo, R., Leblanc, F., Chanteur, G., Johnson, R.E., Luhmann, J.G., 2007. Mars solar wind interaction: Formation of the Martian corona and atmospheric loss to space. J. Geophys. Res. 112, E09009 https://doi.org/10.1029/2007JE002915.

Chaufray, J.-Y., Yelle, R.V., Gonzalez-Galindo, F., Forget, F., Lopez-Valverde, M.A., Leblanc, F., Modolo, R., 2018. Effect of the lateral exospheric transport on the horizontal hydrogen distribution near the exobase of Mars. J. Geophys. Res. Space Physics 123, 2441-2454. https://doi.org/10.1002/2017JA025163.

Chaufray, J.-Y., Gonzalez-Galindo, F., Lopez-Valverde, M.A., Forget, F., Quémérais, E., Bertaux, J.-L., Montmessin, F., Chaffin, M., Schneider, N., Clarke, J.T., Leblanc, F., Modolo, R., Yelle, R., 2021a. Study of the hydrogen escape rate at Mars during Martian years 28 and 29 from comparisons between SPICAM/Mars express observations and GCM-LMD simulations. Icarus 353, 113498. https://doi.org/10.1016/j.icarus.2019.113498.

Chaufray, J.-Y., Mayyasi, M., Chaffin, M., Deighan, J., Bhattacharyya, D., Clarke, J., Jain, S., Schneider, N., Jakosky, B., 2021b. Estimate of the D/H ratio in the Martian upper atmosphere from the low spectral resolution mode of MAVEN/IUVS. J. Geophys. Res.: Planets 126. https://doi.org/10.1029/JE006814 e2020JE006814.

Clarke, J.T., Mayyasi, M., Bhattacharyya, D., Schneider, N.M., McClintock, W.E., Deighan, J.I., Stewart, A.I.F., Chaufray, J.-Y., Chaffin, M.S., Jain, S.K., Stiepen, A., Crismani, M., Holswak, G.M., Montmessin, F., Jakosky, B., 2017. Variability of D and H in the Martian upper atmosphere observed with the MAVEN IUVS echelle channel. J. Geophys. Res. 122, 2336-2344. https://doi.org/10.1002/2016JA023479.

Cravens, T.E., Rahmati, A., Fox, J.L., Lillis, R., Bougher, S., Luhmann, J., Sakai, S., Deighan, J., Lee, Y., Combi, M., Jakosky, B., 2017. J. Geophys. Res. Space Physics 122, 1102-1116. https://doi.org/10.1002/2016JA023461.

Deighan, J., Chaffin, M.S., Chaufray, J.-Y., Stewart, A.I.F., Schneider, N.M., Jain, S.K., Stiepen, A., Crismani, M., McClintock, W.E., Clarke, J.T., Holswak, G.M., Montmessin, F., Eparvier, F.G., Thiemann, E.M.B., Chamberlin, P.C., Jakosky, B.M.,

2015. MAVEN IUVS observation of the hot oxygen corona. Geophys. Res. Lett. 42, 9009-9014. https://doi.org/10.1002/2015GL065487.

Dickeson, Z.I., Davis, J.M., 2020. Martian oceans. Astron. Geophys. 61, 311-317.

Dong, C., Bougher, S.W., Ma, Y., Lee, Y., Toth, G., Nagy, A.F., Fang, X., Luhmann, J., Liemohn, M.W., Halekas, J.S., Tenishev, V., Pawlowski, D.J., Combi, M.R., 2018. Solar wind interaction with the Martian upper atmosphere: roles of the cold thermosphere and hot oxygen corona, (2018). J. Geophys. Res. Space Physics 123, 6639-6654.

Elrod, M.K., Bougher, S., Bell, J., Mahaffy, P., Benna, M., Stone, S., Yelle, R., Jakosky, B., 2017. He bulge revealed: he and CO2 diurnal and seasonal variations in the upper atmosphere of Mars as detected by MAVEN NGIMS. J. Geophys. Res. Space Physics 122, 2564-2573. https://doi.org/10.1002/2016JA023482.

Ergun, R.E., Andersson, L.A., Fowler, C.M., Thaller, S.A., Yelle, R.V., 2021. In situ measurements of electron temperature and density in the Mars' dayside ionosphere. Geophys. Res. Lett. 48 e2021GL093623.

Fedorova, A., Montmessin, F., Korablev, O., Lefèvre, F., Trokhimovskiy, A., Bertaux, J.-L., 2021. Multi-annual monitoring of the water vapor vertical distribution on Mars by SPICAM on Mars express. J. Geophys. Res: Planets 126, e06616. https://doi.org/10.1029/2020JE006616.

Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M., Lewis, S.R., Read, P.L., Huot, J.-P., 1999. Improved general circulation models of the Martian atmosphere from the surface to above 80 km. J. Geophys. Res. 104, 24155-24175.

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

Fox, J., 2015. The chemistry of protonated species in the martian ionosphere. Icarus 252, 366-392.

Gacesa, M., Zhang, P., Kharchenko, V., 2012. Non-thermal escape of molecular hydrogen from Mars. Geophys. Res. Lett. 39, L10203. https://doi.org/10.1029/2012GL050904.

Gilli, G., Gonzalez-Galindo, F., Forget, F., Millour, E., Naar, J., Chaufray, J.-Y., 2022. On the effect of the orbital parameters of Mars to the hydrogen escape and the fate of water in the last millions of years. In: 16th EPSC, Abstract #578. https://doi.org/10.5194/epsc2022-578.

Gonzalez-Galindo, F., Forget, F., Lopez-Valverde, M.A., Coll, M. Angelats i, Millour, E., 2009. A ground-to-exosphere Martian general circulation model: 1 seasonal, diurnal, and solar cycle variation of thermospheric temperatures. J. Geophys. Res. 114, E04001. https://doi.org/10.1029/2008JE003246, 2009.

Gonzalez-Galindo, F., Chaufray, J.-Y., Lopez-Valverde, M.A., Gilli, G., Forget, F., Leblanc, F., Modolo, R., Hess, S., Yagi, M., 2013. Three-dimensional Martian ionosphere model: I. The photochemical ionosphere below 180 km. J. Geophys. Res.: Planets 118, 2105-2123. https://doi.org/10.1002/jgre.20150.

Gonzalez-Galindo, F., Lopez-Valverde, M.A., Forget, F., Garcia-Comas, M., Millour, E., Montabone, L., 2015. Variability of the Martian thermosphere during eight Martian years as simulated by a ground-to-exosphere global circulation model. J. Geophys. Res.: Planets 120, 2020-2035. https://doi.org/10.1002/2015JE004925.

Gregory, B., Elliott, R.D., Deighan, J., Gröller, H., Chaffin, M.S., 2023a. HCO+ dissociative recombination: a significant driver of nonthermal hydrogen loss at Mars. J. Geophys. Res.: Planets 128. https://doi.org/10.1029/2022JE007576 e2022JE007576.

Gregory, B., Chaffin, M.S., Elliott, R.D., Deighan, J., Gröller, H., Cangi, E., 2023b. Nonthermal hydrogen loss at Mars: contribution of photochemical mechanisms to escape and identification of key processes. J. Geophys. Res.: Planets 128. https://doi.org/10.1029/JE007802 e2023JE007802.

Halekas, J.S., 2017. Seasonal variability of the hydrogen exosphere of Mars. J. Geophys. Res.: Planets 122. https://doi.org/10.1002/2017JE005306.

Hanley, K.G., Fadden, J.P., Mitchell, D.L., Fowler, C.M., Stone, S.W., et al., 2021. In situ measurements of thermal ion temperature in the Martian ionosphere. J. Geophys. Res. Space Physics 126.

Hodges, R.R., Johnson, F.S., 1968. Lateral transport in planetary exospheres. J. Geophys. Res. Space Physics 73, 7307-7316.

Ishiguro, E., et al., 1952. On the polarizability of the hydrogen molecule. Proc. Phys. Soc. A 65, 178.

Jakosky, B., Phillips, R.J., 2001. Mars' volatile and climate history. Nature 412, 237-244.

Jakosky, B.M., Grebowsky, J.M., Luhmann, J.G., Connerney, J., Eparvier, F., Ergun, R., et al., 2015. MAVEN observations of the response of Mars to an interplanetary coronal mass ejection, Science, American, Association for the Adv. Sci. 350 (6261) aad0210-1-aad020-7.

Jakosky, B., et al., 2018. Loss of the Martian atmosphere to space: present-day loss rates determined from MAVEN observations and integrated loss through time. Icarus 315, 146-157.

Kallio, E., Chaufray, J.-Y., Modolo, R., Snowden, D., Winglee, R., 2011. Modeling of Venus, Mars and Titan, Space. Sci. Rev. 162, 267-307.

Korolov, I., Plasil, R., Kotrik, T., Dohnal, P., Glosik, J., 2009. Recombination of HCO+ and DCO+ ion with electrons. Int. J. Mass Spect. 280, 144-148. https://doi.org/10.1016/j.ijms.2008.07.023.

Krasnopolsky, V., 2000. On the deuterium abundance on Mars and some related problems. Icarus 148, 597-602.

Krasnopolsky, V., 2002. Mars' upper atmosphere and ionosphere at low, medium, and high solar activities: implications for evolution of water. J. Geophys. Res. 107 (E12), 5128. https://doi.org/10.1029/JE001809.

Krasnopolsky, V., 2010. Solar activity variations of thermospheric temperatures on Mars and a problem of CO in the lower atmosphere. Icarus 207, 638-647.

Krasnopolsky, V., 2015. Variations of the HDO/H₂O ratio in the martian atmosphere and loss of water from Mars. Icarus 257, 377-386.
## References

Anderson, D.E., Hord, C.W., 1971. Mariner 6 and 7 ultraviolet spectrometer experiment: analysis of hydrogen Lyman-Alpha data. J. Geophys. Res. 76, 6666-6673.

Aoki, S., Nakagawa, H., Sagawa, H., Giuranna, M., Sindoni, G., Aronica, A., Kasaba, Y., 2015. Seasonal variation of the HDO/H₂O ratio in the atmosphere of Mars at the middle of northern spring and beginning of northern summer. Icarus 260, 7-22.

Balakrishnan, N., Kharchenko, V., Dalgarno, A., 1998. Slowing of energetic O(3P) atom in collisions with N₂. J. Geophys. Res. 103, 23393-23398.

Belyaev, D.A., Fedorova, A.A., Trokhimovskiy, A., Alday, J., Montmessin, F., Korablev, O.I., et al., 2021. Revealing a high water abundance in the upper mesosphere of Mars with ACS onboard TGO. Geophys. Res. Lett. 48 (10) https://doi.org/10.1029/2021GL093411.

Benna, M., Mahaffy, P.R., Grebowsky, J.M., Fox, J.L., Yelle, R.V., Jakosky, B.M., 2015. First measurements of composition and dynamics of the Martian ionosphere by MAVEN's neutral gas and ion mass spectrometer. Geophys. Res. Lett. 42, 8958-8965. https://doi.org/10.1002/2015GL066146.

Bibring, J.-P.Y., Langevin, J.F. Mustard, Poulet, F., d'Arvidson, R., Gendrin, A., Gondet, B., Mangold, N., Pinet, P., Forget, F., the OMEGA team, 2006. Global mineralogical and aqueous Mars history derived from OMEGA/Mars express data. Science 312, 400-404.

Cangi, E., Chaffin, M., Yelle, R., Gregory, B., Deighan, J., 2023. Fully coupled photochemistry of the deuterated ionosphere of Mars and its effects on escape of H and D. J. Geophys. Res: Planet 128. https://doi.org/10.1029/2022JE007713 e2022JE007713.

Cangi, E., Chaffin, M.S., Yelle, R.V., Gregory, B.S., Deighan, J., 2024. Seasonal Enhancement in Upper Atmospheric D/H at Mars Driven by both Thermospheric Temperature and Mesospheric Water. Geophys. Res. Lett. p. 51.

Carr, M.H., 1996. Water on Mars. Oxford University Press, New York.

Chaffin, M.S., Chaufray, J.-Y., Stewart, I., Montmessin, F., Schneider, N., Bertaux, J.-L., 2014. Unexpected variability of Martian hydrogen escape. Geophys. Res. Lett. 41, 314-320. https://doi.org/10.1002/2013GL058578.

Chaffin, M.S., Chaufray, J.-Y., Deighan, J., Schneider, N.M., Mayyasi, M., Clarke, J.T., Thiemann, E., Jain, S.K., Crismani, M.M.J., Stiepen, A., Eparvier, F.G., McClintock, W.E., Stewart, A.I.F., Holswak, G.M., Montmessin, F., Jakosky, B., 2018. Mars H escape rates derived from MAVEN/IUVS Lyman-α brightness measurements and their dependence on model assumptions. J. Geophys. Res: Planets 123, 2192-2210.

Chaffin, M.S., Kass, D.M., Aoki, S., Fedorova, A.A., Deighan, J., Connour, K., et al., 2021. Regional dust storms on Mars enhance water loss to space. Nat. Astron. 5, 1036-1042.

Chassefière, E., Leblanc, F., 2004. Mars atmospheric escape and evolution: interaction with the solar wind. Planet. Space Sci. 52, 1039-1058 doi:1016/j.pss.204.07.002.

Chaufray, J.-Y., 2021. Departure of the thermal escape rate from the Jeans escape rate for atomic hydrogen at Earth, Mars and Pluto. Planet. Space Sci. 198, 105178 doi: 10.1016/j.pss.2021.105178.

Chaufray, J.-Y., Bertaux, J.-L., Leblanc, F., Quémérais, E., 2008. Observation of the hydrogen corona with SPICAM on Mars express. Icarus 195, 598-613.

Chaufray, J.-Y., Gonzalez-Galindo, F., Forget, F., Lopez-Valverde, M., Leblanc, F., Modolo, R., Hess, S., Yagi, M., Blelie, P.-L., Witasse, O., 2014. Three dimensional Martian ionosphere model: II effect of transport processes due to pressure gradients. J. Geophys. Res: Planets 119, 1614-1636. https://doi.org/10.1002/2013JE004551.

Chaufray, J.-Y., Gonzalez-Galindo, F., Forget, M.A. Lopez-Valverde, Leblanc, F., Modolo, R., Hess, S., 2015. Variability of the hydrogen in the martian upper atmosphere as simulated by a 3D atmosphere-exosphere coupling. Icarus 245, 282-294.

Chaufray, J.-Y., Modolo, R., Leblanc, F., Chanteur, G., Johnson, R.E., Luhmann, J.G., 2007. Mars solar wind interaction: Formation of the Martian corona and atmospheric loss to space. J. Geophys. Res. 112, E09009 https://doi.org/10.1029/2007JE002915.

Chaufray, J.-Y., Yelle, R.V., Gonzalez-Galindo, F., Forget, F., Lopez-Valverde, M.A., Leblanc, F., Modolo, R., 2018. Effect of the lateral exospheric transport on the horizontal hydrogen distribution near the exobase of Mars. J. Geophys. Res. Space Physics 123, 2441-2454. https://doi.org/10.1002/2017JA025163.

Chaufray, J.-Y., Gonzalez-Galindo, F., Lopez-Valverde, M.A., Forget, F., Quémérais, E., Bertaux, J.-L., Montmessin, F., Chaffin, M., Schneider, N., Clarke, J.T., Leblanc, F., Modolo, R., Yelle, R., 2021a. Study of the hydrogen escape rate at Mars during Martian years 28 and 29 from comparisons between SPICAM/Mars express observations and GCM-LMD simulations. Icarus 353, 113498. https://doi.org/10.1016/j.icarus.2019.113498.

Chaufray, J.-Y., Mayyasi, M., Chaffin, M., Deighan, J., Bhattacharyya, D., Clarke, J., Jain, S., Schneider, N., Jakosky, B., 2021b. Estimate of the D/H ratio in the Martian upper atmosphere from the low spectral resolution mode of MAVEN/IUVS. J. Geophys. Res.: Planets 126. https://doi.org/10.1029/JE006814 e2020JE006814.

Clarke, J.T., Mayyasi, M., Bhattacharyya, D., Schneider, N.M., McClintock, W.E., Deighan, J.I., Stewart, A.I.F., Chaufray, J.-Y., Chaffin, M.S., Jain, S.K., Stiepen, A., Crismani, M., Holswak, G.M., Montmessin, F., Jakosky, B., 2017. Variability of D and H in the Martian upper atmosphere observed with the MAVEN IUVS echelle channel. J. Geophys. Res. 122, 2336-2344. https://doi.org/10.1002/2016JA023479.

Cravens, T.E., Rahmati, A., Fox, J.L., Lillis, R., Bougher, S., Luhmann, J., Sakai, S., Deighan, J., Lee, Y., Combi, M., Jakosky, B., 2017. J. Geophys. Res. Space Physics 122, 1102-1116. https://doi.org/10.1002/2016JA023461.

Deighan, J., Chaffin, M.S., Chaufray, J.-Y., Stewart, A.I.F., Schneider, N.M., Jain, S.K., Stiepen, A., Crismani, M., McClintock, W.E., Clarke, J.T., Holswak, G.M., Montmessin, F., Eparvier, F.G., Thiemann, E.M.B., Chamberlin, P.C., Jakosky, B.M.,

2015. MAVEN IUVS observation of the hot oxygen corona. Geophys. Res. Lett. 42, 9009-9014. https://doi.org/10.1002/2015GL065487.

Dickeson, Z.I., Davis, J.M., 2020. Martian oceans. Astron. Geophys. 61, 311-317.

Dong, C., Bougher, S.W., Ma, Y., Lee, Y., Toth, G., Nagy, A.F., Fang, X., Luhmann, J., Liemohn, M.W., Halekas, J.S., Tenishev, V., Pawlowski, D.J., Combi, M.R., 2018. Solar wind interaction with the Martian upper atmosphere: roles of the cold thermosphere and hot oxygen corona, (2018). J. Geophys. Res. Space Physics 123, 6639-6654.

Elrod, M.K., Bougher, S., Bell, J., Mahaffy, P., Benna, M., Stone, S., Yelle, R., Jakosky, B., 2017. He bulge revealed: he and CO2 diurnal and seasonal variations in the upper atmosphere of Mars as detected by MAVEN NGIMS. J. Geophys. Res. Space Physics 122, 2564-2573. https://doi.org/10.1002/2016JA023482.

Ergun, R.E., Andersson, L.A., Fowler, C.M., Thaller, S.A., Yelle, R.V., 2021. In situ measurements of electron temperature and density in the Mars' dayside ionosphere. Geophys. Res. Lett. 48 e2021GL093623.

Fedorova, A., Montmessin, F., Korablev, O., Lefèvre, F., Trokhimovskiy, A., Bertaux, J.-L., 2021. Multi-annual monitoring of the water vapor vertical distribution on Mars by SPICAM on Mars express. J. Geophys. Res: Planets 126, e06616. https://doi.org/10.1029/2020JE006616.

Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M., Lewis, S.R., Read, P.L., Huot, J.-P., 1999. Improved general circulation models of the Martian atmosphere from the surface to above 80 km. J. Geophys. Res. 104, 24155-24175.

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

Fox, J., 2015. The chemistry of protonated species in the martian ionosphere. Icarus 252, 366-392.

Gacesa, M., Zhang, P., Kharchenko, V., 2012. Non-thermal escape of molecular hydrogen from Mars. Geophys. Res. Lett. 39, L10203. https://doi.org/10.1029/2012GL050904.

Gilli, G., Gonzalez-Galindo, F., Forget, F., Millour, E., Naar, J., Chaufray, J.-Y., 2022. On the effect of the orbital parameters of Mars to the hydrogen escape and the fate of water in the last millions of years. In: 16th EPSC, Abstract #578. https://doi.org/10.5194/epsc2022-578.

Gonzalez-Galindo, F., Forget, F., Lopez-Valverde, M.A., Coll, M. Angelats i, Millour, E., 2009. A ground-to-exosphere Martian general circulation model: 1 seasonal, diurnal, and solar cycle variation of thermospheric temperatures. J. Geophys. Res. 114, E04001. https://doi.org/10.1029/2008JE003246, 2009.

Gonzalez-Galindo, F., Chaufray, J.-Y., Lopez-Valverde, M.A., Gilli, G., Forget, F., Leblanc, F., Modolo, R., Hess, S., Yagi, M., 2013. Three-dimensional Martian ionosphere model: I. The photochemical ionosphere below 180 km. J. Geophys. Res.: Planets 118, 2105-2123. https://doi.org/10.1002/jgre.20150.

Gonzalez-Galindo, F., Lopez-Valverde, M.A., Forget, F., Garcia-Comas, M., Millour, E., Montabone, L., 2015. Variability of the Martian thermosphere during eight Martian years as simulated by a ground-to-exosphere global circulation model. J. Geophys. Res.: Planets 120, 2020-2035. https://doi.org/10.1002/2015JE004925.

Gregory, B., Elliott, R.D., Deighan, J., Gröller, H., Chaffin, M.S., 2023a. HCO+ dissociative recombination: a significant driver of nonthermal hydrogen loss at Mars. J. Geophys. Res.: Planets 128. https://doi.org/10.1029/2022JE007576 e2022JE007576.

Gregory, B., Chaffin, M.S., Elliott, R.D., Deighan, J., Gröller, H., Cangi, E., 2023b. Nonthermal hydrogen loss at Mars: contribution of photochemical mechanisms to escape and identification of key processes. J. Geophys. Res.: Planets 128. https://doi.org/10.1029/JE007802 e2023JE007802.

Halekas, J.S., 2017. Seasonal variability of the hydrogen exosphere of Mars. J. Geophys. Res.: Planets 122. https://doi.org/10.1002/2017JE005306.

Hanley, K.G., Fadden, J.P., Mitchell, D.L., Fowler, C.M., Stone, S.W., et al., 2021. In situ measurements of thermal ion temperature in the Martian ionosphere. J. Geophys. Res. Space Physics 126.

Hodges, R.R., Johnson, F.S., 1968. Lateral transport in planetary exospheres. J. Geophys. Res. Space Physics 73, 7307-7316.

Ishiguro, E., et al., 1952. On the polarizability of the hydrogen molecule. Proc. Phys. Soc. A 65, 178.

Jakosky, B., Phillips, R.J., 2001. Mars' volatile and climate history. Nature 412, 237-244.

Jakosky, B.M., Grebowsky, J.M., Luhmann, J.G., Connerney, J., Eparvier, F., Ergun, R., et al., 2015. MAVEN observations of the response of Mars to an interplanetary coronal mass ejection, Science, American, Association for the Adv. Sci. 350 (6261) aad0210-1-aad020-7.

Jakosky, B., et al., 2018. Loss of the Martian atmosphere to space: present-day loss rates determined from MAVEN observations and integrated loss through time. Icarus 315, 146-157.

Kallio, E., Chaufray, J.-Y., Modolo, R., Snowden, D., Winglee, R., 2011. Modeling of Venus, Mars and Titan, Space. Sci. Rev. 162, 267-307.

Korolov, I., Plasil, R., Kotrik, T., Dohnal, P., Glosik, J., 2009. Recombination of HCO+ and DCO+ ion with electrons. Int. J. Mass Spect. 280, 144-148. https://doi.org/10.1016/j.ijms.2008.07.023.

Krasnopolsky, V., 2000. On the deuterium abundance on Mars and some related problems. Icarus 148, 597-602.

Krasnopolsky, V., 2002. Mars' upper atmosphere and ionosphere at low, medium, and high solar activities: implications for evolution of water. J. Geophys. Res. 107 (E12), 5128. https://doi.org/10.1029/JE001809.

Krasnopolsky, V., 2010. Solar activity variations of thermospheric temperatures on Mars and a problem of CO in the lower atmosphere. Icarus 207, 638-647.

Krasnopolsky, V., 2015. Variations of the HDO/H₂O ratio in the martian atmosphere and loss of water from Mars. Icarus 257, 377-386.
Krasnopolsky, V., 2019. Photochemistry of water in the martian thermosphere and its effect on hydrogen escape. Icarus 321, 62–70*.

Krasnopolsky, V., Feldman, P.D., 2001. Detection of molecular hydrogen in the atmosphere of Mars. Science 294, 1914–1917.

Krasnopolsky, V., Bjoraker, G.L., Mumma, M.J., 1997. And D.E. Jennings, high resolution spectroscopy of Mars at 3.7 and 8 µm: a sensitive search of H₂O₂, H₂CO, HCl, and CH₄, and detection of HDO. J. Geophys. Res. 102, 6525–6534.

Laskar, J., Correia, A.C.M., Gastineau, M., Joutel, F., Levrard, B., Robutel, P., 2004. Long term evolution and chaotic diffusion of the insolation quantities of Mars. Icarus 170, 343–364. https://doi.org/10.1016/j.icarus.2004.04.005.

Lasue, J., Mangold, N., Hauber, E., Clifford, S., Feldman, W., Gasnault, O., Grima, C., Maurice, S., Mousis, O., 2013. Quantitative assessments of the Martian hydrosphere. Space Sci. Rev. 174, 155–212.

Leblanc, F., Chaufray, J.-Y., Modolo, R., Leclercq, L., Curry, S., Luhmann, J., Lillis, R., Hara, T., McFadden, J., Halekas, J., Schneider, N., Deighan, J., Mahaffy, P.R., Benna, M., Johnson, R.E., Gonzalez-Galindo, F., Forget, F., Lopez-Valverde, M.A., Eparvier, F.G., Jakosky, B., 2017. On the origins of Mars' exospheric nonthermal oxygen component as observed by MAVEN and modeled by HELIOSARES. J. Geophys. Res.: Planets 122, 2401–2428. https://doi.org/10.1002/2017JE005316.

Leblanc, F., Benna, M., Chaufray, J.-Y., Martinez, A., Lillis, R., Curry, S., Elrod, M.K., Mahaffy, P.R., Modolo, R., Luhmann, J.G., Jakosky, B., 2019. First in-situ evidence of Mars non-thermal exosphere. Geophys. Res. Lett. 46, 4144–4150. https://doi.org/10.1029/2019GL082192.

Leshin, L.A., 2000. Insight into martian water reservoirs from analyses of martian meteorite QUE94201. Geophys. Res. Lett. 27, 2017–2020.

Lewkow, N.R., Kharchenko, V., 2014. Precipitation of energetic neutral atoms and induced non-thermal escape fluxes from the Martian atmosphere. Astrophys. J. 790, 98. https://doi.org/10.1088/0004-637X/790/2/98.

Lillis, R.J., et al., 2017. Photochemical escape of oxygen from Mars: first results from MAVEN in situ data. J. Geophys. Res. Space Physics 122. https://doi.org/10.1002/2016/JA023525.

Mahaffy, P.R., et al., 2015. MSL Science Team, Mars atmosphere: the imprint of atmospheric evolution in the D/H of Hesperian clay minerals on Mars. Science 347, 412–414.

Maltagliati, L., Montmessin, F., Korablev, O., Fedorova, A., Forget, F., Määtänen, A., et al., 2013. Annual survey of water vapor vertical distribution and water-aerosol coupling in the martian atmosphere observed by SPICAM/MEx solar occultations. Icarus 223, 942–962. https://doi.org/10.1016/j.icarus.2012.12.012.

Mayyasi, M., Clarke, J.T., Bhattacharyya, D., Deighan, J., Jain, S., Chaffin, M., Thiemann, E., Schneider, N., Jakosky, B., 2017. The variability of atmospheric deuterium brightness at Mars: evidence for seasonal dependence. J. Geophys. Res. Space Physics 122, 10811–10823.

Mayyasi, M., Clarke, J., Chaufray, J.-Y., Kass, D., Bougher, S., Bhattacharyya, D., Deighan, J., Jain, S., Schneider, N., Villanueva, G.L., Montmessin, F., Benna, M., Mahaffy, P., Jakosky, B., 2023. Solar cycle and seasonal variability of H in the upper atmosphere of Mars. Icarus 393, 115293. https://doi.org/10.1016/j.icarus.2022.115293.

Modolo, R., Hess, S., Mancini, M., Leblanc, F., Chaufray, J.-Y., Brain, D., Leclercq, L., Esteban-Hernandez, R., Chanteur, G., Weill, P., Gonzalez-Galindo, F., Forget, F., Yagi, M., Mazelle, C., 2016. Mars-solar wind interaction: LatHyS, an improved parallel 3D multispecies hybrid model. J. Geophys. Res. Space Physics 121, 6378–6399. https://doi.org/10.1002/2015JA022324.

Montmessin, F., Fouchet, T., Forget, F., 2005. Modeling the annual cycle of HDO in the Martian atmosphere. J. Geophys. Res. 110, E03006 https://doi.org/10.1029/2004JE002357.

Naar, J., 2023. Modélisation du cycle de l'eau actuel et des derniers âges glaciaires de la planète Mars. Phd Thesis. https://doi.org/10.13140/RG.2.2.31506.84162.

Nakai, Y., Shirai, T., Tabata, T., Ito, R., 1987. Cross section for charge transfer of hydrogen atoms and ions colliding with gaseous atoms and molecules. Atom. Data and Nuc. Data Tab. 37, 69–101.

Navarro, R., Madeleine, J.-B., Forget, F., Spiga, A., Millour, E., Montmessin, F., Määtänen, A., 2014. Global climate modeling of the Martian water cycle with improved microphysics and radiatively active water ice clouds. J. Geophys. Res.: Planets 119, 1479–1495. https://doi.org/10.1002/2013JE004550.

Novak, R.E., Mumma, M.J., Villanueva, G.L., 2011. Measurement of the isotopic signatures of water on Mars; implications for studying methane. Planet. Space Sci. 59, 163–168.

Owen, T., Maillard, J.-P., De Bergh, C., Lutz, B.L., 1988. Deuterium on Mars: the abundance of HDO and the value of D/H. Science 240, 1767–1770.

Qin, J., 2021. Solar cycle, seasonal, and dust-storm driven variations of the Mars upper atmospheric state and H escape rate derived from the Lyα emission observed by NASA's MAVEN mission. Astrophys. J. 912, 77. https://doi.org/10.3847/1538-4357/abed4f.

Rahmati, A., et al., 2018. Seasonal variability of neutral escape from Mars as derived from MAVEN pickup ion observations. J. Geophys. Res.: Planets 123, 1192–1202. https://doi.org/10.1029/2018JE005560.

Romanelli, N., Modolo, R., Leblanc, F., Chaufray, J.-Y., Martinez, A., Ma, Y., Lee, C.O., Luhmann, J.G., Halekas, J., Brain, D., DiBraccio, G., Espley, J., McFadden, J., Jakoskky, B., Holmstrom, M., 2018. Responses of the Martian magnetosphere to an interplanetary coronal mass ejection: MAVEN observations and LatHyS results. Geophys. Res. Lett. 45, 7891–7900. https://doi.org/10.1029/2018GL077714.

Rossi, L., Vals, M., Alday, J., Montmessin, F., Fedorova, A., Trokhimovskyi, A., Korablev, O., Lefèvre, F., Gonzalez-Galindo, F., Luginin, M., Bierjon, A., Forget, F., Millour, E., 2022. The HDO cycle on Mars: comparison of ACS observations with GCM simulations. J. Geophys. Res.: Planets 127. https://doi.org/10.1029/2022JE007201.

Schmidt, F., Way, M.J., Costard, F., Bouley, S., Séjourné, A., Aleinov, I., 2022. Circumpolar Ocean stability on Mars 3 Gy ago. PNAS 119. https://doi.org/10.1073/pnas.2112930118.

Shaposhnikov, D.S., Medvedev, A.S., Rodin, A.V., Hartogh, P., 2019. Seasonal water “pump” in the atmosphere of Mars: vertical transport to the thermosphere. Geophys. Res. Lett. 46, 4161–4169.

Shematovich, V.I., 2013. Suprathermal oxygen and hydrogen atoms in the upper martian atmosphere. Sol. Syst. Res. 47, 437–445.

Stone, S.W., Yelle, R.V., Benna, M., Elrod, M.K., Mahaffy, P.R., 2022. Neutral composition and horizontal variations of the Martian upper atmosphere from MAVEN/NGIMS. J. Geophys. Res. 127 https://doi.org/10.1029/2021JE007085 e2021JE007085.

Touma, J., Wisdom, J., 1993. The chaotic obliquity of Mars. Science 259, 1294–1297.

Vals, M., Rossi, L., Montmessin, F., Lefèvre, F., Gonzalez-Galindo, F., Fedorova, A., Luginin, M., Forget, F., Millour, E., Korablev, O., Trokhimovskyi, A., Shakun, A., Bierjon, A., Montabone, L., 2022. Improved modeling of Mars' HDO cycle using a Mars global climate model. J. Geophys. Res.: Planets 127. https://doi.org/10.1029/2022JE007192 e2022JE007192.

Vidal-Madjar, A., Bertaux, J.-L., 1972. A calculated hydrogen distribution in the exosphere. Planet. Space Sci. 20 (8), 1147–1162 doi:10/1016/0032-0633(72)90004-9.

Villanueva, G.L., Mumma, M.J., Novak, R.E., Käufl, H.U., Hartogh, P., Encrenaz, T., Tokugana, A., Khayat, A., Smith, M.D., 2015. Strong water isotopic anomalies in the martian atmosphere: probing current and ancient reservoirs. Science 348, 218–221.

Webster, C.R., Mahaffy, P.R., Flesch, G.J., Niles, P.B., Jones, J.H., Leshin, L.A., Atreya, S. K., Stern, J.C., Christensen, L.E., Owen, T., Franz, H., Pepin, R.O., Steele, A., 2013. Isotope ratios of H, C, and O in CO2 and H2O of the Martian atmosphere. Science 341, 260–263.

Zhang, P., Kharchenko, V., Jamieson, M.J., Dalgarno, A., 2009. Energy relaxation in collisions of hydrogen and deuterium with oxygen atoms. J. Geophys. Res. 114, A07101. https://doi.org/10.1029/JA014055.