Variations in the Ionospheric Peak Altitude at Mars in Response to Dust Storms: 13 Years of Observations from the Mars Express Radar Sounder

TL;DR

This study analyzes 13 years of Mars Express radar data to show that dust storms cause the ionospheric peak altitude to rise and increase atmospheric variability, revealing enhanced coupling between atmospheric layers.

Contribution

It provides the first long-term observational evidence of ionospheric peak altitude responses to dust storms on Mars across multiple Mars Years.

Findings

Peak altitude increased by 10-15 km during dust storms.

Atmospheric variability was significantly higher during dust seasons.

Global dust storms caused the largest increases in variability.

Abstract

Previous observations have shown that, during Martian dust storms, the peak of the ionosphere rises in altitude. Observational studies of this type, however, have been extremely limited. Using 13 years of ionospheric peak altitude data from the Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS) instrument on Mars Express, we study how the peak altitude responded to dust storms during six different Mars Years (MY). The peak altitude increased $\sim$10-15 km during all six events, which include a local dust storm (MY 33), three regional regional dust storms (MYs 27, 29, and 32), and two global dust storms (MYs 28 and 34). The peak altitude's orbit-to-orbit variability was exceptionally large at the apexes of the MY 29 and MY 32 dust seasons, and dramatically increased during the MY 28 and MY 34 global dust storms. We conclude that dust storms significantly increase upper atmospheric variability, which suggests that they enhance dynamical processes that couple the lower and upper atmospheres, such as upward propagating gravity waves or atmospheric tides.

Tables

Table 1. Table 1: Summary of the observed ionospheric response to dust for the six Mars years (MY) analyzed in Section  4 . Columns 2-6 list the solar longitudes (L s , ∘ ), Earth years, dates (mm/dd), solar zenith angles (SZA, ∘ ), and minimum and maximum latitudes (LAT, ∘ ) of the MARSIS observations that were analyzed in each MY. Column 7 lists the EUV irradiance level (EUV, mW/m 2 ). Columns 8-9 summarize each dust storm, including the increase in the globally-averaged dust optical depth ( Δ ​ τ Δ 𝜏 \Delta\tau ), and the maximum equatorial-averaged (between ± plus-or-minus \pm 30 ∘ ) dust optical depth ( τ m ​ a ​ x subscript 𝜏 𝑚 𝑎 𝑥 \tau_{max} ). Columns 10-12 summarize the observed changes in the ionospheric peak, including its increase in altitude ( Δ ​ h m ​ a ​ x Δ subscript ℎ 𝑚 𝑎 𝑥 \Delta h_{max} , km), and its variability before ( σ h m ​ a ​ x i superscript subscript 𝜎 subscript ℎ 𝑚 𝑎 𝑥 𝑖 \sigma_{h_{max}}^{i} , km) and after ( σ h m ​ a ​ x f superscript subscript 𝜎 subscript ℎ 𝑚 𝑎 𝑥 𝑓 \sigma_{h_{max}}^{f} , km) dust storm onset as defined throughout Section  4 .

MY	Ls	Year	Date	SZA	LAT	EUV	$\Delta \tau$	${\tau }_{max}$	$\Delta h_{max}$	${\sigma }_{h_{max}}^i$	${\sigma }_{h_{max}}^f$
27	305-330	2005	10/12 - 11/24	20-50	(-20, -70)	1.3	0.2	0.3	10-15	4	4
28	250-295	2007	06/03 - 08/14	30-40	(-50, -65)	1.6	0.6	1.2	10-15	4	11
29	230-245	2009	03/20 - 04/05	45-55	(-65, -75)	1.2	0.1	0.5	10-20	11	11
32	215-230	2014	10/16 - 11/10	75-80	(40, 50)	2.4	0.1	0.5	10-15	5	7
33	304-310	2017	01/23 - 02/21	70-80	(-20, -40)	1.1	0.0	0.2	$\sim$10	–	–
34	180-195	2018	05/22 - 06/26	60-85	(60, 85)	1.0	0.5	1.5	10-15	–	–

Equations

$$n(h_{max})\sigma H\sec \chi =1$$

$$n_0{e}^{-(h_{max}-h_0)/H}\sigma H\sec \chi =1.0$$

$$n_{0_1}{e}^{-(h_{ma{x}_i}-h_0)/H_i}{H}_i=n_{0_f}{e}^{-(h_{ma{x}_f}-h_0)/H_f}{H}_f=1.0$$

$$\frac{n_{0_f}{H}_f}{n_{0_i}{H}_i}=\exp \left(\frac{\Delta h_{max}}{H_f}\right)$$

$$\frac{P_f}{P_i}=\exp \left(\frac{\Delta h_{max}}{H}\right).$$

Full text

\draftfalse\journalname

arXiv

Department of Physics and Astronomy, University of Iowa, Iowa City, IA, USA Space Science Institute, Boulder, CO, USA Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO, USA Faculty of Mathematics and Physics Charles University, Prague, Czech Republic

\correspondingauthor

Zach Girazianzachary-girazian@uiowa.edu

{keypoints}

We show how the ionospheric peak altitude at Mars varies during six different dust storms.

The peak altitude rises 10-15 km during each dust storm.

Dust storms increase the peak altitude’s variability, suggesting they enhance dynamical processes coupling the lower and upper atmospheres.

1 Introduction

Dust storms have proven to be an important source of atmospheric variability at Mars. Dust particles lifted into the atmosphere are heated by solar radiation, which causes the atmosphere to expand, and global circulation patterns to change (Haberle \BOthers., \APACyear1993; S\BPBIW. Bougher \BOthers., \APACyear1997; Heavens \BOthers., \APACyear2011; Wolkenberg \BOthers., \APACyear2018). Dust also affects the distribution of water vapor, which has consequences for atmospheric photochemistry, hydrogen escape, and climate evolution (Chaffin \BOthers., \APACyear2017; Heavens \BOthers., \APACyear2018; Daerden \BOthers., \APACyear2019; Krasnopolsky, \APACyear2019; Vandaele \BOthers., \APACyear2019). Although the dust itself is mostly confined to altitudes below 80 km (Clancy \BOthers., \APACyear2010), the effects of the dust extend well into the thermosphere ($>$100 km) - even to geographic locations that are far from where the dust originated (Withers \BBA Pratt, \APACyear2013; Liu \BOthers., \APACyear2018).

