FULLY-COUPLED PHOTOCHEMICAL MODELING OF THE DEUTERATED
IONOSPHERE AND NON-THERMAL ESCAPE OF D.
E. M. Cangi, Laboratory for Atmospheric and Space Physics (LASP), University of Colorado Boulder, USA
(eryn.cangi@colorado.edu), M. S. Chaffin, J. Deighan, B. Gregory, LASP, Colorado, USA, R. V. Yelle, Lunar and
Planetary Laboratory, Arizona, USA.

Introduction
Context: The Mars water D/H ratio is 5-6 × the Earth

ratio, and is considered to be an indicator of past water
loss on Mars. Because the isotope D is heavier than H,
thermal escape is less effective for D, and escape of D
is expected to occur primarily via non-thermal escape
driven by ion processes. Neutral H loss from both thermal and non-thermal processes is well-studied at Mars,
but the deuterated ionosphere is relatively unstudied.
This leaves us with a lack of understanding of the true
scale of non-thermal D escape and its variation throughout time.
Aims: To fill this knowledge gap, we have performed
a modeling study of the deuterated ionosphere at Mars.
Our primary goals are to characterize the overall composition of the deuterated ionosphere, quantify the thermal
and non-thermal escape of H and D, and understand associated changes to the exchangeable water inventory
over time.
Methods: To accomplish our goals, we simulated a
long-term equilibrium state of the atmosphere using a
1D photochemical model and atmospheric parameters
representing the present-day mean atmosphere on Mars.
Studying the atmosphere in equilibrium allows us to
establish a mean probable magnitude of non-thermal escape of D. To do this effectively, we have upgraded our
existing 1D photochemical model to be the first such
model to couple the lower and upper atmospheres and
the neutral and ion chemistry–without assuming photochemical equilibrium or fixing a background atmosphere.
Results: We find that DCO+ is the dominant deuterated ion at all altitudes, with D+ , OD+ , HDO+ , and
DCO+
2 also important, similar to their H-bearing analogues. Our density profiles for the H-bearing analogue
ions are generally in agreement with Fox (2015). For
non-thermal escape of D during solar mean conditions,
we find an upper limit of ∼ 1.6 × 104 cm−2 s−1 . A detailed analysis of the contributions of different processes
and implications for water loss will be discussed in the
presentation.

Motivating questions about non-thermal D escape
Many open questions remain about D at Mars. The
MAVEN mission is addressing questions about the upper

atmospheric D/H ratio measured in the atomic species
(Clarke et al., 2021), and past results from many missions have characterized the D/H ratio in water (Alday
et al., 2021; Encrenaz et al., 2018; Fedorova et al., 2020;
Vandaele et al., 2019; Villanueva et al., 2021). D ions
must exist on Mars, but have been neither modeled nor
observed. In this work, we seek to use our 1D model
to simulate the detailed D chemistry and answer the following:
• What are the abundances of deuterated ions in the

martian atmosphere?
• Which chemical reactions dominate D escape?
• How do non-thermal processes affect D/H in the

martian atmosphere?
• What are the relative contributions of thermal and

non-thermal processes to D escape, and how does
this balance differ from H escape?

1D photochemical modeling of deuterium ions
Neutrals and ions behave very differently: they have
vastly different equilibration timescales, different density distributions, and different primary production pathways. 1D photochemical models that include both neutrals and ions must thus find ways to solve large systems
of coupled equations for these very different populations. Most tend to either (1) have a lower boundary in
the middle atmosphere (80-100 km) (V. A. Krasnopolsky, 2019, and references therin), (2) hold many species
in the background neutral atmosphere fixed (Fox, 2015),
or (3) assume that some or all short-lived species (such as
ions) are in photochemical equilibrium (Banaszkiewicz
et al., 2000; Dobrijevic et al., 2016; Fox, 2015). These
approaches improve computation time, but risk losing
fine details, precluding a complete understanding of the
connections between D and H chemistry, escape, and
lower-upper atmospheric coupling on Mars.
To achieve our goals and answer our scientific questions, we have extended the model from our previous
paper (Cangi et al., 2020) to be the first fully coupled,
surface-to-upper-atmosphere model of the neutral atmosphere and ionosphere that does not require any of the
three previously described modeling approaches. The

