A fast bow shock location predictor-estimator from 2D and 3D analytical models: Application to Mars and the MAVEN mission

We present fast algorithms to automatically estimate the statistical position of the bow shock from spacecraft data, using existing analytical two-dimensional (2D) and three-dimensional (3D) models of the shock surface. We derive expressions of the standoff distances in 2D and 3D and of the normal to the bow shock at any given point on it. Two simple bow shock detection algorithms are constructed, one solely based on a geometrical predictor from existing models, the other using this predicted position to further refine it with the help of magnetometer data, an instrument flown on many planetary missions. Both empirical techniques are applicable to any planetary environment with a defined shock structure. Applied to the Martian environment and the NASA/MAVEN mission, the predicted shock position is on average within 0.15 planetary radius $R_p$ of the bow shock crossing. Using the predictor-corrector algorithm, this estimate is further refined to within a few minutes of the true crossing $\approx$0.05 $R_p$). Between 2014 and 2021, we detect 14,929 clear bow shock crossings, predominantly quasi-perpendicular. Thanks to 2D conic and 3D quadratic fits, we investigate the variability of the shock surface with respect to Mars Years (MY), solar longitude (Ls) and solar EUV flux levels. Although asymmetry in $Y$ and $Z$ Mars Solar Orbital coordinates is on average small, we show that for MY32 and MY35, Ls = [135-225$^\circ$] and high solar flux, it can become particularly noticeable, and is superimposed to the usual North-South asymmetry due in part to the presence of crustal magnetic fields.

More advanced physics-based models have also been proposed as a complement to those empirical attempts. A good introduction into analytical models of the bow shock, based on gasdynamic theory and magnetohydrodynamics solutions, is given in Verigin, Slavin, Szabo, Gombosi, et al. (2003) and recently in Kotova et al. (2021). These studies present analytical functions describing the curvature, bluntness and skewing angle of the shock structure, which are arguably better suited to the fitting of the shock flanks; they are applicable to many planetary bow shock conditions. At Mars, due to the sparsity of early data and the non-collisional nature of the shock, the tendency has been to use in priority the simplest fitting model available with least free parameters, that is, an empirical 2D polar equation (Russell, 1977;Slavin et al., 1991;Trotignon et al., 2006;Edberg et al., 2008;Hall et al., 2019). Only recently with the NASA/Mars Atmospheric and Volatile EvolutioN (MAVEN) mission were quadratic fits used to characterise the general structure of the Martian bow shock, with Gruesbeck et al. (2018) providing fits to a careful subset of identified crossings in the first year of operations of the MAVEN mission.
In recent years, many studies have attempted to characterise the Martian shock position and shape and its evolution under various solar wind and EUV conditions  and references therein). Two missions have been used for this goal, the ESA/Mars Express mission and the NASA/MAVEN mission. Mars Express (hereafter MEX for brevity) was launched in 2003 and has been orbiting Mars since 2004, whereas MAVEN was launched ten years later in November 2013 and has been orbiting the planet since 22 September 2014. MAVEN's scientific payload includes among others a fluxgate magnetometer (MAG, Connerney et al., 2015), two ion spectrometers including the Solar Wind Ion Analyzer (SWIA, Halekas et al., 2015) and the Suprathermal and Thermal Ion Composition instrument (STATIC, McFadden et al., 2015), and an electron spectrometer (Solar Wind Electron Analyzer, SWEA, Mitchell et al., 2016). MEX unfortunately does not carry any magnetometer but includes a plasma suite (ions and electrons) as part of the ASPERA-3 package , which was used to investigate the plasma boundaries at Mars (Dubinin et al., 2006). Both missions aim at studying the upper atmosphere and the magnetospheric environment of Mars.
Detections of the bow shock in spacecraft data have relied on manual determinations using as many instruments (including plasma instruments and magnetometer) as available to avoid -3-manuscript submitted to JGR: Space Physics ambiguous detections (Gruesbeck et al., 2018). Recently, Nemec et al. (2020) proposed a region identification scheme based on selected plasma parameters and applied their technique to the MAVEN dataset in order to identify upstream solar wind, magnetosheath and magnetosphere regions crossed by the spacecraft. This method has the advantage of mapping these regions statistically, removing certain biases usually associated with manually picked individual boundary crossings which may be orbit-dependent. However, it requires the reliable knowledge of flow speed, ion density and magnetic field, which may not all always be available. A parallel trend has also been to apply machine-learning techniques to plasma data for the labelling and identification of the regions crossed by a given spacecraft (see Hall et al., 2016;Breuillard et al., 2020, for Mars and Earth, respectively) as part of online databases (Génot et al., 2021). However, these studies require extensive amounts of time and patience, as well as large computer resources to become efficient. Sometimes, a precise determination of the shock position in the data is not of paramount importance and simpler, faster, more straightforward approaches, such as the one presented in this study, can be advantageous. This may be the case in statistical studies where one of those regions needs to be systematically excluded, as in space weather databases monitoring the solar wind. This can also be of interest when areas around the predicted bow shock must be excluded for qualified reasons, as in wave studies focusing on regions outside of foreshock and shock wake structures, or when a first guess of the location and geometry of the shock is needed.
We present in this study new, simple analytical algorithms using two types of historical fitting techniques (2D and 3D), which make it possible to quickly estimate from spacecraft spatial coordinates the statistical geometrical position of the shock in planetary atmospheres.
Special emphasis on the Martian environment and the MAVEN dataset is given throughout, however the method is applicable to any planetary environment and spacecraft dataset. This first crude estimator can be refined further by applying additional criteria, for example on the magnetic field amplitude measured by the MAVEN spacecraft. This provides a fast means to approximately and quite reliably identify the position of the shock so that solar wind and magnetosheath regions can be studied on a statistical level in the data. Moreover, other characteristics of the shock crossing, such as the quasi-parallel (q ) or quasi-perpendicular (q ⊥ ) geometry of the shock can easily be obtained by deriving the perpendicular direction to the shock at any point on the surface.
After a review of 2D and 3D bow shock fitting models at Mars in Section 2, we present in a consistent manner the leading equations behind these models in 2D and 3D, give the an-alytical expressions for standoff distances, and propose a geometric calculation of the normal to the bow shock at any point on the surface. Starting in Section 3, we introduce a predictor algorithm for a fast estimation of the shock position in spacecraft orbital coordinates and its timing (Section 3.1). In Section 3.2, we propose a simple correction on this position and timing with the sole help of magnetometer data (predictor-corrector algorithm). Application to the MAVEN MAG dataset is then given as validation on a few examples and then extended to the whole available dataset. Finally, as a result of the automatic detection proposed here, statistical analytical fits are given for the MAVEN mission between November 2014 and February 2021, with a discussion of the shock's asymmetry based on terminator and standoff distances. Applications for space weather-related databases are also mentioned.

Bow shock models at Mars
In this section, following a survey of existing fitting models at Mars, we present comprehensive formulae for analytical fits in 2D polar coordinates and 3D Cartesian coordinates, with a calculation of bow shock subsolar and terminator standoff distances. We also show how to calculate the normal to the surface at a given point in space, in order to estimate the q ⊥ − q shock conditions.

