Surprising Decrease in the Martian He Bulge During PEDE-2018 and Changes in Upper Atmospheric Circulation

Using the Neutral Gas and Ion Mass Spectrometer (NGIMS) on the Mars Atmosphere Volatile and Evolution spacecraft (MAVEN), we analyzed data from Mars Year (MY) 32, 34, and 35 to examine the He bulge during the northern winter solstice (Ls ∼ 180–240), specifically focusing on the effects from the planet encircling dust event (PEDE-2018). He collects on the dawn/nightside winter polar hemisphere of Mars. The seasonal migration of the He bulge has been observed and modeled (M. Elrod et al., 2017, https://doi. org/10.1002/2016JA023482; Gupta et al., 2021, https://doi.org/10.1029/

The impact of the global dust storm has been shown to have an increase in water ice at higher altitudes than previously measured (Stone et al., 2020) along with compositional changes to the upper atmosphere (M. K. Elrod et al., 2020;Farahat et al., 2021). The compositional changes in the upper atmosphere indicated that the lighter species, specifically the O and N 2 were reduced during the PEDE compared to the main species (CO 2 and Ar). Examination of the nightside winter and dusk data taken by MAVEN/NGIMS during the PEDE revealed that the He density was lower than what had been previously observed for the He bulge during winter solstice observations (see M. Elrod et al., 2017;Gupta et al., 2021 and Figure 1).
Mars has a He bulge (M. Elrod et al., 2017) that peaks in density where the temperatures of the upper atmosphere are relatively cool. MAVEN/NGIMS has obtained in situ observations of He for 5 MYs, MY 32-36. He density peaks are highest around equinox (Ls ∼ 0, or 180) near mid to high latitudes (lat ∼45°-70°) and winter solstice (Ls ∼ 90 or 270) near corresponding high latitudes (lat ∼55-80) on the night or dawn side local solar time (LST) of ∼22-5 hr. Due to the eccentricity of the Martian orbit, there is some observed asymmetry between the northern and southern He bulge (Gupta et al., 2021).
In order for NGIMS to be able to observe the He bulge, because it is an in situ instrument, MAVEN needs to be moving through the He bulge. Therefore, in order to follow the northern winter He bulge, we need to constrain our observations to those when MAVEN is in the northern hemisphere from Ls 180-280. Since the He bulge is constrained from LST 22-6 hr (midnight to dawn) if NGIMS is only observing on the dayside during this period, then a Helium bulge is not likely to be observed. MAVEN was in the northern hemisphere from Ls 180-280 in MY 32, 34, and 35. In all 3 years, NGIMS was observed from midnight to dawn within this Ls range. In MY 33 and 36, MAVEN was in the southern hemisphere from Ls 180-280, making it impossible for NGIMS to detect the winter He bulge at this time. For this reason, we limited our data range to MY 32, 34, and 35 from Ls 180-280 and 0°N-70°N latitude. The latitude versus Ls position with color representing LST for MY 32, 34, and 35 in Figure 1d. The ovals are added to highlight the northern winter season that is the focus of study for this paper.
Helium data from Martian years during which dust storms are only small regional storms (constrained to just the southern hemisphere or are short lived) were compared to the MY 34 global dust storm event. The He densities were found to be substantially lower during the global dust event where the bulge should have been observed (Figures 4 and 5). Mars Global Ionosphere Thermosphere Model (M-GITM)  numerical model simulations demonstrated that under normal circumstances during northern winter (Ls ∼ 180-240) a He bulge should form in the upper latitudes (lat ∼50°-70°) and night to dawn sectors (LST ∼ 20-6 hr) (M. Elrod et al., 2017). MAVEN NGIMS began observations during the PEDE at mid-latitudes and on the dawn-side at mid to low latitudes (∼20°-35°) moving to dusk/night (LST ∼ 15-3 hr; lat 35°-65°).

Data Set Collection
MAVEN has an inclined orbit with a periapsis of ∼150 km from 2014 to 2020; this periapsis was then raised to ∼180-200 km after August 2020 (mid MY 35). NGIMS is an in situ mass spectrometer that takes measurements of the upper atmosphere during periapsis from below ∼500 km (Mahaffy et al., 2015;Stone S. W. et al., 2022). These orbits and measurements allow NGIMS to vary observation between all local times and latitudes throughout the Martian year. Due to the limit of observations, and the fact that dust storm season typically occurs around Ls ∼ 180-290, we needed to select data from the northern hemisphere, latitudes >45° and preferably also after 18 hr and before 6 hr LST. MAVEN/NGIMS had this orientation during MY 32, 34, and 35. In MY 33, the spacecraft was in the southern hemisphere for Ls ∼ 200-360 and unable to observe the He bulge ( Figure 1d).
While the higher altitudes have not been ideal for study of the Martian upper atmosphere, NGIMS has continued to conduct regular science. The higher altitude measurements had an impact on the heavier species (CO 2 , Ar, and O) measurements, but since He is a light element and has a large-scale height (nearly vertical density variation), the higher altitude measurements had little impact. Due to its large scale height, the He densities measured between 175 and 220 km are very similar in the thermosphere. For this reason, we used data from 170-225 km from MY 32, 34, and 35 for coverage and for better statistics. Since the PEDE began at Ls ∼ 180 and extended through Ls ∼ 240, and since the He bulge should be forming in the northern polar region as this is near northern winter, we restricted our data to this range and latitude 0°N-70°N.
The NGIMS density data is obtained from the level 2 closed source neutral files (version 8 revision 6) available on the Planetary Data System (PDS). The NGIMS temperature data are obtained from the scale height-derived Ar density profiles from the level 3 product (version 6 revision 2) also available on the PDS. NGIMS derives the Ar scale heights using the inbound Ar densities from 165 to 195 km this range was modified after the periapsis raised in August of 2020 to 185-220 km and an additional flag was included to indicate when there were insufficient densities, or the scale heights were computed with very high error. We excluded these temperature data from our study. This only impacted the MY 35 temperature data in Figure 1c. We only included the inbound He densities as per recommendation by the NGIMS team due to instrument effects on the outbound densities.

Modeling
A three-dimensional (3D) General Circulation Model is employed to address the structure and dynamics of the upper atmosphere of Mars, with the immediate objective of simulating the helium distribution as it varies with season and dust conditions. Specifically, the existing Mars Global Ionosphere-Thermosphere Model (M-GITM) is utilized to solve the Navier-Stokes equations for energy, momentum, and neutral composition on a regular spherical grid (5 × 5 degree horizontal resolution and 2.5 km vertical resolution over 120 levels). This Mars GITM model now combines the terrestrial GITM framework (e.g., Ridley et al., 2006) with unique Mars fundamental physical parameters, specific ion-neutral chemistry, and essential CO 2 radiative processes in order to simulate the observed global features of the thermal, compositional, and dynamical structure of the Mars atmosphere up to ∼300 km (e.g., Bougher et al., 2015).
The M-GITM code now simulates the conditions of the Martian atmosphere near the traditional exobase all the way to the surface . Thermal and chemical drivers are carefully incorporated into M-GITM. For instance, a state-of-the-art correlated k-radiation code was adapted from the legacy NASA Ames Mars General Circulation Model (MGCM) (Haberle et al., 2013) for incorporation into the M-GITM lower atmosphere below ∼80 km. This provides solar heating (long and short wavelength), aerosol heating (seasonally variable), and CO 2 15-micron cooling in the local thermal equilibrium region of the Mars atmosphere. Conversely, a modern CO 2 non-local thermal equilibrium 15-μm cooling scheme was implemented (e.g., González-Galindo et al., 2013) within M-GITM for the upper atmosphere (∼80-250 km). External solar fluxes are taken from the MAVEN Extreme ultra-violet monitor (EUVM) via the corresponding Flare Irradiance Spectrum Model (Mars) (FISM-M). This empirical model provides 1-195-nm solar fluxes on a daily cadence throughout the MAVEN mission (Theimann et al., 2017). The resulting M-GITM thermospheric heating, dissociation and ionization rates are simulated at each model time step (e.g., Bougher et al., 2015).
A few other parameters, formulations, and needed boundary conditions are important to describe as well. Eddy diffusion in the model is prescribed as a globally uniform pressure-based vertical profile that yields a maximum eddy coefficient (∼1,000 m 2 /s) in the thermosphere and a minimum value (∼500 m 2 /s) in the middle atmosphere and below . Important boundary conditions include those at the top of the model (near the exobase) and at the surface . Topside conditions allow vertical density gradients to be continuous, consistent with molecular diffusion; a non-zero flux boundary condition is thereby maintained. In addition, topside temperatures are specified to be isothermal, in accord with the exobase approximation. At the surface, densities are specified from empirical estimations, along with observed diurnal and seasonal variations for temperatures. In addition, surface thermal inertia and albedo are prescribed from 2-D maps utilized in the NASA MGCM legacy model (e.g., Haberle et al., 2013).
It is clear from both observations and modeling studies that the Mars upper atmosphere is greatly impacted by waves of various spatial and temporal scales propagating upwards from the lower atmosphere (e.g., England et al., 2017;Forbes et al., 2020;Leelavathi et al., 2020;Medvedev & Yiğit, 2012;K. J. Roeten, Bougher, Yiğit, et al., 2022;Yiğit E., 2023;Yiğit et al., 2021). Investigating the interaction of these waves with dust storms is essential for expanding our understanding of the vertical coupling (i.e., energy, momentum, and mass exchange between atmospheric layers) of the lower and upper atmospheres of Mars (Leelavathi et al., 2020;Yiğit E., 2023). For this paper, we address gravity wave (GW) impacts within the M-GITM framework during the PEDE-2018 dust storm.
For this PEDE study, a simulation without GW effects is conducted first. The setup for this simulation is as follows. The period for capturing the PEDE storm is targeted from June through July 2018, which includes the ramp-up, peak, and initial decay phase of the storm (see Kass et al., 2020). Model inputs for time varying solar EUV-UV fluxes are taken from the FISM-M empirical model (Theimann et al., 2017), which includes Mars seasonal variations. Model inputs for time varying dust optical depths and dust mixing ratio distributions are taken from MRO MCS data sets (Kass et al., 2020). Specifically, both zonal and sol averaged values for these parameters are extracted as a function of standard MCS pressure levels (103) and new latitude elements (36). This 2-D grid is chosen to match the latitude grid of M-GITM, while assuming that zonal averaging is a reasonable first step for our simulations. Interpolation in time between even sol intervals is conducted to provide a smooth transition between sol points for the dust distribution. Aerosol heating rates are subsequently computed based upon these specified dust distributions. In addition, a simulation with GW effects is also conducted for this PEDE study. Recently, an existing non-orographic whole atmosphere GW momentum and energy deposition scheme (e.g., Medvedev & Yiğit, E., 2012;Yiğit et al., 2021) was incorporated in the M-GITM code in order to examine the impacts of GWs upon the thermospheric wind, density, and temperature structure (see K. J. Roeten, Bougher, Yiğit, et al., 2022). Adding schemes to parameterize unresolved physical phenomena, such as GWs, is a common approach for relatively low-resolution global models such as M-GITM. Under typical (non-dust storm) conditions, it was previously reported that with the addition of this GW scheme into M-GITM, simulated global winds can be reduced by a factor of two, and temperatures are generally cooled at all latitudes at thermosphere heights, but especially at high latitudes (K. J. K. J. Roeten, Bougher, Yiğit, et al., 2022), compared to simulations without the scheme.
Since the addition of GW effects into M-GITM was found to produce large impacts on the modeled general circulation, an M-GITM simulation during the PEDE which utilizes the same GW scheme was also examined here to see if there are any substantial differences in the He distribution. Previously (e.g., K. J. Roeten, Bougher, Yiğit, et al., 2022), it was found that the key adjustable GW parameters to specify for driving the GW code were the: (a) horizontal wavelength, (b) source momentum flux, and (c) launching level for these GWs near the top of the planetary boundary layer (PBL). The last two of these parameters may change over the course of the evolution of a global dust storm (i.e., yielding reduced GW forcing in the lower atmosphere), and when coupled with more favorable vertical propagation conditions (due to aerosol heating) lead to larger GW activity observed at thermospheric altitudes at the peak of the PEDE storm (e.g., Yiğit et al., 2021). Nevertheless, as a baseline for our investigation, M-GTIM simulations in the current analysis utilize the static GW parameters previously specified in M-GITM simulations (e.g., horizontal wavelength = 300 km; maximum momentum flux at the source level = 0.0025 m 2 /s 2 ; and the source level height at the top of a mean PBL = 8.0 km) (K. K. J. Roeten et al., 2019).

Analysis
In Figure 1a, we plot the He density at ∼200 km from MY 32 (gray diamonds), 34 (red stars), and 35 (dark gray triangles). Due to the high orbit to orbit variability, in order to help guide the eye to pick out the peak, the black line in each data set represents the running 10-orbit average. The blue highlighted oval is the northern winter region from Ls 180-240 that is the focus of interest for this paper. This region is blown up in Figures 1b and 1c. This oval corresponds to the ovals in Figure 1d representing northern winter for each MY. Figure 1b is the NGIMS density data in a 10-orbit average from Ls 180-290 for MY 32 (gray diamonds), MY 34 (red stars), and MY 35 (dark triangles). The error bars are for the standard deviation of the mean, this is higher than the single error of the data in one orbit, though both have been accounted for in these error bars. Figure 1c is the NGIMS temperature in a 10-orbit average for MY 32 (gray diamonds), MY 34 (red stars), and MY 35 (dark triangles). The error bars are for the standard deviation of the orbit-to-orbit variation about the 10-orbit mean. This error accounts for the single orbit variation as well which is less than the orbit-to-orbit variation at this altitude. Figure 1d is the latitude, Ls and LST positional data for NGIMS for MY 32, 34, and 35. The key for this paper is that MY 34 was observing in the northern hemisphere from dawn through the nightside (∼6 to ∼20°hr) during Ls 180-240. MY 32 and 35 both were observing in the same hemisphere and either also had observation on the night or dawn side and thus were able to observe the He bulge during those years.
The 10-orbit average was not used for the data binning in Figures 4 and 5 when comparing the NGIMS data to the model. Each orbit was binned by 5° latitude and 0.5 hr LST for Figures 4 and 5 to compare with the M-GITM data. Model data were binned in the same manner for Figures 4 and 5 and the M-GITM data were reduced to keeping only the matching data from bins that also contained observational data.
There is a distinctly different pattern in the densities for MY 34 during the PEDE (Ls ∼ 240-260), for which the He density is lower where there is a measured increase observed in MY 32 and 35 due to the typical He bulge (Figure 1b). In addition, the MY 34 temperatures are higher on the nightside (LST 22-24 hr) and dawn (LST 2-6 hr) regions compared to the other MY's at the same latitude and LST (Figure 1). The higher MY 35 temperatures are from the lower latitudes and dayside. All densities are from the inbound portion of the orbit segment and are measured at ∼200 km altitude (M. Elrod et al., 2017;Gupta et al., 2021). We examined the impacts on the global circulation by conducting M-GITM simulations for MY 34 during the PEDE from its onset at Ls ∼ 184 to its peak (Ls ∼ 208-210) and through the early decay phase (Ls ∼ 240). We compared this with the observed NGIMS data available through the period. We ran the M-GITM simulation for solar forcing, globally uniform eddy diffusion parameters, and an evolving dust distribution that matches MRO/ MCS dust measurements during the evolution of the PEDE but without gravity waves (Figure 2a). We referred to this version of the model as the M-GITM without gravity waves and then compared the model density results to the data in Figure 5, the temperature results in Figure 6. The M-GITM computed He bulge is larger than expected near the dawn region above ∼50°N latitude. The simulated He bulge after midnight is also too large, with a much bigger difference to the data results. In short, M-GITM meridional and zonal winds appear to be too strong in this vicinity (Miyamoto et al., 2021;K. J. Roeten, Bougher, Yiğit, et al., 2022).
We then revised the M-GITM model to include the addition of the new non-orographic GW momentum and energy deposition scheme previously described (K. J. Roeten, Bougher, Yiğit, et al., 2022). This is referred to as the M-GITM with gravity waves version (Figure 2b). We compare the densities to the NGIMS observations in Figure 4, and the temperatures in Figure 6. Standard GW parameters are used that are proven to be robust but have not included any changes with season, or disruption that may have occurred specifically with this dust event. The major impact from the inclusion of the gravity waves is to slow winds (both zonal and meridional) and cool high latitude temperatures especially. The computed model vertical and meridional winds are shown in Figure 3. The computed He bulge is shown to move further northward (poleward of 50°N at the dawn terminator), providing a reduction in magnitude that reasonably matched available NGIMS measurement as shown in Figure 4 (with gravity waves) and 5 (without gravity waves). However, the mismatch with the post-midnight bulge indicated that the simulated bulge is still much larger than the NGIMS measurements.  and 3d are the computed vertical wind distributions for the M-GITM model without gravity waves (c) and with gravity waves (d). The red is for up direction and the blue is down direction. The bright teal is the zero velocity or calm winds. He is most likely to collect in the teal region further north (>50°) and dawn (0-6 hr).
We acknowledge that gravity waves could be evolving (i.e., non-constant) throughout the PEDE storm (Yiğit et al., 2021), yet our current GW modeling is presently limited to unchanging GW parameters. We plan to make this revision in future efforts to determine how variable gravity waves are impacting the He bulge during the PEDE.

Results
Figures 4a and 4c compare the NGIMS density data with the M-GITM computed densities without gravity waves. The M-GITM outputs are binned by the NGIMS trajectory through the dust event (Ls ∼ 180-240). Figure 4 plots the M-GITM density data that are binned with the NGIMS density data. Figure 4a plots the NGIMS log 10 density (b) plots the M-GITM computed densities and (c) the log 10 of the difference factor of the NGIMS data and the M-GITM model with gravity waves log 10 (|data-model|)/data. We chose to take the log of the difference factor because the difference between the region that matches reasonably well and that which does not was large and this was the best way to be able to see the variations in the difference factor.
Figures 5a-5c compares the NGIMS density with the M-GITM computed densities with gravity waves. The log 10 NGIMS density is plotted in Figure 5a, with the log 10 M-GITM computed densities binned to the NGIMS observed densities in Figure 5b. Figure 5c is the log 10 difference factor for the NGIMS data and the M-GITM model without gravity waves log 10 (|data-model|)/data. We chose to take the log of the difference factor because the difference between the region that matches reasonably well and that which does not was large and this was the best way to be able to see the variations in the difference factor. Both Figures 4c and 5c are on the same scale. While Figures 4 and 5b are not on the same scale. without gravity waves model data binned to match the NGIMS data during the dust storm from Ls 180-240. The He densities are plotted as a log 10 (dens) in order to bring out the variation and features. (c) The difference factor between the observed data and the model (|data-model|)/data. This is a multiplicative factor rather than a percentage as the factor was easier to plot since the numbers were high.
The NGIMS data matches best with the M-GITM model that was revised to include gravity waves. The data and the model agree well for latitudes below 50° latitude and dawn (LST 1-5 hr). The comparison between the model without the gravity waves and NGIMS data matches better for above 50° latitude and night . Because the comparison with the densities is much better with the M-GITM GW, this indicates that gravity waves play a critical role. However, the fact that there is large discrepancy above 50° N latitude and on the night side  indicates that throughout the dust storm, gravity waves possibly need variability in the implementation.
The densities in the computed M-GITM model, which include GW effects, match the NGIMS density data much better on the dawn side (LST ∼ 1-5 hr) and in the mid to high latitudes (∼0°N-50°N lat). But once the spacecraft moves closer to midnight and above 50°N, the mismatch between the model and data becomes much greater. This was also where the model predicted that the presence of the He bulge should be able to be observed. NGIMS did not observe a significant increase in the He data above 50°N latitude on the nightside around midnight (LST ∼ 22-24 hr).
Figures 6a-6f compares the NGIMS measured temperatures (a) with the computed M-GITM modeled temperatures without gravity waves (b) and with gravity waves (c) during the dust event. The M-GITM model data has been binned to correspond with the NGIMS data without gravity waves in (d) and with gravity waves in (e) and the difference between them taken and plotted in each figure respectively. Note that the scales are not the same for both (d) and (e). The temperatures do match better for the GW case for the dawn (LST 1-5 hr) and lat 20°-50° than for the non-gravity wave case. However, the case without gravity waves matches better for the latitudes >50° and night . This indicates that the model upgrade employing gravity waves is a good improvement to the model, but that the GW scheme likely needs to be employed with a variable scheme throughout the dust storm.

Discussion
The reduced He bulge density during the PEDE in 2018 indicates that there is a change to the circulation in the upper atmosphere. While NGIMS wind measurements were obtained at the beginning and during the dust with gravity waves model data binned to match the NGIMS data during the dust storm from Ls 180-240. The He densities are plotted as a log 10 (dens) in order to bring out the variation and features. (c) The difference factor between the observed data and the model (|data-model|)/data. This is a multiplicative factor rather than a percentage as the factor was easier to plot since the numbers were high. storm, the campaign averaged wind speeds were not significantly different in magnitude but were significantly different in direction. NGIMS wind measurements were limited to 2, 10-orbit passes one at the onset of the storm and a second near the end of the storm. This is rather limited data to confirm the M-GITM model given the significant variations that were measured in direction and speed from NGIMS (K. J. . While substantial winds are present throughout the dust storm in the upper atmosphere, large variations of these winds from orbit-to-orbit were also observed, possibly connected to increased GW activity, which was also observed in the thermosphere at the same time (Yiğit et al., 2021).
NGIMS neutral data has been used to catalog GW observations both on the dayside and night-side during non-dust seasons (England et al., 2017). Yiğit et al. (2021) conducted an analysis of the GW activity seen in the upper atmosphere during the PEDE by analyzing perturbations induced in the NGIMS CO 2 densities and found that thermospheric GW activity doubles during the peak of the dust storm. In addition to gravity waves causing short time-scale variability, the momentum deposited by gravity waves also modifies the middle-upper atmospheric circulation, as seen in several Global Circulation Model studies (e.g., Haberle  Medvedev & Yiğit, 2012). A modeling study examining the full impact of variations in the gravity waves on the atmospheric circulation or the He bulge has not yet been completed.
The He bulge forms in the cool and calmer parts of the atmosphere (Figures 2 and 3). Therefore, the highest levels of He are found from ∼2 to 5 hr LST at high latitudes during the winter season, and at mid to high latitudes from ∼2 to 5 hr LST during the equinoxes (M. Elrod et al., 2017;Gupta et al., 2021). The fact that in MY 34 during the dust event, this He bulge did not appear from LST 0-5 hr in the high latitudes indicates a change in the circulation of the upper atmosphere, that the atmosphere was warmer or that gravity waves and winds redistributed the He. Specifically due to the higher cooling effects of the GWs included in the M-GITM model, the computed high latitude temperatures are much cooler than observed but the non-GW temperatures were a slightly better match. This leads us to hypothesize that the gravity waves could vary throughout the dust storm and using a static GW formulation (i.e., constant GW parameters) in the model is likely producing too much cooling.
Overall, our main conclusion thus far is that the GW scheme we have employed in the M-GITM model is a significant improvement for driving a more realistic Mars thermospheric circulation (K. J. Roeten, Bougher, Yiğit, et al., 2022). As a result, an improved picture is emerging describing how He is moving in the upper atmosphere and how the thermospheric circulation is changing during the PEDE. Gravity waves may be reducing the strength of the global thermospheric circulation and thereby improving the dawn region He bulge distribution and magnitude soon after the peak of the PEDE. However, a refined GW scheme may be needed that accounts for the time evolution of GW parameters and the resulting momentum/energy deposition impacting the circulation from the PEDE peak through the decay phase. The presently computed He bulge magnitude is too large (especially above 50°N latitude around midnight) and does not match the NGIMS observed values.
MAVEN continues to operate, and we hope to have additional observations of the upper atmosphere for the next PEDE from all the Mars assets available for the next major dust storm. Limitations with NGIMS observations, given that it is the only in situ upper atmosphere instrument currently delivering He density, is one other factor to be considered in the distribution of the He bulge. While other Mars orbiters such as Mars Orbiter Mission, which has a similar instrument to NGIMS, were also observing the PEDE , it is unfortunately very limited in the light species observations (Bardwaj et al., 2016).