model is built in Julia 1.7 (Bezanson et al., 2017), a language designed to combine the fast numerical modeling
of languages like MATLAB and FORTRAN and the
user-friendliness of Python. The original version of the
model included 23 neutral species, a fixed background
ionosphere of CO+
2 , a fixed water profile, and 117 chemical reactions. Our recent improvements include 13 new
neutral and 43 new ion species, 900 new ion and neutral
reactions, ambipolar diffusion for ion vertical transport,
optional photochemical equilibrium for user-specifiable
species, a water profile that is fixed in the lower atmosphere but freely solved in the middle and upper atmosphere, and of course, calculation of non-thermal escape for H and D. Computation-specific improvements
include updated ODE solver algorithm options, vectorization of repeatedly called functions, and other general
best-practice optimizations. Using the model, we simulate the atmospheric chemistry, photochemistry, and
vertical diffusion in 2 km slices between 0-250 km. We
choose a set of parameters and boundary conditions for
the nominal Mars atmosphere (intended to represent a
baseline). As the model runs, the simulated atmosphere
evolves forward in time until it reaches a state of chemical and diffusive equilibrium, usually within 10 million
years.
Model inputs: The primary model inputs include atmospheric temperatures at the surface, mesosphere, and
base of the exosphere (i.e., roughly the neutral exobase)
and the total atmospheric water content. Many other
model parameters are modifiable, but typically don’t
change once we have arrived at reasonable choices for
each, such as solver-specific options (e.g. number of
timesteps, tolerance limits).
For a baseline, mean atmosphere, we assume a neutral temperature with Tsurface = 216 K, Tmesosphere = 130
K, and Texobase = 205 K (Cangi et al., 2020), with ion
and electron temperature profiles based on the Deep Dip
8 MAVEN data as described by Hanley et al. (2021, Figure 7a). The baseline water content of the atmosphere
is 10 pr µm.
Model outputs: Our model output comprises atmospheric species densities (cm−3 ) at each layer of the
atmosphere at discrete simulation times and for the final
converged state. From this output, we can calculate and
plot quantities relevant to answering our science questions, including:
• Densities of deuterated ions by altitude and species
• Chemical reaction rates by altitude and species
• Final resulting thermal and non-thermal escape

rates of D, either cumulative or broken down by
process
The escape rates are the primary output of interest. While the model does support calculation of hot
H escape, we focus in this work on hot D escape. For

thermal escape, we consider Jeans escape only and assume that all atoms with a velocity falling in the tail of
the Maxwell-Boltzmann distribution (v ≥ vesc ) will escape. Non-thermal escape is a bit more complicated, as
it depends on collisions, and therefore the ever-changing
density distributions of other species.
How we calculate non-thermal escape
Although non-thermal H escape at Mars has been studied
before (Nagy et al., 1990), there are no studies available
in the literature that quantify production of superthermal (sometimes called “hot”) D atoms (those with high
velocities that are not well described by the standard
Maxwell-Boltzmann distribution). V. A. Krasnopolsky
(2002) has provided estimates of certain non-thermal
escape processes affecting this population, but only includes a few processes and excludes a detailed explanation of the calculations. This presents an opportunity,
as non-thermal escape of D is expected to be far more
important than thermal escape (Gacesa et al., 2018; V.
Krasnopolsky, 2000). By applying the non-thermal escape velocity values from V. A. Krasnopolsky (2002)
to their own model output, Cangi et al. (2020) showed
that including non-thermal escape of D in calculations
of the fractionation factor could change it by an order of
magnitude and alter the estimation of integrated water
loss by 10s of m GEL, indicating a need for updated,
targeted modeling of non-thermal D loss.
Work to characterize the hot H production mechanisms in detail is presently underway (see Gregory et
al., 2022, this meeting). In this study, we consider hot D
produced via four main processes: (1) charge exchange
reactions such as D+ + H → Dhot + H+ , (2) photoionization reactions such as H2 O + hν → H+ + O2 , and (3)
dissociative recombination such as DCO+ + e− → CO
+ D, and (4) other exothermic bimolecular ion reactions
that produce a lone atomic D. To calculate how much
escaping hot D is produced, we first calculate (from
the model output) the volume production rate Vproduced
of any and all atomic D produced by exothermic reactions (Fox, 2015), generally falling under one of the four
categories described previously. We then estimate the
probability P of that atomic D to both be hot and to
escape as follows:
P = Ae−aσN (z)

