Annual variations in the Martian bow shock location as observed by the Mars Express mission

The Martian bow shock distance has previously been shown to be anticorrelated with solar wind dynamic pressure but correlated with solar extreme ultraviolet (EUV) irradiance. Since both of these solar parameters reduce with the square of the distance from the Sun, and Mars' orbit about the Sun increases by ∼0.3 AU from perihelion to aphelion, it is not clear how the bow shock location will respond to variations in these solar parameters, if at all, throughout its orbit. In order to characterize such a response, we use more than 5 Martian years of Mars Express Analyser of Space Plasma and EneRgetic Atoms (ASPERA‐3) Electron Spectrometer measurements to automatically identify 11,861 bow shock crossings. We have discovered that the bow shock distance as a function of solar longitude has a minimum of 2.39RM around aphelion and proceeds to a maximum of 2.65RM around perihelion, presenting an overall variation of ∼11% throughout the Martian orbit. We have verified previous findings that the bow shock in southern hemisphere is on average located farther away from Mars than in the northern hemisphere. However, this hemispherical asymmetry is small (total distance variation of ∼2.4%), and the same annual variations occur irrespective of the hemisphere. We have identified that the bow shock location is more sensitive to variations in the solar EUV irradiance than to solar wind dynamic pressure variations. We have proposed possible interaction mechanisms between the solar EUV flux and Martian plasma environment that could explain this annual variation in bow shock location.