The response of the thermosphere to lower atmospheric dust is generally marked by a rapid increase in the neutral density at a fixed altitude, followed by a slow density decay back to nominal levels over several weeks (Keating \BOthers., \APACyear1998; Lillis \BOthers., \APACyear2010; England \BBA Lillis, \APACyear2012; Withers \BBA Pratt, \APACyear2013; Zurek \BOthers., \APACyear2017). During typical dust events, thermospheric neutral densities increase by a factor of $\sim$1.5-3.0 at fixed altitudes (Withers \BBA Pratt, \APACyear2013; Zou \BOthers., \APACyear2016; Liu \BOthers., \APACyear2018), but thermospheric neutral temperatures do not drastically change at fixed pressure levels (McElroy \BOthers., \APACyear1977; Wang \BBA Nielsen, \APACyear2003; Fang \BOthers., \APACyear2019). When thermospheric pressure surfaces rise in response to dust loading, the peak of the ionosphere - which typically forms between 120-150 km and at a fixed pressure level (Withers, \APACyear2009) - rises in altitude.

Elevated ionospheric peak altitudes during dust storms have been observed by radio occultation (RO) experiments on several spacecraft, the first being Mariner 9, which arrived at Mars in 1971 during a global event. Mariner 9 observed that, during the global dust storm, the ionospheric peak altitude was $\sim$20 km higher than usual (Hantsch \BBA Bauer, \APACyear1990), and then slowly decayed back to typical values during the dust storm’s waning stages (Withers \BBA Pratt, \APACyear2013). The RO experiments on the Mars Global Surveyor (MGS) and Mars Atmosphere and Volatile EvolutioN (MAVEN) spacecraft have also observed the ionospheric peak rise during regional dust storms. MGS observed the peak rise $\sim$5 km during a dust storm at Solar Longitude (Ls) 130*∘* in Mars Year (MY) 27 (Withers \BBA Pratt, \APACyear2013; Qin \BOthers., \APACyear2019), and MAVEN observed the peak rise $\sim$10 km during a dust storm at Ls $\sim$305*∘* in MY 33 (Withers \BOthers., \APACyear2018).

Motivated by the limited amount of observations showing how the ionospheric peak responds to dust storms, we utilize a 13-year span of peak altitude measurements from the Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS) instrument on Mars Express (MEX) (Picardi \BOthers., \APACyear2004; Gurnett \BOthers., \APACyear2005). By combining these data with dust optical depth measurements from the same time period (Montabone \BOthers., \APACyear2015), we investigate how dust affected the ionospheric peak altitude during six different MYs. In two of these years, MY 28 and MY 34, there was a global dust storm. Our objectives are (1) to compare the response of the peak altitude during six different dust storms and (2) to determine if dust storms affect the variability of the peak altitude.

The paper is organized as follows. In Section 2, we present the theory that describes the formation of the ionospheric peak, and explain how variations in the ionospheric peak altitude can be used to estimate changes in the thermospheric pressure. In Section 3, we describe the data sets that are used in our analysis. In Section 4, we show how dust affected the peak altitude during several MYs. In Section 5, we discuss our results and present our conclusions.

2 Theory of the Ionospheric Peak Altitude

The main peak of the Martian ionosphere is well-described by Chapman theory (Chapman, \APACyear1931; Schunk \BBA Nagy, \APACyear2009; Withers, \APACyear2009; Girazian \BBA Withers, \APACyear2013; Mendillo \BOthers., \APACyear2017). Chapman theory predicts that, under photochemical equilibrium conditions, the ionospheric peak forms at the altitude where the optical depth of ionizing extreme ultraviolet (EUV) photons is equal to one (Breus \BOthers., \APACyear2004; Withers, \APACyear2009). Mathematically, this is approximated by

[TABLE]

where $n(h_{max})$ is the neutral CO2 density at the peak altitude $h_{max}$, $\sigma$ is the CO2 absorption cross section at EUV wavelengths, $H$ is the neutral scale height, and $\chi$ is solar zenith angle (SZA). Equation 1 requires many simplifying assumptions to be valid (Schunk \BBA Nagy, \APACyear2009; Withers, \APACyear2009) but has proven to be an adequate description of the ionospheric peak at Mars (Withers, \APACyear2009; Girazian \BBA Withers, \APACyear2013; Fallows \BOthers., \APACyear2015; Mendillo \BOthers., \APACyear2017). Using Equation 1, several studies have shown that observed changes in the ionospheric peak altitude can be used to estimate variations in the thermospheric neutral density or pressure (S\BPBIW. Bougher \BOthers., \APACyear2001, \APACyear2004; Withers \BBA Pratt, \APACyear2013; Zou \BOthers., \APACyear2011, \APACyear2016; Qin \BOthers., \APACyear2019).

For our purposes, we are interested in using Equation 1 to quantify how the peak altitude rises or falls in response to changes in the thermospheric pressure. Equation 1 is derived under the assumption of a static, isothermal atmosphere so that it can be recast as

[TABLE]

where $n_{0}$ is the neutral density at some reference altitude $h_{0}$. If $n_{0}$ or $H$ changes, $h_{max}$ must rise or fall such that Equation 2 remains satisfied:

[TABLE]

where the subscripts $i$ and $f$ represent the initial and final states, respectively. In Equation 3, we have removed the SZA dependence by assuming that it is fixed. Next, by setting the reference altitude to $h_{0}=h_{max_{i}}$, Equation 3 becomes

[TABLE]

where $\Delta h_{max}=h_{max_{f}}-h_{max_{i}}$. According to Equation 4, if we assume that $H_{i}=H_{f}$, then a factor of $\sim$3 increase in $n_{0}$ results in the peak altitude increasing by one scale height. Furthermore, Equation 4 states that increases in the neutral density or scale height during dust storm onset will cause the peak altitude to rise, while decreases in the density or scale height during the waning stages of a dust storm will cause the peak altitude to fall.

Equation 4 can also be recast in terms of atmospheric pressure. The pressure is given by $P=\rho gH$, where $P$ is pressure, $\rho$ is mass density, and $g$ is the gravitational acceleration, which is assumed to be constant. This equation for atmospheric pressure is equivalent to the ideal gas law for an isothermal atmosphere with a fixed scale height $H$. The assumption of a fixed $H$ is used throughout this work because observations and modeling suggest that dust storms significantly affect thermospheric densities, but only modestly affect thermospheric temperatures (McElroy \BOthers., \APACyear1977; Wang \BBA Nielsen, \APACyear2003; Liu \BOthers., \APACyear2018; Fang \BOthers., \APACyear2019; Qin \BOthers., \APACyear2019). This assumption does, however, add uncertainty to our analysis.

Substituting the pressure into Equation 4 gives

[TABLE]

