The Effect of the Martian 2018 Global Dust Storm on HDO as Predicted by a Mars Global Climate Model

The deuterium to hydrogen (D/H) ratio is commonly used to investigate the history of water on Mars, yet the mechanisms controlling present‐day HDO behavior are poorly understood. Significant variations of the D/H ratio were first predicted on the basis of a 3D global climate model, which were later confirmed by ground‐based observations. This behavior, consisting of lower HDO/H2O ratios in the colder regions of Mars, is related to the isotopic fractionation occurring at condensation. We leverage this previous effort and present an updated implementation, using the modern version of the model, that remains in agreement with the older version. We explore the impact of the global dust storm (GDS) that occurred during Martian year 34 (MY34) on HDO. Our simulations indicate that HDO is on average 40% more abundant at 100 km during the MY34 GDS year than during a regular year, with likely large consequences for the escape flux of water that year.

HDO and the warmer, condensation-free, regions where D/H is found to be maximum. This led M05 to predict a latitudinal gradient of D/H (with variations greater than a factor of 5) between the warm and moist summer hemisphere and the cold and dry winter hemisphere. This gradient was later confirmed by several Earth ground-based and Mars orbiting assets Encrenaz et al., 2018;Krasnopolsky, 2015;Villanueva et al., 2015). Yet some observations (Khayat et al., 2019;Villanueva et al., 2015) have revealed longitudinal variations of D/H ratios in appearance correlated to topography that were not reproduced by any model.
Since HDO and H 2 O dictate eventually how many D and H atoms will populate the Martian exosphere, studying their escape and extrapolating back in times to infer the amount of water lost to space requires modeling their present-day behavior correctly in the lower atmosphere. It is therefore relevant to address the HDO cycle and its associated fractionation with respect to water in a 3D Mars climate model (GCM) akin to M05. The present work leverages on this past effort, and implements important changes that occurred between the 2005 GCM version and the one currently in use. Upgrading the HDO model enabled the possibility to test the theoretical reaction of HDO during the most extreme meteorologic event that can happen on Mars, namely a global dust storm (GDS). The Martian year 34 (henceforth MY34, corresponding to the 2017-2019 time frame) was the siege of a GDS that started in the middle of the year (L s = 176°) and stopped shortly before perihelion (L s = 235°) (Guzewich et al., 2018). This event profoundly affected the behavior of water vapor in the Martian atmosphere (Fedorova et al., 2020) with hints of a comparable affect for HDO .
The present study describes a 3D climate model projection of the behavior of HDO in the context of a dust annual scenario mimicking the dust seasonal and spatial evolution observed during MY34 which include the occurrence of the GDS. The simulations reported here emphasize the contrast in HDO behavior between a perturbed year such as MY34 and a more traditional MY where the seasonal evolution of dust obeys a recurrent and more quiet, so-called "climatological" scenario.

