Seasonal and Geographical Variability of the Martian Ionosphere From Mars Express Observations

We study the seasonal and geographical variability of the peak electron density and the altitude of the main ionospheric peak at Mars. For this purpose, we use the data obtained by the ESA Mars Express mission, namely by the radar MARSIS and the radio occultation experiment MaRS. The accumulation of data during the long lifetime of Mars Express provides for the first time an almost complete seasonal and geographical coverage. We first remove the dominant variability factors affecting the main ionospheric peak, namely the effect of changes in the solar zenith angle (SZA), and the changes in the solar ultraviolet radiation output at the Sun. When averaging results obtained at all latitudes, we find that the seasonal variation of both the peak density and the peak altitude can be well reproduced by sinusoidal functions with amplitudes about 8%–9% of the annually averaged peak density, and between 8 and 9.5 km for the peak altitude. We also find elevated peak electron densities in the region of strong crustal fields and latitudinal asymmetries in both the peak density and altitude. Comparing the seasonal evolution of the peak altitude during Mars Year 28, a year with a global dust storm, and the rest of the years, we find that the global dust storm raised the altitude of the ionospheric peak by about 10–15 km.

(e.g., Bougher, Brain, et al. 2017), and produced by the interaction between solar extreme ultraviolet photons and the neutral atmosphere. A secondary peak is usually found below, although often is more a shoulder in the profile than a well-defined peak, and it is produced by soft solar X-rays and secondary ionization (Fox & Yeager, 2009). An sporadic peak has been observed sometimes below the secondary peak, at about 90 km altitude, which was first attributed to the deposition of meteoritic material (Pätzold et al., 2005). This interpretation has however been defied by recent observations (Crismani et al., 2017(Crismani et al., , 2019Peter et al., 2021). We will focus here on the study of the seasonal variability of the two parameters defining the main ionospheric peak: the peak electron density and the peak altitude.
The ionosphere is known to be a very variable region affected by changes in the solar illumination conditions, the changes in the solar radiation and the solar wind, the density of the neutral atmosphere, or the presence of crustal magnetic fields, among other factors (e.g., Bougher, Brain, et al. 2017;Withers, 2009). The ionospheric variability induced by changes in some of these factors, namely the solar illumination and the solar UV radiation output, can be well represented by the Chapman theory.
The dominant variability of the main ionospheric peak is the one produced by the change in the solar illumination defined by the solar zenith angle (SZA). When the SZA increases, the optical depth at a constant altitude increases, and as a consequence the peak electron density decreases and the peak altitude increases. The SZA variability of the peak electron density can be represented by (e.g., Withers, 2009): where N e (SZA) is the peak electron density at SZA, N e (SZA = 0) the peak electron density at SZA = 0, and a an exponent. According to Chapman theory, a = 0.5. Fits of the SZA variability of different datasets have provided values ranging from about 0.45 to 0.55, very close to the theoretical value (Fox & Yeager, 2009;Němec et al., 2011;Peter et al., 2014). When approaching the terminator, the approximation of the cosinus loses validity, and cos(SZA) needs to be replaced by a more sophisticated representation, such as the Chapman function (Smith & Smith, 1972).
Regarding the altitude of the peak, according to Chapman theory: where z max (SZA) is the peak altitude at SZA, and H the neutral atmospheric scale height. Again, when approaching the terminator, the cosinus needs to be replaced by the Chapman function. Analysis of different ionospheric datasets are consistent with that expression (Fox & Weber, 2012;Němec et al., 2011;Peter et al., 2014).
The variability of the ionizing UV solar flux with the 11-year solar cycle and the 27-days solar rotation induces also a variation in the values of the peak electron densities (e.g., Sánchez-Cano et al., 2016;Venkateswara Rao et al., 2014). Usually, this variability is represented by: where F 10.7 is a solar activity proxy index (measured at the Earth but usually corrected for the Sun-Mars distance in studies about Mars ionosphere) and c a constant (e.g., Morgan et al., 2008;Němec et al., 2011).
The peak altitude is to first order determined by the temperature structure in the middle atmosphere (Zou et al., 2005), where UV solar photons do not penetrate, and it is thus expected not to exhibit strong changes with the solar UV radiation output. However, Fox and Weber (2012) showed that the altitude of the main ionospheric peak showed a negative slope (quite small, though) with increasing solar activity, due probably to the penetration of X-ray energetic photons to deeper layers.
Mars has a very eccentric orbit, and thus the amount of solar radiation getting to the planet changes significantly with season. So, the Martian ionosphere is expected to present a seasonal variation. Previous studies have identified signatures of this seasonal variability. Zou et al. (2005) found, by analyzing Mars Global Surveyor (MGS) Radio Science (RS) data obtained in the high latitudes of both hemispheres at L s ≈ 140°, that the peak altitudes in the Southern (winter-spring) hemisphere were 7.7 km lower than those in the Northern (summer-fall) hemisphere. A later analysis of MGS RS data obtained in the seasonal range L s = 80-180° showed an increase in the peak altitude (Zou et al., 2011). Morgan et al. (2008) fitted Chapman layers to MARSIS Active Ionospheric Sounding (AIS) mode electron density profiles obtained at different seasons and found significant differences in their derived values of N e (SZA = 0) and z max (SZA = 0) at northern summer and southern summer solstices. Němec et al. (2011) found that the MARSIS AIS peak altitude varied with the Sun-Mars distance. Similarly, Sánchez-Cano et al. (2013) also found an increase in the MARSIS AIS peak electron densities from the Northern Summer solstice to the Southern Summer solstice, but not in the peak altitudes. It is important to note that in some cases the F 10.7 solar proxy index at the actual Sun-Mars distance is used to study the variability of the Martian ionosphere due to changes in the solar UV radiation. The variation of the F 10.7 index at Mars mixes then the changes in the solar radiation output and the variability due to the Sun-Mars changing distance.
In addition to the direct effects of the Sun-Mars distance, Sánchez-Cano et al. (2018) have found an increase in the Total Electron Content between L s = 30-70°, attributed to the seasonal cycle of CO 2 condensation-sublimation in the polar caps.
In this study, we use the observations of the peak ionosphere region by Mars Express to improve the current characterization of the seasonal variability of the Martian ionosphere, isolating it from the changes produced by variations in the SZA and in the solar UV radiation output. Our main advantage with respect to previous studies is that the accumulation of data obtained during the more than 15 (terrestrial) years of life of Mars Express offers for the first time an almost complete seasonal and geographical coverage. We will also take advantage of this good coverage to try to identify possible geographical variations, an aspect less thoroughly studied of the ionosphere. For this purpose, we remove first the dominant variability with SZA and that due to the solar UV radiation change at the Sun. This will also allow us to compare the obtained SZA and solar UV variations with those obtained previously with other smaller datasets. Note that, differently to previous studies, our data set covers more than a full 11-year solar cycle.
Section 2 describes the data sets used in this study. Section 3 shows the results obtained. A discussion about the geographic variations, the effects of a global dust storm, and the sensitivity of the results to some of the assumptions made is presented in Section 4. Finally, a summary of the main results is presented in Section 5

Data
This paper uses data from MARSIS and MaRS on board Mars Express.
The Mars Express Radio Science experiment (MaRS) performs radio occultations at two different frequencies (X-band and S-band) to study the neutral lower atmosphere and the ionosphere . MaRS data obtained until mid-2017 are used here. Details on the processing of the data can be found in Peter et al. (2014). Typical uncertainties in the peak magnitudes are between about ±1% and 5% for the peak density, with an average value of less than ±2%, and between ±1 and 5 km for the peak altitude, with average value around ±2.5 km.
Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS) is a radar designed to investigate the surface and subsurface of Mars and the Martian ionosphere (Orosei et al., 2015). Due to delays in its radar deployment, MARSIS only started operating in August 2005. We use in this analysis the electron density profiles obtained when MARSIS operates in the Active Ionospheric Sounding mode (AIS). In this mode, the instrument successively transmits short pulses at 160 quasilogarithmically spaced frequencies f(δf/f ≈ 2%) from 0.1 to 5.5 MHz. The time delay between the pulse transmission and the detection of its echo reflected vertically from the ionosphere is measured. Additionally, the local electron density can be readily determined from the local plasma oscillations, which allows to calculate electron density profiles spanning from the spacecraft altitude down to the altitude of peak electron density (Morgan et al., 2008(Morgan et al., , 2013. However, due to signal-to-noise limitations, the lowest sounding frequencies (f < ≈0.5 MHz) typically do not produce a detectable echo. This results in an effective gap in the time delay measurements between the local plasma frequency and the lowest frequency available from the remote ionospheric sounding and a need to assume an empirical profile shape in the corresponding range of altitudes. We use a realistic profile shape and a technique described by Němec et al. (2016Němec et al. ( , 2017 to maximize the accuracy of the calculated electron density profiles. Moreover, only the electron density profiles obtained at spacecraft altitudes lower than 500 km, which are deemed to be the most precise as they minimize the altitudinal extent where the interpolation is needed, were used. The expected uncertainty in the peak altitude estimation is lower than 7 km (Němec et al., 2016). Note that the peak electron density can be determined directly from the maximum frequency reflected from the ionosphere, with no need to calculate the electron density profile. The uncertainty of the peak electron density determination is typically lower than about 2% (e.g., Němec et al., 2019). Only data until October 2015 have been inverted to electron density profiles.
For both datasets, we limit our study to observations performed in the dayside, in particular SZA ≤ 88°. In total, we use 278 observations from MaRS and 34,260 from MARSIS. The distribution of the observations used in this study can be seen in Figure 1, where the data obtained at different Mars Years (MY) are represented in different colors. It can be seen that the seasonal and latitudinal coverage is quite good, except in the period L s 100°-180°, when there is a gap in MARSIS observations. This is due to the characteristics of Mars Express orbit. Usually there are eclipse periods at this time, where the instrument is not working. Also Mars Express orbit is close to the poles at this time, where priority is given to observations in the MARSIS SubSurface (SS) mode. In addition, in some of the MYs Mars Express orbit is at the nightside during this particular period, while as stated above we only focus here on the dayside. are represented by the green, blue, and red lines, respectively. Note that MARSIS only obtains measurements above the main ionospheric peak, while MaRS data extend to the lower ionosphere. For each of the measured profiles, we extract the electron density at the peak and its altitude. The peak electron density for all observations with SZA ≤ 88° is shown in the bottom left panel of Figure 3. The anticorrelation with the SZA values is clearly seen at first sight. The corresponding peak altitudes are shown in the bottom right panel of Figure 3. The links with SZA are less clearly seen, although can still be appreciated. For example, peak altitudes in the Southern mid latitudes around L s ∼ 300° are significantly higher for MY29 and MY30, when SZA was larger than 80°, than for MY28, when the SZA was below 50°. It can already be appreciated that the observations in the first half of the Martian year tend to have a lower peak altitude than those obtained around perihelion, a first hint of seasonal variability.
In what follows, we will isolate the seasonal variability of both the peak electron density and the peak altitude. For this purpose, we first remove the dominant variability factors, those due to the changes in SZA and the solar UV radiation output at the Sun. We limit our study to dayside data, in particular to those with SZA ≤ 88°. For the nightside ionosphere, the expressions used to describe the SZA variability of the ionosphere are no longer valid and it is thus not possible to remove it and isolate the seasonal variability.
GONZÁLEZ-GALINDO ET AL.    Figure 4 shows the SZA variability of the peak electron densities measured by both instruments, displaying the well-known decrease of the densities with increasing SZA. Note that the apparent stripes at discrete peak electron densities are due to the discrete nature of the MARSIS radar sounding frequencies, as described above.

Peak Electron Density
In order to remove the SZA variability, we calculate for each observation the peak electron density at SZA = 0 using expression (1), but replacing the inverse of the cosinus by the Chapman function, in order to get a better representation of the behavior close to the terminator. We also fix the exponent a to its theoretical value of 0.5. Namely: For the Chapman function, we use the first term of the asymptotic expansion provided by Huestis (2001).
Comparisons with more rigorous (and much more computationally expensive) implementations of the Chapman function shows that the differences for SZA = 88 are about 0.1%. Figure 4 (bottom panel) shows the SZA-corrected peak electron densities as a function of SZA. It is clear that the decrease of the peak electron density with increasing SZA has been corrected, with the subsolar peak densities clustering around 1.5⋅10 5 cm −3 for all SZA. However, there is still significant variability. For example, it seems that data for MY32 present larger subsolar peak densities than for the rest of the years.
10.1029/2020JE006661 5 of 15  Figure 5 shows the SZA-corrected peak electron densities as a function of the solar radiation output, represented by the F 10.7 value at 1 AU. A clear increase of the SZA-corrected peak electron densities with the solar radiation output is observed. This explains the larger subsolar peak densities obtained during MY32, when solar activity was higher than for the rest of the years. In order to remove this variability, we fit the solar variability to expression (3). The fit is represented by the gray line in Figure 5. The value for the fitted exponent is 0.324 ± 0.002. This value is of the same order than determined from previous analysis of MARSIS AIS data (e.g., 0.30 ± 0.04 for Morgan et al. (2008), 0.388 ± 0.003 for Němec et al. (2011)) and other datasets (0.36 for Hantsch and Bauer (1990) using Mariner and Viking data, or 0.37 ± 0.06 for Breus et al. (2004) using MGS data).
In order to remove the variability of the UV solar radiation output from the SZA-corrected peak electron densities, we calculate the equivalent peak electron density at F 10.7 = 100. The peak electron densities obtained after applying the SZA and solar activity corrections can be seen, as a function of latitude and season, in Figure 6. In contrast with the behavior of the uncorrected data shown in Figure 3, the corrected peak electron densities show a clear seasonal contrast, with lower values in the first half of the year and higher values around the perihelion season.
In order to have a better insight into this seasonal variability, we bin the results in intervals of 5° in L s and average the SZA and solar radiation output corrected peak electron densities in these bins. The binned densities are shown as the black line and error bars in Figure 7. Peak electron densities are found to be highest close to the perihelion season and lowest around the date of aphelion. We fit a sinusoidal function (shown as the gray dashed line in Figure 7) to the seasonal variability of the binned electron densities. The result of the fit is: 14300 800 25 6 (5) The minimum/maximum of the fitted function is obtained at L s = 65/245°, very close to the aphelion/perihelion value at L s = 71/251°. The amplitude of the seasonal variability is about 9% of the yearly averaged value, so that from aphelion to perihelion the peak electron density increases a bit less than 20% due to the eccentricity of the Martian orbit. This rate of increase corresponds well to the 20% change in the Sun-Mars distance, indicating that the peak electron density varies linearly with the inverse of the Sun-Mars distance, in agreement with what is expected from Chapman theory (Mendillo et al., 2003) and from previous observational studies (Sánchez-Cano et al., 2013). Simulations with a Global Climate Model (González-Galindo et al., 2013) have also shown a similar, though slightly lower (≈16%), enhancement in the peak electron densities from aphelion to perihelion.
Note that, behind the seasonal variation described by the sinusoidal function, there are still variations of the peak electron density, such as the increase between L s = 0° and 20° or the oscillations in the L s ≈ 210°-300° range. These variations can be due to factors different to the variation with SZA, solar UV radiation output, and season analyzed until now, known to produce ionospheric variability. For example, variations in the electron temperature (with season, latitude, etc.) can affect the electron densities via the electron-ion chemical recombination rates. Tides and waves are also known to produce variations in the ionospheric parameters (Bougher et al., 2001(Bougher et al., , 2004Mahajan et al., 2007). The presence of crustal magnetic fields in some portions of the Martian surface can also be a factor producing ionospheric variability (Krymskii et al., 2003;Nielsen et al., 2007). These factors can produce a geographical variability of the peak electron densities.
GONZÁLEZ-GALINDO ET AL.  Given that the data coverage is not uniform and changes with season, potential geographical variations can translate into temporal variations such as the ones displayed in Figure 7.
We also note that we do not see any significant increase of the corrected peak electron densities around L s = 30°-70°, as it was observed for the TEC (Sánchez-Cano et al., 2018). This may mean that the increase of the TEC is not due to an increase in the densities at the region of the main peak, but in higher or lower ionospheric regions.

Peak Altitude
We move now to the study of the seasonal variability of the peak altitude. The top panel of Figure 8 shows the SZA variability of the peak altitudes measured by both Mars Express instruments. A very large variability in the peak altitude, from 100 to 200 km, is found. This variability is larger than previously found by the MGS RS data set, from which peak altitudes range between 120 and 165 km. Note however that our data set contains 10 times more observations, and that the MGS RS data set was obtained in quite restricted latitudinal and seasonal corridors, while our data set covers most seasons and latitudes. Despite the strong variability, a tendency to a higher peak altitude can be identified when approaching the terminator. Note also that the variability of the peak altitude tends to increase when approaching the terminator, so that for SZA ≈ 85° the measured peak altitudes range from as low as 100 km to higher than 180 km. Note that early ionospheric measurements by the Mariner and Viking missions also displayed a similar strong variability of the peak altitude close to the terminator, for example, see Figure 3 in Hantsch and Bauer (1990) where peak altitudes around SZA ≈ 80-85° vary between about 115 and 180 km.
As discussed in Section 1, Chapman theory indicates that the SZA variability of the peak altitude can be represented by Equation 2. The SZA variability depends on the scale height, H, that is, on the temperature in the lower thermosphere region. This temperature is known to be variable and nonuniform, and in particular should vary with season and with the 11-year solar cycle (Bougher, Roeten, et al., 2017;Forget et al., 2009;González-Galindo et al., 2009. For this reason, we have not applied a single fit to the whole data set. Instead, we have divided the data set in 12 subsets, one corresponding to each half Mars Year, with the exception of the first half of MY27 and the second half of MY33, where very few data exist and have been merged with the rest of MY27 and MY33, respectively. We will later discuss the sensitivity of the obtained seasonal variability to the assumptions made for the SZA variability. However, caution needs to be taken when interpreting the obtained values for the scale height. First, because the temperature is still expected to present an important variability during half Mars year. Second, because of the inherent limitations of Chapman theory. A good discussion about these limitations and about the caveats of identifying the obtained scale heights with temperatures can be found in Fox and Weber (2012), who found scale heights as low as 3 km when analyzing the MGS radio science measured peak altitudes.
We have then fitted Equation 2 independently to each of the 12 subsets, replacing the cosinus of the SZA by the inverse of the Chapman function.
The fitted values are summarized in Table 1  (see, e.g., Withers, 2009). However, for two of the datasets the values of the fitted scale height are significantly lower, 2.8 and 1.3 km. For the first half of MY31 the very low value of the fitted scale height is probably due to the fact that the SZA coverage is quite restricted, only covering SZA > 70. For the second half of MY28, the results may be affected by the well-known presence of a global dust storm . This is suggested also by the value of the fitted subsolar peak altitude, about 10 km larger than the average for the other datasets spanning the second half of the Mars year. While higher temperatures (and thus higher scale height) could be expected due to the dust storm, the thermospheric response to the dust storm is not directly driven by a local heating, but rather by the modification of the circulation due to the heating of the lower atmosphere, and areas/times of lower thermospheric temperatures are also possible due to the dust storm Medvedev et al., 2013). Note also that, as discussed before, we made a single fit for all the second half of MY28, including data obtained before and during the storm, which could also affect the results. Other apparently counter-intuitive result is the higher scale height during the first half of MY29 than during the second half. All these results caution against the direct identification of the scale heights obtained from the SZA variation of the peak altitude with thermospheric temperatures, as discussed by Fox and Weber (2012).
Then we apply Equation 2 using the fitted scale heights to calculate, for each observation, the equivalent peak altitude at SZA = 0. This is shown, as a function of SZA, in the bottom panel of Figure 8. It can be seen that there is still a considerable amount of scatter in the corrected peak altitudes. As discussed in Section 1, the peak altitudes are not expected to vary with the solar UV radiation output. In order to confirm this, and to check if some of the obtained scatter can be caused by the effect of changing solar UV radiation output, we plot in Figure 9 the variation of the SZA-corrected peak altitudes as a function of the solar radiation output, represented by the F 10.7 solar proxy index at 1 AU. No clear indication of a solar activity effect is seen. The Spearman's rank correlation has been calculated, and the correlation coefficient obtained is −0.134, indicating a very weak negative correlation between the SZA-corrected peak altitudes and the solar UV radiation output, in agreement with previous findings by Fox and Weber (2012). Figure 10 displays the latitudinal and seasonal variability of the SZA-corrected peak altitudes. It can be clearly seen that the peak altitudes during the second half of the year are larger than those in the first half. As for the peak density, we have binned the peak altitudes in intervals of 5° of L s . The result is shown in Figure 11 as the black solid line. Also shown as a gray dashed line is the result of a sinusoidal fit to that variation: So, minimum/maximum peak altitudes are found at L s = 74/254°, in agreement within the uncertainty to the aphelion/perihelion timing at L s = 71/251°. A total increase of about 16 km in the peak altitude from aphelion to perihelion is obtained, in line with the result of a previous analysis of a MARSIS reduced data set (Němec et al., 2011). Note that this result is also in excellent agreement with the prediction of a Global Climate Model of a 16 km seasonal variation of the subsolar ionospheric peak altitude (González-Galindo et al., 2013), due to the inflation-con-GONZÁLEZ-GALINDO ET AL.  traction of the atmosphere produced by the changes in the temperature structure in the lower atmosphere.
Besides the seasonal variability represented by the sinusoidal function, there are still some variations, as in the case of the electron density. The most notable one is the decrease in the peak altitude around L s = 90°. This is due to a period of very low peak altitudes consistently observed in MY31 during the L s = 92-97 range at all the observed latitudinal range (15S-50S). We have checked that these observations with very low peak altitudes have a distribution in SZA, latitude, longitude, and solar radiation output similar to the other observations in the first half of MY31. We have also checked that there is no influence on the spacecraft altitude on these observations. The temperatures measured in the lower atmosphere by PFS during MY31 show no particular or unusual feature at this L s /latitude range (Giuranna et al., 2021). So, we currently have no explanation for this period of very low peak altitudes. A dedicated in-depth study may be needed to clarify this issue.

Geographic Variations
The unprecedented good seasonal and latitudinal coverage of the Mars Express datasets allows the exploration of the possibility of latitudinal differences in the seasonal behavior of the ionosphere. For example, Figure 6 suggests that corrected peak densities are larger in the Northern hemisphere than in the Southern hemisphere around L s 210°-270°. We have binned the SZA, F 10.7 corrected peak electron densities in bins of 5° of L s and in three latitudinal ranges: low latitudes (Lat < 30) and mid latitudes (30 < Lat < 60) in both hemispheres. The results are represented by the blue, orange, and green lines and error bars in Figure 7. There are hints of slightly larger densities in the corresponding winter hemispheres (green line during the first half of the year, orange line during the second half of the year). However, the obtained variability for the summer/winter hemispheres are within the obtained error bars for the average of all latitudes. We have performed a Kolmogorov-Smirnov test in each L s bin to check the statistical significance of the differences between the mid latitudes peak electron density distributions in both hemispheres. With the exception of the bins L s = 15°-20° and L s = 75°-80°, the obtained p-values are always lower than 2% (and in most of the bins lower than 0.1%), meaning that the differences between the mean densities in both hemispheres are statistically significant.
What potential physical mechanisms could be at the origin of this hemispherical asymmetry? One possibility is differences in the electron temperatures between the summer and the winter hemisphere. Given that the rate of electron recombination with  2 O , the dominant loss mechanism for the dayside ionosphere, depends inversely on the electron temperatures, a direct relation between electron densities and electron temperatures is to be expected. Thus, larger electron densities in the winter hemisphere would imply larger electron temperatures there, which is somehow counter-intuitive. However, the variation of neutral temperature in the mesopause/lower thermosphere region, close to the altitude of the ionospheric peak, is affected by processes such as dynamics and variations in the amount of atomic oxygen, which can produce departures from a purely radiative-dominated situation Medvedev et al., 2015). Also, MAVEN measurements (Ergun et al., 2015;Vogt et al., 2017) have found an anticorrelation between electron temperatures and electron densities, instead of the expected direct correlation. Other variations in the background neutral atmosphere, such as changes in the CO 2 density, can also produce variations in the electron densities GONZÁLEZ-GALINDO ET AL.   Figure 9. Variability of the SZA-corrected peak altitudes with the F 10.7 solar proxy index at 1 AU, different colors represent data obtained at different MYs. MY, Mars Years; SZA, Solar Zenith Angle. (Mendillo et al., 2017). But again, larger CO 2 densities are expected in the summer hemisphere. A dedicated comparison with a GCM is clearly needed to get more insight into the potential origin of this hemispheric asymmetry.
As for the peak electron density, we explore possible latitudinal differences in the seasonal variability of the ionospheric peak altitude. The average peak altitude in different latitudinal bins (from −30° to 30°, from 30N to 60N, and from 30S to 60S) is shown as blue, orange, and green lines, respectively, in Figure 11. At the solstices, peak altitudes tend to be higher in the summer hemisphere. A Kolmogorov-Smirnov test in each L s bin to the peak altitudes between 30N and 60N and between 30S and 60S produces p-values lower than 1% in all bins, except between L s = 255°-265°. This indicates that the latitudinal differences are statistically significant. A similar hemispheric asymmetry in the peak altitudes has been found in the analysis of MGS observations (Zou et al., 2005), with peak altitudes around L s = 140° 7.7 km higher on average on the summer hemisphere. This asymmetry is due to the expected differences in the temperatures in the lower atmosphere in the summer and the winter hemispheres. Note that in the second half of the year another factor can contribute to establish these latitudinal differences. The presence of the global dust storm in the second half of MY28 can raise the peak altitudes in the Southern hemisphere, where the measurements for this period are located. We discuss in more details the effects of the MY28 global dust storm below.
In order to gain more insight into possible geographical variations in the peak electron density and the peak altitude, we correct for the seasonal variability of both magnitudes, already corrected for SZA and solar radiation output, by substracting to each observation the sinusoidal term in Equations 5 and 6, respectively. Then we bin the resulting magnitudes in bins of 5° in latitude and 5° in longitude. For the peak altitude, to avoid the effects of the global dust storm in MY28, we exclude the results for this MY. The resulting average peak electron densities and peak altitudes in each bin are shown in Figure 12.
Elevated peak electron densities are found in the mid and high latitudes of the Southern hemisphere in the longitude sector around lon = 180°. This is in good agreement with the location of the strongest crustal fields in the Martian surface (Connerney et al., 2005). Previous works have found that the presence of crustal magnetic fields affects the ionosphere. Oblique echoes in the MARSIS data set have been interpreted as a result of an inflated topside ionosphere over regions of strong crustal field Gurnett et al., 2005). Focusing on the ionospheric peak region, (Krymskii et al., 2003) found in the MGS RS data set enhanced peak electron densities within mini-magnetospheres formed in the crustal fields region. Nielsen et al. (2007) found significant increases in the peak electron densities over regions of strong and vertical crustal field in a set of three MARSIS orbits, not linked to increases in the solar flux or the particle precipitation rate. This increase was attributed to a heating of the plasma, produced by the excitation of plasma waves due to the interaction of the solar wind with the crustal field, and leading to reduced electron recombination rate. However, the analysis of MAVEN observations during deep dip campaigns did not show any effect of the crustal fields on the peak density , which could be due to the small size of the data set used in that study. Our result, obtained when analyzing a large data set, confirms the presence of the density increases over the regions of crustal fields. The long-term nature of our data set points to an stable mechanism at the origin of this increase, instead of transient phenomena such as solar flares or episodes of enhanced particle precipitation, in agreement with the result in Nielsen et al. (2007).
Note that other geographic variations in the peak electron density, apart from the increase in the crustal field region, are also present. For exam-  ple, the low latitudes of the Southern hemisphere display an alternance in longitude of regions with relatively high and relatively low peak electron densities. This can be the signature of nonmigrating tides (Mahajan et al., 2007). A detailed analysis of the longitudinal variations and the potential effects of tides and waves in this data set is left for a future study.
Regarding the peak altitude, we see no apparent effect of the crustal fields. Previous works have found oblique echoes in MARSIS observations over regions of strong vertical crustal field, attributed to an inflation of the ionosphere (Gurnett et al., 2005). Our result suggests that this process does not affect the ionospheric peak region, and thus must be acting only in the topside ionosphere. Apart from potential effects of the crustal fields, an increase of the peak altitude toward the high latitudes of both hemispheres is obtained. This could be due to the presence of polar warmings and temperature inversions in the high latitude middle atmosphere (Forget et al., 1999;McCleese et al., 2010). As for the peak density, longitudinal variations are found, particularly in the low latitude region, which again can be the signature of nonmigrating tides.

Effects of the MY28 Global Dust Storm
As discussed above, the presence of a global dust storm in MY28 may be affecting some of the obtained results, in particular producing a subsolar ionospheric peak altitude significantly higher (10 km) than usual for the season. A previous, independent analysis of the MARSIS AIS data set has concluded that the presence of the MY28 global dust storm raised the ionospheric peak in about 10-15 km and induced also significant variability of the peak altitude (Girazian et al., 2020). To distinguish between the variability due to the global dust storm and that expected from the seasonal evolution of the ionosphere, we represent in Figure 13 the latitudinally averaged seasonal variability of the SZA-corrected peak altitudes for MY28 and for the average of all the other years covered by Mars Express observations. As before, we bin the results in bins of 5° of L s .
Before the onset of the dust storm at around L s = 265° the peak altitudes are slightly larger (though within the standard deviation) for MY28 than for the average of the other years. This may be a result of the latitudinal asymmetry in the peak altitudes described above, as observations for MY28 were collected in the mid and high latitudes of the Southern (summer) hemisphere, with most other observations at that season obtained in the Northern (winter) hemisphere. After L s = 265° the MY28 peak altitudes display clearly a different seasonal variability to that of the rest of the years. While the peak altitude averaged for all MYs except MY28 decreases with time, the MY28 peak altitude sharply rises by about 10 km and remains elevated until the end of the measurement period around L s = 300°, so that by L s = 280° the MY28 peak altitudes are about 15 km higher than those for other MYs. This altitude increase is consistent with the one obtained by a different, independent analysis of the same data set (Girazian et al., 2020). While the data set used is the same in both studies, different methods are used to examine the effects of the storm. In particular, while we use the SZA-corrected peak altitudes, Girazian et al. (2020)   ing the recent MY34 dust storm (Felici et al., 2020). Simulations with a multifluid magnetohydrodynamic model show also an increase of 15 km in the peak altitude during the global dust storm in 1971-1972(Fang et al., 2020. The reason of this altitude increase is the expansion of the atmosphere due to the strong heating in the lower atmosphere produced by the storm.

Sensitivity to the Assumptions About the SZA Variability
As mentioned in Section 3.2, the expected seasonal and solar cycle variability of the scale height has made us to perform different fits of the SZA variability of the peak altitude for each half Mars year. Here we test if the obtained seasonal variability is affected by this choice. We have tested three other different fitting strategies to correct for the SZA variability of the peak altitude: a single fit for the whole data set (fit 1); one fit for all the data in the first half of the year (all MYs together) and another fit for all the data in the second half of the year (fit 2); and one fit for all the data in each MY (fit 3). We have then subtracted the SZA effect to the observed peak altitudes using each of these fits, and repeated the sinusoidal fit to the obtained seasonal variability of the SZA-corrected peak altitudes. The obtained fits obtained in each case are: While the assumption made for the SZA variability does not change the big picture of the obtained seasonal variability, the precise numbers are sensitive to it. In particular, the amplitude of the fitted seasonal variability can change by about 1.5 km, that is, about 20% of its value, depending on the strategy used to fit the SZA variation of the peak altitude.

Summary and Conclusions
We have used the data provided by two instruments on board the Mars Express mission to isolate the seasonal variability of the main ionospheric peak induced by the eccentricity of the Martian orbit. The main conclusions we obtain are: 1. The seasonal variability of the peak electron density can be reproduced by a sinusoidal function with an amplitude of about 8%-9% of the annually averaged peak electron density. This amplitude is in line with the expectations of Chapman theory (Mendillo et al., 2003) and with predictions by a Global Climate Model (González-Galindo et al., 2013) 2. We do not find an increase of the peak electron density around L s = 30°-60°, implying that the increase in the TEC in Sánchez-Cano et al. (2018) is produced higher or lower in the ionosphere 3. The seasonal variability of the peak altitude can also be reproduced by a sinusoidal function with an amplitude between 8 and 9.5 km, depending on the assumptions made to remove the SZA variability. This amplitude is in excellent agreement with the seasonal variability of the peak altitude predicted by a Global Climate Model and originated by the changes in the lower atmospheric temperatures (González-Galindo et al., 2013) 4. At solstices, we find latitudinal differences in the peak electron densities and the peak altitudes, with larger electron densities and lower peak altitudes in the winter hemispheres. A similar asymmetry in the peak altitude, due to the temperature contrast in the lower atmosphere in both hemispheres, was found in the MGS RS data set (Zou et al., 2005). Potential explanations for the latitudinal variation of peak densities are variations in the electron temperatures and/or the background neutral atmosphere, but a dedicated study with a GCM is needed 5. We find that the peak densities, once corrected for the SZA, solar radiation output and seasonal variations, are larger over the crustal field regions, confirming previous results obtained with much smaller datasets (Nielsen et al., 2007) 6. The global dust storm in MY28 raised the altitude of the ionospheric peak by about 10-15 km, in agreement with a previous independent analysis of the MARSIS AIS data set (Girazian et al., 2020), with studies of the ionosphere during the MY34 dust storm (Felici et al., 2020), and with modeling results (Fang et al., 2020) Data Availability Statement MARSIS and MaRS data are available in the ESA PSA archive (ftp://psa.esac.esa.int/pub/mirror/MARS-EX-PRESS/MARSIS/ and ftp://psa.esac.esa.int/pub/mirror/MARS-EXPRESS/MRS/). All the data used in this study, including also derived magnitudes, have been made publicly available at Zenodo (https://doi. org/10.5281/zenodo.4320793).