Here, $P_{i}$ and Pf are the atmospheric pressures at the ionospheric peak and $H$ has been assumed to not change. Equation 5 states that, at a fixed SZA, the peak altitude forms at a constant atmospheric pressure level. Equation 5 was used in the study by Withers \BBA Pratt (\APACyear2013) to estimate changes in the thermospheric pressure during the waning stages of the Mariner 9 dust storm.

In addition to allowing one to quantify how the peak altitude responds to changes in the thermospheric density, scale height, or pressure, Eqs. 1, 4, and 5 summarize the conditions that control variations in the ionospheric peak altitude. Equation 1 can be inverted to show that the peak altitude increases with increasing SZA proportional to $ln(\sec\chi)$, indicating a steep rise in the peak altitude near the day-night terminator (Withers, \APACyear2009; Fallows \BOthers., \APACyear2015). Equation 5 describes how the peak altitude can vary with latitude and local time due to diurnal pressure gradients in the thermosphere (Zou \BOthers., \APACyear2011; S\BPBIW. Bougher, Pawlowski\BCBL \BOthers., \APACyear2015), and further shows that any physical process that alters the thermospheric pressure will also alter the peak altitude. These processes include, but are not limited to (1) the annual variation in solar irradiance at Mars due to its large orbital eccentricity; (2) atmospheric tides and waves; (3) solar extreme ultraviolet (EUV) heating; and (4) atmospheric circulation patterns. All of these processes must be considered when attempting to identify variations in the peak altitude caused solely by dust storms.

3 Data

We use three types of data: ionospheric peak altitudes from the MARSIS radar sounder, dust optical depths compiled from several instruments, and solar EUV irradiance from an empirical model. In this section, we describe each data set and our processing techniques.

3.1 Peak Altitudes

Ionospheric peak altitudes are derived from MARSIS radar sounding observations. MARSIS sounds the ionosphere during the periapsis segment of the $\sim$7-hour, near-polar orbit of MEX. When MARSIS sounds the ionosphere, it transmits radio pulses and records the return echoes from pulses that are reflected off the ionosphere. The transmitter sweeps through 160 qausi-logarithmically spaced frequencies between 0.1 and 5.4 MHz over 1.257 seconds. Each sweep produces an ionogram - the echo intensity as a function of frequency and time delay. MARSIS makes frequency sweeps every 7.54 seconds while MEX is below $\sim$1500 km, returning several hundred ionograms during each periapsis pass.

Figure 1a shows an ionogram from 16 May 2018. The vertical stripes at low frequencies are a common feature; they are produced by electron plasma oscillations induced by the radio transmitter during ionospheric soundings. The spacings between the stripes are used to determine the local electron density at the spacecraft location (Duru \BOthers., \APACyear2008; Gurnett \BOthers., \APACyear2005, \APACyear2008). The thin horizontal stripe between $\sim$1.0-3.0 MHz is the ionospheric echo, which represents radar pulses that were reflected off the ionosphere below $\sim$200 km (Gurnett \BOthers., \APACyear2005, \APACyear2008). The highest frequency of the ionospheric echo, near 3 MHz in this example, is the reflection from the ionospheric peak. The time delay at this frequency in the ionospheric trace is related to the altitude of the peak. The horizontal stripe at frequencies greater than 3 MHz, with a time delay of $\sim$4 ms, is the return echo from the surface of Mars.

If an ionogram provides a local electron density measurement, and also has a clear ionospheric trace, then it can be inverted into an altitude profile of the electron density. The inversion procedure requires one to make an assumption about the shape of the electron density profile between the altitude of MEX, where the local density is measured, and the altitude of the first electron density measurement in the ionospheric trace. In this work, we adopt the inversion technique described in Němec \BOthers. (\APACyear2016), which assumes that the shape of the density profile within the measurement gap is characterized by a lower and upper topside scale height, with a smooth transition between them. Figure 1b shows the electron density profile that was derived by applying this technique to the ionogram shown in Figure 1a. Once inverted into an electron density profile, the peak altitude is easily extracted as marked in Figure 1b.

We use MARSIS ionograms obtained between 11 July 2005 and 14 July 2018. All of the ionograms obtained up to 22 May 2016 were inverted into electron density profiles. Only a subset of the ionograms obtained after 22 May 2016 were used, because they have not yet been inverted into electron density profiles due to the time consuming and hands-on processing techniques that are required (Gurnett \BOthers., \APACyear2008; Morgan \BOthers., \APACyear2008). In lieu of this, we have inverted two subsets of ionograms from after this date specifically for this study. The first subset of ionograms is from the month of Jan. 2017. This month was targeted because elevated peak altitudes were observed by the MAVEN RO experiment during a dust storm (Withers \BOthers., \APACyear2018). The second subset is from 14 May 2018 through 14 July 2018. This period was targeted because it covers a significant portion of the 2018 global dust storm (Guzewich \BOthers., \APACyear2019; Vandaele \BOthers., \APACyear2019).

After inverting the ionograms into electron density profiles and extracting their peak altitudes, we filter the data set based on several criteria. First, we keep density profiles only if they monotonically decrease with increasing altitude, which is a requirement of the Němec \BOthers. (\APACyear2016) inversion technique. Second, we remove any peak altitudes below 80 km or above 220 km, well outside the expected peak altitude and likely the result of bad inversions (Withers, \APACyear2009; Fox \BBA Weber, \APACyear2012; Němec \BOthers., \APACyear2016; Vogt \BOthers., \APACyear2017). Third, we limit the profiles to times when MEX was below 1000 km. After applying these criteria, the complete data set includes more than 180,000 peak altitude measurements from 2401 MEX orbits.

Uncertainties in the peak altitudes are, at best, $\pm$7 km, as determined by the intrinsic 91.4 $\mu$s time resolution of the MARSIS receiver. The Němec \BOthers. (\APACyear2016) inversion technique also adds to this uncertainty, given that it relies on an assumption about the shape of the electron density profile within the measurement gap. In Section 4, we analyze orbit-averaged peak altitudes (from a limited SZA range), and assign uncertainties based strictly on the spread of the observed peak altitudes during each orbit. In particular, for each orbit, we define the uncertainty in the orbit-averaged peak altitude as the standard deviation of the peak altitudes used to calculate the average.

3.2 Dust Optical Depths

Dust optical depth maps are derived by combining measurements from several instruments, as described in Montabone \BOthers. (\APACyear2015). The maps provide a continuous measure of the dust content in the lower atmosphere, with good coverage in latitude, longitude, and Ls from 1999 until present. The dust data for the time period considered here are synthesized from infrared observations by the Mars Odyssey Thermal Emission Imaging System (THEMIS) (Christensen \BOthers., \APACyear2004) and the Mars Reconnaissance Orbiter Mars Climate Sounder (MCS) (McCleese \BOthers., \APACyear2007). The maps provide the optical depth at 9.3 $\mu$m absorption, normalized to an atmospheric pressure level of 610 Pascals. An example of a dust optical depth map is shown in Figure 1c.