Coordinate systems and solar wind flow aberration
All spacecraft coordinates in this study are in Mars Solar Orbital coordinates (MSO) for simplicity, in accordance with most previous studies. In the MSO system, identical to the Sunstate coordinate system, the +X MSO axis points towards the Sun from the planet's centre, +Z MSO is towards Mars' North pole and perpendicular to the orbital plane defined as the X MSO -Y MSO plane passing through the centre of Mars, and Y MSO completes the orthogonal system.
Because of the orbital motion of Mars with respect to the average direction of the solar wind, the apparent direction of the solar wind in the rest frame of the planet deviates from the anti-sunward direction. As a result, an anti-clockwise rotation by an angle α around the Z axis must be applied so that the bow shock's major axis is aligned with respect to the X axis along the solar wind flow. This aberration, first seen in cometary tails and at the origin of the hypothesis by Biermann of a stellar wind (Biermann, 1951), is taken into account in the socalled aberrated MSO coordinates, denoted X MSO , Y MSO and Z MSO (although Z is left unchanged by the transformation). To unclutter notations, the 'MSO' subscript is now dropped. Follow- ing Formisano et al. (1979), Slavin and Holzer (1981) define the angle α as α = tan −1 V p / sw where V p is Mars' orbital velocity and sw is the solar wind velocity, for example expressed in km s −1 . Mars' average orbital velocity is V p = 24.07 26.5 22.0 km s −1 . For the maximum value (26.5 km s −1 ), the angle is α = 3.8 • for a typical solar wind speed of 400 km s −1 . The angle assumed by all studies except those of Slavin and Holzer (1981) and Slavin et al. (1991) is 4 • . In Slavin and Holzer (1981), the aberration angle was chosen to be varying with solar wind speed conditions. In Slavin et al. (1991), α = 3.2 • . Figure 1 shows the aberration angle with respect to orbital velocity (abscissa) and to solar wind velocity (colour code). Bow shock models at Mars have been proposed since the end of the 1970s, including (but not limited to): Russell (1977); Slavin and Holzer (1981); Slavin et al. (1991); Schwingenschuh et al. (1990); Trotignon et al. (1991);   Verigin et al. (1999); Vignes et al. (2000); Trotignon et al. (2006); Edberg et al. (2008Edberg et al. ( , 2010.