(1)

Where σ is a mean collisional cross section (Zhang
et al., 2009), N (z) is the column density above altitude z ,
and A and a are parameters determined using by fitting to
Monte Carlo simulations of hot H escape (see Gregory et
al., this meeting). The volume production rate Vproduced
multiplied by this probability gives a volume production
rate of escape-capable hot H and D:
Vesc [cm−3 s−1 ] = P Vproduced

(2)

Figure 1: Volume production rates of escaping hot D from the
top 5 chemical reactions. DCO+ dissociative recombination is
the most important mechanism up to 180 km, at which point
charge exchange of D+ with CO2 briefly takes over before
giving way to O+ charge exchange with HD.

We perform this calculation for each chemical reaction, and the results can also be totaled. The expressions
are vectors in altitude; by integrating them (summing
over all altitudes), we arrive at a total non-thermal escape flux (cm−2 s−1 )that we use to define an additional
flux boundary condition at the upper boundary of the
model. This flux is a function of the atmospheric conditions, so it constantly updates as the model runs, and its
value at the end of the simulation reflects non-thermal
loss in an equilibrium state.

Figure 2: Volume production rates of escaping hot D from the
top 5 chemical reactions. DCO+ dissociative recombination is
the most important mechanism up to 180 km, at which point
charge exchange of D+ with CO2 briefly takes over before
giving way to O+ charge exchange with HD.

Results
Density profiles for D-bearing ions in our model are
shown in Figure 1. For our mean atmosphere, we find
that DCO+ is the dominant deuterated ion at all altitudes, with D+ , OD+ , and HDO+ the most common
deuterated ions near the base of the exosphere. DCO+
2
is also present in non-negligible concentrations at most
altitudes. The order of these ions’ appearances is similar
to their H-bearing analogues, and the presented density
profiles are useful for guiding possible future detection
attempts.

The dominant chemical reactions contributing to
non-thermal D escape are shown in Figure 2, and the
total production of escaping non-thermal H and D in Figure 3. We find an upper limit for non-thermal escape of
D to be ∼ 1.7×104 cm−2 s−1 , double the 2×102 cm−2 s−1
thermal escape reported by Cangi et al. (2020) for similar atmospheric conditions (Figure S5). The solar cycle
conditions are for solar mean (Woods et al., 2019); it

Figure 3: Volume production rate of escaping hot H and D
from all contributing chemical reactions.

REFERENCES