Throughout this work, we use the observation-only dust maps that sometimes have incomplete coverage in latitude and longitude (Montabone \BOthers., \APACyear2015). We also use dust maps that were developed specifically for the MY 34 global dust storm and this special issue. These MY 34 dust maps use estimated column dust opacities from MCS as described in detail in Montabone \BOthers. (\APACyear2019, submitted to this issue). We use the v2.5 version of the maps.

We assign a “local” average and global average dust optical depth to each MARSIS electron density profile. The local average accounts for our expectation that the ionospheric peak will respond to dust storms that are nearby. The global average accounts for our expectation that the ionospheric peak might also respond to dust storms that are far away, due to their effects on atmospheric circulation (Bell \BOthers., \APACyear2007; Withers \BBA Pratt, \APACyear2013; Withers \BOthers., \APACyear2018). To assign the global average, we match each MARSIS electron density profile with the dust optical depth map from the closest date, and then average the dust optical depth map over all latitudes and longitudes. Since the dust maps have a resolution of $\sim 1^{\circ}$ in Ls, the global average dust optical depth is constant during each MEX orbit. We assign the local average in a similar way but only after restricting the dust data to the latitude range of the MARSIS observations during each orbit.

3.3 Solar EUV Irradiances

Solar EUV irradiances are from the Flare Irradiance Spectral Model-Mars (FISM-M), which is an empirical model that provides daily-averaged solar EUV spectra at Mars (Thiemann \BOthers., \APACyear2017). The spectra have 1 nm resolution and cover wavelengths between 0.5-189.5 nm. An example FISM-M spectrum, from 16 May 2018, is shown in Figure 1d.

We assign a single number of the EUV irradiance (W m*-2*) to each MARSIS electron density profile. We do this by first matching each electron density profile with the solar EUV spectrum from the closest date, then integrating the matched EUV spectrum over wavelengths between 0.5-92.5 nm. The cutoff wavelength of 92 nm is chosen because it is the longest wavelength photon that can ionize O and CO2, which are the most abundant neutral species in the thermosphere of Mars (Schunk \BBA Nagy, \APACyear2009; Girazian \BBA Withers, \APACyear2013, \APACyear2015; Mahaffy \BOthers., \APACyear2015).

3.4 Overview of the Data

Figures 2a-c show the 13-year time series of ionospheric peak altitudes along with the SZAs and geographic latitudes of the MARSIS observations. Data gaps are present throughout the time series because MARSIS frequently toggles between ionospheric and subsurface sounding modes, and does not make observations during eclipse seasons when there is insufficient sunlight to recharge the MEX solar panels. The error bars in Panels b-c show the $\sim$80*∘* in latitude and the $\sim$30*∘* in SZA that MARSIS covers during a typical periapsis pass. We also note that, from one periapsis pass to the next, the observational SZA and latitude are nearly identical, but the longitude is shifted by $\sim$100*∘*.

Figure 2d shows the solar EUV irradiance and the inverse-square of the Mars-Sun distance during the MARSIS observations. The latter is representative of the solar insolation at Mars, which varies annually due to the planet’s eccentric orbit around the Sun. The 13-year period covers more than a complete solar cycle, starting with the declining phase of Solar Cycle 23 in 2005, and continuing through the declining phase of Solar Cycle 24 in 2018. The EUV irradiance varies over the 11-year solar cycle, reaching a minimum in 2008, a maximum in 2014, and also annually due to the varying Mars-Sun distance.

Figure 2e shows the global average dust optical depth from the same time period. The repeated peaks in the optical depth during each Martian Year are the signature of the well known annual dust cycle at Mars (Montabone \BOthers., \APACyear2015; Fang \BOthers., \APACyear2019). During most years, the optical depth reaches a maximum between Ls 210*∘-240∘, and then exhibits a smaller, secondary peak between Ls* 300*∘-340∘. Exceptions to the typical annual dust cycle are seen in MYs 28 and 34, when global dust storms caused large spikes in the dust optical depth at atypical Ls* values (Montabone \BOthers., \APACyear2015, \APACyear2019, submitted to this issue; Wolkenberg \BOthers., \APACyear2018; Guzewich \BOthers., \APACyear2019; Sánchez-Lavega \BOthers., \APACyear2019).

Now that we have provided an overview of the data used in our study, we are ready to focus on our science objectives, which are (1) to compare the response of the peak altitude during six different dust storms and (2) to determine if dust storms affect the variability of the peak altitude. In the next section, we will focus on these objectives by closely examining how the ionospheric peak altitude responded to dust storms during six different MYs. For each MY, we analyze a subset of MARSIS peak altitudes from a limited range of Ls and SZA. With the exception of MY 33, the Ls range is chosen to capture the onset of a dust storm, and the SZA range is chosen to capture the smallest SZAs observed by MARSIS during that period.

To provide global context, Figure 3 shows the complete dust maps from the six MYs we perform our analysis. Each panel in Figure 3 is marked with a rectangle that shows the Ls and latitude coverage of the MARSIS peak altitudes that we analyze during that MY. Figure 3 also highlights that, with the exception of MY 27, the dust maps have full latitudinal coverage during the time periods considered.

The observations in MYs 27-29 provide the most favorable conditions because MARSIS observations cover dayside SZAs $<$ 55*∘* during times when the global dust content significantly increases. The MY 29 observations at southern latitudes, and the MY 32 observations at northern latitudes, are from similar dust conditions. This allows us to compare the effects of dust in two different hemispheres. The MY 30 and MY 31 observations are not presented because dayside observations were not obtained during the onset of dusty periods. The MY 33 observations are specifically targeted to compare with MAVEN observations from the same period (Withers \BOthers., \APACyear2018). Finally, the MY 34 observations, which are from northern polar latitudes, cover the onset of the MY 34 global dust storm.

Based on previous studies, we expect that high local dust content will elevate the peak altitude, and that high global dust content may elevate the peak altitude due to its effects on atmospheric circulation patterns (Hantsch \BBA Bauer, \APACyear1990; Bell \BOthers., \APACyear2007; Withers \BBA Pratt, \APACyear2013; Zou \BOthers., \APACyear2016; Liu \BOthers., \APACyear2018; Fang \BOthers., \APACyear2019; Qin \BOthers., \APACyear2019). Additionally, increases in the peak altitude might be less noticeable during low solar EUV levels when the upper atmosphere is intrinsically more variable due to increased gravity wave activity (England \BOthers., \APACyear2017; Terada \BOthers., \APACyear2017; Harada \BOthers., \APACyear2018; Siddle \BOthers., \APACyear2019).