Parametric models
These studies were performed with several spacecraft including Viking, Mars Global Surveyor (MGS) and Phobos-2, and for varying solar conditions. Gringauz et al. (1976) quoted in Russell (1977 reported 11 crossings for the Russian Mars-2,3,5 satellites, but Slavin and Holzer (1981) later reanalysed the datasets and found 14 crossings in total. Slavin et al. (1991) reported 94 crossings for Phobos-2, upped to 126 by Trotignon et al. (1993). Using the same dataset, Verigin et al. (1993) reported a weak dependence of terminator bow shock position on solar wind dynamic pressure P sw ; additionally in an analogy with magnetised planets regarding the dependence of the magnetopause thickness to P sw , these authors anticipated the discovery of crustal magnetic field sources at Mars, later vindicated with MGS measurements (Acuña et al., 1998).
In contrast with what was found at Venus, Slavin and Holzer (1981) and later Vignes et al. (2000later Vignes et al. ( , 2002 (using MGS and data from previous missions) suggested that the mean bow shock standoff distance was likely independent of the solar activity. Slavin et al. (1991) showed that the terminator distance, which is a marker of the swelling of the cavity flanks, varied by as much as 11% between the Mars-2,3,5 observations (low activity) and the Phobos-2 observations (high activity), although the number of crossings for each mission largely differed. Mars EXpress (MEX), with its very long activity spanning the end of solar cycle 23 and cycle 24 up to now (towards the beginning of new cycle 25), had the best chance to conclusively solve this aspect: Hall et al. (2019) found that for the years 2004-2017, the terminator distance underwent variations up to ∼ 7%, in agreement with the preliminary results of Trotignon et al. (1993). Of note, Mazelle et al. (2004) made a review of all available measurements before MEX started observing, and discussed the solar cycle variations and the differences observed with Venus (for which up to 35% increases of the bow shock location in the terminator plane with increasing activity have been reported, see Russell et al., 1988;Zhang et al., 2008;Edberg et al., 2010).
Of great import, Edberg et al. (2009Edberg et al. ( , 2010 used MGS and MEX data in combination with ACE data extrapolated to Mars to study the dependence of the bow shock location to solar EUV flux and magnetosonic Mach number (noted M ms ). They pointed out that the shape of the magnetosonic shock wave depends on the ratio of the solar wind speed to the magnetosonic speed. Later, for the entire period 2004 identified 11, 861 crossings in the MEX database, using electron spectrometer data with machine-learning algorithms.
This study was extended by Hall et al. (2019) up to 2017 with a total of 13, 585 crossings, totalling 13 years of operations of MEX. Both works used a standard 2D conic fit depending on the Mars Year (MY), with observed variations up to a few percent in terms of standoff bow shock distance. In addition to solar wind cycle variations and hemispherical changes from a 2D perspective, these studies confirmed the dependence of the shock position with P sw and, most drastically, with solar EUV flux. Correspondingly, Ramstad et al. (2017) studied with coincident electron and ion data a subset of only 1, 083 inbound and outbound MEX orbits for the period 2005 − 2016. They evaluated the dependence of the Induced Magnetospheric Boundary (IMB) and bow shock (BS) to EUV flux and the solar wind's lowest moments (density, bulk velocity), showing that the BS mostly expands and contracts with the IMB. However, they also found that the BS swelling in the flank due to increased EUV fluxes cannot be solely explained by a corresponding swelling of the IMB. Simultaneously with MAVEN (both magnetometer and ion measurements), Halekas et al. (2017) investigated how the Martian magnetosphere and bow shock responded to EUV flux, M ms and P sw between October 2014 and May 2016 (0.85 Mars year). In agreement with previous studies, they showed that the shock inflates with increasing EUV flux and contracts with increasing dynamic pressure and M ms ; this in turn leads to EUV flux and dynamic pressure competing against one another because of their common 1/d 2 h dependence on heliocentric distance d h .
Recently, Nemec et al. (2020) used MAVEN plasma and magnetometer data to construct maps of the solar wind/magnetosheath regions. Their modelled bow shock locations, explicitly dependent on P sw , EUV flux and crustal field intensities, were in good agreement with the average fits of Trotignon et al. (2006), with appreciable differences in terminator extensions as compared to the results of Ramstad et al. (2017).
Broadly speaking, two approaches fitting the shape of bow shocks have historically been employed, one using a simple 2D polar form (e.g. Slavin & Holzer, 1981), the other the 3D general Cartesian conic form (e.g. Formisano et al., 1979;Formisano, 1979). Semi-empirical models based on gasdynamic and MHD predictions such as those of Verigin et al. (1993Verigin et al. ( , 1999 are not discussed in the following.
This range of angles depends on the nature of the conic section: for a parabola ( = 1) θ ∈ ] − π, π[ (borders excluded), for an ellipse ( < 1) θ ∈] − π, π], and for a hyperbola ( > 1) . Typical 2D conic bow shock shape in the aberrated (X , √ Y 2 + Z 2 ) coordinate system. For a point M on the shock surface, ρ is the Euclidean distance to the shock from the centre of the planet of radius R p , and r is the distance to the shock surface from the focus x F of the conic with semilatus rectum L and making an angle θ with the X direction, so that Equation (1) holds. ϑ is the usual polar angle, with respect to the centre of Mars. R ss and R td are the standoff subsolar and terminator distances.
(Cartesian) form of this equation is (Trotignon et al., 2006): In this representation, Trotignon et al. (2006) derived two additional useful quantities, the standoff shock distance along the X axis, R ss (also called subsolar aerocentric distance in Trotignon et al., 2006) and the standoff terminator distance R td along the cylindrical coordinate (which is none other than the diameter of the tail at X = 0 divided by 2, or as in Russell, 1977, the "dawn radius"): We can derive another parameter of interest, that is, the aperture of the Mach cone related to the shock structure -the limiting Mach cone angle (Verigin, Slavin, Szabo, Kotova, & Gombosi, 2003). In the gasdynamic approach, it is defined as = arcsin (1/M s ), where M s = SW / s is the sonic Mach number. The sonic speed is s = γP/ρ, with ρ the solar wind ion mass density, γ = 5/3 the ratio of specific heats, and P = n sw k B (T e + T i ) the solar wind thermal pressure, with T e and T i the electron and ion temperatures, respectively. For a hyperbola ( > 1), the limiting Mach cone angle is exactly the angle made by the asymptotes of the hyperbola (Slavin et al., 1984): with : as uncertainty. In a canonical form for the hyperbola, with a the distance from the nose to the intersection of the asymptotes on the X axis, and b that from the shock nose to the asymptote on the Y axis, arctan (b/a). Since, by definition = 1 + b 2 /a 2 , the substitution readily yields expression (7). It is noteworthy to remark that for close to 1, the uncertainty increases to infinity; any determination of is thus unreliable for quasi-parabolic curves.
For the fit, Slavin and Holzer (1981) and Slavin et al. (1991) rewrote Equation (1) as y = ax + b (posing y = 1/r, x = cos θ, a = /L and b = 1/L) and performed simple linear regressions for a range of foci locations. As pointed out by Vignes et al. (2000), this may result in fitting biases when observations are widely disparate in their location: in this case direct fitting methods to Equation (1) should be preferred. With a direct polar fit to MGS data, Edberg et al. (2008) gave for example the following fitted values: = 1.05±0.04, L = 2.10± 0.09, x F = 0.55 ± 0.12. However, to match the results plotted in Figure 1 of Edberg et al. (2008), the value of must be modified down to = 1.03, a marginal difference likely due to rounding errors. For comparison, the corresponding values derived by Hall et al. (2019), for MEX data but with a larger sample, are = 0.998 ± 0.003, L = 1.802 ± 0.002 and x F = 0.76. For MGS and MEX data (Edberg et al., 2008;Hall et al., 2016Hall et al., , 2019, the subsolar standoff distance is R ss = 1.63 ± 0.04 R p , whereas the terminator standoff distance is R td = 2.50 ± 0.09 R p . By comparison, Halekas et al. (2017) found large variations of the bow shock standoff distances in the early MAVEN data, with R ss ∼ 1.6-1.9 R p and R td ∼ 2.5-3.1 R p depending on EUV flux levels and combined with either M ms or the solar wind dynamic pressure.   (7) in the case of a hyperbolic shape. The mean value for each mission is also given, with Mars 2-3-5 and Mariner 4 (Slavin et al., 1991), Phobos 2 (Trotignon et al., 1993), Mars Global Surveyor (MGS) (Edberg et al., 2008) and Mars Express (MEX) .
The planetary radius of Mars is by definition R p = 3389.5 km.
Russell ( a Because the Mars 2, 3 and 5 measurements reported by Gringauz et al. (1976) (in total 11 crossings) did not specify local times, aberration angle α was assumed to be zero. b These authors use the full definition of the aberration angle, resulting in α = arctan V p / sw , in contrast to the more recent studies. See Section 2.1. c Mariner 4, and Mars 2, 3, 5 data only here. Uncertainties on fitted values assumed to be 1%. d "Direct fit" method with all 3 parameters varying simultaneously. e "Slavin's method", using a linear regression in (1/r,cos θ) space. f Note that = 1.03 matches better with Figure 1 of Edberg et al. (2008), for which the Mach cone aperture would instead be = 13.9. g Fits were performed on 2D-gridded density data, co-depending on M ms and EUV flux levels on the one hand, and solar wind dynamic pressure and EUV flux levels on the other. The coordinate system adopted by Halekas et al. (2017) was the Mars Solar Electric (MSE) system, with the X axis lying anti-parallel to the solar wind flow. h Ramstad et al. (2017) use the following rectangular function (required to be cylindrically symmetric with respect to the solar wind direction): with ρ the radial distance to the bow shock on the Y axis from the centre of Mars, R ss the subsolar standoff bow shock distance on the X axis, and ς the so-called scale length. This function is valid ∀x R ss since ρ(y = 0) = R ss . By definition, ς is a constant equal to 33.54 R p derived in Ramstad et al. (2017) from the bow shock model values for R ss , R td and of Vignes et al. (2000) and can be calculated as ς = − 1 2Rss R 2 ss − R 2 td /( 2 − 1) . The original values of R ss and in their study were fitted to a function a n b sw ( sw /100) c + d; we have assumed here nominal conditions (n sw , sw ) = (2 cm −3 , 400 km s −1 ) for simplicity. We calculate the semilatus rectum as L = R 2 td − ( 2 − 1)R 2 ss /(2R ss ) from formulae (5) and (6). Uncertainties on fitted values assumed to be 1%. i Here we only recall the results for all MYs (MY27 − 33). Individual MYs have eccentricities below 1 (ellipse), except for MY28 − 29 (hyperbola). j That is, Mariner 4 and Mars 2-3-5 (Slavin et al., 1991), Phobos 2 (Trotignon et al., 1993), MGS (Edberg et al., 2008), MEX  and MAVEN (Halekas et al., 2017). The listed uncertainties are the standard deviations of the series. Accordingly mean angles are calculated only for 3 values and are only given for for completeness here. k Although this is a hyperbola with cone angle = 8.1 • , the large eccentricity uncertainty leads to a cone angle uncertainty of 44 • , hence no value is provided here. l Calculated from formula (5). m Calculated from formula (6).
We present in Table 2 the fitted conic parameters of the main quoted references in Table 1, in chronological order. It is interesting to remark that most shapes fitted are stricto sensu hyperbolic ( ≥ 1), but in practice can be considered quasi-parabolic as eccentricity ∼ 1, which makes it possible to calculate the limiting Mach cone angle. For example, when fitting bow shocks from different MYs, Hall et al. (2019) showed that eccentricities varied around = 1 by less than 5% between MY27 and MY33, with a marked tendency towards ellipsoidal shapes (only 2 consecutive years, MY28 and MY29, had eccentricities above 1).

3D Cartesian form
The more general way of characterising the bow shock shape does not assume any symmetry with respect to any axis. A 3D shape model can be constructed in the form of a quadratic equation (for example, Formisano et al., 1979;Gruesbeck et al., 2018;Simon Wedlund et al., 2017, for Earth, Mars and comets): Here and for clarity in the equations, ( Although quadratic surfaces are not necessarily centered on the planet nor is their main axis directed along the X MSO axis (see Appendix A), a simple estimator of the shock's position in the subsolar and terminator directions can be of interest. We propose here such an expression, based on Equation (9). We calculate thus the subsolar standoff distance along the X MSO axis at coordinates (x, y = 0, z = 0) by finding the positive root of the simplified quadratic equation (i.e. the intersection of the surface with the X MSO axis): whereas the terminator standoff distances in the Y−Z plane (non-aberrated MSO coordinates) are similarly given by: Because of the small angles involved at Mars, non-aberrated coordinates are rather accurate for the subsolar standoff distance. That said, rotating the MSO coordinate system by a certain small angle α about the Z axis does impact the terminator distances by a few 0.01 R p . For the parameters given above, R ss = 1.56 R p , R td,z = 2.50 R p and R td,y = 2.62 R p . Because the shock is a 3D object, the exact position of the tip of the ellipsoid may vary with respect to the values taken at the origin.

Comparison of historical models
A comparison of a representative selection of historical bow shock models (some of them as listed in Vignes et al., 2000;Trotignon et al., 2006) in the X −Y /Z plane is given in Figure 3. The 3D quadratic model of Gruesbeck et al. (2018) was rotated anticlockwise by 4 • around the Z axis to ease the comparison. It is noteworthy to remark that although the fit is not valid for X −0.5 R p , the figure displays the fits for X > −2.5 R p to illustrate the differences in shock surface swelling.
All models are in excellent agreement around the subsolar point, with a mean subsolar standoff distance value of R ss = 1.59 ± 0.05 R p . The terminator standoff distance is also in very good agreement -however, for X 0 R p , the difference between fits becomes substantial, especially (i) between the recent MEX investigations of Hall et al. (2016Hall et al. ( , 2019 and all other fits on the one hand, and (ii) between the MAVEN fits and the other fits on the other.
For MAVEN, this is due, as previously mentioned, to the lack of orbital coverage by the spacecraft for X < 0 R p . In this sense, MEX has a much better antisolar spatial coverage. Of note, the North-South asymmetry in the 3D fits of Gruesbeck et al. (2018), coinciding with the presence of crustal fields in the Southern hemisphere of Mars, can easily be seen, a characteristic which no axisymmetric model directly quantifies.

Quasi-perpendicular or quasi-parallel shock?
A collisionless shock may have different behaviours depending on the upstream solar wind magnetic field (the IMF), which conditions how the solar wind is losing its energy to the magnetosheath. Two main cases are conveniently studied for their varying properties: q and q ⊥ shocks. Additional important physical quantities driving the shock structure and dynamics are the magnetosonic Mach number (which defines the shock's criticality) and the plasma- It is useful to recall that a q shock condition is defined so that the background IMF lines are intersecting normally the shock surface, whereas a q ⊥ shock describes an IMF that is in effect in the tangent plane to the surface shock. Thus, the angle of importance is the angle between the average IMF vector upstream of the shock and the shock normal. This angle is in the literature almost always named θ Bn , which is kept here for convenience. The geometry of the shock is defined as follows: Starting in the magnetic field compression region in the solar wind, q ⊥ shocks have structures, from the point of view of B-fields, almost always characterised by (i) a foot, (ii) a fast ramp, and (iii) a wider overshoot followed by a more gradual undershoot (see Kennel et al., 1985, Figure 11). This classic picture is a first approximation as fine electron-scale structures in the foreshock, foot and ramp can be seen with high-cadence magnetic field measurements. Q ⊥ shocks reflect particles back upstream to satisfy the shock conditions and are on average diffusive. Magnetic structures trapping particles such as mirror modes are observed to predominantly take place in the magnetosheath behind a q ⊥ shock (Gary, 1992 In order to determine the q or q ⊥ geometry of the shock crossed by a spacecraft, the normal direction to the shock surface needs to be first estimated. For a single spacecraft, this can be done either with methods that take advantage of upstream and downstream magnetic field measurements (coplanarity method as in Schwartz, 1998;Horbury et al., 2002, although prone to rather large uncertainties) or through geometrical considerations only, as we propose in Section 2.3.2. The accuracy of the geometrical method presented here is linked to the assumption that the shock surface is smooth and does not possess any kinks or local structures where the current curls on itself. In practice, we do not expect such a smooth surface as the shock may assume a more rippled shape which depends on the upstream solar wind condition and the turbulence at the boundary (Moullard et al., 2006). However, our geometric determination may still be a useful first approximation of the geometry of the shock.

Determination of the shock normal
The normal to the shock surface at point (r 0 , θ 0 ) in polar coordinates, (x 0 , y 0 , z 0 ) in Cartesian coordinates or (r 0 , ϑ 0 , ϕ 0 ) in spherical coordinates, is simply defined as the gradient vector of the (assumed) smooth surface f at that point. Mathematically we can express this condition as: where v is a vector tangential to the surface at that point. This leads to the following expressions in the 2D and 3D cases.
2D case. For the 2D polar coordinate case, let f be equal to f (r, θ) = r−L/(1+ cos θ) following Equation (1) where θ is the angle from the focus x F on the X axis. The gradient of f depends on the two variables (r, θ): At point (r 0 , θ 0 ) vector ∇ f = (R 0 , Θ 0 ) is perpendicular to the surface. Note that because of the peculiarity of a conic, values in x must always be corrected by the focus distance x F , because the typical polar angle ϑ is not strictly the same as the conic angle θ used in Equation (1) (see Figure 2).
In Cartesian coordinates, using Equation (4), the gradient will be against directions along X and Y and equal to: This expression circumvents the ambiguity on the angle direction of the polar formula, and as such should be preferred when calculating the normal direction. Figure 4a displays our estimates of the normal direction to several points on the shock surface applied to the 2D bow shock polar fit of Edberg et al. (2008), and converted to Cartesian coordinates.
3D case. For the 3D Cartesian quadratic equation, f is simply equal to the left member of equation (9).
The gradient of f is then: is perpendicular to the surface. Additionally, the equation of the tangent plane to the smooth surface at that point is of the general form: where subscript "0" in the gradient components denotes the gradient taken at points (x 0 , y 0 , z 0 ) for brevity. Thus, knowing the 3D position of the spacecraft at the expected bow shock position, one can calculate the transverse and tangent directions to the bow shock surface. With the knowledge of the normal direction to the shock at the spacecraft location, it becomes possible, from the 2D and 3D model cases, to calculate the angle θ Bn from the average direction of the magnetic field at the spacecraft position near the shock, noted B bg : For added robustness, the function arctan2 is recommended for the calculation of the inverse tangent, as it returns a value corresponding to the correct quadrant of the Euclidean plane. Because of the inherent 3D nature of a spacecraft orbit and of the local magnetic field, we prefer the 3D calculation with Equation (19) over the corresponding 2D case (Equation 18). It is however important to recall here that the shock's local shape may be assuming that of a "corrugated iron" section, as evidenced for example at Earth with the Cluster quartet of spacecraft (Moullard et al., 2006). No method is foolproof in estimating θ Bn : in case studies, the local normal to the shock can be more carefully checked, for which several complementary methods exist, such as the magnetic coplanarity method mentioned earlier.
When applied to the MAVEN dataset, the geometric calculation of θ Bn , assuming a smooth surface, is expected to reach an uncertainty of about ±5 • depending on the upstream field determination. We obtained this estimate for a small sample of crossings, by extending over a few minutes the time spans used to calculate the upstream magnetic field direction.

Detecting the bow shock in spacecraft orbits
Estimating the bow shock position from spacecraft spatial coordinates can be achieved either empirically or theoretically, depending on the precision needed. Semi-empirical but computationally intensive techniques using machine-learning imaging algorithms are currently attempted to detect automatically and precisely the exact position of the shock. Such technique makes use of the full plasma instrumental payload on board planetary missions, when available. Other techniques like that of Nemec et al. (2020) use ion and magnetic field data in combination to statistically identify the regions crossed by the spacecraft. However, a faster approach, based on a simple geometrical estimator using a static analytical bow shock model (see Section 2), may still prove valuable for statistical studies or for new datasets. We present such an approach and its possible refinements in 2D and 3D coordinate systems based on magnetic field-only measurements. Because in the following we do not take into account plasma measurements, and because the true signature of the shock may be difficult to detect with magnetic field data only, additional criteria to identify solar wind and magnetosheath regions are required to mitigate this ambiguity. These criteria are presented in Section 3.2.

Predictor algorithm for the shock position from existing analytical models
The predictor algorithm is based on the calculation of polar (2D, θ) or spherical (3D, [ϑ, ϕ]) angles in the corresponding frame of reference centred on the planet in MSO coordinates. These angles unequivocally define the predicted distance to the bow shock at the position of the spacecraft. Comparing this bow shock distance with the Euclidean spacecraft dis--18-manuscript submitted to JGR: Space Physics tance to the centre of the planet gives access to the region in which the spacecraft is located either in the solar wind or in the magnetosheath.
2D case. Our algorithm is (see Figure 2 for definitions of angle and distances): • Choose conic model r(θ), with eccentricity , semilatus rectum L and focus' position (x F , 0, 0), • Calculate the spacecraft's Euclidean distance r sc from the chosen conic model focus x F , in aberrated MSO coordinates, so that: r sc = (X sc − x F ) 2 + Y 2 sc + Z 2 sc , • Calculate the angle θ sc at the position of the spacecraft: θ sc = arctan2 Y 2 sc + Z 2 sc , X sc − x F , • Calculate the predicted bow shock distance R bs at the corresponding spacecraft θ sc from the focus x F : R bs = r(θ sc ) = L/(1 + cos θ sc ) following Equation (1), • Calculate ∆R = r sc −R bs . If ∆R goes from negative to positive (respectively, from positive to negative) values, the spacecraft is expected to move from the magnetosheath to the solar wind (respectively, from solar wind to magnetosheath). At the temporal precision of the spacecraft ephemerides, the closest value to ∆R = 0 defines the shock position (X bs , Y bs , Z bs ) and crossing time t bs .
This purely geometrical approach was tested for polar coordinate models (such as Edberg et al., 2008;Hall et al., 2016Hall et al., , 2019, provided that all spacecraft coordinates are first rotated 4 • into the aberrated MSO system. 3D case. For 3D quadric models (Gruesbeck et al., 2018), the aberration is already taken into account and there is no need to correct the spacecraft coordinates for the position of the focus of the conic. Thus we only need generalise the approach above to spherical coordinates (ρ, ϑ, ϕ), where ρ = √ X 2 + Y 2 + Z 2 is the planetocentric distance, whereas ϑ = arctan (Y/X) and ϕ = arctan Z/ √ X 2 + Y 2 represent azimuth and elevation (measured from the X-Y plane) by convention. To compensate for the inherent ambiguity on azimuth depending on the quadrant and gain robustness, the function arctan2 is preferred throughout for simplicity, in a programming sense. Equation (9) becomes a second-degree equation of the form: with: a = A cos 2 ϕ cos 2 ϑ + B cos 2 ϕ sin 2 ϑ + C sin 2 ϕ + D 2 cos 2 ϕ sin 2ϑ + E 2 sin 2ϕ sin ϑ + F 2 sin 2ϕ cos ϑ and: b = G cos ϕ cos ϑ + H cos ϕ sin ϑ + I sin ϕ In the case of an ellipsoid of revolution (as it is the case for the parametrisation of Gruesbeck et al., 2018), the bow shock distance R bs at the angles (ϑ,ϕ) corresponds to the positive root of this equation: • Calculate the spacecraft's Euclidean distance in non-aberrated coordinates, R sc = X 2 sc + Y 2 sc + Z 2 sc , • Calculate (azimuth, elevation) angles (ϑ, ϕ) at position of the spacecraft: ϑ sc = arctan2 (Y sc , X sc ) and ϕ sc = arctan2 Z sc , X 2 sc + Y 2 sc , • Calculate the bow shock distance R bs at the corresponding spacecraft spherical angles with Equation (23), • Calculate ∆R = R sc − R bs . If ∆R goes from negative to positive (respectively, from positive to negative) values, the spacecraft is expected to move from the magnetosheath to the solar wind (respectively, from solar wind to magnetosheath). At the temporal pre-cision of the spacecraft ephemerides, the closest value to ∆R = 0 defines the shock position (X bs , Y bs , Z bs ) and crossing time t bs .
Application to MAVEN orbits. Using 1-min-averaged MAVEN orbits, we perform the automatic detection of the bow shock as shown in Figure 5. Because of the relatively poor temporal resolution of this dataset, as well as the fast approach in the early stages of the orbit insertion, some points in the orbit yield false positive detections which disappear when increasing the orbital data resolution to 1 s.
Thanks to the simple algorithms presented above, we can statistically predict bow shock crossings in a given spacecraft orbit. To help identify the solar wind region, we can distinguish between trajectories moving from the magnetosheath to the solar wind region, and vice-versa.
For each orbit intersecting the bow shock model, two points per orbit will be identified. At 1 s resolution, we predicted a total of 16, 515 bow shock crossings using the 3D analytical model of Gruesbeck et al. (2018) for the MAVEN dataset between November 2014 and February 2021, including 8, 256 crossings from the sheath to the solar wind and 8, 259 crossings from the solar wind to the sheath.
Actual crossings may in practice be different than the predictions and the algorithms may fail to pinpoint the location of the shock sometimes by several tens of minutes, mostly because of solar wind varying conditions, or due to the geometry of the shock at the point of passage for an individual orbit (Halekas et al., 2017). We estimate thus the precision of these automatic estimates to be of the order of ±0.08 R p (±270 km) around the 'true' bow shock location, from a representative subset of data. Because of variable shock position from orbit to orbit and the geometric average nature of the detection, some orbits that may have experienced shock crossings but lie inside the average shock position will not be tested for potential detection. It is estimated that only a few hundred potential crossings were ignored in the process. Consequently, the true shock structure location should be checked directly in the magnetic field and plasma data. Moreover, since the bow shock is a dynamical object, it may experience fast forward and backward motions, crossing the spacecraft trajectory several times more per orbit. This can be seen for example in Figure 6c, where the total magnetic field undergoes sharp intermittent jumps in the foreshock area (23 June 2018, around 01:00 UT). As in Gruesbeck et al. (2018), these multiple crossings, usually the first one in the temporal sequence for crossings into the sheath, and the last one for crossings into the solar wind, are identified as one with our algorithm.
We successfully applied this 3D algorithm to the retrieval of undisturbed solar wind density and velocity moments in the MAVEN/SWIA data with 1 min resolution, as part of the He-lio4Cast solar wind in-situ data catalogue, enabling the statistical study of interplanetary coronal mass ejections and high speed streams (Möstl et al., 2020).
3.2 Refining the position of the shock: a predictor-corrector algorithm Because of variations in the shock position, the automatic detection may give inaccurate predictions. We present here a fast method to correct to a certain extent for these discrepancies. It makes use of the magnitude of B to identify the position of the shock structure, either when crossing from the magnetosheath into the solar wind or vice-versa. As before, the main assumption is that the shock is crossed twice per orbit at maximum, although in practice the shock structure may be crossed several times due to the fast motion of the boundary across the spacecraft trajectory (Halekas et al., 2017). As we are interested in the statistical position of the shock, this assumption nonetheless provides a valuable estimate of the average position of the shock during those times. All magnetic field data are assumed here to be of the order of 1-s resolution.
From the point of view of a single spacecraft's magnetic field measurements across the shock boundary, the total magnetic field intensity increases sharply, from typically 5-10 nT in the solar wind at Mars to about twice that level on average when moving into the magnetosheath. Additionally, fluctuations increase, going from small standard deviations around the mean in the solar wind, to comparatively larger fluctuations in the magnetosheath. Typical crossings of the Martian shock illustrating these behaviours are shown for example in Figure 6 which presents examples of bow shock crossings around Mars as seen with the MAG instrument on board MAVEN throughout the mission. As discussed in Section 2.3, q ⊥ shocks usually display distinct features (foot, sharp ramp, overshoot) as in Figure 6a,b,d), whereas q || geometries have less clear signatures in magnetic field data, making it difficult to detect reliably (Figure 6c). Consequently, the corrector method presented here is biased towards the detection of q ⊥ crossings.
Our predictor-corrector detection algorithm attempts to consistently identify solar wind undisturbed regions in the chosen dataset, characterised by relatively low magnetic field intensities in combination with small fluctuations. For each orbit of a spacecraft, it proceeds as follows, with user-defined values (∆t bs , B th , T , ς th ) discussed afterwards: 1. Calculate predictor estimate of the shock crossing time t bs with the automatic algorithm in 2D aberrated polar coordinates or in 3D (as in Section 3.1).
2. If t bs exists in the considered orbit, choose a small time interval symmetric around the estimated shock so that [t bs − ∆t bs , t bs + ∆t bs ] (user-defined), 3. Calculate a robust estimate of the average magnetic field, e.g. the median of the magnetic field in half of this interval, noted |B sw | 1/2 . By definition |B sw | 1/2 corresponds to the assumed solar wind region magnetic field. For crossings from the solar wind into the sheath, the first half-interval is selected. For crossings from the sheath to the solar wind, the second half-interval is selected. There are two possibilities at this junction: (a) If |B sw | 1/2 > B th , where B th is the solar wind-to-magnetosheath threshold (user-defined), the position of the solar wind region is difficult to assess. In that case: i. Make the time interval float around the estimated location of the shock, by increments of ∆t bs /3 in one direction or the other, so that t bs = t bs ± ∆t bs /3, ii. Repeat interval shift until |B sw | 1/2 ≤ B th or until a maximum shift of 2∆t bs is effected from the estimated shock timing in either direction. If |B sw | 1/2 > B th still, either the spacecraft is always in the sheath or in an usually high solar wind B-field region, in which case the crossing is altogether ignored and removed from the database. If not, go to step (b): (b) If |B sw | 1/2 ≤ B th , this half-interval is a good candidate for undisturbed solar wind conditions. In that case: i. Calculate the running Median Absolute Deviation (MAD) ς mad,B of the total Bfield signal over a temporal window of duration T (user-defined) in the chosen half-interval, and smooth further the result with a running median over a time span enclosing the shock structure in its entirety (for example 2T or 3T ). This also helps remove potentially abrupt but temporally isolated changes in the signal. Note that this particular choice of ς mad,B is somewhat arbitrary. After several tests including running standard deviations, normalised or not to the "solar wind" signal, the choice of a smoothed running MAD was empirically found to work consistently well with the MAVEN dataset at 1 s resolution. For all times t at which the total magnetic field B tot is measured over a running interval [t i , t i + T ]: ii. Compare ς mad,B to threshold value ς th (user-defined). A jump above a certain threshold ς th indicates a transition between a less turbulent region (ς mad,B < ς th ) to a more turbulent one (ς mad,B > ς th ), an indicator of the presence of a shocklike structure. If this threshold is reached, take the first (respectively, last) time this happens in the chosen interval for solar wind-to-sheath (respectively, sheathto-solar wind) crossings, and correct the original timing t bs of step 1 to new t bs,corr .
If not, discard crossing.
(c) Repeat for each orbit.
At Mars, we tested step 1 (predictor) in the previous section with the MAVEN dataset: on average, the detected shock was within ±0.08 R p (±270 km) of the true shock crossing, corresponding to about ±30 min of data along the orbit. We thus set the interval of study for the corrector algorithm to ∆t bs = 30 min. This value depends on the orbit inclination with respect to the bow shock surface and is thus mission-dependent.
We determined the typical user-defined threshold values for the MAVEN mission manually, for simplicity, on a reduced dataset. We found a good compromise by trial and error with B th = 11 nT, because the undisturbed solar wind magnetic field in Mars' vicinity is of the order of 2-6 nT on average (Slavin & Holzer, 1981), but can reach up to about 10 nT or more when solar transient effects such as coronal mass ejections or co-rotating interaction regions are involved (Liu et al., 2021). In future studies, a more dynamic criterion in step (3a) may be preferred, i.e. where the amplitude of the magnetic field is normalised to the assumed upstream solar wind value. The criterion for being in the magnetosheath could for example become γ = |B|/|B sw | 1/2 > γ th , where γ th is an adequately chosen threshold (γ ∼ 1.5 for a clear increase of magnetic field when moving into the magnetosheath, with nominal solar wind levels γ ∼ 1).
Because the shock appears in measurements as a turbulent structure whereas the solar wind is on average less so, step (3b-i) calculates a measure of the variability of the magnetic field in the vicinity of the shock. The value of T is also mission-and instrument-dependent; for MAVEN/MAG data at a resolution of 1 s, a value T = 120 s was chosen, which adequately captures the mean magnetic field variations. When in the solar wind, we found that the intrin-sic variability of the signal ς mad,B , calculated over a running window of duration T = 120 s, was on average less than 0.5, which is adopted as the threshold ς th . This makes it possible to detect the very first perturbations in the solar wind leading to the creation of the shock structure, a point which is identified here as the position of the shock proper at time t bs . Again, in future studies, the threshold can be normalised to the magnetic field level, as in Halekas et al. (2017) where, together with constraints on plasma parameters, the normalised root-sum-squared value of the magnetic field was chosen so that RSS(B)/|B| < 0.15 to identify undisturbed solar wind intervals.
Application to MAVEN dataset. Applied to the MAVEN dataset at 1-s orbital and magnetic field resolution, the corrector algorithm reaches an accuracy of ±0.02 R p (±70 km) around the "true" shock (manually picked on a reduced dataset for comparison), a factor 4 increase in accuracy with respect to step 1. In the temporal datasets, this corresponds to only a few minutes of continuous data along MAVEN's orbit. This is epitomised in Figure 6, The final corrected timings for the detections yield with this algorithm a lower estimate of the total actual number of crossings encountered by a spacecraft throughout its mission. Events occurring when |B sw | 1/2 > B th , even after shifting the temporal window significantly, or when ς mad,B < ς th , were discarded in the final selection as can be seen in Figure 6. They may indicate that the magnetic field was either too turbulent or too complex in its structure (e.g. multiple crossings as is regularly the case with q crossings) for the corrector algorithm to capture. Moreover, as discussed previously, the analytical approximation model used for the determination in step 1 (either 2D or 3D) is likely to underestimate the true number of crossings due to the planetary bow shock variability (Halekas et al., 2017). Because our study is primarily interested in the statistical position of the bow shock throughout the mission, this loss of potential detections may be compensated by the large number of orbits of the considered spacecraft. only ∼ 1% of the q shocks are highly parallel ones (θ Bn ≤ 10 • , 30 crossings). This is shown in Figure 7, where the highest number of crossings occurs for θ Bn ∼ 80 • . This is also in qualitative agreement with the results of Vignes et al. (2002) for the MGS mission, when they investigated a proportion of 93 q ⊥ shocks for only 23 q shocks.  (Thiemann et al., 2017). The median of the EUV flux in the 2014-2021 period is 0.0028 W/m 2 and defines two EUV flux levels, one "high" for fluxes above that limit, one "low" for fluxes below.

Application to the study of the Martian bow shock variability
In Section 4.1, we perform 2D and 3D fits to the found bow shock positions using the spacecraft ephemerides. A discussion of these fits and what they imply is given in Section 4.2.

Statistical position of the Martian bow shock
2D case. In 2D MSO aberrated coordinates, the polar equation, Equation (1), can be rewritten in the linear form y = ax + b (see Trotignon et al., 2006): with a linear regression in the (r, X − x F ) space performed for a chosen focus location x F .
First we chose a focus location randomly between 0 and 1 R p , and for each linear fit performed, the residuals are calculated. The adopted focus point is the one that minimises the residuals.
Because MAVEN's orbits are not suited to bow shock detections for X < −0.5 R p , additional constraints on the tail distributions are necessary to obtain a more realistic conic fit. This can fits are essentially the same around the nose of the shock downstream to about −1 R p where patent differences start to appear. Because of the added cloud of ghost points, this is expected and thus 2D fits presented below are only valid in practice in the range [−0.5 R p , R ss ]. Incidentally, differences on the terminator and subsolar standoff distances are less than < 2% in each case. It is noteworthy to add that the determination of the nature of the conic section found from the fits can be significantly altered when using tail models either from E08 (hyperbola) or from H19 (fits' nature given by these authors, ranging from ellipse to hyperbola, depend on the Mars year considered): in that case, the fit's nature will naturally be biased towards matching that of their respective parent tail model.   (5) and (6). For hyperbolae, the Mach cone aperture is also given as calculated by Equation (7). For each fit, the coefficient of determination R 2 gives a measure of the goodness of the linear regression. Due to the large data spread, uncertainties on R ss and R td are of the order of 5% and of the order of 2% for the other quantities. 3D case. In 3D, we perform quadric fits using the method first put forward by Taubin (1991), adapted to the quadratic surface of Equation (9). This fitting method constructs scatter matrices from local gradients S of tested model T and finds the diagonal matrix of the generalised eigenvalue problem so that Tv = λSv, where v is the generalised eigenvector of T and S, and λ are the eigenvalues. Because of the scatter of points in the database, uncertainties on the found parameters A to I are of the order of 1%, in a least-squares sense. Table 4 collects all 3D fit parameters for each case; all fitted surfaces are ellipsoids of revolution. For completeness, we present and give the physical interpretation of these parameters in Appendix A, in terms of principal axes, their direction and lengths and the centering of the ellipsoids. Figure 9 shows the corresponding fits and their sections onto the X MSO −Y MSO Table 4. Martian bow shock 3D conic parameters from quadric surface fits applied to the MAVEN orbits and magnetic field data (predictor-corrector algorithm). See Equation (9) for the definition of parameters A to I and Equations (11), (12) and (13) for those of the subsolar standoff distance along the X MSO axis and the terminator standoff distances along the Y MSO and Z MSO axes. Uncertainties on the parameters are of the order of 1% in a least squares sense. All quadrics below are ellipsoids. The domain of validity for each fit is shown in Figure 9: fits are valid for X MSO ≥ −0.5 R p on average. The number of fitting points used for each case is the same as for the 2D fits, see Table 3 (last column). Also, see Appendix A for a physical interpretation of the tabulated parameters.
Case  Because of the spacecraft changing orbits during the mission, some of the ellipsoid fits appear anomalous in their orientation. This is especially obvious for MY35 when MAVEN, as of Mars 2020, decreased its apogee to ∼ 4500 km and hence its revolution period to 3.5 h to accommodate Mars 2020 rover operations on the ground. Consequently, MAVEN only seldom explored regions below X < 0 R p for half of MY35: this makes it difficult to constrain the fit, and we end up with an ellipsoid having its longest principal axis unphysically tilted almost 90 • in the X-Z plane (panel a, purple curves of Figure 9). A similar issue is found for Ls = 135-225 • (panel b, yellow, Figure 9), which has the lowest number of detections among Ls ranges and for which the orbit was never favourable for detections due to orbit precession.
Thus, in these two cases, no physical interpretation should be drawn from the axes orientations of the ellipsoid and the fit should only be valid for near-subsolar crossings. More robust physics-based analytical models could be used to overcome these fitting issues (Kotova et al., 2021).

Discussion
Our 2D and 3D fits give some insight on how the Martian bow shock is moving globally for different conditions of Mars Year, Ls and EUV flux and complements previous studies with the MAVEN mission (Halekas et al., 2017;Gruesbeck et al., 2018;Nemec et al., 2020).
As the bow shock position is connected to the balance between thermal pressure from the plasma in the ionosphere and the dynamic pressure from the solar wind, any variation of these two quantities will have repercussions on the position of the shock.
When assuming axisymmetry around the aberrated axis X in the 2D polar rectangular coordinates case, Table 2 and the average subsolar and terminator distances can be a first guide for our interpretation. Our new results with MAVEN (all points, Table 3) agree rather well with past measurements (Table 2) considering the data spread and estimated uncertainties: R ss = 1.74 ± 0.09 compared to 1.61 ± 0.08 R p and R td = 2.46 ± 0.13 compared to 2.56 ± 0.20 R p .
More specifically, for: • Mars Years: the subsolar standoff distance decreases by as much as 10% between MY32 and MY33-MY34, from 1.87 to 1.69 R p , although some of this variation may be stem- conditions, the bow shock appears to expand in the subsolar direction by about 7−13% from its summer and autumn position (R ss ≥ 1.72 R p compared to R ss ≈ 1.61 R p ). Simultaneously, the area encompassed by the bow shock conic is also increased during those two instances. One possible driving factor behind these changes may be in turn linked to changes in Mars' dayside upper atmosphere and extended exosphere, and how they expand and contract with seasons, increasing or decreasing the size of the obstacle to the solar wind flow (Hall et al., 2016, and references therein). globally increasing the ionisation rates in the ionosphere and through photoionisation of the extended exosphere as well as heating up and expanding the neutral atmosphereexosphere system (Forbes et al., 2008;Edberg et al., 2009;Hall et al., 2016). Photoionisation of exospheric neutrals creates newly born ions that are picked up by the solar wind convective electric field, resulting in mass-loading and slowing down of the solar wind flow (Yamauchi et al., 2015, with the presence of pickup ions in the foreshock region). Such combined effects have been shown to expand the bow shock in the solar wind direction (Mazelle et al., 2004). The two fits we present here, one for higher and one for lower EUV fluxes (more than 6500 points each), display the expected behaviour, with a larger standoff distance by 7% and a noticeably larger flaring of the fitted conic for the higher EUV fluxes (terminator distances increasing from 2.41 to 2.54 R p , i.e. 5%).
• Shock conditions: q and q ⊥ bow shock crossings are related to the average interplanetary magnetic field's (IMF) direction and the spacecraft's orbit (more precisely, the spherical quadrant in which the spacecraft emerges into the solar wind). Because the predictorcorrector algorithm favours q ⊥ detections (Section 3.2), the statistics between the two cases is heavily unbalanced (see Vignes et al., 2002, for a similar result). On average, we find no significant difference between the two conditions, with the shock surface slightly contracting and flaring up in q ⊥ conditions with respect to q conditions (|∆R ss,td | ∼ 4%).
Such a tendency is marginal considering that these percentages are at the precision limit obtained with the fits.
Let us now look at our 3D fit results. Figure 9 clearly shows several asymmetries depending on the Mars year, Ls, EUV flux and shock condition. The usual pronounced North-South asymmetry (X MSO − Z MSO plane, second column, and also Y MSO − Z MSO plane, third column), mostly ascribed to the presence of crustal magnetic fields in the southern hemisphere (Gruesbeck et al., 2018), is clearly seen for all cases with the standoff subsolar distances being skewed towards that hemisphere. This is shown by crosses representing the tip of the pro- A comparison between all standoff distances, subsolar and terminator alike, and calculated by our 2D and 3D algorithms is shown in Figure 10, and based on Tables 3 and 4. The standoff distances calculated from the 3D fits show the same tendencies as their 2D counterparts, although since the X and Y coordinates are not solar-wind aberrated and hence no axisymmetry is considered, the comparison between the 2D and 3D cases can only be that of general trends. From MY32 to MY35, a general decrease of standoff distances can be seen.
This result is arguably in contrast to those presented in Vignes et al. (2002) although our statistics with MAVEN is much larger than in their study. In an identical way to the 2D fits, larger EUV fluxes result in a bow shock surface significantly expanding in the solar wind towards the subsolar direction. With respect to the geometry of the shock, the subsolar standoff distances R ss appear to marginally increase from q ⊥ to q conditions, although the inverse trend is seen for the terminator distances. Again, these differences are slight, which may reflect in part the bias against q conditions of our bow shock estimator (thus yielding a low amount of q shock detections).
As a preliminary conclusion, we note that: • The X MSO -Y MSO and X MSO -Z MSO asymmetry seems particularly marked for Ls = [135 • -225 • ] (labelled "Ls3" on the figure), MY32, MY35 and higher EUV fluxes: it can readily be seen by comparing the length of the blue and red bars. As explained earlier, the number of points used for fits for MY32 is the lowest of all the cases because MAVEN arrived at Mars late in MY32. This is in qualitative agreement with the conclusions of Gruesbeck et al. (2018). On average and outside of those special cases, the shock's shape stays rather symmetric about the X MSO axis: the terminator distances R tdz (3D fits) and R td (2D fits) indeed seem to match rather well most of the time. This axisymmetric tendency can be further amplified by aligning the X MSO -Y MSO plane with the solar wind aberration system, rotating the 3D quadric surface 4 degrees anticlockwise around the Z MSO axis; new standoff distances for the 3D fits (R tdz and R ss ) differ by less than 5% with their corresponding 2D fits values (not shown).
• Although the 3D and 2D conic fits retain strong similarities in their behaviour, the 3D fits (seemingly paradoxically) appear more robust and less affected by external assumptions. It is recalled here that not only do the 2D fits assume axisymmetry around X MSO , but certain 2D fits had to also be constrained at larger Euclidean distances from the centre of the planet due to the poor coverage of MAVEN for X MSO < 0.5 R p . This superiority of the 3D fitting algorithm is due to: (i) the number of fitting variables (A to I, allowing more flexibility despite risking over-determination of the linear system of equations), (ii) the natural asymmetry of the shock (albeit small), and (iii) the fitting points being statistically better distributed over a larger space (both in X MSO -Y MSO and X MSO -Z MSO planes instead of a single polar plane) and thus optimising the fits.
-34-manuscript submitted to JGR: Space Physics • Because the Martian seasons (monitored by Ls ranges) to a degree and the EUV flux both depend on Mars' heliocentric distance, correlations between these fits are to be expected. For example, similar fits for low Ls values (< 135 • ) and low EUV flux can be seen in Figure 8 (panels c and e, orange curves) and 9 (panels b and c, red curves).
• Because the solar cycle is a continuous underlying driver of the shock's position regardless of the binnings adopted here , we expect also correlations between EUV flux and Mars year results. This effect is most clearly exemplified with the shock standoff decrease when going from the declining phase of solar cycle 24 (MY32 and MY33, high fluxes) to the next solar minimum (MY34, low fluxes).

Conclusions
In this study, we presented a fast method to estimate automatically the position of the bow shock in real spacecraft orbits, as well as analytical expressions for the normal direction to the shock surface at any point in its close vicinity. After a survey of existing analytical smooth models of the bow shock surface at the planet Mars based on 2D and 3D fits, we used these models as a first prediction of the shock position in the data and refined this prediction further with a predictor-corrector algorithm based on the median absolute deviation of the magnetic field around the predicted shock. This method, biased towards the detection of q ⊥ shocks but not entirely limited to them, does not substitute for a detailed analysis of the crossing or for machine-learning techniques currently developed for space missions. It however finds a useful application when it is necessary to quickly determine the position of the spacecraft, or at least an estimate thereof, with respect to the bow shock.
As part of the solar wind and space weather database Helio4Cast, our technique was successfully used to retrieve solar wind undisturbed parameters from the MAVEN mission (Möstl et al., 2020). We also successfully applied the predictor-corrector method to the MAVEN orbit and magnetic field data between November 2014 and February 2021 (see list compiled in Simon Wedlund et al., 2021), and performed a series of fits, in 2D and in 3D, to test our method and investigate statistically the shape of the shock depending on Mars Year, solar longitude Ls, and two solar EUV flux levels. The 3D fitting has obvious advantages over the 2D polar axisymmetric geometry usually used to describe the shock structure, namely, a more accurate estimate of asymmetries in the global structure, and taking full advantage of the 3D distribution of bow shock detections in space. This is especially important for bodies such as Mars with large orbital eccentricities and axial tilts to the ecliptic, and for which the heliocentric distance is a strong driver of the EUV flux input and seasonal changes on the planet.
Expectedly, we found the Martian shock to be highly asymmetric with respect to the North-South hemispheres, in agreement with previous studies (see for example Hall et al., 2016;Halekas et al., 2017;Gruesbeck et al., 2018). Such an asymmetry is in part linked to the presence of crustal magnetic fields at Mars; however, a specific study taking into account the planet's rotation and the location of these crustal magnetic sources on the nightside or on the dayside is left for the future. Bow shock fits for quasi-perpendicular and parallel shock conditions were, to the precision of our approach, almost identical. In addition, the shock appeared noticeably asymmetric with respect to Y MSO and Z MSO directions in specific conditions, namely, for MY32 and MY35, Ls = [135 • -225 • ] and larger EUV fluxes. Despite this observed asymmetry, solaraberrated axisymmetric models may still provide a worthy first approximation of the shock's shape and position.
To investigate further the conditions of the shock's asymmetry throughout different solar cycles, solar drivers and internal drivers such as crustal magnetic fields, and isolate their respective contribution, a full exploitation of MAVEN's continuously growing datasets is warranted; likewise, a reanalysis of past encounters at Mars using 3D quadric fits would be a welcome addition. These are left for future studies. Applications of these methods, especially in 3D, to other bodies with large orbital eccentricities (such as Mercury) may also prove of interest.

Appendix A Characteristics of a quadric surface
The 3D planetary bow shock in this paper is approximated as a quadratic surface described by the Cartesian equation (9). Mathematically, 17 different quadrics can exist. However, here only 3 are physically acceptable for the approximation of a bow shock surface. These are the 'real' ellipsoid, the elliptic paraboloid, and the hyperboloid of two sheets. From coefficients A to I defining the quadric's surface equation, it is possible to extract some more 'physical' quantities of these surfaces such as the centre of the surface, the direction of the principal axes, the typical length, or the 'nose' of the surface. This requires the analysis of one particular matrix M given by: Determinant det M yields useful pieces of information on the considered surface. If det M < 0, the surface is an ellipsoid or an hyperboloid of two sheets. If det M = 0, it is an elliptic paraboloid.
The coordinates of the centre of the surface is given by: if M −1 exists. In the case of an elliptic paraboloid, there is an infinite number of centres placed along the intersections of the two planes of symmetry.
Finally, one may be interested in the position of the tip (tail-like direction) or of the nose (subsolar direction) of the surface. These extremum points are at a distance L 1 , L 2 , and L 3 from the centre in the direction ±V i . Therefore they are given by: and P ±,ext = P centre 1 1 1 where the columns of P ±,ext are the locations of the extrema. For a hyperboloid of two sheets, only the real solution associated with λ 3 should be considered: this gives the position of the noses or tips of both sheets. The sunward-most position of the ellipsoid (its nose) is referred to as P nose .
The tip or nose of the ellipsoid is in our context along the direction of the eigenvectors with the largest X (in absolute value) component.
The volume of the ellipsoid is: Table A1 presents the length of each of the three principal axes of the quadric L i , the ellipsoid's volume, eigenvectors V i and the coordinates of the centre and sunward nose of the surface for each ellipsoid in Table 4.

Appendix B Subsolar tip of the trace of an ellipsoid surface in Cartesian coordinates
The subsolar point of the projection of a 3D ellipsoid in 2D planes, as shown in Figure 9 (crosses), can be obtained by finding the roots of the corresponding 2D conic in the plane considered. For z = 0, Equation (9) becomes a second order equation: Ax 2 + By 2 + Dxy + Gx + Hy − 1 = 0.
Fixing variable y, the equation can be put in quadratic form with the following positive root: Finding the maximum of this function is equivalent to finding a y value that maximises this function. Posing ξ = Dy + G, its derivative has the form: Solving ∂x M ∂y = 0 for y and using that result in Equation (B2) makes it possible to calculate the final (x, y) coordinates of the projected ellipsoid's tip in the corresponding x − y plane.
An identical reasoning can be made for the x − z plane.
This tip in a plane is however not necessarily the farthest subsolar point of the ellipsoid's surface. Its position in 3D is by contrast given by Equation (A2) in Appendix A. Table A1.
Characteristics of the 3D Martian bow shock as derived from the MAVEN dataset, see Table 4 for the parameters of the 3D surfaces considered. All quadrics are ellipsoids. Mars' radius is R p = 3, 389.5 km.   where their corresponding yearly fits are plotted. Candidate detections points for each case are also drawn as filled circles of varying colours, with the opacity giving a measure of the density of points in that area, giving more or less weight to the fitting method.   Figure 10. Comparison of standoff distances, both at the subsolar point and at the terminator, calculated from the 2D and 3D fits, and for each case as in Tables 3 and 4