The Martian Bow Shock Over Solar Cycle 23–24 as Observed by the Mars Express Mission

The Martian bow shock position is known to be correlated with solar extreme ultraviolet irradiance. Since this parameter is also correlated with the evolution of the solar cycle, it is expected that the Martian bow shock position should also vary over such a period. However, previous reports on this topic have often proved contradictory. Using 13 years of observations of the Martian bow shock by the Mars Express mission over the period 2004 to 2017, we report that the Martian bow shock position does vary over the solar cycle. Over this period, our analysis shows the bow shock position to increase on average by 7% between the solar minimum and maximum phases of solar cycle 23–24, which could be even larger for more extreme previous solar cycles. We show that both annual and solar cycle variations play major roles in the location of the bow shock at Mars.


Introduction
The Martian bow shock develops to slow the supermagnetosonic flowing solar wind to subsonic speeds such that it can be diverted about the Martian plasma system. Unlike celestial bodies that are enclosed within an intrinsic global magnetic field (e.g., the magnetosphere of Earth), Mars' ionosphere and extended exosphere instead acts as an obstacle to the solar wind flow. In this sense, Mars and Venus are very similar planets, sharing an ionosphere as well as an exosphere extending outside the bow shock and so, contributing to the obstacle to the solar wind plasma flow.
The location of the bow shock depends on several drivers that can be either external, such as the solar extreme ultraviolet (EUV) flux, the solar cycle phases, the solar wind dynamic pressure, the interplanetary magnetic field (IMF) orientation, or the magnetosonic Mach number; or internal, such as the crustal magnetic field sources location rotating with the planet. Both types of drivers produce different types of bow shock variability, which can manifest on short and long timescales. More details can be found in, for example, Mazelle et al. (2004) or Hall et al. (2016).
Focusing on the EUV parameter as a driver, the solar cycle variation is well known to play a major role in the location of the bow shock at Venus. Alexander and Russell (1985) indicated that this behavior was caused by neutral atmosphere variations with the solar cycle and its consequent effect on the mass loading of the solar wind. However, solar cycle variations in the Martian bow shock location have yielded contradictory results to date. The reason could be the different relative distances to the Sun, as previous missions made measurements at different radial distances and levels of solar EUV flux.
Recent studies have found that the bow shock position at Mars is sensitive to EUV radiation (e.g., Edberg et al., 2009;Hall et al., 2016;Halekas et al., 2017;Gruesbeck et al., 2018). Therefore, we may initially expect some form of solar cycle period variability in the average bow shock position. Solar cycle studies of the Martian bow shock have been attempted by comparing the observations of bow shock location by the variety of the older missions (see a summary in Table 2). The Mariner-4 (pre-maximum phase of cycle 20) and Mars-2,3,5 (post-maximum phase of cycle 20) missions recorded a total of 24 bow shock crossings which occurred during a period of moderate solar activity with an average sunspot number, SSN, of 59 (Slavin et al., 1991). The Phobos-2 mission recorded up to 126 bow shock crossings (e.g., Trotignon et al., 1993), 94 reported by Slavin et al. (1991) occurring during a period of higher solar activity (near the maximum phase of cycle 22)   (Slavin et al., 1991). By comparing statistically computed model shock surfaces of the observations at each phase of the solar cycle, Slavin et al. (1991) reported no apparent solar cycle variation in the bow shock position. Modolo et al. (2006) reported a similar finding based on hybrid simulations at the maximum and minimum of solar activity. However, using very similar data sets and methods, Russell et al. (1992) and Trotignon et al. (1993) reported the opposite. Vignes et al. (2000Vignes et al. ( , 2002 extended the previous work by including the 450 bow shock crossings recorded by the MGS mission. These observations occurred during a period where the SSN ranged between 30 and 90 (premaximum phase of cycle 23), and as with the Slavin et al. (1991) study, Vignes et al. (2000Vignes et al. ( , 2002 also reported no obvious solar cycle variability in the Martian bow shock position. In order to determine whether the Martian bow shock position varies with the solar cycle, we require nearcontinuous observations of the bow shock over a solar cycle period by the same instrumentation. The MEX mission (Chicarro et al., 2004) has been sampling the Martian plasma system since December 2003 in an elliptical (~10,000-km apoapsis,~350-km periapsis), 7-hr period orbit that precesses around Mars with time.
Moreover, due to the MEX trajectory, the Martian dawn hemisphere at apoapsis altitudes is sampled more than the dusk hemisphere. Therefore, there are much more bow shock crossings over the dawn region. Since this orbit typically results in MEX crossing the Martian bow shock several times a day and MEX carries appropriate plasma instrumentation to detect this boundary, it is the appropriate mission for studying the entire solar cycle variability of this aspect of the Martian plasma system.
The objective of this paper is to assess for the first time with the same data set how the Martian bow shock varies across the solar cycle and to what extent. This is an important question because we know that the Martian plasma system reacts differently for similar levels of EUV flux at different phases of the solar cycle (e.g., Sánchez-Cano et al., 2015. We present an investigation using the entire MEX data set and present an analysis of the bow shock position in relation to annual EUV variations, as well as different solar activity phases of the solar cycle. We do not focus on the dynamic pressure parameter, because using the same data set of this study, Hall et al. (2016) recently reported that the bow shock location is more sensitive to variations in the solar EUV irradiance than to solar wind dynamic pressure variations.
The remainder of this paper will detail the data sets and methods used (section 3) to determine whether or not the Martian bow shock position does vary over the period of solar cycle 23-24 (sections 4-5).

Data Set and Analysis
The main data set used in this study originates from the Hall et al. (2016) report on the annual (Martian year) variation in the Martian bow shock position. Hall et al. (2016) devised an automatic algorithm to scan the full (at the time) MEX Analyzer of Space Plasmas and Energetic Atoms Electron Spectrometer (ASPERA-3 ELS; Barabash et al., 2006) data set for bow shock crossings. Across the period of January 2004 to May 2015, the authors identified 12,091 bow shock crossings. In addition to this Hall et al. (2016) bow shock crossing list, henceforth referred to as the HALL16 data set, we incorporate additional recently released MEX ASPERA-3 ELS observations to extend the HALL16 crossing list to the end of 2017 (and a total of 13,585 crossings). The additional bow shock observations from this extended data set were extracted using the HALL16 automated identification algorithm. We note that only MEX data have been used in this study in order to avoid wrong statistical comparisons with data sets from other missions that were not crosscalibrated with respect to MEX observations, such as the crossings shown in Table 2 and the different data sets from the Mars Atmosphere and Volatile EvolutioN (MAVEN) mission. Moreover, all of these missions only cover a small portion of the solar cycle. Therefore, in order to perform a coherent statistical study, only observations from the same data set are compared.
To describe the different phases (or levels of activity) of the Sun over the solar cycle, three supplementary data sets from spacecraft that observe the Sun from 1 AU have been used: (1) The 10.7-cm solar radio flux, F 10.7 , from the OMNI-2 data set (King & Papitashvili, 2005). This index is an excellent indicator of solar activity, and it is measured daily at the Penticton Radio Observatory in British Columbia, Canada.
(2) The SSN also from the OMNI-2 data set, which takes the data from the WDC-SILSO, Royal Observatory of Belgium, Brussels.
(3) Solar EUV flux, I EUV , measured by the Thermosphere, Ionosphere, and Mesosphere Energetics and Dynamics (TIMED) Solar EUV Experiment (SEE; Woods et al., 1998). The solar EUV flux has also been extrapolated from Earth to Mars to account for the angular and radial separation between the two planets with respect to the solar surface (e.g., Hall et al., 2016;Yamauchi et al., 2015).
To study the spatial evolution of the Martian bow shock position over the period of a Martian year, Hall et al. (2016) computed a 2-D statistical model of the Martian bow shock surface in an axisymmetric plane (about Sun-oriented axis) of the solar wind aberrated Mars-Centric Solar Orbital (MSO) coordinate system. To study the spatial effects with respect to temporal changes in the bow shock position, the authors then used this model surface to extrapolate each individual crossing to a common reference point, that is, the terminator plane. In this study, we provide a recalculated and more robust model from the entirety of the extended HALL16 crossing data set such that we can use it to complete the same extrapolation. For the most part, we used the same fitting procedure as the Hall et al. (2016) study, which included the following: 1. Rotating the crossings into the solar wind aberrated MSO system (4°rotation about the Z MSO axis), thus accounting for the average direction of solar wind with respect to Mars' orbital motion. 2. Transforming the Cartesian aberrated MSO coordinates into polar coordinates as defined by equations (1) and (2), where r and θ are the polar coordinates of each crossing with respect to a conic's focus located at (x 0 , 0, 0). 3. Computing a regression analysis of all crossings with respect to the linearized form of the conic section equation (equation (3) in the form of r −1 ∝ cos θ), where L is the semilatus rectum and ϵ is the conic's eccentricity.
The Hall et al. (2016) study repeated steps 2 and 3 of this method for many different x 0 parameters, determining the best final set of conic parameters (x 0 , L, and ϵ) from the regression results giving the maximum coefficient of determination (i.e., the R 2 correlation statistic). Here we have improved upon the selection of the x 0 parameter by adopting basic machine learning techniques such as train-test and k-fold cross-validation (CV) schema (e.g., Fushiki, 2011, and references in there). For each input x 0 parameter, this new process adhered to the following procedure: 1. Split the full data set into 10 random subsets or folds (note: the same random subsets are used for each x 0 ). 2. Select one of the folds as a validation/testing data set and set aside. Perform regression analysis on the data in the remaining nine folds (or training subsets). Next, apply resultant model on reserved validation fold to measure its predictive performance (i.e., R 2 statistic). 3. Repeat step 2 until all folds have been used as a training or validation/testing subset, and average the fit statistics of each run to give a performance metric for each input x 0 parameter.
The final x 0 fit parameter was chosen from the CV run giving the largest R 2 output. Finally, this optimal x 0 was used in the fitting procedure applied to the full data set, giving the full set of conic section fit parameters defined earlier. This CV technique gives a more statistically robust estimate of x 0 than the Hall et al. (2016) method that effectively used a one fold CV. Using this more robust fitting procedure, we obtain the following for the whole data set (including the new additional crossings; Table 3 To extrapolate each bow shock crossing to the terminator plane, each set of aberrated MSO coordinates and the final x 0 and ϵ conic fit parameters are substituted into equations (1) to (3) to determine an L parameter for each crossing. All three conic parameters are then substituted into equation (4) to obtain the extrapolated bow shock terminator distance, R TD .

10.1029/2018JA026404
Journal of Geophysical Research: Space Physics This is the same method as used in the Hall et al. (2016) study. In doing this, a proxy for the bow shock positions at a common reference point can be tracked over time.
The fit parameters and model terminator distance for the updated bow shock model from the entirety of this study's crossing list and for the same procedure applied to each Martian year of observations are given in Table 3. The uncertainties given in this table are from the standard error on the fit parameters and from their consequent use in standard error combination formulae (Hughes & Hase, 2010) for any computed values (e.g., model R TD ). These models will be referred to throughout the remainder of the paper.  Hall et al. (2016) highlighted that this factor drives spatial variations in the average Martian bow shock position over the Martian year. However, we note that despite the MY-period variation, the solar cycle variability is clearly visible in the EUV flux extrapolated to Mars, especially when comparing the EUV flux at same heliocentric distances for different MY (e.g., each EUV minimum and maximum that correspond to farther and closer Mars-Sun positions, respectively). Figure 1c shows each individual crossing extrapolated to the terminator plane (grey dots) and a median filtered (27-day temporal window size) profile (black line) to show its average variability within a solar rotation time period. The blue envelope about this profile represents the median absolute deviation of the extrapolated R TD within each filter box. The absence of a filtered profile at the start of the data set is due to the limit set, which is set on the filtering of least 54 crossings per filter window, corresponding to a minimum of two bow shock crossings per day. For nominal MEX orbital trajectories, two crossings per orbit are expected, or approximately six crossings per day. We have relaxed our filter-window observation limit from this in order Journal of Geophysical Research: Space Physics to present as much data as possible, while minimizing spurious filter results. Finally, the median filtered R TD value at each timestamp is calculated from the preceding 27 days of observations.

Results
Compared with the median R TD across the entire unfiltered data set ( g R TD ¼ 2:48 R M , dashed red line in Figure 1c), we see two periodic variations in the filtered R TD profile (black line). The first corresponds to the period of a MY and the annual modulation of the EUV flux at Mars (clearest in MY29 and 30), as noted previously by Hall et al. (2016). The second is a longer-term trend where the filtered R TD is mostly below g R TD during the solar minimum phase (MY28-30), then mostly above g R TD during the solar maximum phase (MY31-32), before reducing again to around g R TD values as the solar activity declines again around MY33. This longer-term variation is more relatable to the evolution of the solar activity proxies over a period of the solar cycle. Figure 2 highlights both of these variations in more detail, in particular with a focus on quantifying the apparent solar cycle variation. Figure 2a presents the solar EUV flux at Mars (blue), the median-filtered profile of the extended HALL16 with a pass band of 27 days in red to remove the solar rotation period, and the median-filtered profile of the extended HALL16 with a pass band of 687 days in green to remove the variation caused by a Martian year. The three profiles in Figure 2a are on the same ordinate scale. The scale of these parameters are in standardized units (denoted R′ TD and I′ EUV ), where the data set population mean, μ, is subtracted from each observation before normalizing by the population standard deviation. This scaling allows us to directly compare multiple quantities that have different magnitudes, means, and variances (i.e., in standardized form they both have μ = 0 and σ = 1). A line at the standardized value of 0 (R TD − μ (R TD ) = 0σ) is superimposed to aid the description. The MY, a color-coded season (in terms of the solar longitude, L S , in the key below panel), has also been provided above this panel for more context.

Journal of Geophysical Research: Space Physics
The EUV trend and the 27-day filtered R′ TD trends clearly vary in almost unison across the period of the solar cycle, although the bow shock total variation between Mars' aphelion and perihelion is larger after 2011 (higher solar activity phase) than before (lower solar activity phase). There is a notable exception at the aphelion of MY32 where the bow shock variability when compared to perihelion is smaller than for other years. Moreover, during the "through perihelion" passages of Mars orbit (red color-coded season sectors), R′ TD,27day can peak by a larger amount than the EUV flux (e.g., see MY29 and 30 for clear examples). We do not know the origin of this extra variability, which seems not related to EUV fluxes, but we note that at that time, there was a notable rise on the soft X-ray background radiation, which could explain these larger values from 2011 to mid-2015 (see e.g. Sanchez-Cano et al., 2015. X-ray fluxes are known to be the ionization source below the main peak of the ionosphere (~130 km; e.g. Mendillo et al., 2006). Therefore, if the X-ray flux is more intense during a particular period, the lower ionosphere is more robust and can contribute to an enhanced thermal pressure of the ionosphere. Consequently, the Mars' plasma obstacle to the solar wind is enhanced, and the bow shock location could be found further from Mars. The MY R′ TD profile (green) also had a minimum observation limit for the filter window set to two observations per day (or 1,374 observations per filter window). Consequently, this filtered profile starts later into the data set during the perihelion sector of MY28. Nevertheless, this filtered profile clearly shows a change from the average bow shock terminator position being around or below −1 standardized units (R TD − μ(R TD ) = −1σ) during solar minimum (MY28 and 29; see Figure 1a), before increasing to~1 standardized units R TD − μ(R TD ) = 1σ) around solar maximum (MY31 and 32). Moreover, as the solar activity then begins to decrease again (e.g., MY33), this profile starts to decrease back to the average amount of 0 standardized units. The R′ TD,687day profile (green) clearly shows that despite any other variation caused by Mars' orbit about the Sun, there is a clear solar cycle trend on the bow shock position. This confirms the finding of previous studies such as Trotignon et al. (1993).
While the time series shown in Figure 2a shows both clear annual and solar cycle variations in the Martian bow shock position, to extrapolate each crossing to the terminator plane, we have assumed that the model Martian bow shock remains the same across the entire period. To investigate this assumption, a model for the bow shock (using the method outlined earlier) was calculated for each MY of observations across the period (annual average). It was not possible to perform this annual fit for same Sun-Mars distance conditions due to scarce data coverage for all the periods at a same distance. The results of this are shown across Figures 2b-2i and in Table 3.  Table 3) have been doubled for visibility. This panel also shows reference lines at the mean (μ = 2.43R M , solid black) and 1σ (σ = 0.06R M , dashed black) levels of the R TD,model population.
Figures 2b-h shows the variation in the spatial coverage and observation frequency of the Martian bow shock by MEX between each MY. For example, MY27 (panel 2b) had fewer observations than any other period but were nevertheless well distributed along the shock surface in this reference frame. All other MYs had many more observations (also see Table 3) but sometimes had significant gaps in spatial uniformity about the terminator (e.g., MY29-31, panels 2d-e). Fortunately, such MYs also had significant numbers of observations either side of the terminator, resulting in the fit being well constrained about the terminator and giving an R TD,model estimate expected to be representative of the period. Figure 2a (and 1c earlier) shows that each MY of crossings are not necessarily well distributed throughout all seasons of that year. For example, MY28 and 29 (panels 2b and 2c) both have reduced crossing frequency during Mars' aphelion passage (aqua color-coded season sector). Conversely, MY28 and 32 have limited numbers of crossings during the perihelion phase (red sector). As reported previously by Hall et al. (2016), the Martian bow shock exhibits a significant variation in average position between these two seasonal sectors. It is due to these above factors that we often opt to use the entire data set, where all seasonal periods are covered multiple times and the crossings present an equal spatial uniformity along the shock surface, to extrapolate each observation to a common reference point. Despite the spatial and temporal coverage issues noted previously, the R TD,model -MY time series in Figure 2i demonstrates the same solar cycle period variation in the average Martian bow shock position that was indicated earlier by the extrapolation method (Figure 2a). Figures 1c and 2b-2h are most probably caused by short timescale external variations of the incident solar wind plasma such as dynamic pressure, IMF direction (convection electric field influence and dawn to dusk asymmetry), or Mach numbers, but also possibly due to intrinsic microscopic variations of the collisionless shock structure even during steady external conditions (nonstationary character). All these variations occur on an orbit-to-orbit basis and do not influence the average bow shock location on a much larger timescale when considered statistically.

The observed large dispersion in
In summary, the average bow shock position at the Martian terminator is closer to Mars at low solar activity (e.g., MY28-29) than it is during the higher solar activity of solar maximum (MY31-32). Quantitatively, this is an~0.17R M (or~7%) increase in the terminator position between the low and high solar activity periods of solar cycle 23/24. This 7% solar cycle variation (which could be larger for other more extreme solar cycles; see next section) occurs in addition to the overall 11% variation throughout the Martian year . We note that the 11% annual variation was obtained after averaging over all levels of solar activities. Due to the nonuniformity of observations across every season and every year, it is not possible to get a robust estimation of the annual bow shock location at each MY. Nevertheless, we show for the first time with the same data set that both annual and solar cycle variations play major roles in the variation of the bow shock average location at Mars on a long-term basis.

Discussion
This study has aimed to answer the long-standing question of whether the Martian bow shock position varies over the period of a solar cycle. We have used the Hall et al. (2016) Martian bow shock crossing list and extended it with a newly released data set to give~13 years of MEX observations of the Martian bow shock crossings over the period of 2004 and 2017. This period includes~7 Martian years occurring over solar cycles 23 and 24. By modeling the average bow shock surface from the entire data set ( Figure 1) and from each Martian year of observations (Figure 2), we have been able to demonstrate the bow shock position varying over the period of a Martian year (first reported by Hall et al., 2016) and solar cycle. We have also demonstrated that the magnitudes of both of these spatial and temporal variations are similar in magnitude to those seen in the solar EUV flux over the period of the solar cycle. This suggests that along with driving the annual variability in the bow shock position , the solar EUV flux could also be a main driver of the solar cycle period variations.
Earlier studies of this topic did not converge on a clear consensus regarding the role of the solar EUV or Mars season regarding bow shock position (e.g., Russell et al., 1992;Trotignon et al., 1993;compared to Slavin et al., 1991;Vignes et al., 2000Vignes et al., , 2002. This was likely due to the previous studies being limited by low numbers of bow shock observations that often had both low sampling coverage across space and the entirety of even a single MY (e.g., maximum of 450 by MGS mission over sub-MY period). This and the Hall et al. (2016) study have both demonstrated the importance of considering how Mars' position through its orbit of the Sun can significantly impact the average location of the bow shock position. Unfortunately, the solar cycle period covered by MEX consisted of both the lowest minimum and maximum levels of the recent solar cycle history (e.g., McComas et al., 2013). Therefore, we expect that the percentage variation between minimum and maximum at previous solar cycles were larger, proportional to their EUV fluxes. Trotignon et al. (1993) compared the Mars-2,3,5 (low activity, solar cycle 20) and Phobos-2 (high activity, solar cycle 22) bow shock crossings (see Tables 1 and 2) and reported a solar cycle period variability with a total variation of~11% in position at the terminator. This is similar to the 7% we report here (e.g., comparing models from MY32 as high activity and MY28 as low activity; Figures 2c, 2g, and 2i), which might be slightly lower due to the general lower activity of solar cycle 23/24 as compared to preceding cycles. Accordingly, one would expect the bow shock position closer to Mars during the very low minimum of solar cycle 23/24 than at previous solar minima because of a much lower EUV-flux level. However, it is not possible to conclude anything in this respect because the number of observations from former missions is not sufficient to make a proper statistical comparison. We can only say that R TD at minimum of solar cycle 20 and maximum of solar cycle 22 are reasonably comparable to the MEX estimates of solar cycle 23/24 (see Table 2). In addition, previous missions have observed bow shock crossings over different regions over Mars. Only the near-continuous MEX and MAVEN (since 2014) observations are able to provide a more coherent sample to assess this issue, which will be 10.1029/2018JA026404 Journal of Geophysical Research: Space Physics improved in the coming years when two continuous solar cycles will be sampled, and in situ solar wind and solar radiation observations will be available for several Martian years.
The orbits belonging to single-spacecraft observations of the Martian bow shock can impact the results of studies such as this in two main ways. If the spacecraft has an orbital apoapsis below the bow shock, it will simply not be observed. Similarly, if any hemispherical asymmetries exist in the bow shock's location about Mars (e.g., Hall et al., 2016) and a single spacecraft observes the bow shock in these hemispheres by different amounts, a bias can quickly develop. Although this is undoubtedly a source of error, we expect it to impact the results by a very small amount when compared to the whole MEX data set.
While the limitations of the earlier studies can explain the differences to our newfound understanding reported here, recent attempts at theoretically modeling the entire Martian plasma system during differing levels of solar activity have also reported no obvious solar cycle variations in the Martian bow shock position (e.g., Modolo et al., 2005Modolo et al., , 2006Modolo et al., , 2016. Here we only comment that this difference is surprising and that future work should attempt to seek out why theoretical modeling attempts show no obvious solar cycle variability. In Figure 2a of section 4, we noted how the standardized solar EUV flux (extrapolated to Mars) and extrapolated bow shock terminator distance over the period of the solar cycle both varied by similar extents over the Martian year and over the solar cycle periods. Consequently, we agree with the mechanisms proposed by the Hall et al. (2016) study where varying levels of solar EUV flux can modulate the way the solar wind interacts with the Martian exosphere (exospheric ion pickup) and ionosphere (enhanced ionization and thermal pressure), resulting in a change in location of the bow shock position. However, we noted in section 4 that the extrapolated bow shock terminator distance peaks more during perihelion than the solar EUV flux does during the same period (see Figure 2a blue and red profiles). While we have no direct explanation for this, it could be related to other more energetic radiation, such as the soft X-ray background irradiance flux, which ionized the ionosphere below the main peak. This is a clear future study line to be done with MAVEN and MEX conjoined observations. In addition, another possible explanation could reside in an annual (MY) variability in the average solar wind magnetosonic Mach number at Mars. Edberg et al. (2010) provided preliminary evidence that the Martian bow shock terminator distance is linearly anticorrelated to the magnetosonic Mach number. According to the Edberg et al. (2010) study, the expected lower average Mach number during Mars' perihelion (cf. aphelion) could drive the bow shock further from the planet than what the solar EUV can do alone. This could explain the higher peaking of R TD around perihelion (compared to the solar EUV) noted here in our results (Figure 2a). This agrees with the recent MAVEN results that have found that at the terminator distance, the bow shock position moves around 0.2R M closer to the planet during high magnetosonic Mach number time periods, being 2.41-2.66R M (high-low Mach number) over the North pole and 2.53-2.83R M over the South (Gruesbeck et al., 2018). Moreover, the magnetosonic Mach number also anticorrelated with the solar cycle. Extrapolating once again from the Edberg et al. (2009) results, an on average lower Mach number at solar maximum should result in the Martian bow shock being on average further away from the planet, giving a similar effect on the bow shock position as described above at Mars perihelion. While a MY filtered EUV profile is not shown in Figure 2, we have checked and found that the 687-day filtered R′ TD during solar maximum does exceed the solar EUV flux filtered to the same period, matching our observations at perihelion. Thus, the impact of the magnetospheric Mach number on the average bow shock position with solar cycle is likely an important driver that needs better parameterization. Moreover, as previously mentioned, high magnetosonic Mach number also favors rapid oscillations of the bow shock front location (nonstationarity). The role of the Mach number in this type of study could be done in a future study together with MAVEN data sets. Currently, the MAVEN data set covers less than two Martian years and cannot be used today for a similar purpose. However, one can take advantage of the MY32-33 overlapping of both missions and compare similar data sets on both spacecraft, with the advantage that MAVEN can provide all the necessary drivers.
While this study has provided clear evidence of the average Martian bow shock position varying over both the Martian year and solar cycle, there are several routes that future work could follow. Such topics include a full parameterization of statistical Martian bow shock models (e.g., including dynamic pressure, IMF direction, Mach numbers, and intrinsic microscopic variations of the shock structure), which may aid in explaining the highly variable bow shock position on short timescales. In addition it is a comparison with similar (e.g., Venus, comets) and dissimilar (e.g., Earth, Mercury) plasma systems as well as a utilization of multispacecraft observations of the bow shock and inner Martian plasma system (e.g., induced magnetosonic boundary and ionopause) in order to explain the dynamics of the full plasma system as a whole. Finally, an assimilation of observational data sets with theoretical models could be done to explain any differences.

Conclusion
In this paper we have used~13 years of near-continuous observations of the Martian bow shock by the MEX mission to report that its average position is sensitive to solar activity and thus varies over the solar cycle. We have provided evidence of this by producing statistical models of the bow shock over long (entire bow shock crossing data set) and short (each of the 7 Martian years across solar cycle 23-24) time periods, consequently using these results to estimate the extent of the variability. At the Martian terminator, we have found that between a period of low and high solar activity, the average Martian bow shock position increases bỹ 0.17R M , or alternatively by an increase in~7%. This is similar in magnitude and in addition to the annual variability noted by Hall et al. (2016), which reported an increase in the average Martian bow shock position between Mars' aphelion and perihelion by~0.17R M (~11%). Therefore, future studies into the spatial variability of the Martian bow shock should take these aspects into account, and with the continued support of both the MEX and MAVEN missions into the coming years, future works will be able to provide a much more detailed understanding of the variability reported here and of the entire Martian plasma system.