4 Peak Altitude Variations During Six Dust Storms

4.1 Mars Year 27

Figure 4 summarizes the observations from late in the MY 27 dust season between Ls 305*∘-330∘. In Figure 4a, the MARSIS observational coverage in latitude and SZA is plotted on top of the dust optical depth map. The MARSIS data from this time period are restricted to SZAs between 20∘-50∘* to rule out SZA being a significant factor in any observed peak altitude variations. The dust map shows a rapid increase in the atmospheric dust content starting at Ls 310*∘* at most latitudes. Dust optical depth measurements at latitudes greater than $+$45*∘* were unavailable during this period (Montabone \BOthers., \APACyear2015). Figure 4b shows the global and local average dust optical depths (Section 3.2) for each MEX orbit during which MARSIS obtained at least 10 peak altitude measurements. The optical depths increase sharply at Ls $\sim$305*∘* and then slowly decline through Ls $\sim$320*∘*.

Figure 4c shows the corresponding orbit-averaged ionospheric peak altitudes, with error bars representing the spread in the peak altitude measurements (1$\sigma$) from each MEX orbit. The peak altitude rapidly rises 10-15 km between Ls 307*∘-315∘, and then slowly descends back to pre-dust storm values over $\sim 15^{\circ}$ of Ls*. The peak altitude’s rapid rise and slow descent is consistent with previous reports of the peak altitude’s response to dust storms (Withers \BBA Pratt, \APACyear2013; Withers \BOthers., \APACyear2018).

Using the observed peak altitudes and Equation 5, we can estimate how the thermospheric pressure changed during the dust storm (Withers \BBA Pratt, \APACyear2013; Qin \BOthers., \APACyear2019). In Equation 5, we set $h_{max_{1}}$ to a constant reference altitude equal to the average peak altitude prior to dust storm onset (L${}_{s}<307^{\circ}$). Then, we set the neutral scale height, $H$, to a fixed value of 12 km, which is a typical value for the dayside thermosphere (Withers, \APACyear2006; Mahaffy \BOthers., \APACyear2015; Zurek \BOthers., \APACyear2017). This choice introduces some uncertainty in the derived pressure changes because $H$ can vary between 8-16 km. With $h_{max_{1}}$ and $H$ fixed, Equation 5 gives the pressure ratio, $P_{f}/P_{i}$, at every point where there is an orbit-averaged peak altitude. This method was used by Withers \BBA Pratt (\APACyear2013) to estimate thermospheric pressure changes during the Mariner 9 global dust storm. We adopt their terminology by calling $P_{f}/P_{i}$ the “relative pressure” ($P_{rel}$).

The peak in the relative pressure, shown in Figure 4d, coincides with the maximum in the dust optical depth at L${}_{s}\sim$315*∘*, and reaches a value of $P_{rel}=2.4\ (\pm 2.4)$. The standard deviation of $P_{rel}$ does not take into account the range of possible $H$ values, but is large due to the highly variable peak altitudes observed during this orbit. The maximum pressure increase of 2.4 is on the same order as the value derived during the waning stages of the Mariner 9 global dust storm (Withers \BBA Pratt, \APACyear2013).

Next, we test if the variability of the peak altitude is affected by the dust storm. To accomplish this, we define a variability metric, $\sigma_{h_{max}}$, and compare its value before and after dust storm onset. The metric is a measure of the peak altitude’s orbit-to-orbit variability. We determine $\sigma_{h_{max}}$ by first calculating the absolute differences between adjacent orbit-averaged peak altitudes, and then setting $\sigma_{h_{max}}$ equal to the average of the differences within a specified Ls range. Using this procedure, we find that $\sigma_{h_{max}}$ = 4 km both before (L${}_{s}<312^{\circ}$) and after (L${}_{s}>312^{\circ}$) the peak of the storm. Thus, the dust storm did not significantly affect the variability of the peak altitude.

To conclude, the peak altitude increased by 10-15 km during the MY 27 dust storm, but the variability of the peak altitude was unaffected. This MY 27 case study is summarized in Table 1, which includes a description of the MARSIS observing conditions, the dust storm, and the response of peak altitude. Table 1 also includes a summary of the five other case studies presented throughout this section.

4.2 Mars Year 28

Figure 5 summarizes the observations from the MY 28 global dust storm between Ls 250*∘-300∘. During the global dust storm, the lower atmospheric dust content was atypically large, and wide-spread around the planet (Montabone \BOthers., \APACyear2015; Wolkenberg \BOthers., \APACyear2018). As Figure 5a shows, the MARSIS data from this period are from SZAs 30∘-40∘, and from southern latitudes between -50∘* and -65*∘. Figures 5a-5b show that the local and global dust optical depths begin to increase at Ls* $\sim$265*∘* and reach maxima between Ls $\sim$275*∘-$\sim$280∘. The MARSIS observations do not cover the entire event because dust levels remained elevated through Ls* $\sim$320*∘* (Montabone \BOthers., \APACyear2015; Wolkenberg \BOthers., \APACyear2018).

The ionospheric peak altitude, shown in Figure 5c, begins to rise at the onset of the dust storm, and continues to rise throughout the observational period. Furthermore, the variability of the peak altitude significantly increases after the peak of the dust storm (Ls $\sim$280*∘*). From before dust storm onset (L${}_{s}<265^{\circ}$) to after the dust storm peak (L${}_{s}>280^{\circ}$), $\sigma_{h_{max}}$ more than doubles, increasing from 4 km to 10 km. (Although not shown, we also note that the highly variable peak altitudes have no longitudinal trend).

Such variability was not observed in MY 27 (Figure 4), during which the MARSIS observations covered similar latitudes and SZAs, but during which the dust optical depth was a factor of $\sim$2 smaller (Figure 4). One possible explanation is that the atypically high dust levels during the MY 28 global dust storm increased upper atmospheric variability, perhaps by enhancing upward propagating waves, atmospheric circulation, or atmospheric tides (Bell \BOthers., \APACyear2007; Medvedev \BOthers., \APACyear2013; S\BPBIW. Bougher, Cravens\BCBL \BOthers., \APACyear2015, and references therein). The relative pressure, shown in Figure 5d, is extremely variable after the peak of the dust storm, changing by as much as a factor of 15 across the time period, implying that the upper atmosphere is perturbed both spatially and temporally.