Model Overview
The simulations presented here were conducted with the Martian Global Climate Model (MGCM) developed at the Laboratoire de Météorologie Dynamique (LMD) in collaboration with several European teams (LATMOS, IAA Granada, University of Oxford, The Open University), with support from ESA and CNES. A general description of the model is given in Forget et al. (1999). A first implementation of the water cycle was presented in Montmessin et al. (2004) (used in M05) and was then improved in Madeleine et al. (2011) and subsequently in Navarro et al. (2014), with the implementation of the radiative effect of clouds and the microphysics accounting for the fine processes of cloud formation as nucleation, ice growth, scavenging of dust by condensing ice, and supersaturation. However, those were not activated for the simulations conducted in this study. This is motivated by the desire to track changes added to the GCM since M05 before increasing model complexity. The simpler cloud approach of M05 employed here is nevertheless augmented by the predictive dust representation of Madeleine et al. (2011), which better constrains the availability of dust nuclei for water ice clouds to condense onto, even if the details of cloud-dust interactions (nucleation, scavenging, and supersaturation as mentioned above) are not yet introduced. Observations from Earth (Khayat et al., 2019;Novak et al., 2011;Villanueva et al., 2015) are in agreement with predictions of M05 regarding the latitudinal variability of the D/H ratio over the column of atmosphere, supporting the idea that a simple cloud model carries enough physics to capture the effect of condensation-induced fractionation on HDO. We also mention that since Colaïtis et al. (2013) the model integrates a thermal plume model to better represent the sub-grid scale convection occurring in the planetary boundary layer (PBL). This parametrization has been shown to significantly improve the wind and temperature structures near to the surface in comparison to observations and LES studies, which has a direct impact on the transport of water and dust above the PBL.
In this upgraded version of the HDO model, the vertical distribution of dust obeys a semi-interactive scheme (Madeleine et al., 2011) based on an order two moment representation of the dust tracer. Dust is injected everywhere and at every numerical time-step (15 min in our simulations) at the surface of the model by a constant lifting rate. Both the mass and the number of dust particles are transported, enabling the interactive derivation of the dust particle size in every model cell to better constrain dust sedimentation. While the vertical profile of the dust evolves freely, sustaining changes due to winds and gravitational settling, the resulting opacity profile is multiplied by a scaling factor so that the column dust opacity of the model matches the observed column dust opacity as compiled by Montabone et al. (2015Montabone et al. ( , 2020. The latter being the main parameter controlling the radiative balance of the atmosphere and thus the background thermal state. It also linearly scales the number of dust nuclei available for water ice cloud formation and thus ice crystal size and subsequent gravitational fall.
The MGCM used for our HDO simulations relies on a 64 × 48 longitude-latitude grid with 32 vertical levels, covering the atmospheric column from the surface up to 120 km (10 −3 Pa). An extension to the thermosphere up to the exobase around 250 km (Angelats i Coll et al., 2005;González-Galindo et al., 2009) and a photochemical module with 15 species (Lefèvre et al., 2004(Lefèvre et al., , 2008 are available in the model. However, these modules have not been activated in the simulations used for the present paper, since their integration will be part of a future study analyzing the whole HDO and D/H cycle.

Simulation of HDO
HDO is tracked in both its vapor and ice phases in the atmosphere and also on the surface as ice deposit. It is treated as a tracer separate from H 2 O, but undergoes the same condensation and transport processes. Isotopic fractionation is introduced during condensation following the approach described by Fouchet and Lellouch (2000), Bertaux and Montmessin (2001) and later by M05.
The fractionation factor α has been traditionally derived from the formula established by Merlivat and Nief (1967). This factor describes the relative enrichment of HDO in the ice with respect to the surrounding vapor phase at equilibrium (i.e., no net flux between the two phases), for a temperature T, and can be expressed as: We note that Equation 1 is based on experiments performed at temperatures above 233 K, warmer than typical temperatures encountered on Mars, leading M05 to extrapolate Equation 1 down to temperatures as low as 100 K. However, recent experimental measurements from Lamb et al. (2017) were conducted between 234 and 194 K. In doing so, they updated the fractionation law and proposed a new formula: The difference between the two formulas increases exponentially with decreasing temperature, from 3% at 200 K, to 8% at 150 K and 21% at 100 K ( Figure 1). The Lamb formula yields a fractionation factor systematically smaller than the classical Merlivat formula at Martian temperatures, which suggests the studies conducted so far on Martian HDO have slightly overestimated the fractionation factor. Lamb's formula is used as a reference for this study.
We assume that the condensation flux is in isotopic equilibrium with the vapor phase. This is justified by the fact that diffusivity of HDO inside ice is too slow to homogenize HDO inside the ice particle. Therefore, the newly added icy layer is independent of the existing ones. In that context, the condensing mass dM D of HDO can be expressed as: where M indicates the mass in the vapor phase, dM the mass condensing, with subscript H referring to water and D to HDO.
In the model, the flux of water during cloud condensation is computed assuming that any supersaturated excess of water vapor is turned into atmospheric water ice . At the surface, the flux of water is computed with the surface turbulent flux (Forget et al., 1999;Montmessin et al., 2004). In both cases, once the flux of water is known, the flux of HDO is computed using Equation 3. The temperature used to compute α(T) is that of the atmospheric layer condensing (for the particular case of the surface flux, it corresponds to the temperature of the first layer above the surface).
A significant difference with M05 is the implementation of the "tracer genealogy" scheme developed by Risi et al. (2010) to study the HDO atmospheric cycle in the LMDZ Earth GCM. M05 had an older version of the dynamical core of the GCM that did not include or require this specific transport of isotopes in the tracer scheme. In this scheme, the isotopic ratio is transported by the dynamical transport scheme instead of HDO itself (more details about this transport scheme can be found in the appendix of Risi et al., 2010, and in the Supplementary Information). The implementation of the "tracer genealogy" scheme has significantly improved the GCM results by eliminating the unphysical isotopic ratios initially observed.

Dust Scenarios
We employ two dust scenarios for this study. We are calling "dust scenario" the prescription in the model of the seasonal evolution and the spatial variation of the dust column opacity. The scenarios for each MY have been compiled by Montabone et al. (2015), and the dust scenario relevant to MY34 that encompasses the 2018 GDS was recently added by Montabone et al. (2020) to the existing database.
The first one is the climatology scenario, which corresponds to an average of non-GDS Martian years, and that is used as our reference run. The other one is the dust scenario of MY34, which includes the dust storm occurrence between L s = 176° and L s = 235°.

Initial Conditions
The initial conditions of the simulation impose a constant D/H ratio of 5 everywhere in the atmosphere and in the surface ice. The perennial ice in the northern polar cap acts as an infinite reservoir of water ice, which is also prescribed to have a D/H ratio of 5.
ROSSI ET AL.  Merlivat and Nief (1967). (Right) Change on the zonally averaged seasonal variation of the D/H ratio. The absolute difference is showed between the use of the formula from Lamb et al. (2017) and that of Merlivat and Nief (1967). Contour lines indicate where the D/H ratio reaches 3 and 5 in the Lamb simulation.
In order to reproduce the column of water vapor, in particular the peak of the subliming north pole in northern spring and summer, we had to adjust the albedo of the perennial surface ice, as described in Navarro et al. (2014), to a value of 0.32 yielding a maximum column abundance of 70 pr.-μm similar to observations (Trokhimovskiy et al., 2015).
We start using the climatology scenario and run the model for a couple of Mars years, allowing us to reach a stable state. The resulting model state, after this two-year spin-up phase, supplies the new initial conditions of two following simulations: one based on the "climatology" dust scenario, and one based on the MY34 dust scenario.

Reference Run: The "Climatological" Scenario
Figure 2 (center column) shows the latitudinal and seasonal distribution of the zonally averaged column abundance of HDO vapor along with the corresponding D/H ratio. In essence, HDO vapor and ice reproduce the global cycle of water vapor and ice with a peak of vapor and atmospheric ice abundance in the Aphelion period where the north polar cap, fully exposed to the Sun, rejects massive amounts of water and heavy water in the atmosphere. At the same time, water ice clouds are building up in the northern tropics where air masses ascend and adiabatically cool as a consequence of a global convergence of the atmosphere in the lower atmosphere, causing the recurrent apparition of the "Aphelion" cloud belt.
In summary, the polar night in both hemispheres, which encompasses the L s = 180° to L s = 360° period in the north and the L s = 0° to L s = 180° period in the south, extends equatorward up to 45° near solstices. These regions host the coldest atmospheric temperature overall and as a result the most active water ice condensation mechanism. Likewise, HDO condenses a lot in these regions and is subsequently heavily fractionated, even more because of the exponential increase of the fractionation factor with decreasing temperature (α doubles between 150 and 100 K). This results in a 5:1 gradient between the polar night regions and the regions elsewhere in both hemispheres. Remarkably, the edge of the polar night corresponds to a sharp gradient in D/H, indicating fractionation is very effective there and rather ineffective in regions equatorward.
ROSSI ET AL.
10.1029/2020GL090962 5 of 10 Therefore, the latitudinal gradient suggested by M05 between the cold regions and the warm regions, and whose existence has been confirmed later by several observers, is also present in our simulation.
Clearly, an event such as the Aphelion Cloud belt has no strong D/H signature in a column-integrated perspective since (1) the cloud belt results from moderately cold temperature (180 K) condensation and (2) forms above 10-15 km where atmospheric density has already dropped by a factor 3 to 10, fractionating HDO over a small portion of the column. Altogether, this explains why the global effect of fractionation on Mars is so much dominated by the polar regions.

The Martian Year 34 Run: Impact of the Global Dust Storm
Figure 2 (left column) shows the latitudinal and seasonal evolution of the HDO vapor column and D/H ratio, this time for the MY34 dust scenario. The difference with the climatology scenario is shown in the right panels. One can note that the D/H ratios over the column are quite similar, suggesting that the dust storm does not affect the isotopic ratio of the atmosphere as a whole.
In order to assess the effect that the dust storm may potentially have on the vertical transport of H 2 O and HDO and its evolution with time, we analyze zonally averaged meridional profiles of these quantities (Figure 3). This figure shows a comparison of the computed meridional profiles between the climatology and the MY34 scenarios. The influence of the dust storm is evident on the temperature field with a strong (>20 K) increase of the temperatures in the whole atmosphere due to heating by dust.
Increased dust content has a direct effect on fractionation: higher temperatures restrain condensation, and the subsequent cloud formation in the middle atmosphere exhibited in the Climatological run is pushed to higher altitudes in the MY34 run. Condensation being the main source of isotopic fractionation for water, the deuteropause, that is the level above which D/H declines due to condensation-induced fractionation (Bertaux & Montmessin, 2001), is found around 100 Pa in the climatology run throughout the same L s period as the GDS run. In the MY34 run, clouds are forming above 10 Pa, moving the deuteropause higher in altitude, leaving the D/H ratio mostly unchanged vertically below that level with values of 4.5 reached at 10 Pa, instead of 3.5 in the climatology run. This implies that the vertical barrier of deuterium becomes 30% more porous in a MY34 configuration, likely letting the same excess of deuterium atoms accessing the upper atmosphere where they can escape to space.
The GDS therefore affects the D/H profile in two ways: (1) higher temperature suppresses the vertical confinement imposed by condensation with associated fractionation, (2) intensified atmospheric circulation conveys higher amount of water and its isotope to much higher altitude. This effect on D/H, predicted here by our HDO model, has been already highlighted for water by concomitant observations from the instruments ACS and NOMAD onboard TGO Fedorova et al., 2020;Vandaele et al., 2019) and also reproduced by a GCM (Neary et al., 2020). . The effect of the GDS is quite visible on the water and HDO vapor with a sharp increase in the mixing ratio at most altitudes above 100 Pa in the case of the MY34 scenario. The amount of water is increased in the upper atmosphere by about 10%, but the effect is even stronger for HDO with an increase by 40% which directly scales with the escaping flux of deuterium and can be attributed to the reduced efficiency of the deuteropause. On average over the dust-storm season, the D/H at high altitude is greater by 25% compared to a regular year represented by the climatology run.

Discussion
Our scenario of the MY34 GDS is not fully representative of the observed conditions. In particular, regarding the dust profile and the altitude of the hygropause. In our simulations, the dust remains mostly below 70 km, while observations from MCS show that the dust reached altitudes up to 80 km during the storm . And the height of the hygropause is limited to 60-65 km in our simulation, while observations show that it extended up to 75-80 km during the GDS (Heavens et al., 2019).
These discrepancies may be attributed to missing physical processes in the model. First, we ignore the radiative effect of clouds and the microphysics of cloud formation. These affect the temperatures and enhance the global circulation, as discussed in Madeleine et al. (2011) andNavarro et al. (2014). This could result in an overestimation of the ice content and a reduced transport of water vapor. In addition, without taking microphysics into account, it is not possible to reproduce the important supersaturation observed by ACS during the GDS, and which is thought to help water vapor reaching higher altitudes (Fedorova et al., 2020).
Second, Madeleine et al. (2011) particularly mention the limits of uniformly lifting dust in the case of strong dust storm events, as it does not account for the strong positive feedback between atmospheric dust heating and lifting occurring then. However, other approaches using finer parametrizations of dust lifting have ROSSI ET AL.
10.1029/2020GL090962 7 of 10 similar difficulties reproducing the dust vertical distribution in the GDS conditions (Neary et al., 2020). Finally, the model doesn't account for the convective events suspected to be responsible for the presence of the detached dust layers (Daerden et al., 2015;Heavens et al., 2018;Spiga et al., 2013;Wang et al., 2018). As discussed by Neary et al. (2020), an underestimation of the amount of dust above 40 km leads to insufficient heating, leaving the formation of ice particles possible and therefore restraining the extent of water vapor in the upper atmosphere.
For these reasons, we suspect that our simulation could in fact overestimate the efficiency of the deuteropause. It is thus possible that the increased amount of HDO at high altitudes caused by the GDS as in our simulation, is in fact a lower estimate of the real effect. While we are not able to represent the full extent of the GDS of MY34, our results feature the essential processes occurring in such situation, in particular the reduction of the condensation-induced fractionation, which defines the altitude of the deuteropause, and the increase in atmospheric circulation.

Conclusions
We have presented an updated model of the Martian HDO cycle inherited from Montmessin et al. (2005) and leveraging on the work done since then in particular by Madeleine et al. (2011). This model, which represents the 3D advection of HDO along with the physical processes controlling the fate of HDO in the Martian atmosphere, has been employed to explore for the first time the theoretical impact of a GDSon HDO. Of particular interest is the access of HDO to the upper atmospheric layers where it is subsequently chemically decomposed into escaping deuterium atoms. In doing so, we find that the presence in increased amounts of dust in the atmosphere at times of the GDS restrains condensation and cloud formation, mostly removing the effect of condensation-induced fractionation. The deuteropause, that is otherwise predicted to remain at an approximate altitude of 30 km in a regular year is pushed up to 50 km and appears more porous than in our climatology run representative of average climatic conditions on Mars. The GDS is predicted to enhance the presence of deuterated water at 100 km by 40%, implying that a similar excess of deuterium atoms might escape from Mars in a year of a GDS.
The detailed effects of GDS on the amount of HDO in the upper atmosphere will be further explored including detailed cloud microphysics and cloud radiative feedbacks, which are known to affect the cloud formation (  details of the photochemistry and processes occurring in the thermosphere are missing from the present simulation, but are essential to represent the differential photolysis of H 2 O and HDO and their subsequent escape. These processes will be the object of future improvements of the LMD MGCM, in particular in the context of detailed comparison with available spacecraft observations from ACS onboard TGO. They will help us better understand the global cycle of deuterium from the troposphere to the exosphere, and supply a self-consistent framework to investigate its relation with the escape of water from the planet.

Data Availability Statement
The data used to produce the figures presented in this article can be obtained on the following link: https:// doi.org/10.14768/7f37269b-7e29-42de-a365-139e703a5fc4.