has been predicted (V. A. Krasnopolsky et al., 1998)
that non-thermal escape should dominate for D at solar
minimum only. Our model does not account for hot H
and D that receive less than all the excess heat produced
in exothermic reactions. With this thermal branching
ratio included, we expect our calculation of non-thermal
escape of D to be less, and likely more comparable to
thermal escape at solar mean as predicted.
Ongoing work
A more thorough discussion of the results and methods
will occur at the in-person presentation at the meeting.
Ongoing and future work will include similar modeling generated for solar minimum and solar maximum
conditions, as well as seasonal cycling.
References
Alday, Juan et al. (2021). “Isotopic fractionation of water and
its photolytic products in the atmosphere of Mars”. In: Nature Astronomy 5, pp. 943–950. DOI: 10.1038/s41550021-01389-x.
Banaszkiewicz, M. et al. (2000). “A Coupled Model of Titan’s
Atmosphere and Ionosphere”. In: Icarus 147.2, pp. 386–
404. DOI: 10.1006/icar.2000.6448.
Bezanson, Jeff et al. (2017). “Julia: A fresh approach to numerical computing”. In: SIAM review 59.1, pp. 65–98. URL:
https://doi.org/10.1137/141000671.
Cangi, E. M., M. S. Chaffin, and J. Deighan (2020). “Higher
Martian Atmospheric Temperatures at All Altitudes Increase
the D/H Fractionation Factor and Water Loss”. In: Journal of Geophysical Research (Planets) 125.12. DOI: 10 .
1029/2020JE006626.
Clarke, J. T. et al. (2021). “The Martian D/H Ratio In Escaping
Atoms and Evidence for Primordial Water”. In: AGU Fall
Meeting Abstracts. Vol. 2021.
Dobrijevic, M. et al. (2016). “1D-coupled photochemical model
of neutrals, cations and anions in the atmosphere of Titan”. In: Icarus 268, pp. 313–339. DOI: 10 . 1016 / j .
icarus.2015.12.045.
Encrenaz, T. et al. (2018). “New measurements of D/H on
Mars using EXES aboard SOFIA”. In: A&A 612, A112.
DOI : 10.1051/0004-6361/201732367.

Fedorova, Anna A. et al. (2020). “Stormy water on Mars: The
distribution and saturation of atmospheric water during the
dusty season”. In: Science 367.6475, pp. 297–300. DOI:
10.1126/science.aay9522.
Fox, Jane L. (2015). “The chemistry of protonated species in
the martian ionosphere”. In: Icarus 252, pp. 366–392. DOI:
10.1016/j.icarus.2015.01.010.
Gacesa, Marko et al. (2018). “Non-thermal escape rates of
light species from Mars using MAVEN in-situ measurements”. In: European Planetary Science Congress, EPSC2018–
604.
Hanley, K. G. et al. (2021). “In Situ Measurements of Thermal
Ion Temperature in the Martian Ionosphere”. In: Journal of
Geophysical Research (Space Physics) 126.12. DOI: 10 .
1029/2021JA029531.
Krasnopolsky, Vladimir (2000). “On the Deuterium Abundance
on Mars and Some Related Problems”. In: Icarus 148.
Krasnopolsky, Vladimir A. (2002). “Mars’ upper atmosphere
and ionosphere at low, medium, and high solar activities:
Implications for evolution of water”. In: Journal of Geophysical Research (Planets) 107.E12, p. 5128. DOI: 10 .
1029/2001JE001809.
– (2019). “Photochemistry of water in the martian thermosphere and its effect on hydrogen escape”. In: Icarus 321,
pp. 62–70. DOI: 10.1016/j.icarus.2018.10.033.
Krasnopolsky, Vladimir A., Michael J. Mumma, and G. Randall Gladstone (1998). “Detection of Atomic Deuterium in
the Upper Atmosphere of Mars”. In: Science 280, p. 1576.
DOI : 10.1126/science.280.5369.1576.
Nagy, A. F., J. Kim, and T. E. Cravens (1990). “Hot hydrogen
and oxygen atoms in the upper atmospheres of Venus and
Mars.” In: Annales Geophysicae 8, pp. 251–256.
Vandaele, Ann Carine et al. (2019). “Martian dust storm impact on atmospheric H2 O and D/H observed by ExoMars
Trace Gas Orbiter”. In: Nature 568.7753, pp. 521–525. DOI:
10.1038/s41586-019-1097-3.
Villanueva, Geronimo L. et al. (2021). “Water heavily fractionated as it ascends on Mars as revealed by ExoMars/NOMAD”.
In: Science Advances 7.7. DOI: 10 . 1126 / sciadv .
abc8843.
Woods, Thomas N. et al. (2019). LISIRD (LASP Interactive
Solar Irradiance Datacenter). Version 3.30.0. URL: http:
//lasp.colorado.edu/lisird/.
Zhang, P. et al. (2009). “Energy relaxation in collisions of hydrogen and deuterium with oxygen atoms”. In: Journal of
Geophysical Research (Space Physics) 114.A7, A07101,
A07101. DOI: 10.1029/2009JA014055.