In summary, the peak altitude increased by 10-15 km and became highly variable during the MY 28 global dust storm.

4.3 Mars Year 29

Figure 6 summarizes a small number of observations from the peak of the MY 29 dust season between Ls 230*∘-245∘. As Figure 6a shows, the MARSIS observations during this period are from SZAs 45∘-55∘* and from southern latitudes between -65*∘* and -75*∘. Figures 5a-5b show that, at the start of observing period, the dust optical depth is already higher than usual because the MY 29 dust season began several months earlier at Ls* 170*∘* (Figure 3). Nonetheless, the local and global dust optical depths increase by 0.1 between 235*∘-245∘*.

The peak altitude, shown in Figure 6c, potentially increases between Ls 234*∘-238∘*, but the trend is weak owing to the peak altitude’s remarkable variability. The variability metric is $\sigma_{h_{max}}$ = 11 km throughout the entire period and does not change after the global average dust optical depth increases (L${}_{s}>238^{\circ}$). The magnitude of the variability metric is comparable to that derived for the MY 28 global dust storm (Section 4.2). The relative pressure, shown in Figure 6d, is also remarkably variable and changes by as much as a factor of 10.

Comparing the MY 29 and MY 27 (Figure 4) case studies reveals some stark differences. MY 27 is marked by low variability ($\sigma_{h_{max}}$ = 4 km) and the peak altitude exhibits a clear increase in tandem with the dust optical depth. Meanwhile, MY 29 is marked by such high variability ($\sigma_{h_{max}}$ = 11 km) that the increasing peak altitude trend is comparable to the orbit-to-orbit variations.

Differences in dust content during the MY 27 and MY 29 observational periods may explain these differences. The global dust optical depths in MY 27 rapidly increase by a factor of three (0.1-0.3) over 7*∘* in Ls, while the optical depths in MY 29 increase by only a factor of 1.5 (0.25-0.35) over the same time period. Another consideration is that, in MY 29 the dust season was well underway when the observing period began, which may explain the exceptional variability. To contrast this, the MY 27 dust storm was weaker as whole (Figure 3), and the global dust content was smaller before the event began (it happened later in the year during the second peak of the annual dust cycle).

To conclude, the peak altitude may have increased 10-20 km during the MY 29 dust storm, but the increase is comparable to the exceptionally high variability observed at that time.

4.4 Mars Year 32

Figure 7 summarizes observations from the MY 32 dust season between Ls 215*∘-230∘. The observing conditions are very similar to those in MY 29 (Figure 6). In both cases, the dust level is already elevated when the observing period begins, and the dust optical depths increase by a factor of $\sim$1.5 over $\sim$10∘* of Ls. In MY 32, however, the MARSIS observations are now in the northern hemisphere instead of the southern hemisphere, and they cover a higher SZA range between 75*∘-80∘*.

The peak altitude, shown in Figure 7c, increases 10-15 km between Ls 221*∘-222∘* and is somewhat elevated through Ls 225*∘*. However, similar to MY 29, is is difficult to conclusively determine if this trend is statistically significant because the peak altitude exhibits substantial variability. The variability metric is $\sigma_{h_{max}}$ = 5 km at L${}_{s}<221^{\circ}$ and $\sigma_{h_{max}}$ = 7 km at L${}_{s}>223^{\circ}$. The relative pressure, shown in Figure 7d, also exhibits substantial variability; it changes by as much as a factor of eight.

Although the variability is high in MY 32, it is not as high as in MY 29 ($\sigma_{h_{max}}$ $\simeq$ 7 km vs. $\sigma_{h_{max}}$ = 11 km), which might be a consequence of the different latitudes of the observations. Both are from the southern spring season, but the MY 29 data are from the southern hemisphere, and the MY 32 data are from the northern hemisphere.

Another consideration is that the solar EUV irradiance was a factor of two smaller during the MY 29 than during the MY 32 period (2.4 mW/m2 vs. 1.2 mW/m2, Figure 2d). The higher variability observed in MY 29, then, might be explained by increased atmospheric gravity wave activity. Gravity wave activity is anticorrelated with the thermospheric temperature (England \BOthers., \APACyear2017; Terada \BOthers., \APACyear2017; Harada \BOthers., \APACyear2018). Given that the solar EUV heating rate was significantly smaller in MY 29, thermospheric temperatures were likely lower (S\BPBIW. Bougher, Pawlowski\BCBL \BOthers., \APACyear2015; Thiemann \BOthers., \APACyear2018), and strong gravity wave perturbations were likely present. These perturbations in the thermospheric pressure would drive significant variations in the ionospheric peak altitude (Equation 5), making the peak altitudes more variable in MY 29 than in MY 32.

If this hypothesis is correct, we might also expect higher variability in MY 27 when the EUV irradiance was only 1.3 mW/m2. The peak altitude, however, is less variable in MY 27 ($\sigma_{h_{max}}$ = 4 km) than in both MY 29 ($\sigma_{h_{max}}$ = 11 km) and MY 32 ($\sigma_{h_{max}}$ $\simeq$ 7 km). This once again points to differences in the magnitudes of the dust storms. During the weak MY 27 dust storm, we observe a stable peak altitude that clearly increases in altitude, but during the relatively stronger MY 29 and MY 32 dust storms, we observe highly variable peak altitudes and less significant increases. This comparison suggests dust storms increase upper atmospheric variability, and that the magnitude of the variability is proportional to the strength of the dust storm.

In summary, the peak altitude may have increased 10-15 km during the MY 32 dust storm, but, like MY 29, the increasing trend is comparable to the orbit-to-orbit variability.

4.5 Mars Year 33

Figure 8 shows a small subset of MARSIS observations from MY 33 between Ls 304*∘-310∘. The MARSIS observations, shown in Figure 8a, are from SZAs 70∘-80∘* and southern latitudes between -20*∘* and -40*∘. These observations are specifically targeted in our study because Withers \BOthers. (\APACyear2018), using observations from the MAVEN Radio Occultation Science Experiment (ROSE), reported that the peak altitude increased $\sim$10 km between Ls* 305*∘-307∘* during a small, localized dust storm near the south pole. The ROSE observations were obtained at similar SZAs ($\sim$75*∘), but at higher latitudes (53∘*N, Figure 8a).

Figures 8a-8b show that, although the dust optical depths are slightly elevated compared to non-dust season they were constant throughout this period (Figures 2 and 3). Nonetheless, the peak altitude shown in Figure 8c, increases from $\sim$160 km to 175 km between Ls 304*∘-306∘, then falls sharply to 140 km by Ls* 309*∘. The Ls* in which the peak altitude rises coincides with the Ls range reported by Withers \BOthers. (\APACyear2018) for the MAVEN ROSE observations, as marked by the red line in Figure 8c.