Introduction
Planetary bow shocks mark the region where the interplanetary solar wind flow starts to be perturbed by the presence of a downstream obstacle such as a planetary magnetosphere or atmosphere. Mars is an unmagnetized body with an atmosphere, which leads to its ionosphere and induced magnetosphere being the main obstacle to the solar wind. Additionally, the relatively small size, mass, and therefore gravity of Mars lead to an extended exosphere which also interacts directly with the solar wind [Mazelle et al., 2004]. Variations in these regions of the Martian plasma environment are expected to play a role in subsequent variations in the bow shock boundary location.
The bow shock shape has been studied in detail [e.g., Russell, 1977;Slavin and Holzer, 1981;Schwingenschuh et al., 1990;Slavin et al., 1991;Trotignon et al., 1991aTrotignon et al., , 1991bTrotignon et al., , 1993Trotignon et al., , 2006Vignes et al., 2000;Edberg et al., 2008Edberg et al., , 2009a, with modeling of the bow shock tending to involve least squares (LSQ) fitting of an axisymmetric conic section to crossings in the solar wind-aberrated cylindrically symmetric coordinate system (as done in this study and described in more detail in section 3.3). Despite individual bow shock crossings showing a large variability in position, and the LSQ fitting being completed on data sets of differing size, spacecraft used, and time periods (e.g., 11 crossings by the Mars series in Russell [1977] and 993 crossings by MGS in Edberg et al. [2008]), the resultant models have presented a surprisingly stable average shape well defined by a conic section. The average subsolar and terminator (i.e., solar zenith angle of 90 ∘ ) distances of the bow shock from all these models are R SS ∼ 1.58 R M and R TD ∼ 2.6 R M (Martian radius, R M ∼ 3390 km), respectively [Bertucci et al., 2011].
The high variability in location between individual crossings has been studied with respect to many possible driving factors. Typically, studies of bow shock variability have involved mapping crossings to the terminator plane in order to remove solar zenith angle (SZA) dependence such that each crossing can be directly compared with another. The bow shock location has been shown to reduce in altitude with increasing solar wind dynamic pressure, presenting a weak exponential relationship [e.g., Schwingenschuh et al., 1992;Verigin et al., 1993;Edberg et al., 2009b]. Furthermore, the bow shock location has been found to increase in altitude with increasing solar extreme ultraviolet (EUV) flux, presenting a strong exponential relationship [Edberg et al., 2009b]. Two ways in which solar EUV can impact the variability of the bow shock are through modulation of the ionization rate in the Martian ionosphere (i.e., main obstacle to solar wind flow) and by creating new ions within the Martian extended exosphere that can further interact with the solar wind. During periods of enhanced solar EUV irradiation, the Martian ionospheric total plasma density and therefore thermal pressure both increase , making the ionosphere a more effective obstacle to the solar wind flow. Newly created ions within the Martian extended exosphere are picked up and accelerated by the electromagnetic fields carried by the solar wind. The solar wind is mass loaded and slowed by these picked up exospheric ions, thereby potentially affecting the location of the bow shock [Mazelle et al., 2004;Yamauchi et al., 2015].
Other bow shock variability studies have included driving factors such as the phase of the solar cycle [e.g., Russell et al., 1992;Vignes et al., 2000Vignes et al., , 2002, the interplanetary magnetic field (IMF) and convective electric field direction [e.g. Dubinin et al., 1998;Vignes et al., 2002Vignes et al., , 2009b, the magnetosonic Mach number [e.g., Edberg et al., 2010], and the presence of the intense localized Martian crustal magnetic fields [e.g., Slavin et al., 1991;Schwingenschuh et al., 1992;Vignes et al., 2002;Edberg et al., 2008]. The impact of solar cycle phase on bow shock variation is not well understood due to the lack of continuous bow shock observations across a full solar cycle inhibiting such studies. The IMF direction has been found to raise the bow shock altitude when having an approximate perpendicular alignment to the shock front normal (i.e., quasi-perpendicular shocks). Similarly, the bow shock location has also been found to increase in altitude as the magnetosonic Mach number of the solar wind decreases. The convective electric field direction and presence of crustal magnetic fields have both been found to lead to asymmetries in bow shock location between different hemispheres, with the former being attributed to enhanced pickup of exospheric oxygen by the solar wind.
Aside from the aforementioned drivers of variations in the Martian bow shock location, one particular area that is yet to be described are any variations that may be linked with Mars' orbit of the Sun. Mars has an elliptical orbit about the Sun (eccentricity of ∼0.09), leading to its farthest distance from the Sun being ∼20% larger than its closest distance (aphelion ∼1.67 AU, perihelion ∼1.38 AU). Solar parameters, such as solar wind density (therefore dynamic pressure), interplanetary magnetic field magnitude, and solar irradiation, are all expected to reduce with the square of the distance from the Sun. Since these parameters impact the bow shock location in different senses (i.e., bow shock moves toward Mars at higher solar wind density and pressure but moves away when the solar EUV is higher), it is important to characterize which is the dominant factor throughout the Martian year. The study presented in this paper investigates this by using more than 5 Martian years (∼10 Earth years) of bow shock crossings identified through MEX electron plasma measurements.

Instrumentation and Data
The MEX mission (late 2003(late -present, Chicarro et al. [2004) is in an elliptical polar orbit of Mars. Its periapse and apoapse altitudes are at ∼300 km and ∼10,100 km, respectively, defining an orbital period of ∼7.5 h.
To obtain a global coverage of Mars, the orbit was designed such that the latitudinal periapse location would shift over time. The high apoapse altitude allows MEX to typically pass through the Martian bow shock twice per orbit (inbound and outbound passage). However, due to the temporal evolution of the MEX orbital geometry, at times MEX may not enter the solar wind at all, thus limiting the observation of a bow shock crossing. Alternatively, the orbital geometry can be such that MEX is traveling quasi-parallel to the bow shock surface, which in turn could manifest as multiple short lasting bow shock crossings (i.e., as MEX passes from solar wind to magnetosheath and back). In general, the MEX orbit has provided excellent (near symmetrical) coverage in the northern and southern hemispheres of Mars. However, for the majority of the mission to date, MEX periapse has been located in the dusk hemisphere, leading to the bow shock being sampled by MEX mostly within the dawn hemisphere of Mars.
The longevity of the MEX mission has provided well calibrated, long-term plasma measurements at Mars spanning a period of more than 5 Martian years. The Analyser of Space Plasma and EneRgetic Atoms (ASPERA-3) package on board MEX has been in operation since the start of the MEX mission and contains instruments for measuring charged particles (electrons and ions) from regions such as the solar wind down to the Martian upper ionosphere . The ELectron Spectrometer (ELS) measures superthermal electrons within the energy range of 0.001-20 keV, has an energy resolution of 8%, and a field of view (FOV) of 4 ∘ (polar) × 360 ∘ (azimuth). Although ELS is designed to have many operational modes  and Hall et al.
[2016] describe a few of these modes), for the majority of the mission it has operated in what is known as survey mode. In this mode ELS samples electrons across the full energy range separated into 128 log-spaced steps at a cadence of 4 s. The Ion Mass Analyzer (IMA) initially measured ions in the energy range 30 eV/q-30 keV/q, but a new energy table uploaded in May 2007 brought the lower limit to 10 eV/q. Since November 2009 the energy table has been updated again and is now divided across two intervals, first 20 keV/q-20 eV/q in 76 log-spaced steps and then 20 eV/q-1 eV/q in 10 linearly spaced steps (R. Ramstad, private communication, 2016). The energy resolution has remained at 7%, and IMA can resolve ions into masses of 1, 2, 4, 16, 32, and 44 amu/q, corresponding to the main ion species seen at Mars. The FOV of IMA is 90 ∘ (polar) × 360 ∘ (azimuth), with the polar FOV broken into 16 sequential elevation steps (−45 ∘ to 45 ∘ ). A full energy scan per elevation takes 12 s, resulting in the full FOV taking 192 s.
MEX also carries a nadir-looking pulse limited radar sounder known as the Mars Advanced Radar for Subsurface and Ionosphere Sounding (MARSIS) [Picardi et al., 2004], which has been in operation since mid-June 2005 [Orosei et al., 2015]. MARSIS operates across two modes and five transmission frequencies. In Active Ionospheric Sounding (AIS) mode [Gurnett et al., 2005], MARSIS samples the Martian topside ionosphere, while in the SubSurface (SS) mode [Safaeinili et al., 2007], it analyzes the Martian surface and subsurface using synthetic aperture techniques. Both modes, although through different means, can characterize the ionosphere (the effective obstacle to the solar wind flow) through its total electron content (TEC). The TEC is a measurement of the number of free electrons in an atmospheric column with vertical extent set by two altitudes and a cross section of 1 m 2 [e.g., Sánchez-Cano et al., 2015]. The topside ionospheric TEC (via MARSIS-AIS) is obtained through numerical integration of the vertical electron density profiles with altitude [e.g., Sánchez-Cano et al., 2012. The TEC of the entire atmosphere is retrieved as a by-product of MARSIS-SS mode radar signal distortions caused by the ionosphere [Safaeinili et al., 2007;Mouginot et al., 2008;Cartacci et al., 2013;Sánchez-Cano et al., 2015. Our study uses daily averaged TEC measurements of the entire atmosphere from this latter method at a solar zenith angle (SZA) of 85 ∘ . This data product is obtained via the Cartacci et al. [2013] algorithm, from which TEC daily averages are calculated from at least three orbits per day. In order to improve the accuracy of the TEC estimation, only data that have a signal-to-noise ratio >25 dB are used [Sánchez-Cano et al., 2015. We note that the near-terminator TEC is typically more variable than at lower zeniths. A near-terminator TEC (e.g., SZA = 85 ∘ ) is used in this study for two reasons.
(1) The MARSIS-SS mode was designed to probe the Martian surface for water and thus operates best when the signal is not distorted by the ionosphere (i.e., on the nightside of Mars where the TEC rapidly reduces to zero due to electron and ion recombination). Since the TEC over the full atmosphere is obtained from the signal distortion by the ionosphere, we require SZAs when the ionosphere is still present. Unfortunately, at lower SZAs (i.e., toward subsolar point) the signal to noise becomes too low for the TEC of the full atmosphere to be determined reliably. This limits our choice of TEC estimates to higher SZAs on the dayside.
(2) Our later analysis of the variability in the bow shock location (section 4 onward) involves extrapolating each bow shock crossing to the terminator plane (see section 3.3). Consequently, the TEC at near-terminator SZAs is more appropriate for our study.
The Martian ionosphere is primarily produced by photochemical reactions driven by the solar EUV flux impinging on the neutral upper atmosphere [e.g., McElroy et al., 1976;Hanson et al., 1977;Mantas and Hanson, 1979]. Thus, to understand how the ionosphere may vary we need a measure of how the solar EUV flux varies. Using the daily averaged measurements of the solar EUV spectrum by the Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics (TIMED) Solar EUV Experiment (SEE) [Woods et al., 1998] at Earth (∼1 AU), we have derived a proxy for the solar EUV irradiance impinging on the Martian upper atmosphere. To derive our proxy measurement we first integrate the TIMED-SEE spectrum across the wavelength range of 10-124 nm to get a single irradiance (I, units of W m −2 ). Next, we extrapolate each daily averaged value to Mars via (1) reducing the magnitude by the square of the distance between Mars and 1 AU, and (2) time-shifting the time stamp of each measurement by the amount of time taken for the same part of the solar surface to rotate from Earth to Mars [see equation (3) of Vennerstrom et al. [2003] and Edberg et al. [2010]. This proxy measurement of the solar EUV flux at Mars is used for two reasons: (1) to be consistent with the proxy solar EUV measurements used in recent studies concerning the Martian plasma environment [e.g., Yamauchi et al., 2015], and (2) the wide range of wavelengths caters for ionization of the neutral Martian atmosphere across a wide range of altitudes.
Regarding point 2, in general, smaller wavelength solar radiation penetrates to lower altitudes of the atmosphere before interacting with, and ionizing a neutral component. Thus, this proxy measurement of the solar EUV irradiance is suitable for describing the main contributing factor to variations in the Martian ionosphere and exosphere, both of which are likely to impact the bow shock location.
From the instruments described above, we are able to (1) identify Martian bow shock crossings by MEX with ASPERA-3 ELS; (2) determine solar wind plasma moments (i.e., density, velocity, and dynamic pressure) at Mars using ASPERA-3 IMA (see Fränz et al. [2007] for a method of determining plasma moments from particle distributions and Ramstad et al. [2015] for more details with respect to IMA); (3) identify the state of the Martian atmosphere/ionosphere using MARSIS TEC measurements; and (4) track the solar EUV irradiation, a major driving factor of ionization in the Martian plasma system, using the extrapolated TIMED-SEE measurements.
In their entirety, each of these data sets covers different time periods ( MEX orbits of Mars (starting from the first science orbit), although only 10,295 of these orbits had the ASPERA-3 ELS instrument in operation. Moreover, a number of these orbits may not contain valid data for identifying the Martian bow shock (e.g., due to orbital geometry of MEX, etc.). The following section describes our automated method used for detecting a crossing of the Martian bow shock by MEX and includes a description of how it was determined whether a MEX orbit would allow for a bow shock detection. Additionally, for all of the above data sets we have only included measurements with corresponding good data flags.

MEX ASPERA-3 ELS Observations of Martian Bow Shock
We identify the Martian bow shock using ASPERA-3 ELS measurements alone since it has both been more consistently in operation and has a much higher measurement cadence than the ASPERA-3 IMA instrument (4 s; cf. IMA's 192 s for full scan). For the MEX orbit #14352 (25 April 2015 with periapsis at ∼01:52 UT), we present the ELS electron energy-time spectrogram in terms of the differential number flux (DNF, quantified in the colorbar to the right of the panel) in Figure 1a, a derived proxy parameter describing the DNF sampled between the energies of 20-200 eV (f , green profile) in Figure 1b, and the MEX ephemerides of altitude (h MEX , left ordinate, blue) and total spacecraft velocity (|v MEX |, right ordinate, magenta) in Figure 1c, with all three panels sharing a common universal time (UT) scale (abscissa of Figure 1c). The DNF presented here is an average calculated across all ELS anodes that look away from the spacecraft body, with no further corrections being made. The aforementioned electron DNF proxy (f , Figure 1b) was calculated by integrating the DNF across the energy range of 20-200 eV. This f proxy was calculated in a similar way to Hall et al. [2016] but with two improvements made. First, we determine which ELS anodes look across the spacecraft body and exclude them from the calculation of the average electron DNF across the entire ELS sensor. Second, when integrating across the energy range of 20-200 eV, instead of summing the electron DNF across the energy range and The crossings are presented in the axisymmetric solar wind abberated MSO coordinate frame. The green curve represents the passage of MEX orbit #14352, with the inbound (black) and outbound (red) crossings identified by colored circles. The solid blue curve represents the best fit conic section to all crossings, and the orange crosses are located at the average subsolar and terminator distances from all previous models [Bertucci et al., 2011]. (e) The same as Figure 1d but with crossings grouped and counted into 0.1 × 0.1 R M spatial bins. Best fit conic model parameters are given in Table 1. multiplying by the total width of the energy range (i.e., 180 eV), we instead compute individual integrations of the electron DNF at each energy step between 20 and 200 eV and then sum each integration. Both methods give the same temporal variations, but these additional consideration will be more effective at removing any contamination effects due to the instrument looking across the spacecraft body. The f parameter is useful in compounding the majority of the information seen in the entire energy-time spectrogram as a single value. We now use these panels to describe a typical bow shock crossing as observed by MEX ASPERA-3 ELS.
As MEX crosses the Martian bow shock the ELS instrument typically registers a sudden increase in flux of electrons across a wide range of energies (typically up to a few hundred eV). In Figure 1a, up until ∼01:00 UT MEX is traveling through the solar wind with peak DNF around 10 eV. Between 01:00 UT and 01:05 UT an order of magnitude increase in the electron DNF across a wide range of energies occurs (up to ∼200 eV, with positive spacecraft potential partly responsible for peaking of the DNF at energies below 20 eV). After this initial enhancement the DNF remains high, signifying that MEX has passed from the solar wind into the Martian magnetosheath via the bow shock. Eventually, the DNF reduces in magnitude and energy, correlating with the MEX passage through the Martian magnetosheath and into the induced magnetosphere. Approaching the MEX passage in reverse (i.e., backward in time from the end of the spectra), we see a similar pattern that corresponds to the outbound pass of the Martian magnetosheath and bow shock (∼02:50 UT). At the time of both of the MEX bow shock crossings described above, the f parameter ( Figure 1b) more clearly shows an enhancement of more than an order of magnitude that occurs over a rapid time span of less than 5 min.
At this point it is worth noting how the period of time it takes MEX to cross the bow shock region compares to the typical size of the region. Typically, MEX is traveling around 2 km s −1 as it passes through the Martian bow shock. In this example orbit (Figures 1a-1c), MEX is traveling slightly faster as it passes the bow shock, with a total velocity of |v MEX | ∼ 2.5 km s −1 . Previous studies [e.g., Bertucci et al., 2011] have suggested that the size of the Martian bow shock (at the terminator) is of the order of the Larmor radius of solar wind protons at Mars. Using typical solar wind velocities and IMF magnitudes at Mars of v SW ∼ 400 km s −1 and B IMF | ∼ 3 nT , we arrive at a solar wind proton Larmor radius of r L,p ∼1400 km. Across a compressive bow shock the magnetic field magnitude (and solar wind density) can increase by up to a factor of 4 [Burgess, 1995], thus r L,p at the inner edge of the bow shock will be reduced by up to a factor of 4. Consequently, to a first order approximation we expect the bow shock width to be a few hundred kilometers in size. Thus, if MEX is traveling at ∼2 km s −1 , it will cross the entirety of the shock in less than 5 min. Of course, this is assuming MEX is traveling directly across (i.e., perpendicular) to the bow shock surface. If MEX crosses the bow shock at a more quasi-parallel incidence to the shock surface, the period of time over which a significant flux enhancement will occur will become significantly longer and difficult to identify through automated or visual means. In terms of the actual electron flux data, we see that the inbound bow shock crossing by MEX in Figures 1a and 1b occurs between 01:00 UT and 01:05 UT, which corresponds to a total altitude variation ∼0.13 R M (minor tick marks in Figure 1c are placed at every 0.25 R M or km s −1 ). This corresponds to an altitude change of 440 km which is in agreement with the bow shock width approximated above. Thus, it is likely that MEX was traveling somewhat quasi-perpendicular to the bow shock surface during this orbit. To identify the MEX crossings of the Martian bow shock, we have designed an automated algorithm (see below) that is sensitive to bow shock crossings represented in Figures 1a, and 1b.

Automatic Identification Algorithm
The formal description of our automated method for identifying crossings of the Martian bow shock by MEX within the f data set is as follows. First, we note that we would normally expect two bow shock crossings per MEX orbit (one per half orbit, e.g., see Figure 1a), so we first separate the flux proxy data within each MEX orbit into two segments (or half orbits) worth of data either side of the MEX periapsis time (at ∼01:52 UT of the example orbit in Figures 1a and 1c). To aid description, we label these segments of data as PRE and POST periapsis data. Bow shock crossings in the PRE (inbound) direction are typically a crossing from solar wind into Martian magnetosheath plasma and manifest as an enhancement in the flux proxy, whereas the opposite is typically true in the POST (outbound) direction. To make our algorithm simpler such that we only have to look for an enhancement in the flux proxy, we run our algorithm backward in time for POST periapsis style bow shock crossings. Thus, all identified crossings are detected in the frame of MEX transiting from solar wind to magnetosheath plasma. Due to this, we require that for each MEX half orbit, the ASPERA-3 ELS data must start (PRE) or end (POST) in the solar wind. A lack of any ASPERA-3 ELS data in the solar wind can be due to either periods when the MEX orbital trajectory has its apoapsis passage entirely within the Martian magnetosheath (i.e., no bow shock crossings expected) or the instrument is simply not in operation during the MEX passage through the solar wind plasma. Visual identification of ELS data (energy-time spectra) in the solar wind would be relatively simple, but, since we use an automatic method and with the a posteriori knowledge that the flux proxy within the solar wind can have a similar magnitude to that in the induced magnetosphere (IM) and ionosphere, we use the Edberg et al. [2008] bow shock model to help determine whether a MEX half orbit is expected to have any solar wind flux proxy data. If ASPERA-3 ELS data are found outside of the Edberg et al.
[2008] bow shock model, we use our fully automatic (FA) algorithm, otherwise we use our semiautomatic (SA) algorithm, the details of which will now be described.
To aid in describing our algorithm, in Figure 2 we have presented each step of the FA algorithm (described below) as it progresses throughout the first half of the MEX orbit shown in Figure 1 up to the point of identifying the inbound bow shock crossing. In all of the panels of this figure, the f proxy is on the ordinate, and UT is on the abscissa.
The first stage of the FA algorithm (Stage 1, Figure 2a) starts by applying a rolling median filter (temporal box width of 30 s) to the flux proxy data set to remove spurious high-frequency variations while minimizing blurring the exact location of a rapid enhancement. The impact of this is shown in Figure 2a where the raw f proxy profile is given in black, and the median filtered profile is shown in green. Just after 01:00 UT MEX crosses the bow shock, and we can clearly see that the filtered profile has a significant enhancement in f that is colocated at a very similar time to the raw profile. Additionally, as noted earlier we typically expect MEX to pass the entirety of the bow shock region within a few minutes, which is significantly longer than the temporal box width used in the filtering. Thus, this method of filtering the data to remove high-frequency variations is deemed suitable for this study.
After preparing the f data, we then start checking for whether solar wind data exist. As noted above the Edberg et al. [2008] bow shock and magnetic pileup boundary (MPB, inner Martian plasma boundary separating magnetosheath and upper atmospheric plasma) models are used to help identify the presence of solar wind data. This is a two step process which starts by preliminarily identifying solar wind f data as everything outside of the Edberg et al. [2008] bow shock model shifted 0.1 R M toward Mars (solid vertical black line at ∼01:05 UT in Figure 2a). We call this our solar wind f distribution and then calculate the median value, defining it as the average value of the solar wind flux proxy, f SW,avg (solid horizontal blue line from ∼22:25 UT to ∼01:05 UT in Figure 2a). It is important to note that it is entirely possible that this method will lead to some of the solar wind data distribution including f values representative of bow shock or magnetosheath region fluxes. This is true in Figure 2a where the solar wind data are expected to be up until ∼01:05 UT, but data between 01:00 UT and 01:05 UT are significantly larger than the data prior to 01:00 UT. Thus, we should consider the impact this could have on our value of f SW,avg . In the limit that the bow shock/magnetosheath data points only contribute a small proportion of the solar wind data distribution, our use of the median to calculate f SW,avg will prevent it being skewed to the higher values. This is the case in Figure 2a where f SW,avg is very similar to the f values prior to 01:00 UT (i.e., the actual solar wind data points). We would only expect the median value to be strongly influenced by additional bow shock/magnetosheath data points if they contribute to the majority of the distribution. In this case we will have very little true solar wind data leading to the bow shock to be unidentifiable via our automated method and also very difficult to identify visually. Nevertheless, a second test is also performed to check that this f SW,avg is representative of solar wind fluxes. This additional test compares the maximum flux proxy value (f max , red circle at ∼01:05 UT in Figure 2a)  topped red hatched box in Figure 2a). If true, we consider that this half orbit includes an enhancement that is large enough that it could describe a bow shock crossing and, therefore, that the f SW,avg value is suitable for describing the average f within the solar wind. An increase in 0.6 orders of magnitude was chosen since it represents a factor of 4 increase in f in log space. As noted earlier, across a compressional bow shock the IMF and solar wind density can increase by up to a factor of 4 [Burgess, 1995]; thus, we used such an increase to add a more physical description to the enhancement we expect to see across the bow shock.
At this stage of the algorithm, we have identified that within a MEX half orbit there are two regimes of f that are significantly separated in magnitude from each other. While this could represent a bow shock crossing, we are yet to evaluate where the enhancement occurs, whether the enhancement occurs on a rapid timescale (i.e., within 5 min), and if the initial enhancement continues. These questions are evaluated in next stage of the algorithm by iterating through each time step of the flux proxy data (f (t iter )) and comparing how the f evolves over the following 5 min of data (f (t box )). This stage of the algorithm is shown in Figure 2b-2g, where in each panel the black profile is the median filtered f (equivalent to smaller segments of the green profile shown in Figure 2a), the solid horizontal blue line is the f SW,avg , the light blue circles represents f (t iter ), and the transparent green box represents the 5 min of data following f (t iter ). At each f (t iter ) we evaluate the following three criteria.
1. Is f (t iter ) similar to solar wind data? The order of magnitude difference between f SW,avg and f (t iter ) (red circle in Figures 2b and 2e) must satisfy −0.3 ≤ log 10 (f SW,avg ∕f iter ) ≤ 0.1 (i.e., the red circle must be within the red hatched box of Figures 2b and 2e). 2. Is there a large enhancement from f (t iter ) within 5 min? The order of magnitude difference between the maximum value within the 5 min box f (t box,max ) (red circle in Figures 2c and 2f ), and f (t iter ) must satisfy log 10 (f box,max ∕f iter ) > 0.6 (i.e., the red circle must be within the red hatched box of Figures 2c and 2f ). 3. Does the f continue to increase? In the 5 min following f (t box,max ) (solid horizontal line in Figure 2g and transparent red box in Figure 2g), there must be more f values larger than f (t box,max ) than there is below it. The first time step of f is being evaluated by the above criteria, finding them all to be false and thus correctly identifying that no rapid enhancements occurs at this stage of the orbit. We note that the first criterion (i.e., is our starting data like solar wind data?) is incorrectly evaluated to false. This is by design as we found that having a smaller lower threshold caused the initial identification of a bow shock like enhancement to occur very close to the base of the increase in flux. Criteria 2 and 3 then evaluate if the enhancement is actually a bow shock or otherwise (e.g., a short lasting increase in flux that was not removed by the initial median filtering). In this evaluation of the criteria, criterion 3 is not evaluated since no significant increase in flux is identified. These three criteria will continually be evaluated as we sequentially step forward in time, until all are evaluated to be true. Figures 2e-2g (Stage 2: True) show an instance of this, with criterion 3 now being evaluated.
Since criterion 3 finds that the f has continued to increase beyond 0.6 orders of magnitude from the initial value over a period longer than 5 min, the algorithm now considers that this enhancement is indeed a bow shock crossing.
The final stage (Stage 3, Figure 2h) of the algorithm is implemented to define the bow shock crossing as a single point. To define this single point we first take the following 10 min of f data after the time step in which all three criteria above proved true (i.e., the transparent green box in Figure 2g that starts from the position of the blue circles in Figures 2e-2g). Next, we calculate the median f value within this 10 min of data (horizontal dashed red line in Figure 2h) and identify the bow shock location as the first point that exceeds this median f value (i.e., the vertical solid orange line in Figure 2h). As seen in the Figure 2h, the 10 min box is composed of both solar wind (low f ) and magnetosheath (high f ) data leading to the median value being situated at a point within the f enhancement (i.e., as MEX passes through the bow shock rather than at the inner and outer edges of the bow shock). An attempt was made to define the width of the bow shock by identifying when f within the 10 min time period exceeded two further thresholds. These thresholds were calculated from the 30th and 70th percentile values of f within the 10 min box, each representing the outer (toward solar wind) and inner (toward Mars) edge of the bow shock region, respectively. We have represented these points as two orange circles in Figure 2h, with the outer edge located near the start of the enhancement and the inner edge near the peak of the enhancement. While these extra locations could be used to define an uncertainty (or bow shock width) on each crossing, as is described later in section 3.3, we have opted to only use the more central bow shock location (i.e., the 50th percentile threshold) in our statistical analysis such that our results can be 10.1002/2016JA023316 compared more directly to other studies that have also opted to describe the bow shock location as a single point. Since we utilize the same algorithm for every MEX orbit, each crossing will be identified with the same criteria, making a single point identification suitable for a large number statistical analysis. In Figure 2h we have also superimposed the Edberg et al. [2008] bow shock model location (vertical dashed black line) which is almost colocated at the same time as our automatically identified bow shock location. Since the Edberg et al. [2008] method was based on visual identifications, this gives us confidence that our algorithm is identifying the bow shock in an appropriate location.
After identifying the time of the MEX bow shock crossing, the time stamp is then converted to a spatial coordinate using the NASA SPICE system [Acton, 1996]. Next, the algorithm skips ahead 5 min in time and continues running the above criteria to determine if there are any further bow shock crossings. Multiple crossings of the bow shock could be due to a temporal (e.g., boundary motion past spacecraft) or a spatial (e.g., MEX crossing the bow shock in multiple places due to orbital trajectory) variation. Unfortunately, as with all single-spacecraft observations, it is difficult, if not impossible to separate the two. The algorithm terminates when all data outside of the Edberg et al. [2008] MPB location shifted 0.1 R M away from Mars have been exhausted. This then sets an extreme inner limit to the bow shock location of which we rarely, if at all, expect the bow shock to reach. As mentioned earlier, the outbound bow shock crossings will be detected in the same way as described above by stepping through the data in the reverse direction, thus leading to the crossing once again appearing as an enhancement in f . Earlier we noted the use of a fully automatic (FA) and semiautomatic (SA) algorithm to detect the bow shock location. If the initial use of the Edberg et al. [2008] bow shock model suggests that no ASPERA-3 ELS data exists within the solar wind, we next check for the availability of any data outside of the Edberg et al. [2008] model MPB location. If any data are found, we run our SA algorithm. The SA algorithm runs exactly the same steps of the FA algorithm described in Figure 2, but we visually inspect each half orbit for the existence of solar wind data. If it is visually determined that a minimum of 5 min of solar wind ASPERA-3 ELS data exists, we input the time stamps for the range of data in which to run the bow shock identification algorithm (i.e., where solar wind data start and an endpoint either within the magnetosheath or induced magnetosphere) and an approximate time stamp at which solar wind data start to become magnetosheath like (used to calculate average solar wind flux proxy, i.e., Stage 1 of Figure 2). When running this algorithm on the entire ELS data set, the FA algorithm determined that there were 802 MEX half orbits that contained no ASPERA-3 ELS data outside the Edberg et al.
[2008] MPB (i.e., data entirely within induced magnetosphere), which were automatically excluded. Moreover, the FA algorithm determined that 2920 half orbits were expected to not have any data within the solar wind but did have data within the magnetosheath. When these 2920 MEX half orbits were visually inspected, it was found that 570 of them did in fact have viable data for running the SA algorithm.
The thresholds for the criteria used in the above algorithms were determined to be optimum via visual inspection of the results. A quantitative analysis of this is given in section 5, in which we compare a subset of automatically identified crossings to that of visually identified crossings. Although not shown, we have performed a sensitivity analysis on the threshold involved in criterion 2 by increasing and decreasing it by 0.05 and reapplying the algorithm to the full ELS data set. The outcomes of this sensitivity analysis are discussed in section 5, but, in general, our algorithm and results were found to be robust to small changes in the criterion threshold. The results presented throughout the remainder of this study are representative of the algorithm, criteria, and thresholds described throughout this section.

Bow Shock Model
To describe the Martian bow shock location and its variability, we describe each crossing in the aberrated Mars-centric Solar Orbital (MSO) coordinate system. In the traditional MSO system, +X MSO is directed toward the Sun, +Y MSO is in the opposite direction to the Mars orbital velocity vector, and +Z MSO is in the northern hemisphere, completing the right-hand cartesian set. The aberrated MSO system, (X ′ MSO , Y ′ MSO , Z ′ MSO ), is a 4 ∘ rotation about the Z MSO axis to account for the relative motion of Mars to the average solar wind direction. The average shape and location of the bow shock is then determined by fitting a conic section to the crossings in a similar way to previous studies [e.g., Russell, 1977;Slavin and Holzer, 1981;Slavin et al., 1991;Trotignon et al., 1991aTrotignon et al., , 1991bTrotignon et al., , 1993Trotignon et al., , 2006Vignes et al., 2000;Edberg et al., 2008], thereby producing a conic section of polar form r = L(1 + cos ) −1 , where r and are polar coordinates of each crossing with an origin at (x 0 , 0), and L, , and x 0 are the fit parameters corresponding to the conic sections semilatus rectum, eccentricity, and focus location along the X ′ MSO axis, respectively. The subsolar and terminator distance of the bow shock with such a model is determined from R SS = x 0 + L(1 + ) −1 and R TD = √ L 2 + ( 2 − 1)x 2 0 + 2 Lx 0 , respectively. As described in section 3.2, instead of using an intrinsic uncertainty of each crossing location (i.e., determined from a shock width estimation), we choose to define the uncertainty in our fitted parameters via a similar method used by Edberg et al. [2008]. This method varied each fit parameter individually from the best fit value until a 5% change in the root-mean-square deviation occurred, finally defining the uncertainty on each parameter as the difference between the new fit parameter value and the best fit value. In our estimates of the uncertainty, unlike Edberg et al. [2008], we increment and decrement upon the best fit parameter values to give us a negative and positive deviation. By using the same method our fitted model can be directly compared to the work of Edberg et al. [2008].
By applying our automatic algorithm to the entire ELS data set available (January 2004 until May 2015), we initially identified 12,091 bow shock crossings (11,309 and 782 from FA and SA algorithms, respectively). A total of 11,861 (11,098 and 763 from the FA and SA algorithms, respectively) of these crossings are shown in Figures 1d and 1e in the aberrated cylindrically symmetric MSO coordinate system, with X ′ MSO on the abscissa MSO on the ordinate. We show each individual crossing as a grey dot in Figure 1d and spatially grouped into bin sizes of 0.1 × 0.1 R M in Figure 1e to illustrate where the majority of events occur. The spatial binning in Figure 1e was used to remove false detections of the bow shock due to instrumental effects by neglecting any crossings in a bin with less than three crossings, resulting in 230 crossings being removed from the entire data set to produce the aforementioned 11,861 crossings. Additionally, the crossings shown in this paper (Figures 1d and 1e) include cases of multiple crossings per half orbit. Although not shown here, it was found that inclusion of multiple crossings per half orbit made at most a difference of ∼1.2% to the subsolar and terminator distances calculated from the resultant models. Due to this we choose to include them in the fitting of an overall average model. In both Figures 1d and 1e we have overlaid our resultant fitted conic section as a blue curve, with the best fit parameters and associated uncertainties given in Table 1. In section 3.2 we noted that the algorithm actually identified the bow shock as three points (two of which were estimates of the inner and outer edges of the bow shock). By fitting conic sections to the inner and outer edges of the bow shock (not shown), we verified that these resultant fits fell within the uncertainty of that presented in Table 1. Thus, including the width of the bow shock region was not needed in the remainder of this study. For comparison with previous studies, we have included the average subsolar (R SS =1.58 R M ) and terminator (R TD =2.60 R M ) bow shock distances of all previous models [Bertucci et al., 2011] as orange crosses in Figure 1d. In Figure 1d we have also included the MEX orbital trajectory (green curve) and bow shock crossings (inbound: black circle, outbound: red circle) of MEX for the electron spectra presented in Figure 1a.
Inspecting the individual crossings shown in Figure 1d, a large amount of variability in location is observed. Toward the subsolar point (i.e., low SZA) the bow shock is located closer to Mars and has smaller variability in position. As the SZA increases (i.e., moving around the flanks of Mars), the bow shock location both moves farther away from Mars, and its variability in location increases. The smaller variability at the subsolar point could be partially due to the limits of the MEX orbital geometry. This bow shock location dependence on SZA position can be removed by using our statistically defined model of the average bow shock location and extrapolating each crossing into the subsolar and terminator plane. This is done by keeping the best fit values of and x 0 constant but then varying L until the conic section passes through each individual crossing [Farris and Russell, 1994]. After this, each crossing will be defined by a set of three conic section parameters, allowing R SS and R TD to be calculated and compared. This method was used by Vignes et al. [2000], Crider et al. [2002], and Edberg et al. [2008and Edberg et al. [ , 2009band Edberg et al. [ , 2010. As with these previous studies, the MEX orbit limits the total number of crossings near the subsolar point; and therefore, we choose to use the mapped terminator distance in the remainder of this study of the variability of bow shock location.

Bow Shock Terminator Distance and Solar Parameter Variations With Solar Longitude
We present, as a function of the Mars solar longitude, L s , the mapped to terminator bow shock distances (R TD ) in Figure 3a; solar wind dynamic pressure determined from ASPERA-3 IMA solar wind density and velocity moments (P dyn ) in Figure 3b; the integrated solar EUV irradiance (I, 10-124 nm) measured by TIMED-SEE and extrapolated to Mars in Figure 3c; the TEC of the full atmosphere determined by MARSIS-SS measurements at a solar zenith angle of 85 ∘ in Figure 3d; and finally the Martian heliocentric distance (R h ) in Figure  The occurrence distributions of each quantity presented in Figure 3 show that N cross∕obs typically exceeds 100 observations per L s bin, allowing us to have confidence in our statistics. The large numbers of observations within each L s bin also contribute to the relatively small standard errors calculated for each average quantity. We note that more bow shock crossings have been identified within the range of L s =70 ∘ -210 ∘ , with the lowest number of crossings just after this range, and remaining low (but still with N cross >100) throughout perihelion. The average value of the mapped terminator distance is R TD =2.49 R M ± 0.11% (∼1% larger than our model value but comfortably within the uncertainty, see Table 1). Dividing the R TD trend (Figure 3a) into the two halves of the Martian year (i.e., northern spring/summer, L s < 180 ∘ , and northern autumn/winter, L s ≥ 180 ∘ ), we see that, aside from regions where the trend is similar to the full data set average, the first half of the year tends to have reductions in R TD , whereas the second half of the year tends to have enhancements in R TD . There is a prolonged reduction in R TD occurring around Mars' aphelion L s = 60 ∘ -120 ∘ , with a minimum value reached of R TD = 2.39 ± 0.01 R M (bin L s = 90 ∘ -100 ∘ ), a ∼4% reduction from the full data set average. There is a deeper minimum prior to this feature (R TD = 2.38 ± 0.02 R M , bin L s = 20 ∘ -30 ∘ ), but it is not prolonged and has a larger error on the mean value. For the enhancements in R TD , the largest and most prolonged feature occurs around Mars' perihelion (L s = 240 ∘ -300 ∘ ), with a maximum value reached of R TD = 2.65 ± 0.02 R M (bins L s = 260 ∘ -280 ∘ ), an ∼8% increase from the full data set average. From the minimum value in the reduction to the maximum value in the enhancement there is an overall increase in R TD of ∼11%. According to a Students t test, the minimum and maximum mean values of R TD per L s bin are deemed to be significantly different from the full data set average at the 99% confidence interval.
We also note in Figure 3a that there appears to be a general trend toward R TD increasing from L s ∼ 190 ∘ onward, but with a sudden large reduction toward the full data set average in bin L s = 230 ∘ -240 ∘ . Although not shown here, this has been identified to be due to the crossings within this bin originating from only two different Martian years (the other bins contain crossings from a minimum of four different Martian years) that are at vastly different points of the solar cycle (around solar cycle 23/24 minimum and 24 ascending/maximum in Martian years 29 and 31 [Clancy et al., 2000;Cantor et al., 2010], respectively). Removal of the crossings from Martian year 29 increased R TD for bin L s =230 ∘ -240 ∘ (without significantly impacting the rest of the trend), but resulted in this bin containing crossings from only a single Martian year that could be biased by as of yet to be investigated solar cycle effects.
The solar wind dynamic pressure is given by P dyn = SW v 2 SW , where SW = m SW n SW is the solar wind mass density (m SW and n SW are the average solar wind particle mass and the number density, respectively) and v SW is the solar wind velocity. Although m SW and v SW are not expected to vary with distance from the Sun, the density is expected to reduce by an inverse square relationship, and thus, on average, we expect P dyn to be larger the closer Mars is to the Sun. In general, for solar longitudes in which Mars passes through its aphelion passage, P dyn (see Figure 3b) tends to be lower than the full data set average of P dyn = 0.31 nPa (±1.22%). Whereas, for solar longitudes where Mars is approaching and passing through perihelion, P dyn tends to be above the full data set average. The maximum P dyn value past L s = 180 ∘ is ∼70% larger than the minimum value prior to L s =180 ∘ . However, the standard error on each average value of P dyn is somewhat larger than the other data sets in this figure. If we take the upper error limit of the minimum P dyn reached (P dyn = 0.23 ± 0.02 nPa), and the lower error limit of the maximum P dyn reached (P dyn = 0.39 ± 0.04 nPa), the percentage variation reduces from 70% to 40%. Clearly, the overall increase in P dyn as Mars moves from its aphelion to perihelion position shows HALL ET AL.
VARIATIONS IN MARS' BOW SHOCK DISTANCE 11,485 the expected variation in P dyn with distance to the Sun. However, the large differences in the total variation due to the uncertainties alone show that P dyn is highly variable throughout all periods of the Martian orbit.
The solar EUV irradiance, I (Figure 3c), and TEC within the Martian ionosphere (Figure 3d) both follow a sinusoidal variation, demonstrating a clear modulation by the heliocentric distance of Mars (Figure 3e). The solar EUV irradiance is known to modulate the TEC within the Martian ionosphere by driving ionization of the atmosphere, and thus, this TEC profile is as expected. As Mars approaches aphelion, both the EUV and TEC reduce toward their minimum values (I = 3.79 ± 0.04 mW m −2 , TEC = 0.36 ± 0.01 TECU (total electron content unit; 1 TECU = 10 16 el m −2 ), within bins L s = 80 ∘ and 90 ∘ , respectively), before gradually increasing as Mars moves toward perihelion (maxima of I = 5.56 ± 0.08 mW m −2 , TEC = 0.58 ± 0.02 TECU, within bins L s = 230 ∘ and 250 ∘ , respectively). The minimum to maximum percentage increase of the EUV and TEC are 46% and 61%, respectively. Edberg et al. [2008] noted that the bow shock was on average located farther out over the southern hemisphere of Mars, and that this could be due to the presence of the intense crustal magnetic fields. Thus, to test whether the variations seen in Figure 3a could be due to bias of crossings over a single hemisphere, in Figure 4 we present the R TD against L s trend for crossings occurring in the northern hemisphere (blue, planetographic latitudes > = 0 ∘ ) and southern hemisphere (red, planetographic latitudes <0 ∘ ). Figure 4a is presented in the same way as Figure 3a. The occurrence distributions in Figure 4a show us that while there are typically more than 100 crossings per L s bin irrespective of the hemisphere, there tends to be more crossings identified in the northern hemisphere for L s < 180 ∘ , while the opposite is true for L s ≥ 180 ∘ . The average southern hemisphere terminator distance (red dashed line, R TD =2.52 R M ± 0.14%) is ∼2.4% larger than that of the northern hemisphere (blue dashed line, R TD =2.46 R M ± 0.15%), presenting a small overall asymmetry between the hemispheres. In Figure 4b we present the fractional change in R TD with respect to the full data set averages for each hemisphere (i.e., the blue/red profile in Figure 4a divided by the corresponding blue/red horizontal dashed line). Figure 4b then allows us to see how R TD in each hemisphere varies on the same scale. In general, we observe that the variation in R TD with L s is similar irrespective of the hemisphere the crossing occurred in. Although, when comparing the trends from each hemisphere, the southern hemisphere terminator distance reduces to smaller values in the first part of Mars' orbit (L s < 180 ∘ ), whereas the northern hemisphere terminator distance reaches larger values earlier in the second part of Mars' orbit (L s < 180 ∘ ).

Influence of P dyn and Solar EUV on Bow Shock R TD
In general, from all panels in Figure 3, we see that the mapped R TD increases and decreases in line with similar deviations in the quantities P dyn , I, and TEC. Additionally, Figure 4 shows that the variations in R TD are not due HALL ET AL. VARIATIONS IN MARS' BOW SHOCK DISTANCE 11,486 Figure 5. Response of Martian bow shock terminator distance with solar parameters. (a) Terminator distance against solar EUV irradiance and further subdivided into two regimes of low (blue) and high (red) solar wind dynamic pressure.
(b) Terminator distance against solar wind dynamic pressure and further subdivided into two regimes of low (blue) and high (red) solar EUV irradiance. The solid lines/curves in both panels show the corresponding models fitted to the data, with the functional form and best fit parameters given in the panels and in Table 2. to the bow shock crossings being biased by a north/south hemisphere asymmetry. However, since the solar parameters are varying in the same way (i.e., reducing toward aphelion and increasing toward perihelion), and that we expect the bow shock to move toward Mars at higher P dyn but to also move away from Mars at higher solar EUV, it is not clear which of these is the dominant factor. To investigate which out of the solar EUV or dynamic pressure has the largest impact on the bow shock terminator distance, we matched (where possible) each of the 11861 bow shock crossings to a corresponding ASPERA-3 IMA dynamic pressure, P dyn,IMA , and TIMED-SEE extrapolated solar EUV irradiance I EUV . We used the requirement that, with respect to the time stamp of each crossing, the P dyn,IMA measurement must be within 1 h, and then the I EUV must be within 12 h (the large time threshold is due to daily averaged extrapolated values used). Finally, to reduce any biasing of the results by multiple crossings being matched to the same P dyn,IMA measurement, multiple crossings on a given orbit were reduced to a single R TD value by taking the mean. Due to the limitations of each data set, this process reduced the total number of crossings available to 5023, which is still significantly larger than any previous study.
Using these crossings, we present the variation in the average terminator distance, <R TD >, with I EUV (Figure 5a), and P dyn,IMA (Figure 5b). In both panels <R TD > is on the ordinate, with the value calculated from the mean R TD of all crossings within bins of size 0.05 W m −2 or nPa. The scatter points show the <R TD > and standard error on the mean within each I EUV /P dyn,IMA bin, with only bins with greater than 20 crossings used to reduce noisy data points. In an attempt to isolate the individual effects of solar EUV irradiance, and solar wind dynamic pressure on the bow shock terminator distance, the data have been further subdivided into two different regimes of P dyn,IMA (for <R TD > against I EUV , Figure 5a) and I EUV (for <R TD > against P dyn,IMA , Figure 5b). The regimes are defined as low dynamic pressures or EUVs (blue data points) for P dyn,IMA <0.34 nPa and I EUV <4.62 mW m −2 , respectively, with high dynamic pressures and EUVs (red data points) at larger values than these limits. The limits of the regimes are defined by the mean value of P dyn,IMA and I EUV , and although calculated from this reduced data set of crossings, the average value is similar to those calculated from the entire data sets of the respective values (see section 4.1 and Figure 3).
The response of R TD to solar EUV ( Figure 5a) presents a linear relationship of the form R TD = aI EUV + b. A model of this form was fitted to these data, and the best fit parameters along with 1 errors are given in Table 2.
In general, both pressure regimes are well described by the linear relationship, and within error limits, R TD increases at the same rate with solar EUV irrespective of the dynamic pressure. However, larger dynamic pressures do indeed tend to dampen the overall distance of R TD across all EUV values presented (i.e., the b fit parameter is smaller for larger dynamic pressures).
For the response of R TD to the solar wind dynamic pressure (Figure 5b), a power law relationship of the form R TD = P dyn,IMA was more suitable (as used by Verigin et al. [1993]). An exponential form similar to that used HALL ET AL. VARIATIONS IN MARS' BOW SHOCK DISTANCE 11,487 Low: 0 < P dyn,IMA (nPa) < 0.34 0.11 ± 10.5% 2.02 ± 2.7% 0.69 High: P dyn,IMA (nPa) ≥ 0.34 0.12 ± 6.0% 1.92 ± 1.7% 0.84 Solar Wind Dynamic Pressure: R TD = P dyn,IMA Low: 0 < I EUV (mW m −2 ) < 4.62 2.39 ± 0.3% −0.02 ± 11.3% 0.69 High: I EUV (mW m −2 ) ≥ 4.62 2.52 ± 0.4% −0.03 ± 11.7% 0.88 a r 2 is the calculated coefficient of determination for each fit. Unless otherwise stated all fit parameters are given to 2 decimal places (d.p.) by Edberg et al. [2009b] resulted in fits where the fit parameters had 1 errors larger than the actual best fit values. The fitted models to this data are also given in Table 2. As expected, we see that the R TD reduces as the solar wind dynamic pressure increases, with the rate of reduction being largest at very small pressures. Above P dyn,IMA ∼ 0.2 nPa the rate of decrease is much more steady. A clear division is seen between the solar EUV regimes, with periods of higher solar EUV (red) having resultant terminator distances that are always larger, and almost completely isolated from those at periods of lower solar EUV (blue). The fitted models also suggest that, when the solar wind dynamic pressure is very low, the bow shock terminator distance is more sensitive to increases in the dynamic pressure during periods of higher solar EUV. However, since the solar wind dynamic pressure and solar EUV tend to vary in phase with each other (i.e., both share an inverse square law relationship with distance from the Sun, we do not expect there to be many periods in which the dynamic pressure is exceptionally low and at the same time the solar EUV is high. This could be interpreted as a good assumption from the errors on the terminator distances at very low pressures for periods of high EUV being somewhat large. In general, the data and model fits presented in Figure 5 and Table 2 show that the solar EUV irradiation plays a much larger role in modulating the bow shock terminator distance, with the solar wind dynamic pressure acting to dampen its impact.

Discussion
The key outcomes of this study are as follows.
1. A Mars' bow shock model has been determined from 11861 crossings automatically identified across a period exceeding 5 Martian years using proxy measurements of the electron flux by the ASPERA-3 ELS instrument on board the MEX spacecraft (see sections 3 and 4, and Figures 1 and 2). 2. For the first time we have characterized the average variation of the bow shock position in terms of the Mars' orbit of the Sun, discovering that on average, the bow shock is closer to Mars around aphelion but farther away from Mars around perihelion (see sections 4.1 and 4.2, and Figures 3 and 4). 3. The impact, and main driving factor out of the solar wind dynamic pressure, and solar EUV irradiance on the bow shock terminator distance have been described, confirming that the solar EUV plays a larger role in the overall variations (see section 4.3 and Figure 5).
These three keypoints will now be discussed in turn. Since we have used an automated algorithm (see section 3) to identify our boundary crossings, it is important to note that there is the possibility that a number of crossings have simply been missed rather than not being present. Beyond visually inspecting our automatically identified results, we also attempted to quantify how effective our algorithm was compared to a manual identification. To do this we took a sample of 340 MEX orbits and visually identified the bow shock crossings (single and multiple per half orbit) as rapid enhancements in the electron flux. These manual results were then directly compared to crossings identified using our automatic algorithm (section 3.2) on the same subset of MEX orbits. The manual method resulted in 663 crossings, whereas the automatic method resulted in 530 crossings. Although it appears that the automatic method only identified 80% of the manual identifications (which is still a large proportion), when calculating the average number of crossings per unique orbit, the manual identification resulted in ∼2.1 crossings per orbit, whereas the automatic method had ∼1.9 crossings per orbit. Typically, we expect two crossings per orbit for the majority of the sample period, so we can see that rather than the automatic method missing 20% of the manually identified crossings, it was actually less sensitive to identifying multiple crossings within the same half orbit. Therefore, we believe that our automatic method has been effective in identifying the Martian bow shock, and although it will miss valid crossings, this is less than 20% when compared to identification via a manual method. We also note that our algorithm is less sensitive to quasi-parallel crossings of the bow shock surface, which typically manifests as slow increases in electron flux with time. Nevertheless, we have identified the largest number of bow shock crossings from a single spacecraft to date, allowing us to easily complete the analysis described throughout the paper.
As noted at the end of section 3.2, a sensitivity analysis on the threshold used in criterion 2 of our bow shock identification algorithm was completed. This criterion described how large an enhancement in the electron flux (within 5 min) had to be to trigger as a potential bow shock crossing. We found that as the threshold was increased the algorithm became stricter and started to miss valid bow shock crossings. Reducing the threshold initially lead to more crossings identified (∼10% increase in FA algorithm detected crossings). However, these extra crossings were found to be additional crossings identified within the same orbits. Thus, these crossings represent an increase in false detections. Nevertheless, the entire analysis presented from section 3.3 onward was also completed for the crossings identified using the slightly smaller and larger criterion 2 threshold. We have not presented the specific results in this paper, but no significant differences to those already shown throughout this study were identified. Due to this, the remainder of our discussion and the resulting conclusions will remain the same irrespective of the criterion 2 threshold used.
Compared to previous studies (see Tables 1 and 2 of Trotignon et al. [2006] and Edberg et al. [2008], respectively), our bow shock model best fit gives a different set of conic section parameters (see Figure 1d and Table 1), but the overall shape remains similar. To get an idea of how our model compares to those that have come before, it is easier to inspect the resultant subsolar and terminator distances. We find that our subsolar distance is ∼4% larger (1.65 R M cf. 1.58 R M ), and our terminator distance is ∼5% smaller (2.46 R M cf. 2.60 R M ) than the average of previous models (orange crosses in Figure 1d [Bertucci et al., 2011]), but these values remain within the error limits (see Table 1). Nevertheless, differences between models could be expected due to (1) the selection methods of each crossing, (2) the data sets used to identify the crossings, and (3) the time periods in which the crossings come from. Previous methods have typically opted for a visual identification of the bow shock crossings and will carry an inherent observer bias that is difficult to quantify. Our crossings were identified via an automated method (described in section 3) and will mitigate any observer bias. Previous models have also used bow shock crossings identified primarily using magnetic field observations (i.e., enhancement and increased variability in magnetic field magnitude [e.g., Bertucci et al., 2011]). Due to a lack of a magnetometer on MEX we have opted to use a proxy of the electron flux, and even though a bow shock crossing manifests in a similar way to that observed in magnetic field data, differences in location could be possible. Finally, our much larger data set provides crossings from a longer time span and all from the same spacecraft, thereby reducing and/or removing bias due to possible solar cycle modulation and cross-calibration errors when multiple spacecraft observations are used. We thus expect our model of the bow shock to be a much more holistic view of its average location. Now relating to the second and third key points, in Figure 3 we presented that the Mars' bow shock terminator distance varied throughout the orbit of the Sun, being situated closer to Mars around aphelion and farther away from Mars around perihelion. To help explain this we also presented the solar parameters of dynamic pressure and EUV irradiance in the same scheme, and both were seen to vary in the same way as R TD . An analysis of where the bow shock crossings occurred with respect to the northern and southern hemispheres of Mars was completed to eliminate any concern of the bow shock position at certain period of the Mars' orbit being biased by any intrinsic hemispherical bias (i.e., impact of the crustal magnetic fields in the southern hemisphere, see section 4.2 and Figure 4). Our results demonstrated that the southern bow shock terminator distance is located on average 2.4% higher than the northern bow shock terminator distance. This difference agrees with the findings of Edberg et al. [2009b]. Despite the authors noting their low coverage of bow shock crossings across the southern hemisphere of Mars, they used dayside crossings of the bow shock by MEX from February 2004 to January 2009 to find that the southern hemisphere bow shock terminator distance was on average ∼1.6% larger than the same quantity in the northern hemisphere. Aside from this preliminary agreement in these two sets of results, we leave a full analysis of this asymmetry to a later study and simply note that we observed no obvious hemispherical bias leading to the overall observed variation in R TD with solar 10.1002/2016JA023316 longitude. We next continued our study and investigated the individual impact of the solar EUV irradiation and solar wind dynamic pressure on the bow shock terminator distance.
Previous studies have shown that the Martian bow shock location varies with solar wind dynamic pressure Verigin et al., 1993;Edberg et al., 2009b] and EUV flux [Edberg et al., 2009b;Russell et al., 1992;Vignes et al., 2000Vignes et al., , 2002. In particular, Edberg et al. [2009b] modeled these variations and found that the R TD decreases via a weak exponential relationship with P dyn , whereas it increases via a strong exponential relationship with increasing EUV flux. Using our larger data set, we produced our own modeled relationships of the response of R TD to changing levels of solar EUV and P dyn (see section 4.3 and Figure 5). Instead of an exponential relationship between the bow shock location and solar EUV, we determined that a linear model was much more suited ( Figure 5a and Table 2). By further dividing the data into multiple pressure regimes, we identified that high dynamic pressures only act to weakly dampen any response of the bow shock due to changes in solar EUV. The linear relationship that we present is likely different to that of the Edberg et al.
[2009b] study due to the nature of their data, in which they use the F 10.7 radio flux measured at Earth and then extrapolated to Mars as a solar EUV proxy, whereas we use measured solar EUV flux (10-124 nm) at Earth and extrapolated to Mars. Our use of this range of solar EUV irradiance is more appropriate since the wide range of wavelengths caters for ionization of the neutral Martian atmosphere across a wide range of altitudes (in general, smaller wavelengths penetrate to lower altitudes of the atmosphere before interacting with a neutral). An additional reason to the difference between our functional forms could arise from the Edberg et al. [2009b] results being derived from crossings identified within a limited period of the solar cycle (declining and minima phase of solar cycle 23/24), which could bias their results toward lower EUV proxies. Now, when comparing our modeled variation in R TD with P dyn (Figure 5b and Table 2), we also see a difference in the form of the fit to that of the Edberg et al. [2009b] results. We found that a power law form (e.g., alike Verigin et al. [1993]) was more suitable than an exponential form. Despite this difference, the overall inference of the bow shock position moving towards Mars as the solar wind dynamic pressure increases remains the same. However, what is much clearer from our results is how differing levels of solar EUV irradiation have a much larger impact on where the bow shock will be located. We do note that our modeled responses will only give the gross response of the bow shock terminator distance to varying solar conditions, and, by themselves, cannot describe the entirety of the physical systems occurring within the Martian plasma system to produce the variations we see as Mars orbits the Sun.
Due to the lack of semicontinuous measurements of upstream IMF conditions by MEX, we have not been able to consider the impact the IMF field direction, magnitude, and solar wind magnetosonic Mach number (M MS ) could have on our observations. We expect that variations in the IMF field direction will have no correlation with the Mars orbit about the Sun and thus are highly unlikely to influence our findings. However, the IMF field strength varies with distance in a similar way to solar wind density (inverse square relationship) and solar EUV irradiance. The M MS value of the solar wind takes these values into account and typically increases with distance from the Sun and thus could be another factor that modulates yearly variations in the bow shock location. Edberg et al. [2010] (and references therein) studied this and found the Martian bow shock to behave similarly to that of Venus, that is, the bow shock would flare outward (R TD increases) as the Mach number decreases. Edberg et al. [2010] derived an inverse linear relationship between R TD and M MS , and using their relationship we would require a value of M MS = 9.1 for our smallest R TD just after aphelion (bin L s = 90 ∘ ) and M MS = 6.5 for our largest R TD (just after perihelion, bins L s = 260°− 280°). According to the distribution of M MS values expected at Mars during the Edberg et al. [2010] study (see Figure 2 of their paper), our required M MS values are from the tail end of their observed distribution. Although we do believe that M MS will play a part in the modulation of the bow shock location, it is likely unfeasible for such extrema of M MS to be a typical value consistently reached solely due to the position of Mars throughout its orbit of the Sun. It is more likely that such values of M MS could occur at any point throughout the Mars orbit, simply due to the intrinsic variability of the solar wind and IMF. In addition to our results presented in section 4.2 (Figure 4), since the location of the Martian crustal magnetic fields rotate with the planet, any impact on the bow shock location will be averaged, and consequently we do not consider them to drive our observations.
From the results in section 4 and the discussion of other studies above, it is likely that our observations are largely driven by the external conditions of solar EUV flux. This solar parameter can modulate both the internal condition of the ionosphere (TEC varies in same way as solar EUV, see Figures 3c and 3d) and have an impact on the amount of ionized particles within the Martian extended exosphere. At perihelion, the EUV flux is larger, leading to more ionization within the ionosphere and larger internal pressures , which could impact the pressures in outer layers of the Martian system (i.e., magnetosheath) acting to raise the bow shock location. Also during higher levels of EUV flux, we might expect both a more heated and expanded Martian exosphere, as well increased ionization of exospheric neutrals. The newborn exospheric ions are electromagnetically susceptive to the solar wind, which could then lead to a higher rate of ion pickup and mass loading of the solar wind causing its flow to slow and consequently raising the location of the bow shock. A study by Yamauchi et al. [2015] looked at seasonal variations of Martian pickup ions using MEX ASPERA-3 IMA measurements. They identified evidence of the exosphere breathing over a Martian year, with a higher frequency of pickup ions observed when Mars is closer to the Sun. These results are in agreement with our observations. We also note that the enhancement of R TD at perihelion from the full data set average is much larger than the reduction seen at aphelion, and the solar EUV ( Figure 3c) follows a similar pattern. Finally, although we note that the largest increase in R TD is observed just after perihelion, similar to the EUV trend we actually see a general upward trend in R TD from L s bin 90 ∘ , which is not consistently similar (i.e., several successive L s bins) to the full data set average until L s bin 290 ∘ . This range of L s is similar to that of a period typically known as the Martian dust storm season [e.g., Farrell and Desch, 2001]. In this season, as Mars approaches the Sun, increased solar irradiance and temperatures lead to a warmer and more circulated Martian atmosphere, which then lifts dust off the surface of Mars and causes the so-called dust storms. These dust storms have been previously shown to interact with the upper atmosphere and ionosphere of Mars in various ways, including the possibility of modulating the neutral densities and TEC [e.g., Wang and Nielsen, 2004;Mouginot et al., 2008;Liemohn et al., 2012;Withers, 2009;Withers and Pratt, 2013;Withers et al., 2015;Xu et al., 2014;Haider et al., 2015;Němec et al., 2015]. However, we do not draw any further conclusions on how the dust storms could directly impact the location of the Martian bow shock and leave such an investigation to a future study. It is perhaps likely that it is not any single one of these mechanisms that are leading to our observations, but rather a combined effect of all of them, and at this point we cannot exclude any of them.
One method of validating the above hypothesis is via comparison with numerical simulations of the Martian plasma system. Modolo et al. [2005Modolo et al. [ , 2006Modolo et al. [ , 2016 presented results from a series of 3-D hybrid simulations (ions treated kinetically, electrons as a massless charge-neutralizing fluid) of the Martian plasma environment at two different levels of solar EUV flux (descriptive of solar minimum and maximum). The authors found the bow shock location to either be invariant with respect to solar EUV, or, in the case of Modolo et al. [2006], that the difference in the subsolar and terminator distance at the two different levels of solar EUV was smaller than their simulation cell size of 0.09 R M . Our results and those of Edberg et al. [2009b] both clearly show that the bow shock position at the terminator varies with solar EUV conditions. This could suggest that the inherent assumptions and simplifications that must be made in order to simulate such a complex plasma system in 3-D and at high resolution may lead to these contradictory results. Nevertheless, continued efforts should be made to understand how the Martian plasma system responds to these changing external conditions.

Conclusion
We have presented a statistical analysis of the variation in the Martian bow shock location throughout the Mars' orbit of the Sun. A conic section was fitted to 11,861 bow shock crossings identified using measurements from the MEX ASPERA-3 ELS instrument. The totality of the crossings come from a period exceeding 5 Martian years allowing for multiple chances to sample the bow shock location at each stage of the Martian year. Also, since this period spans approximately a full solar cycle, which inherently leads to different solar conditions, this study has allowed for a general view of how the Martian bow shock varies throughout the period. Each crossing has been mapped to the terminator plane of Mars in order to remove the bow shock location dependence on zenith angle, thus being able to directly compare each crossing to each other. We have presented this mapped terminator distance along with measurements of external (dynamic pressure and solar EUV) and internal (EUV and TEC) factors expected to drive bow shock variability, all as a function of solar longitude. Clear prolonged reductions (i.e., bow shock moving toward Mars) and enhancements (i.e., bow shock moving away from Mars) have been observed around the periods of aphelion and perihelion, respectively. Such excursions have been discussed in terms of the average solar conditions throughout a Martian year, and models of how the bow shock terminator distance varies with solar EUV irradiance and solar wind dynamic pressure have been produced to aid the explanation of these observations. At the terminator, the bow shock location increases linearly with increasing solar EUV irradiance, whereas it reduces via a power law relationship with solar wind dynamic pressure. We find that the variation in solar EUV flux impacting the Martian plasma environment with heliocentric distance is likely the major driving factor in this phenomenon. At closer proximity to the Sun, enhanced levels of solar EUV can increase the number of exospheric ions, which through pickup processes, can act to slow the solar wind flow, and can raise the bow shock location. Additionally, higher EUV leads to enhanced TEC and internal pressure within the ionosphere, potentially forming a more effective barrier to the solar wind flow. Since the dust storm season is also commonly thought to be related to higher levels of solar irradiation (which is in turn related to proximity to Sun), the period of larger mapped terminator distances (approximately Mars perihelion) could be correlated with the Martian dust storm season. However, it is unclear whether these are two related or independent phenomena. It is clear from this study that the Martian bow shock and its plasma system warrant further study in terms of its variability with external and internal factors. In particular, this study has completed an important step toward further work into identifying any long-term variations in the Martian bow shock location resulting from the ∼11 year solar cycle. In addition to the long-term statistics now achievable due to the success of the continuing MEX mission, the comparatively young MAVEN mission now affords us with continuous measurements of the magnetic conditions at Mars allowing for further studies of variability of the bow shock with magnetosonic Mach number and IMF conditions.