Enhanced water loss from the martian atmosphere during a regional-scale dust storm and implications for long-term water loss

Enhanced water loss from the martian atmosphere during a regional-scale dust storm and implications for long-term water loss Journal Item How to cite: Holmes, J. A.; Lewis, S. R.; Patel, M. R.; Chaffin, M. S.; Cangi, E. M.; Deighan, J.; Schneider, N. M.; Aoki, S.; Fedorova, A. A.; Kass, D. M. and Vandaele, A. C. (2021). Enhanced water loss from the martian atmosphere during a regional-scale dust storm and implications for long-term water loss. Earth and Planetary Science Letters, 571, article no. 117109.


Introduction
The past climate of Mars is now generally considered to have been much wetter than the present (Carr and Head, 2010;Ehlmann and Edwards, 2014, and references therein). The precise value of increased water abundance in the early martian atmosphere is still however largely unconstrained, in part due to the lack of observations quantifying the rates of water loss from the atmosphere over time. Observations from the SPICAM (Spectroscopy for the Investigation of the Characteristics of the Atmosphere of Mars) instrument on Mars Express (Chaffin et al., 2014), the Hub-* Corresponding author.
ble Space Telescope (Clarke et al., 2014) and more recently by the MAVEN/IUVS (Mars Atmosphere and Volatile EvolutioN Imaging Ultraviolet Spectrograph instrument (Chaffin et al., 2021) have advanced our knowledge on the unexpected variability of martian hydrogen escape; but upper atmosphere observations alone are unable to diagnose effectively the global scale evolution of the complex water cycle and associated transport processes. To understand how much water has been lost to space over time and how much it varies seasonally requires an understanding of the processes by which hydrogen escape is coupled to the lower atmosphere water cycle. Recent 1-D modelling efforts indicate that atmospheric escape of hydrogen is increased with the introduction of high-altitude water (Chaffin et al., 2017;Krasnopolsky, 2019). Based on this work, the paradigm of martian atmospheric water loss is shifting focus to high-altitude water being the domi- nant pathway over diffusive transport of molecular hydrogen from the lower atmosphere (Anderson and Hord, 1971;Parkinson and Hunten, 1972;McElroy and Donahue, 1972;Krasnopolsky, 2002). Efforts to calculate integrated water loss over time have been conducted most recently by Jakosky et al. (2018) with an extrapolation back in time of recent MAVEN/IUVS observations. The calculated value however is sensitive to multiple currently unknown factors, including the extent to which the seasonal water and dust cycles can affect the composition of the upper atmosphere and therefore the escape rate of hydrogen.
Advances in modelling techniques and the number of observations now available from several spacecraft currently in orbit around Mars allow for a more complete investigation of the processes driving hydrogen escape. The ability to combine observations with a global circulation model (GCM) through data assimilation techniques (Lewis et al., 2007(Lewis et al., , 2016Holmes et al., 2018Holmes et al., , 2019Holmes et al., , 2020Streeter et al., 2020) is a critical step in providing the most realistic global atmospheric state. While providing a detailed exploration of the parameter space, Chaffin et al. (2017) note that future modelling efforts should look into the response of the atmosphere to realistic variability in upper atmospheric water, of which data assimilation provides the most effective tool for investigation. Previous modelling studies have investigated vertical transport of water during the Mars Year (MY) 28 and MY 34 global dust storms (Shaposhnikov et al., 2019;Neary et al., 2020), where MY follows the designation of Clancy et al. (2000), but neither were constrained directly by observations. In this investigation we focus on the MY 34 regional dust storm that occurred from L S = 320.6 • -336.5 • , hereafter called the C storm based on classification by Kass et al. (2016).
Vertical profiles of water vapour have been retrieved from the Nadir and Occultation for MArs Discovery (NOMAD) and Atmospheric Chemistry Suite (ACS) instruments on the ExoMars Trace Gas Orbiter (TGO) during the initiation and decay of the MY 34 C storm (Aoki et al., 2019;Fedorova et al., 2020). Observations of water could not be made at the time of peak activity observed by MAVEN/IUVS during the time period that covers the MY 34 C storm, a key gap in understanding of lower atmosphere water/hydrogen activity. During the time period of the MY 34 C storm unobserved by ExoMars TGO, however, we can still constrain model simulations of the water/hydrogen activity using observations of the temperature and dust distribution from the Mars Climate Sounder (MCS) aboard the Mars Reconnaissance Orbiter (MRO) spacecraft (Kleinböhl et al., 2017) that cover the entire MY 34 C storm time period, a powerful advantage of multi-spacecraft data assimilation. This multi-spacecraft assimilation can then be analysed and compared to upper atmosphere Lyman alpha brightness measurements (a proxy for hydrogen escape) from MAVEN/IUVS that display an increase in the Lyman alpha brightness (Chaffin et al., 2021) during the MY 34 C storm. If the hydrogen escape flux calculated from a simulation of the atmosphere constrained by observations matches the observed Lyman alpha brightness trends, we will have provided a more complete picture of the processes that lead to variations in hydrogen escape.
In this study we investigate the global distribution of lower atmosphere water and hydrogen and coupling to the upper atmosphere using data assimilation covering the time period leading up to and during the MY 34 C storm. The data combined with the Open University modelling group Mars GCM includes observations of water vapour from NOMAD and ACS (that constrain the initial global distribution of water), temperature profiles from ACS and MCS and dust column from MCS. The results are also compared to an assimilation during the MY 30 C storm to identify interannual variations in hydrogen escape.

Methods
This section details the GCM and observational data used for this study, followed by a description of the different simulations conducted to investigate water and hydrogen abundance and spatial variations during the MY 34 C storm.

Model and data assimilation
For the global simulations, we use the Open University (OU) modelling group Mars GCM, hereafter MGCM, which has been developed in a collaboration between the Laboratoire de Météorologie Dynamique (LMD), the OU, the University of Oxford and the Instituto de Astrofísica de Andalućia. This model combines physical parameterisations Forget et al. (1999) of multiple processes (such as radiative transfer, surface processes and subgrid-scale dynamics) and a photochemical module (Lefèvre et al., 2004(Lefèvre et al., , 2008 shared with the LMD Mars GCM, that are coupled to a spectral dynamical core and semi-Lagrangian advection scheme (Newman et al., 2002) to transport tracers. Tracers such as water vapour and hydrogen are transported by the semi-Lagrangian advection scheme with mass conservation. The advection scheme uses wind fields updated by the dynamical core to determine the transport of each chemical species at each model grid point every 15 minutes.
The MGCM includes several sub-models that encompass modelling of the planetary boundary layer and water and dust cycles. A thermal plume model is used to better represent turbulent structures in the planetary boundary layer (Colaïtis et al., 2013). Regarding the martian water cycle, the most recent cloud microphysics package is included (Navarro et al., 2014) which also accounts for the effects of radiatively active water ice clouds. A 'semi-interactive' two-moment scheme is used to freely transport dust in the model (Madeleine et al., 2011), although the dust column optical depth at each grid point is scaled to match the observed dust distribution maps created by Montabone et al. (2015), updated to include MY 34 in Montabone et al. (2020), using an interpolation of numerous sets of observations from orbiters and landers. The model is truncated at wavenumber 31 resulting in a 5 • longitude-latitude grid for physical variables, with 70 vertical levels extending to an altitude of ∼100 km (chosen to provide a set of pressure levels that correspond roughly to every 2 km for the majority of the middle atmosphere).
To perform the multi-spacecraft data assimilation, we use the Analysis Correction (AC) scheme (Lorenc et al., 1991) with necessary parameters adapted to martian conditions (Lewis and Read, 1995). The AC scheme is a form of successive corrections in which analysis steps are interleaved with each model dynamical time step and increments (observation -model) are first calculated and incorporated in a vertical assimilation step and then spread horizontally to other model grid points (Lewis et al., 2007). The AC scheme has previously been used for multiple studies covering thermal tides (Lewis and Barker, 2005), the dust cycle (Montabone et al., 2005(Montabone et al., , 2006, surface warming during the MY 34 global dust storm (Streeter et al., 2020) and several chemical cycles including water (Steele et al., 2014b,a), ozone (Holmes et al., 2018) and carbon monoxide (Holmes et al., 2019). The AC scheme has also been utilised to create the Open access to Mars Assimilated Remote Soundings (OpenMARS) dataset , a publicly available global record of martian weather from 1999 to 2015.

Temperature and water vapour profiles
The temperature profiles used in this study are from the MCS and ACS instruments on the MRO and ExoMars TGO spacecraft respectively, with version 5.3.2 MCS temperature profiles used  and the ACS temperature profiles described in Fedorova et al. (2020). The differing orbits of the MRO and Ex-oMars TGO spacecraft mean that while MCS temperature profiles are available at 5 km vertical resolution in 12 strips of data per sol and at local times of 3 a.m. and 3 p.m. away from the poles, the ACS temperature profiles span a much wider range of local times related to the terminator, with good agreement between both datasets when coincident measurements occur (Fedorova et al., 2020). Therefore, combining both datasets through data assimilation provides more complete constraints on the diurnal temperature variations than possible when using either dataset alone. For the assimilation of temperature profiles from both the MCS and ACS instruments, we use the same method that has been used previously for studies that included Thermal Emission Spectrometer data (Steele et al., 2014b) or MCS data (Steele et al., 2014a;Holmes et al., 2019;Streeter et al., 2020) and compared the model to retrievals in the form of layer thicknesses (Lewis et al., 2007).
Water vapour profiles are assimilated from both NOMAD and ACS instruments on the ExoMars TGO spacecraft, described in Aoki et al. (2019) and Fedorova et al. (2020) respectively, that are retrieved from solar occultation measurements and therefore restricted to the terminator. The water vapour profiles have a vertical resolution of 1 to 3 km and can span altitudes from 5 to 100 km depending on atmospheric conditions at the time of observation. For the assimilation of water vapour profiles, a similar method to the assimilation of MCS ice opacity detailed in Steele et al. (2014a) is used. Water vapour number density in each profile are converted from observation levels to model levels (assuming the number density varies linearly with the natural log of pressure). Background vertical correlations are approximated using the same Gaussian function and parameter values in Steele et al. (2014a), chosen to allow spreading of data to at most two model levels outside the bounds of the vertical profile coverage, since water vapour can rapidly vary based on local saturation conditions. Regarding water vapour, there are three sources of information; NOMAD data, ACS data and the model forecast, that are analysed by the assimilation scheme to provide the best fit atmospheric state of water vapour. This provides a more robust investigation than pure model studies that simply compare to observed water, by ensuring that multiple parameters (i.e. water vapour, dust opacity, temperature) are all realistically constrained and physically consistent at the same time.

Simulations
The primary simulation conducted for this study is an assimilation of temperature and water vapour profiles covering the MY 34 C storm. Initial conditions were obtained from a 10 year model spin-up of the water cycle, with the dust mass and mixing number ratios in each layer scaled to the MY 34 dust scenario  for the final year. The simulation is initiated at L S = 270 • MY 34 during the decay stage of the MY 34 global dust storm  to provide an assimilation spin-up of 82 sols before the C storm initiates, with results displayed covering the time period of L S = 315 • -345 • MY 34 to focus on the impact of the C storm on the lower atmosphere water and hydrogen distribution.
To investigate whether the exceptional dust circumstances in MY 34 were representative of hydrogen loss during a standard Mars year, for sections 4 and 5 we perform a simulation that assimilates MCS temperature profiles from MY 30 to provide a comparison with another MY. This particular Mars year was chosen as it had good coverage of temperature observations and a similar initiation time of the C storm, albeit with much less dust present (Montabone et al., 2015;Kass et al., 2016), allowing us to investigate the extreme cases expected during the C storm period (the MY 34 C storm is amongst the strongest of this type of event whereas the MY 30 C storm is close to the minimum). This sim-ulation has identical initial conditions to the primary simulation (the interannual variability in water at this time of year is small (Smith, 2004)) but from L S = 310 • the simulation switches to assimilate MY 30 temperature profiles and the dust column is scaled to match the MY 30 dust scenario (Montabone et al., 2015) with no water vapour profiles assimilated.

Water distribution during the C storm
The water vapour vertical distribution from the C storm time period in MY 34 is shown in Fig. 1. ACS observations assimilated before the ExoMars TGO observational gap and NOMAD observations after the gap are predominantly located at latitudes between 30 • S-45 • S and 60 • S-90 • S respectively, and the water vapour abundance in the assimilation are in good agreement with those presented in Chaffin et al. (2021). The distribution shown in the simulated global water vapour distribution and the ExoMars TGO observations are not directly comparable because the NO-MAD and ACS retrievals are plotted for specific observation occurrences whereas the assimilated water vapour distribution shown is a time-zonal-mean. At the initiation of the C storm (L S = 320.6 • ) the hygropause (defined as where the water vapour volume mixing ratio drops below 70 ppmv) is at an altitude of around 50 km for the majority of latitudes, with a sharp drop-off northward of 60 • N as a result of northern winter conditions. During the peak of the MY 34 C storm (Fig. 1b), the hygropause altitude at southern latitudes increases to 80 km and is comparable to the hygropause altitude reported during the MY 28 global dust storm through indirect measurements (Heavens et al., 2018) and more directly through investigation of water vapour profiles retrieved by the SPICAM instrument (Fedorova et al., 2018).
The increase in high-altitude water abundance coincides broadly with the ExoMars TGO gap in observations (L S = 326.1 • -333.5 • ) and covers all latitudes albeit with varying levels in peak abundance. Even if NOMAD/ACS were able to continue observing, orbital constraints mean that the instruments can only capture a glimpse of the complex spatial structures apparent throughout the global atmosphere at this time due to the latitudinal sampling of the solar occultations. While it is unfortunate that there are no water vapour retrievals to constrain the lower atmosphere during the peak of the MY 34 C storm, the initial conditions for this time period come from an assimilation of water vapour profiles and during the peak the dynamics of the atmosphere remain constrained through continued assimilation of MCS temperature profiles and dust column, a huge benefit of multi-spacecraft assimilation over a model-only GCM simulation.
The results found in Fig. 1 can largely be explained through the mean meridional circulation patterns that occur during the MY 34 C storm. The time periods are so as to reasonably cover the initiation, peak and decay of the C storm and the assimilation before, during and after the gap in water vapour profile data from NO-MAD and ACS. The majority of water vapour in the atmosphere is located over the southern pole after perihelion, and the subsolar point is moving away from the south pole and towards the northern hemisphere. At the onset of the C storm (Fig. 1a), the anticlockwise cell situated over southern high latitudes is relatively weak, with a more dominant clockwise cell covering 30 • S-60 • N, and no water vapour in excess of 20 ppmv above 80 km altitude. This is in contrast to the expanded water vapour distribution during the peak of the C storm shown in Fig. 1b, with high-altitude water vapour exceeding 70 ppmv across all latitudes above 70 km. During this time period (i.e. in the ExoMars TGO observation gap and during peak brightness in MAVEN/IUVS observations), the subsolar point has shifted northward and the strength of both the anti-clockwise and clockwise cell have been increased due to the heating of the increased atmospheric dust. The meridional cells The strengthened cells during the C storm provide an influx of water vapour to higher altitudes and also to more northerly latitudes throughout the lower atmosphere through the dusty deep convection mechanism (Heavens et al., 2018(Heavens et al., , 2019. The upwelling branch of circulation that separates the anti-clockwise and clockwise cells at this time is located at around 35 • S, and so highaltitude water abundance is higher at nearby latitudes. Above around 65 km, the increased latitudinal extent of the clockwise cell results in water vapour in excess of 80 ppmv reaching the majority of the northern hemisphere. After the decay of the C storm (Fig. 1c), the water distribution reverts to a similar pattern to before the C storm (Fig. 1a), although the subsolar point and hence upwelling branch of circulation are shifting further northward before transitioning from a cross-hemispheric clockwise cell towards split cells in each hemisphere by northern spring equinox. Chaffin et al. (2017) indicated high-altitude water vapour as a possible explanation for seasonally variable rates in the loss of hydrogen from the upper atmosphere of Mars. The water vapour that reaches high altitudes produces hydrogen through photodissociation, providing a direct source for increased hydrogen and enhanced escape. Our simulations include the LMD photochemical module (Lefèvre et al., 2004(Lefèvre et al., , 2008 and hence hydrogen is also simulated, with its spatial distribution constrained by the photochemical module and additionally by the NOMAD/ACS water vapour profiles (when available) and the MCS temperature profiles throughout the assimilation. The following results sections investigate the high-altitude hydrogen distribution, followed by a calculation of the hydrogen escape through coupling our lower atmosphere assimilation with the 1-D modelling work of Chaffin et al. (2017) that has recently been updated and adapted by Cangi et al. (2020) to investigate D/H fractionation and water loss. These results can also be compared to 1-D photochemical modelling by Krasnopolsky (2019). Through this lower-upper atmosphere coupling, we can estimate the hydrogen loss rate expected during the C storm time period across the entire globe, rather than in a region or at a specific point such as is the case with observations alone. Modelling without constraints imposed by available retrievals can deviate from the true state of the atmosphere, but using data assimilation we can provide robust estimates of the hydrogen escape rate across the globe. Fig. 2b,d,f shows the zonally averaged water vapour and hydrogen abundance at 80 km and the column dust optical depth during the C storm time period of MY 34 respectively. Before the C storm time period, water vapour levels are low after peak southern summer sublimation and hydrogen levels are maximum around the southern polar region associated with Mars moving away from the perihelion season. Krasnopolsky (2019) indicate hydrogen abundance of around 6 ppmv at 80 km during typical near-perihelion conditions using a 1-D model, while the simulations here have a slightly higher hydrogen abundance that can be accounted for by an expected shorter timescale for transport from the lower to upper atmosphere in a fully 3-D simulation. During peak MY 34 C storm activity, water vapour and hydrogen levels reach around 80 ppmv and 45 ppmv respectively, with hydrogen abundance increased across all latitudes compared to the low abundance at L S = 322.5 • during the C storm initiation (Fig. 2f).

Hydrogen distribution during the C storm
To determine whether the high-altitude hydrogen abundance during the MY 34 C storm is representative of a typical Mars year (the C storm in MY 34 was a particularly strong dust event) and hence can be used as a basis for the extrapolation of hydrogen escape during C storm events under the present obliquity configuration, the results need to be compared against another Mars year to identify if the C storm response is the same. Fig. 2a,c,e shows the zonally averaged water vapour and hydrogen abundance at 80 km and the column dust optical depth during the C storm time period of MY 30 respectively. As previously stated, the MY 30 assimilation has initial conditions at L S = 310 • identical to the MY 34 simulation and the only difference being that MCS temperature profiles and dust column observations are assimilated from MY 30.
Before the initiation of the C storm, water vapour and hydrogen abundance and distribution in the MY 30 simulation (Fig. 2c) is largely similar to the MY 34 simulation displayed in Fig. 2d, indicating any increase in hydrogen found associated with the C storm can be decoupled from the seasonal trends on Mars that occur each Mars year. During the time period of the weaker C storm in MY 30, water vapour and hydrogen abundance peak at values of around 20 ppmv and 15 ppmv meaning that water vapour and hydrogen abundance during the MY 34 C storm were at least a factor of 3 times higher, and increased further during the decay phase of the C storm.
These results indicate that weaker C storm events do not have a similar effect on upper atmosphere hydrogen abundance as stronger C storm events such as the MY 34 C storm. Therefore differences in the strength of C storm events, that lead to differential heating and expansion or contraction of the water vapour distribution, are likely to lead to annual variations in the hydrogen escape from Mars during the C storm time period. This hypothesis will be investigated in the next section.

Hydrogen escape during the C storm
Number densities at 80 km for each simulated chemical species in the lower atmosphere global assimilation are provided as initial lower boundary conditions for the 1-D photochemical model that spans the altitude range from 80 km to 200 km, with the top altitude representing the exobase at which the hydrogen escape flux can be calculated. A fixed temperature profile (Fig. 3a) is used that is consistent with upper atmosphere simulations in the Mars Climate Database version 5.3 (Millour et al., 2018) during  the investigated time period. Chemical number densities at 80 km are updated daily from the global assimilation to mimic the lower atmosphere variations simulated for MY 34 and MY 30, with the initial atmosphere above 80 km for the 1-D simulations calculated using the GCM boundary conditions and an example shown in Fig. 3b. Fig. 4 displays the MAVEN retrievals of Lyman alpha brightness during the MY 34 C storm and the hydrogen escape flux calculated during the C storm time period for both the MY 34 and MY 30 assimilation, with the hydrogen escape flux calculated across the globe on a 30 • by 25 • longitude-latitude grid. The sampling was chosen to attempt to represent global variations in hydrogen escape that are likely to occur but are not explicitly modelled as the 1-D model implementation results in a lack of horizontal transport above 80 km. While a quantitative assessment is not possible since the mapping of Lyman alpha brightness to hydrogen escape flux is not 1-to-1, analysis of the general trend in both sets of data and timing can be made. The time required for transport from the lower to upper atmosphere is expected to be shorter in a fully 3-D simulation rather than the current coupled 3-D lower atmosphere and 1-D upper atmosphere set-up and so the calculated rates can be interpreted as a lower bound on the hydrogen escape rate. For the study of hydrogen escape over longer timescales (e.g. annual and multi-annual) it would be beneficial to run a 3-D simulation throughout the entire atmosphere, but this is beyond the scope of the current study. In the initiation phase of the MY 34 C storm, the hydrogen escape flux across the globe range from 2.8 − 4.2 × 10 8 cm −2 s −1 depending on the exact spatial location, consistent with present day escape rates at nominal conditions (Chaufray et al., 2008;Chaffin et al., 2014). Once the hydrogen variations at 80 km during the MY 34 C storm have propagated to higher altitudes, the hydrogen escape flux increases to a peak of around 0.9 − 1.5 × 10 9 cm −2 s −1 , a value that falls within the derived 1 − 5 × 10 9 cm −2 s −1 during a previous global dust storm from SPICAM measurements (Chaffin et al., 2014;Heavens et al., 2018). This means that strong regional dust storms can enhance water loss rates on Mars to those levels observed during global-scale dust storms. These values are also in good agreement with hydrogen escape values calculated for typical conditions near perihelion by Krasnopolsky (2019). Under the assumption that the hydrogen escape flux is at a nominal value of 1.5 − 2.0 × 10 8 cm −2 s −1 (Bhattacharyya et al., 2015;Krasnopolsky, 2019) for a large majority of a Mars Year (and increased at peak southern summer conditions), the conditions during the MY 34 C storm time period would contribute a loss of 15% of the total annual water loss during only 5% of the year.
The time period over which the hydrogen escape flux increases and general trend is consistent with the MAVEN/IUVS Lyman alpha observations in Fig. 4a, with a similar week lag from the lower atmosphere expansion of water vapour as a result of the C storm ( Fig. 2b and d). The decrease in hydrogen escape flux from around L S = 332 • onwards is also consistent with the reduction in Lyman alpha brightness observed by MAVEN/IUVS.
While the hydrogen escape flux calculated during the MY 34 C storm shows a distinct increase, the MY 30 assimilation is in stark contrast with similar hydrogen escape rates of 4 − 6 × 10 8 cm −2 s −1 across the globe at the initiation of the MY 30 C storm followed by a steady decline in hydrogen escape as the MY 30 C storm progresses and decays. These results indicate that C storm strength and associated variations in hydrogen escape need to be taken into account when extrapolating loss of water back in time.

Conclusions
A multi-spacecraft assimilation of ExoMars TGO and Mars Reconnaissance Orbiter retrievals into the OU modelling group GCM has been performed to investigate the lower atmosphere distribution of water vapour and hydrogen during the MY 34 C storm.
These 4-D simulations have then been coupled to an upper atmosphere 1-D photochemical model to calculate global hydrogen escape rates associated with the MY 34 C storm, a particularly intense regional dust event.
The inclusion of MCS and ACS temperature profiles and for the first time water vapour profiles from NOMAD and ACS in the assimilation process provides the best possible representation of the actual temperature structure and circulation of the lower atmosphere during the MY 34 C storm and associated evolution of water vapour (and hydrogen). During the peak of the MY 34 C storm the dynamical structure was constrained by continuing MCS observations.
Water vapour retrievals from ExoMars TGO were unavailable during the peak hydrogen activity and elevation of the hygropause altitude associated with the MY 34 C storm that caused an expansion of the lower atmosphere and increased abundance of water vapour above 80 km during L S = 326.1 • -333.5 • . The increase in water vapour during this time period led to an associated increase in hydrogen abundance, which after propagating to the exobase led to an increase in the globally averaged hydrogen escape rate from 3.5 × 10 8 cm −2 s −1 before the MY 34 C storm to around 1.0 × 10 9 cm −2 s −1 at its peak, which means strong regional dust storms can enhance water loss rates on Mars to those levels observed during global-scale dust storms. A MY 30 C storm assimilation showed no associated increase in the hydrogen escape flux, indicating the influence of the MY 34 C storm was particularly intense, and that the hydrogen escape rates for any given C storm can be highly variable. These results indicate that interannual variations in the C storm strength need to be taken into account to provide a robust estimate of the integrated loss of water that has occurred on Mars.

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
The simulation data used in this study are publicly available via the Open Research Data Online (ORDO) data repository at the following DOI: https://doi .org /10 .21954 /ou .rd .13622699.
ST/S00145X/1. The NOMAD experiment is led by the Royal Belgian Institute for Space Aeronomy (IASB-BIRA), assisted by Co-PI teams from Spain (IAA-CSIC), Italy (INAF-IAPS), and the United Kingdom (Open University). This project acknowledges funding by the Belgian Science Policy Office (BELSPO), with the financial and contractual coordination by the ESA Prodex Office (PEA 4000103401, 4000121493), by Spanish Ministry of Science and Innovation (MCIU) and by European funds under grants PGC2018-101836-B-I00 and ESP2017-87143-R (MINECO/FEDER), and Italian Space Agency through grant 2018-2-HH.0. This work was supported by the Belgian Fonds de la Recherche Scientifique -FNRS under grant number 30442502 (ET_HOME). The IAA/CSIC team acknowledges financial support from the State Agency for Research of the Spanish MCIU through the 'Center of Excellence Severo Ochoa' award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). US investigators were supported by the National Aeronautics and Space Administration. The MAVEN mission is supported by NASA through the Mars Exploration Program in association with the University of Colorado and NASA's Goddard Space Flight Center. Mike Chaffin was supported by NASA through the MAVEN Project. AF was supported by the Ministry of Science and Higher Education of the Russian Federation. Work at the Jet Propulsion Laboratory, California Institute of Technology, was performed under a contract with the National Aeronautics and Space Administration. US Government sponsorship is acknowledged.