Although our results are consistent with Withers \BOthers. (\APACyear2018), it is interesting that the dust optical depths used in our work do not show any dust increase during this time. Figure 4b in Withers \BOthers. (\APACyear2018) clearly shows that a small, localized dust storm near the south pole was observed by the Mars Climate Sounder at Ls 305*∘*. This inconsistency can likely be explained by considering two major differences in the dust data used. First, Withers \BOthers. (\APACyear2018) used MCS dust opacities at the specific pressure level of 50 Pa, while we use MCS estimated column dust optical depths, derived by integrating the MCS dust opacity over all available levels, including over the extrapolated part of the dust opacity profile down to the ground, assuming the dust is well mixed (Montabone \BOthers., \APACyear2015). Second, recent analysis of MY 34 MCS column dust optical depths have highlighted large differences between dayside and nightside values during dust storms, with dayside values generally being larger (Montabone \BOthers., \APACyear2019, submitted to this issue). While these differences are still a topic of ongoing research, the MY 33 reconstructed maps of column dust optical depth used in our work were constructed primarily from nightside MCS observations (Montabone \BOthers., \APACyear2015). Hence, it is possible that the optical depths used in our work are low-biased such that the small, localized dust storm in MY 33 is missing from the maps.

For completeness, Figure 8d shows the MCS dust extinction map at Ls 275*∘-330∘. We produced the map by bin-averaging the MCS dust extinction data from the 50 Pa level (McCleese \BOthers., \APACyear2007) over 3∘* in latitude and 1*∘* in Ls. The map is nearly equivalent to Figure 4b in Withers \BOthers. (\APACyear2018) and confirms the small, localized dust storm near the south pole that started at Ls $\sim$ 305*∘*.

It is also interesting to consider that the MAVEN ROSE observations are from northern latitudes near $+50^{\circ}$ while the MARSIS observations are from southern latitudes near $-30^{\circ}$ (Fig. 8a). Despite the large separation between them, both instruments observe the peak altitude rise at nearly the same time. This implies that the dust storm affected the upper atmosphere over large spatial distances on short timescales. Furthermore, compared to MARSIS, the MAVEN ROSE peak altitude descends more slowly after being elevated by the dust storm (see Figure 4 in (Withers \BOthers., \APACyear2018)). This implies that the effects of the dust storm lasted longer in the northern hemisphere than in the southern hemisphere, despite the dust storm being localized near the south pole.

To summarize, the peak altitude increased $\sim$10 km during the MY 33 local dust storm.

4.6 Mars Year 34

Figure 9 summarizes the MARSIS observations from the MY 34 global dust storm between Ls 180*∘-195∘. The observations are somewhat unfavorable because the MEX periapsis segment covered high SZAs and was rapidly evolving towards the nightside. Consequently, the observational period covers from before the onset of the dust storm (L${}_{s}\simeq 185^{\circ}$) to after it became a global event (L${}_{s}\simeq 190^{\circ}$) (Montabone \BOthers., \APACyear2019, submitted to this issue). It does not cover the entirety of the storm as dust levels remained highly elevated through $\sim$Ls* 240*∘* (Figure 3). Given the constraints of the observations, we slightly modify our analysis for this event.

The MARSIS observations, shown in Figure 8a, are from SZAs 60*∘-85∘* and northern polar latitudes between +59*∘* and +86*∘. The global average dust optical depth, shown in Figure 9b, begins to rapidly increase at Ls* $\sim$185*∘, from $\sim$0.15 before dust storm onset, to 0.6 by Ls* 197*∘. The local average dust optical depth, however, is relatively constant, indicating that the lower atmospheric dust content at northern polar latitudes, where the MARSIS observations are from, did not significantly increase. In Figure 8b we also show a third type of averaged optical depth, which we derive by averaging the dust data from a confined latitude range between +0∘-45∘. This “northern hemisphere” average optical depth is a better indicator of dust at locations near the peak altitude observations. Figure 8b shows that the northern hemisphere dust optical depth increases more rapidly than the local or global dust optical depths, a consequence of the dust storm having originated in the northern hemisphere near +30∘* (Sánchez-Lavega \BOthers., \APACyear2019).

Figure 9c highlights the limitations of the SZA and latitude coverage during this period. The beginning of the period contains measurements from the southernmost latitudes ($<$ 65*∘) that are closest to where the dust storm originated. Also, after Ls* 187*∘, the SZA range is skewed towards higher and higher values with each passing orbit. Since the peak altitude is strongly dependent on SZA near the terminator (Eq 1), we must be careful when analyzing the peak altitudes from this period. Therefore, we calculate orbit-averaged peak altitudes from two different SZA ranges: 60∘-85∘* and 60*∘-70∘*.

The peak altitude at SZAs between 65*∘-85∘, shown in Figure 9d, rises 10-15 km between L$s$ 185∘-188∘. The increase occurs at the onset of the dust storm when the northern hemisphere dust optical depth increases from 0.15 to 0.35. Following its increase, the peak altitude becomes highly variable. Between Ls* 188*∘-190∘*, the peak altitude decreases from 153 ($\pm$ 14) km to 125 ($\pm$ 6) km, then increases to 148 ($\pm$ 6) km the next orbit. Using Equation 5 with $h_{max_{i}}=137$ km, these peak altitude changes correspond to relative pressures of 4.0 ($\pm$ 4.0), 0.4 ($\pm$ 0.2), and 3.0 ($\pm$ 1.0), respectively. The variability metric also increases from $\sigma_{h_{max}}$ = 5 km (L${}_{s}<185^{\circ}$) to $\sigma_{h_{max}}$ = 11 km (L${}_{s}>189^{\circ}$), but this increase is not necessarily meaningful because the metric is being compared across two periods that have vastly different SZA coverage (Figure 9c).

The peak altitude at SZAs between 60*∘-70∘* is shown in Figure 9e. In this case, we require five or more peak altitude observations to calculate an orbit-average (as opposed to 10 or more in Figure 9d). Again, the peak altitude rises 10-15 km between 185*∘-188∘, concurrent with the increase in the northern hemisphere dust optical depth. The two different SZA ranges (Figures 9c-d) provide ample evidence that the peak altitude increased at the onset of the MY 34 global dust storm between Ls* 185*∘-188∘*.

MAVEN ROSE also observed peak altitudes during the MY 34 global dust storm, but at a later time period between Ls 195*∘-270∘* (Felici \BOthers., \APACyear2019, submitted to this issue). The observed peak altitudes in both the northern (+50*∘) and southern (-20∘) hemispheres were somewhat elevated relative to their expected values. However, in the northern hemisphere, ROSE did not observe the peak altitude abruptly increase as one would expect during a dust storm. Felici \BOthers. (\APACyear2019, submitted to this issue) suggested that the peak altitude in the northern hemisphere might have already been elevated before the ROSE observations began. Our results our consistent with this scenario since we observe the peak rise between Ls* 185*∘-188∘, before the first ROSE observation at Ls* 195*∘*.

In summary, the peak altitude increased $\sim$10-15 km during the MY 34 global dust storm, and its variability increased after shortly after it was elevated.

5 Discussion and Conclusions

Table 1 provides a description of the MARSIS observing conditions, the dust storm, and the response of the ionospheric peak altitude for the six case studies that were analyzed in Section 4. Elevated peak altitudes were observed during each of the case studies, although in MY 29 and 32 the increases were less significant because the peak altitudes were highly variable.

The local dust storm in MY 33, the regional dust storm in MY 27, and the global dust storms in MYs 28 and 34 provided the clearest examples of elevated peak altitudes during dust storms. During the MY 27 regional dust storm, the peak altitude sharply increased by 10-15 km over $\sim 5^{\circ}$ of Ls, and then slowly decreased back to pre-storm altitudes over $\sim 15^{\circ}$ of Ls. This rapid rise and slow descent is consistent with previously reported observations of the thermosphere and ionosphere during dust storms (Keating \BOthers., \APACyear1998; Lillis \BOthers., \APACyear2010; England \BBA Lillis, \APACyear2012; Withers \BBA Pratt, \APACyear2013; Zurek \BOthers., \APACyear2017). During the MY 28 global dust storm, the peak altitude also increased by 10-15 km, but the response was more gradual, occurring over $\sim 10^{\circ}$ of Ls. The elevated peak altitudes during these storms were caused by thermospheric pressure levels rising in response to solar heating of dust, and the subsequent expansion of the lower atmosphere (Hantsch \BBA Bauer, \APACyear1990; Haberle \BOthers., \APACyear1993; S\BPBIW. Bougher \BOthers., \APACyear1997; Keating \BOthers., \APACyear1998; Bell \BOthers., \APACyear2007; Heavens \BOthers., \APACyear2011; Withers \BBA Pratt, \APACyear2013; Wolkenberg \BOthers., \APACyear2018; Fang \BOthers., \APACyear2019; Qin \BOthers., \APACyear2019).

During MY 33, the peak altitude increased by 10-15 km at Ls 305*∘* and at southern latitudes between -20*∘* and -40*∘. Withers \BOthers. (\APACyear2018) reported MAVEN ROSE observations showing a similar increase in the peak altitude at the same time, but at northern latitudes near +50∘. They attributed the elevated peak altitude to a small, localized dust storm near the south pole. Since the MARSIS and MAVEN observations were separated by 80∘* in latitude, these concurrent observations suggest that this small dust storm spread quickly and affected the upper atmosphere across large distances.

Another interesting aspect of these two observations is that, after being elevated by the dust storm, the peak altitude observed by MARSIS descended more rapidly than the peak altitude observed by MAVEN ROSE. Hence, even though the dust storm was localized near the south pole, its effect on the upper atmosphere lasted longer in the northern hemisphere than in the southern hemisphere. This surprising result points to the important role of interhemispheric circulation, which allows localized dust storms to affect the upper atmosphere across the planet (Bell \BOthers., \APACyear2007; Withers \BBA Pratt, \APACyear2013). As shown in the model simulations (Bell \BOthers., \APACyear2007; Medvedev \BOthers., \APACyear2013), meridional circulation – which transfers energy from the summer hemisphere to the winter hemisphere – is enhanced by increased dust levels. The MARSIS and MAVEN ROSE observations support this result by showing that a dust storm located near the south pole during southern winter can significantly affect the upper atmosphere in both hemispheres.

Finally, from the six case studies we conclude that dust storms significantly enhance upper atmospheric variability. The clearest example being the MY 28 global dust storm, during which the peak altitude variability metric more than doubled. Substantial variability was also observed throughout the entire MY 29 and MY 32 periods, which covered the peak of the annual dust cycle near Ls 230*∘. Less variability was observed throughout MY 27 near Ls* 315*∘*, when dust levels were significantly lower. During the MY 34 global dust storm, the peak altitude varied by more than 20 km immediately after rising in altitude. Observations of increased upper atmospheric variability during dust storms is not unprecedented. The Mars Global Surveyor’s accelerometer observed a more than 100% increase in the orbit-to-orbit variability of thermospheric densities during the 1997 Noachis regional dust storm (Keating \BOthers., \APACyear1998; S. Bougher \BOthers., \APACyear1999).

The increased variability strongly suggests that dust storms enhance dynamical processes that couple the lower atmosphere to the upper atmosphere, such as upward propagating gravity waves or atmospheric tides (S. Bougher \BOthers., \APACyear1999; Medvedev \BOthers., \APACyear2013; England \BOthers., \APACyear2017). Our results point to the importance of including these processes in global atmospheric models, which can successfully reproduce the neutral density and peak altitude increases observed during dust storms (S. Bougher \BOthers., \APACyear1999; Bell \BOthers., \APACyear2007; González-Galindo \BOthers., \APACyear2010; Medvedev \BOthers., \APACyear2013), but cannot currently reproduce the observed orbit-to-orbit variability. Improvements to these global models will improve our understanding of how dust storms affect the dynamical processes that link the lower and upper atmospheres, and drive substantial thermospheric perturbations.

Acknowledgements.

We thank Joe Groene and Chris Piker for helping us acquire MARSIS data using the das2py software package (https://das2.org/das2py/). We also thank an anonymous reviewer whose suggestions greatly improved the structure and content of the paper. Z. Girazian gratefully acknowledges Paul Withers, Marianna Felici, and Shane Stone for their valuable discussions. This research was supported by NASA through Contract No. 1560641 with the Jet Propulsion Laboratory. F. Němec was supported by MSMT grant LTAUSA17070. The MARSIS, EUV, and dust data compiled for this study can be downloaded at https://www.researchgate.net/publication/335104578_MARSIS_EUV_DUST_DATA and are publicly available. The complete dust maps can be downloaded at www-mars.lmd.jussieu.fr/mars/dust_climatology/. The MCS data used in Figure 8d is publicly available and can be downloaded at https://atmos.nmsu.edu/data_and_services/atmospheres_data/MARS/aerosols.html.