A 3D parametric Martian magnetic pileup boundary model with the effects of solar wind density, velocity, and IMF

Using global magnetohydrodynamic simulations, we construct a 3D parametric model of the Martian magnetic pileup boundary (MPB). This model employs a modiﬁed parabola function deﬁned by four parameters. The effects of the solar wind dynamic pressure, the solar wind densities and velocities, and the intensity and orientation of the interplanetary magnetic ﬁeld (IMF) are examined using 267 simulation cases. The results from our parametric model show that (1) the MPB moves closer to Mars when the upstream solar wind dynamic pressure ( P d ) increases, the subsolar standoff distance decreases and the ﬂaring degree of the Martian MPB increases with an increasing P d according to the power-law relations. For the same P d , a higher solar wind velocity (a lower density) leads to a farther location of the MPB from Mars, along with a larger ﬂaring degree, which is explained by the higher solar wind convection electric ﬁelds and a stronger magnetic pileup process under these conditions. (2) Larger Y or Z components of the IMF, B Y or B Z , result in a thicker pileup region and a farther MPB location from Mars, as well as a decrease in the ﬂaring degree. The radial IMF component, B X , has little effect on the geometry of the MPB. (3) In most of the simulations used to derive the current parametric model, the strongest Martian crustal magnetic ﬁeld is located on the dayside. However, for a larger value of the southward IMF, the Martian MPB is located farther away in the northern hemisphere instead of the southern hemisphere. The north-south asymmetry of the Martian MPB with the southern hemisphere being farther away is observed for other IMF directions. We suggest that the magnetic reconnection of the southward IMF with the crustal ﬁeld that occurs at middle latitudes of the southern hemisphere results in differ- ent magnetic ﬁeld topologies and the closer location of the MPB under these conditions. Our model results show a relatively good agreement with the previous empirical and theoretical models.


Introduction
The interaction of the solar wind with Mars, which is a weakly magnetized planet with an atmosphere, is complicated (e.g., Nagy et al. 2004). Generally, due to the mass loading and the magnetic pileup processes, the interplanetary magnetic field (IMF) can pile up outside the Martian ionosphere; the compressed magnetic field balances the solar wind dynamic pressure and results in the formation of the Martian induced magnetosphere (e.g., Holmberg et al. 2019). The Martian magnetosphere has an important boundary that is similar to the magnetopause in the terrestrial magnetosphere. This boundary is referred to either as the magnetic pileup boundary (MPB) when it is defined based on the strength and fluctuation of the magnetic field (along with the suprathermal electron fluxes) from the Mars Global Surveyor (MGS) mission (e.g., Vignes et al. 2000;Edberg et al. 2009), or as the induced magnetosphere boundary (IMB) based on the particle characteristics for the Mars Express (MEX) spacecraft (e.g., Lundin et al. 2004;Dubinin et al. 2006). Usually, the MPB and IMB are considered to be the same boundary (e.g., Matsunaga et al. 2017;Holmberg et al. 2019), and this assumption is used in the present work as well.
The geometry of the Martian MPB has been investigated on a number of different occasions (e.g., Vignes et al. 2000;Edberg et al. 2008;Němec et al. 2020). The solar wind dynamic pressure (P d ) is recognized as one of the main factors controlling its position (e.g., Crider et al. 2003;Dubinin et al. 2007). Generally, it is reported that the Martian MPB is located closer to Mars during higher P d intervals (e.g., Brain et al. 2005;Matsunaga et al. 2017). Specifically, using data from Phobos 2, Verigin et al. (1993) found a power-law dependence of the asymptotic magnetotail diameter of the IMB on P d . Similarly, Crider et al. (2003) presented a power-law correlation of the terminator distance of the MPB with P d based on MGS data. However, based on the observations of MGS and MEX, Edberg et al. (2009) reported an exponential relation between P d and the extrapolated terminator distance of the IMB. Even more recently, also using MEX data, Ramstad et al. (2017) argued that the subsolar standoff distance of the IMB decreases with P d according to a power-law function, instead of an exponential one. In addition, using the Mars Atmosphere and Volatile Evolution (MAVEN) spacecraft data, Lentz et al. (2021) further reported that, compared with the magnetic pressure of the magnetic pileup region and the thermal pressure of the ionosphere, the effect of P d on the standoff distances of the IMB was relatively weak during the interplanetary coronal mass ejection event of September 2017. The effects of P d corresponding to different solar wind densities and velocities have also been studied. For the Earth, using global magnetohydrodynamic A&A 664, A74 (2022) (MHD) simulations, Samsonov et al. (2020) reported that under the same P d , the change in the subsolar standoff distance of the magnetopause is greater when velocity increases compared to the change when density is enhanced, especially for a southward IMF. In contrast, for Mars, using MAVEN data, Ramstad et al. (2017) suggested that changing solar wind density or velocity individually for the same P d had little effect on the location of the Martian IMB. However, in their observational study, it was hard to remove the influence of other factors completely, such as the IMF and the crustal field. Recently, using a 3D multispecies MHD model, Wang et al. (2021) confirmed the power-law relation between the subsolar standoff distance of the Martian MPB and P d . They further indicated that for the same P d , a higher solar wind velocity (lower density) can result in a larger subsolar standoff distance, which is explained by the stronger magnetic pileup process. However, Wang et al. (2021) only focused on the subsolar point of the Martian MPB, and did not study the dependence of the whole geometry of the Martian MPB on the individual solar wind density and velocity.
Another factor that affects the location of the Martian MPB is the IMF. Although several publications have reported that the Martian MPB location is sensitive to the IMF direction (e.g., Brain et al. 2003Brain et al. , 2005Dubinin et al. 2006), the details of the influences of the IMF on the geometry of the MPB remain controversial. Furthermore, the magnetic field topology of the Martian induced magnetosphere is also dependent on the upstream IMF (e.g., DiBraccio et al. 2018;Weber et al. 2020;Xu et al. 2020). Specifically, using MGS observations, Brain et al. (2005) indicated that the MPB altitude in the northern hemisphere is partly controlled by the IMF draping direction. Later, Dubinin et al. (2006) showed a strong "dawn-dusk" asymmetry of the magnetospheric boundary in the IMF coordinates that is associated with the IMF draping asymmetry. Based on MEX and MGS data, Edberg et al. (2009) found that the IMF direction influence on the MPB is weaker than P d , but still significant. That is, in the hemisphere of locally upward solar wind electric field (E SW = −V SW × B IMF ), the MPB generally moves outward, while it is the opposite in the southern hemisphere where the MPB actually moves inward. However, using MAVEN observations, Matsunaga et al. (2017) argued that the IMB is always pushed to high altitudes in the downward E SW hemisphere. Moreover, Ramstad et al. (2017) further indicated that by contributing additional magnetic pressure and inducing stronger ionospheric currents, an enhanced IMF can affect the IMB location. However, the details of this effect were not discussed in their work.
The influence of the Martian crustal magnetic field on the MPB has also been studied. The intrinsic crustal magnetic field of Mars is complex and the most intense crustal sources are mainly concentrated in the area of the heavily cratered southern highlands in the southern hemisphere (e.g., Acuna et al. 1999;Connerney et al. 1999). As a result, it was reported that, on average, the MPB is located farther away in the southern hemisphere than that in the northern hemisphere (e.g., Crider et al. 2002;Brain et al. 2003;Edberg et al. 2008). The average difference of the MPB location between the two hemispheres is estimated to be ∼10% (Edberg et al. 2009). Moreover, Matsunaga et al. (2017) suggested that this north-south asymmetry is permanent, regardless of the intense crustal fields locations, while the boundary locations become higher when the intense crustal fields are located on the dayside. Using MHD simulations, Fang et al. (2017) confirmed the crustal field distribution effect on the IMB location and further indicated that the IMB locations cannot be simply parameterized only by the subsolar longitude; the subsolar latitude information is equally important.
The solar extreme ultraviolet (EUV) radiation also affects the background Martian atmosphere and ionosphere conditions (e.g., Kim et al. 1998;Bougher et al. 2000;Lundin et al. 2013) and, therefore, is expected to influence the MPB location. However, some early studies did not find obvious relations between them (Vignes et al. 2000;Modolo et al. 2006). More recently, based on the MEX and MGS data, Edberg et al. (2009) showed that the MPB radius decreases linearly with increasing EUV flux. However, Edberg et al. (2009) did not simultaneously constrain solar wind parameters (Ramstad et al. 2017), which might have affected their results. Ramstad et al. (2017), also using MEX data, indicated that the IMB expands significantly with increasing EUV radiation under the constrained nominal solar wind conditions. In addition, Lentz et al. (2021) further confirmed that, during the interplanetary coronal mass ejection event of September 2017, the subsolar standoff distance of the IMB increased with the EUV flux.
As various factors influencing the Martian MPB location have been investigated, efforts have also been made to model its size and shape. Initially, static models for the MPB location were developed to represent its position under typical conditions (e.g., Vignes et al. 2000;Dubinin et al. 2006;Trotignon et al. 2006;Edberg et al. 2008). This was largely due to the limited number of MPB crossings at Mars. Usually, the conic section function, r = L/(1 + ε cos θ), was employed. Here, r and θ are the polar coordinates with origin at X 0 , which is displaced from the Mars center. L and ε are the semilatus rectum and the eccentricity, respectively. Generally, the values of X 0 were in the range of 0.6−0.9 R M , and L was 0.9−1.1 R M , giving the subsolar standoff distance r 0 of 1.25−1.33 R M . The eccentricity, ε, was typically slightly less than 1.0, meaning that the MPB shape was closed on the nightside. However, the locations of the Martian MPB are dynamically changing in response to the upstream and downstream conditions, which could lead to a remarkable variability. Therefore, a dynamic model with the corresponding control factors is needed. Using MEX data, Ramstad et al. (2017) constructed the first dynamic Martian IMB model depending on the upstream solar wind density and velocity, which showed that the standoff distance and the eccentricity (flaring degree) of the Martian IMB decrease with P d . However, due to the limited number of available crossings, their model only employed the solar wind density and velocity as controlling parameters, while other factors such as the intensity and orientation of the IMF, and the EUV flux were not considered. Recently, based on an automated region identification approach and the MAVEN data, Němec et al. (2020) derived another dynamic empirical model of the Martian MPB, which includes the dependence on the solar wind dynamic pressure, the solar ionizing flux, and the crustal magnetic fields. However, due to the sparse data coverage, this MPB model was deemed unreliable beyond the terminator (X ≤ −0.5 R M ). In addition, the models of Ramstad et al. (2017) and Němec et al. (2020), are both axisymmetric, which means they do not describe the north-south and dawn-dusk asymmetries of the Martian MPB.
Using MHD simulations, we constructed a 3D parametric Martian MPB model, which uses a modified parabola function depending on six parameters. In this paper, the effects of different solar wind densities and velocities, as well as the intensity and orientation of the IMF, are examined.

Simulation model
In this work, we employed the 3D multispecies MHD model originally developed by Ma et al. (2004) to study the interaction between the solar wind and the Martian magnetosphere (Fang et al. 2010(Fang et al. , 2015(Fang et al. , 2017Ma et al. 2014a,b). This model is implemented within the Block Adaptive-Tree Solar wind Roe-type Upwind Scheme (BATSRUS) in the Space Weather Modeling Framework (SWMF; Tóth et al. 2012). It is worth mentioning that, besides the single fluid multispecies MHD model developed by Ma et al. (2004), BATSRUS also has other different versions of the Martian plasma models, such as the multifluid MHD model (e.g., Najib et al. 2011;Dong et al. 2015) and the multifluid electron pressure MHD model (e.g., Ma et al. 2013). Also, hybrid simulations have been employed to study the interaction of the solar wind with the Martian upper atmosphere (e.g., Modolo et al. 2016;Jarvinen et al. 2018). On the one hand, the model comparison by Egan et al. (2018) shows that, in general, the plasma boundaries in all these models vary only slightly in extent and appear symmetric overall. On the other hand, the single fluid multispecies MHD model of Ma et al. (2004) obtains the more accurate plasma boundaries, while the multifluid MHD model and the multifluid electron pressure MHD model get a more extended shock region, and the hybrid simulations show a more compressed shock region (Egan et al. 2018). In addition, we used the model of Ma et al. (2004) to construct a 3D parametric Martian bow shock model (Wang et al. 2020a) and study the effect of solar wind density and velocity on the subsolar standoff distance of the Martian MPB (Wang et al. 2021). Our parametric bow shock model also shows the north-south and dawn-dusk asymmetries, which are related to the IMF condition and the location of the intense crustal field, and it is in good agreement with the previous empirical and theoretical models. Hence, the single fluid multispecies MHD model developed by Ma et al. (2004) is consistently used in this work.
In this model, the input solar wind parameters include the components of the solar wind velocity (V X , V Y , V Z ), the solar wind number density (n), the solar wind temperature (T ), and the components of the IMF (B X , B Y , B Z ) in the Mars-centered solar orbital (MSO) coordinate system. Our model uses the MSO coordinate system in which the origin is located at the Mars center, the X-axis points to the Sun, the Z-axis is normal to the Martian orbital plane, and the Y-axis completes the right-handed system.
The computational domain of a typical simulation is defined km is the radius of Mars, and the inner boundary is taken to be 100 km above the Martian surface. This model applies a spherical grid near Mars in which the resolution in longitude and latitude is uniform at 3 • , while the radial resolution varies from 10 km at the lower boundary to 630 km at the outer boundary. In the radial direction, the grid size is <0.01 R M for 1.20 R M ≤ r ≤ 1.40 R M . This is the region where the subsolar part of the Martian MPB is generally located. The other details of the model are given in Ma et al. (2004).
The Martian crustal field is described using the 60-order spherical harmonic model developed by Arkani-Hamed (2001). For most of the simulations analyzed in this work, following Wang et al. (2020a) and Wang et al. (2021), the subsolar location was fixed to 180 • west longitude and 0 • north latitude, which meant that the strongest Martian crustal magnetic field was located in the dayside region (e.g., Crider et al. 2002;Fang et al. 2017). As a result, the Martian season in the simulations was assumed to be the spring or autumn equinoxes. Preliminary results for the strongest crustal fields located on the nightside are discussed in Sect. 4.3.
Solar extreme ultraviolet (EUV) radiation, which varies with solar cycle, affects the background Martian atmosphere and ionosphere conditions (e.g., Kim et al. 1998;Bougher et al. 2000;Lundin et al. 2013). In this work, most of the simulations correspond to typical solar maximum conditions that are identical to Case 1 of Ma et al. (2004). Two simulations based on solar minimum conditions are briefly discussed in Sect. 4.3. In all cases, the density profiles of the neutral species in the Martian atmosphere are assumed to be spherically symmetric.
The other solar wind parameters were set as follows: the solar wind velocities V Y and V Z were chosen to be 0, the upstream solar wind ion temperature T i = 5 × 10 4 K, and the electron temperature T e = 3 × 10 5 K, representing typical solar wind conditions upstream of Mars (e.g., Halekas et al. 2017;Liu et al. 2021). We simulated a total of 267 cases with the different solar wind dynamic pressure P d corresponding to different number density n and velocity V, along with the different magnitudes and orientations of the IMF. Specifically, in this work P d varied from 0.15 to 9.0 nPa (n ∈ [1, 16] cm −3 , V X ∈ [−250, −1334] km s −1 ), and the IMF strength range was from 0 to 12 nT. The detailed information for all the simulations is listed in the supplementary document, which is available online 1 .
In this work, the IMF cone angle, µ, is defined as the angle between the IMF direction and the X-axis in the MSO coordinates, For the B X = B Y = B Z = 0 case, µ was set to 0. The IMF clock angle λ is defined as the angle between the Z-axis and the projection of the IMF on the Y−Z plane in the MSO coordinates, as follows: For B Y = B Z = 0 cases, λ was set to π/2.

3D parametric model function of the magnetic pileup boundary
We employ a modified parabola to describe the 3D shape of the Martian MPB, as follows: where r = √ X 2 + Y 2 + Z 2 is the radial distance from the Mars center, θ = arccos(X/r) is the angle between the Mars-Sun line (the X-axis in the MSO coordinate system) and the radial direc- is the azimuthal angle measured from the positive Z-axis of the MSO coordinate system in polar coordinates. This function is based on the shape of the terrestrial magnetopause proposed by Shue et al. (1997), and has been modified and extensively used in numerous other models (e.g., Lin et al. 2010;Liu et al. 2015). The parameters in the function can directly represent the shape and size of the Martian MPB, such as the subsolar standoff distance and the flaring degree, which we will discuss in detail next. The azimuth angle φ is introduced to describe the north-south and dawn-dusk asymmetries, because of its different trigonometric values on the Y−Z plane.
The shape of the MPB given by Eq. (3) is open on the nightside as long as the value of the exponent is positive. However, for the values of the exponent less than 0.5, the MPB asymptotically approaches the negative X axis. Thus when considered only for relatively small values of X on the nightside, the shape of MPB may appear to be closed (see Fig. 1 of Shue et al. 1997). There are currently not enough observations to constrain the MPB far in the tail, and therefore we do not consider our model to be applicable to X < −4 R M .
In Eq.
(3), four parameters, r 0 , α 0 , α 1 , and α 2 are used to describe the size and shape of the MPB. The parameter r 0 represents the MPB subsolar standoff distance, which is the distance from the Mars center to the intersection of the MPB and the Mars-Sun line. Here, α 0 describes the flaring degree of the MPB tail, which plays the same role as in other magnetopause models (Shue et al. 1997;Lin et al. 2010;Liu et al. 2015). The parameter α 1 introduces the north-south asymmetry of the MPB; for α 1 < 0 the flaring degree in the southern hemisphere is larger than that in the northern hemisphere, and for α 1 > 0 the reverse is true. Similarly, α 2 denotes the dawn-dusk asymmetry; for α 2 < 0 the flaring angle in the eastern hemisphere is smaller than that in the western hemisphere. Thus, we use a set of the parameters of r 0 , α 0 , α 1 , and α 2 to describe the shape and size of the Martian MPB with the north-south and dawn-dusk asymmetries.

Identification of the magnetic pileup boundary
In this section, we discuss the identification of the location of the MPB in the simulation data. The Martian MPB, the interface between the solar wind and the Mars-induced magnetosphere, is created by the mass loading and magnetic pileup processes. As a result, some of the physical properties change dramatically across the MPB. For instance, inside the MPB, the magnitude of the magnetic field and the density of the predominant planetary heavy ions (O + and O + 2 ) are expected to increase suddenly while the solar wind protons (H + ) density decreases rapidly. At the same time, the fluctuations of the magnetic fields and the electron fluxes are expected to decrease as well (e.g., Trotignon et al. 2006;Dubinin et al. 2007;Matsunaga et al. 2017). Therefore, several methods were used to define the location of the MPB in previous studies. The observational work employed methods based on increase in the magnetic field strength and the absence of strong magnetic fluctuations (Matsunaga et al. 2015), the ratio of (n O + + n O + 2 )/n H + (Matsunaga et al. 2017), and the ratio of (P th + P d )/P B (Xu et al. 2016, here, P th , P d , and P B represent the plasma thermal pressure, the dynamic pressure, and the magnetic pressure, respectively). In numerical studies, the MPB location was calculated based on the balance between the magnetic pressure and the solar wind thermal pressure (Ma et al. 2014a), the gradient of the magnetic field magnitude (Fang et al. 2015;Wang et al. 2021), and the flow speed gradient (Fang et al. 2015). In general, the methods based on the pressure balance and the magnetic field magnitude gradient are more accurate in the dayside region, while the methods using the flow speed gradient work better on the nightside (Fang et al. 2015).
In this work, we combined the methods of the pressure balance and the flow speed gradient to determine the location of the MPB. First of all, we divided the azimuth angle φ into 32 sections and found the MPB locations along the radial direction for each θ resolution in each φ section, which was expected to be asymmetric in 3D. In each φ section, we determined the MPB locations using two different methods, depending on the solar zenith angle (SZA). (1) In the region of SZA ≤ 45 • , the MPB was identified according to the pressure balance ratio of (P th _sw + P d _sw)/(P B + P th _ion) ≤ 1.0. Here, P th _sw, P d _sw, P B , and P th _ion are the thermal pressure and dynamic pressure of solar wind, the magnetic pressure, and the thermal pressure of the predominant planetary heavy ions (O + , O + 2 , and CO + 2 ), respectively. Although P d _sw and P th _ion are almost negligible near the MPB (especially at the subsolar region), taking them into account allowed us to determine the location of the MPB more accurately. (2) In the region of SZA ≥ 90 • , the MPB locations are defined by the extreme value of the flow speed gradient when r < (r BS _model − 0.1). Here, r BS _model is the corresponding radial distance of the bow shock from the Mars center, calculated by the Martian bow shock model of Wang et al. (2020a). Therefore, the extreme value of the flow speed gradient corresponding to the bow shock location is excluded.
(3) For SZA ∈ (45 • , 90 • ), we employed these two methods to find the two locations of the MPB, which refer to P ratio and dUt, and the MPB locations (X , Y , Z ) were calculated according to the relations as follows: ]. This technique makes the method of determining MPB gradually change from the pressure balance ratio of SZA ≤ 45 • to the flow speed gradient of SZA ≥ 90 • . The final method combines two definitions of the MPB in a smooth manner. Figure 1 shows an example of the identified locations of the Martian MPB from simulation data. From the figure, we can see that the determined MPB locations agree well with the place where the solar wind total pressure and the vectors of the solar wind velocity change sharply. Moreover, Fig. 1 also displays some features of the Martian MPB (such as the north-south and dawn-dusk asymmetries), which we will discuss in detail in the next section.
After determining the 3D MPB locations, we fit them with the model Eq.
(3) for all 267 simulation cases. The method of nonlinear least squares, which employs the algorithm of Trustregion with 95% confidence bounds, was used to do the fitting. For the case in Fig. 1, the fitting results are r 0 = 1.270 R M , α 0 = 0.3874, α 1 = −0.03509, and α 2 = 0.01443. The determination coefficient of the fitting R 2 = 0.9437 (R 2 = 1 corresponds to a perfect fitting), which denotes a relatively high goodness of fit and demonstrates the good performance of our MPB shape function. For the entire set of the 267 simulations used in this study, the average value of R 2 was 0.9268. The detailed information on all the fitting results is listed in the supplementary document.

Model results
In this section, we discuss the specific effects of the solar wind parameters on the size and shape of the Martian MPB. First, we compare the results of fitting the dependence of r 0 and α 0 on P d with the exponential and the power-law functions. The powerlaw function gives a better fit than the exponential one, which is consistent with Ramstad et al. (2017) and Wang et al. (2021).  solar wind dynamic pressure P d calculated with different solar wind number densities n and velocities V. Panel a of Fig. 2 shows that the dependence of r 0 on P d can be presented by the power-law relations; r 0 decreases with increasing P d . For the same P d , a higher solar wind density results in a smaller r 0 (the differences can reach ∼0.05 R M ), which is in agreement with the study by Wang et al. (2021). Panel b shows that α 0 is enhanced with increasing P d as a power-law; a higher solar wind density results in a smaller α 0 under the same P d , and the differences can reach ∼0.1. These differences indicate that r 0 and α 0 can be better described by fitting n and V individually than by just fitting P d . This conclusion is in agreement with Wang et al. (2021). Panels c and d indicate that for increasing P d , α 1 and α 2 increase. Their relations on P d can be approximated well by the hyperbolic tangent functions. These two parameters are mainly P d -dependent; changing n and V individually in such a way that the dynamic pressure remains the same and has little effect. In addition to analyzing the effects of n and V on the MPB shape for these IMF conditions (B X = −1.6776 nT, B Y = 2.4871 nT, and B Z = 0 nT), we also studied them for other IMF conditions and the results are qualitatively similar (see the supplementary document). Figure 3 shows the dependence of the MPB parameters on the IMF components. It is worth mentioning that the IMF conditions of Fig. 3 were set with one nonzero component while the other two components were null. From Fig. 3 we can see that  B X barely affects all the MPB parameters, except that α 1 slightly decreases with increasing B X . In contrast, B Y and B Z have a significant influence on the MPB parameters. Panels a and b show that when the absolute values of B Y and B Z increase, r 0 increases while α 0 decreases. The dependence of α 1 on the IMF is a bit complex. Specifically, for B Y and a positive B Z condition, α 1 is negative (which means that the MPB in the southern hemisphere is located farther than in the northern one). Also, the value of α 1 increases with the increasing absolute value of B Y and decreases with positive B Z . For negative B Z , α 1 is positive and its value increases with the increasing absolute value of B Z , which implies a more distant northern hemisphere and increasing north-south asymmetry. We suggest that such different effects of positive and negative B Z are caused by different magnetic field topologies on the dayside, which we will discuss in detail next. Panel d demonstrates that the parameter α 2 is sensitive to the IMF direction. That is to say, for negative B Y and B Z , α 2 is a bit less than 0, while α 2 ≈ 0.02 for positive B Y and it is close to 0.01 for positive B Z . The effects of the IMF components on the MPB parameters for other values of P d are also included in the supplementary document.
Next, we examine the effects of the IMF clock angle λ on the MPB parameters when B X = 0 nT, which are shown in Fig. 4. The parameter r 0 follows a trigonometric relation (mostly the cosine function) with λ. A similar trigonometric relation seems to exist between λ and α 0 . The parameter α 1 generally increases with an increasing absolute value of λ, while α 2 displays no distinct dependence on λ. Figure 5 shows the effects of the IMF cone angle µ on the MPB parameters. The parameter r 0 can be approximated by quadratic functions of µ. When µ increases from 0 to π/2, r 0 increases; when µ ∈ [π/2, π], r 0 decreases with µ. The parameter α 0 also appears to follow a quadratic relation with µ, while its change is opposite to that of r 0 . The parameter α 1 is generally positive for some negative B Z cases when µ > π/2, and α 1 < 0 for the other IMF conditions. This also agrees with the result from panel c of Fig. 3 and will be discussed in detail next. The dependency of α 2 on µ is a bit complex. For µ that is changed by B X and B Z , its dependence on α 2 is weak. For µ varying with B X and B Y , α 2 is generally positive for the positive B Y cases, and α 2 appears to have a quadratic relation with µ varying with negative B Y . In general, the effect of B Y on α 2 is larger than that of B Z . Further results concerning the µ effects for other values of P d and the IMF strengths are in the supplementary document.
In addition, according to Figs. 3 and 5, we suggest that the radial IMF B X effect on the geometry of the Martian MPB is relatively small and that the IMF cone angle influence is also mainly due to the B Y and B Z components. As a result, we employed the B t = B 2 Y + B 2 Z and IMF clock angle λ as the IMF factors in the MPB model. Also, to minimize the impact of the pure B X cases in this work, we excluded the cases of B t = 0 (23 in 267) to fit the MPB model. Using the results of the selected set of 244 simulations, we obtained the following best fit for the MPB shape parameters: r 0 = p 1 n p 2 |V| p 3 (p 4 B t + 1.0)[p 5 cos(p 6 λ) + 1.0], α 0 = p 7 n p 8 |V| p 9 (p 10 B t +1.0)[p 11 cos(p 12 λ)+p 13 sin(p 12 λ)+1.0], The numerical values of the 32 coefficients in these equations are given in Table 1. Here, it is assumed that n is measured in cm −3 , V in km s −1 , P d in nPa, B t in nT, λ in radians, and r 0 in R M . The fitting determination coefficients R 2 for Eqs. (4)-(7), are 0.9212, 0.8211, 0.6427, and 0.5589, respectively.
Next, for all 267 cases, we further examined the R 2 model and the root mean squared error (RMSE) as shown in Fig. 6. These two coefficients were used to measure the errors between the identified MPB locations and the final model functions (Eqs. (4)-(7) and   a larger R 2 and a better fitting. Our results show that, although R 2 model is generally smaller than R 2 (which is directly fitted from identified simulation data), the mean value of R 2 model is 0.9079 (the black dashed line in panel a), which indicates a good fitting and a relatively high accuracy of the model result. Moreover, the mean value of RMSE is 0.1572 (the black dashed line in panel b), which indicates that the average quadratic mean value of the differences between the model results and the identified MPB locations is 0.1572 R M . However, for a very few cases, R 2 model and RMSE are far from the mean value, which may be due to an error in the simulation or the fitting. In general, both R 2 model and RMSE demonstrate the high goodness of fit of the model result. Figure 7 illustrates the effects of the solar wind number density n and velocity V on the shape and size of the Martian MPB. Panels a and b are the influences of n and V; the left, middle, and right panels show the modeled MPB in the X−Y, X−Z, and Y−Z planes in the MSO coordinates, respectively. For panel a, the solar wind dynamic pressures are 1 nPa (the number density n = 2.4 cm −3 , the red lines), 2 nPa (n = 5.0 cm −3 , the green lines), and 3 nPa (n = 7.2 cm −3 , the blue lines). The solar wind , π/6 (the red lines), π/3 (the blue lines), and π/2 (the green lines) with B X = 0 nT and B t = 5 nT. The solar wind dynamic pressure is P d = 1 nPa (n = 4.0 cm −3 , V X = −400.0 km s −1 ) for all these cases.
From Fig. 8 and Eqs. (4)- (7), we can see that as the absolute value of B Y or B Z increases, the subsolar standoff distance increases and the flaring degree of the Martian MPB decreases. The degree of the north-south asymmetry decreases with increasing |B Y | and increases with |B Z |. Moreover, for southward IMF, the southern hemisphere is generally located closer to Mars than the northern hemisphere (α 1 > 0) (for northward B Z , α 1 < 0). In other words, the well-known north-south asymmetry is reduced by the southward B Z when the strongest crustal field is on the dayside. Finally, the degree of the dawn-dusk asymmetry also slightly changes with B Y . Panel c and Eqs.  indicate that the subsolar standoff distance is proportional to a cosine function of λ, and the flaring degree, the degrees of the north-south asymmetry, and the dawn-dusk asymmetry are approximated by the Fourier functions of λ. For the range of the control parameters, P d (n, V), B t , and λ, considered, the range of the Martian MPB model parameters, r 0 , α 0 , α 1 , and α 2 are 0.38 R M (from ∼1.10 to ∼1.48 R M ), 0.16 (from ∼0.27 to ∼0.43), 0.08 (from ∼−0.058 to ∼0.024), and 0.06 (from ∼−0.020 to ∼0.039), respectively.

Mechanism of solar wind effects on the Martian MPB
In general, a planetary magnetopause forms in the area where inside pressure approximately equals the outside pressure (e.g., Chapman & Ferraro 1931;Shue & Chao 2013;Lu et al. 2015). Although Mars is an unmagnetized planet, the inside pressure is mainly the magnetic pressure provided by the crustal fields (e.g., Holmberg et al. 2019). The outside pressure consists of the dynamic, thermal, and magnetic pressures (e.g., Brain et al. 2010;Wei et al. 2012;Holmberg et al. 2019;Zhang et al. 2019). As a result, when the upstream P d increases, it pushes the MPB closer to the surface of Mars. In addition, the dayside region usually moves deeper than the flank region, and hence the flaring degree is enhanced with an increasing P d (e.g., Shue et al. 1997). For a higher P d when the MPB moves to a lower altitude and the asymmetry caused by the magnetic pileup process is expected to be more pronounced. Therefore, the dawn-dusk asymmetry increases with P d . However, the north-south asymmetry caused by the intrinsic magnetic field is compensated by the stronger magnetic pileup process, and so it decreases with P d . The differences between the effects of the solar wind density and the velocity on the MPB can be explained by the magnetic pileup process. As the planetary ion pickup and solar wind mass loading processes are driven by the solar wind convection electric field (E SW = −V SW × B IMF ), the IMF can pile up outside the ionosphere, in which the unmagnetized planet can prevent the penetration of the solar wind magnetic field by a diamagnetic current at the ionopause (e.g., Nagy et al. 2004;Fang et al. 2018). As a result, for a higher solar wind velocity under the same P d , the solar wind convection electric field is higher for the same IMF conditions. The mass loading effect due to the pickup process also becomes stronger, which would cause a stronger magnetic pileup region (e.g., Modolo et al. 2006;Vaisberg et al. 2017;Chang et al. 2020). In this vein, for a higher solar wind velocity, the magnetic pressure at the Martian MPB should be higher (Wang et al. 2021). Moreover, at a higher magnetic field strength, the plasmas near the plasma depletion region (the inner part of the Martian magnetosheath where the magnetosheath magnetic fields are enhanced and the magnetosheath plasmas are partially depleted) are more easily squeezed out along the pileup magnetic field, which leads to a decrease in the plasma density, plasma β, and the ion temperature, especially for the hotter plasmas (e.g., Lee & Lee 2020;Wang et al. 2020b). Thus, the stronger magnetic pileup process for the higher solar wind velocity causes the magnetic pressure to increase and the thermal pressure to decrease at the Martian MPB. In addition, in the magnetosheath the thermal pressure dominates, while in the magnetosphere it is the magnetic pressure (Brain et al. 2010;Cui et al. 2018;Sánchez-Cano et al. 2020). As a result, for the same P d , a higher solar wind velocity (a lower density) leads the MPB located farther away from Mars and having a higher flaring degree.
The mechanism by which the IMF affects the geometry of the Martian MPB is also due to the magnetic pileup process. As we mentioned above, the magnetic pileup process is caused by the interaction of the Martian ionosphere and the solar wind, in which the unmagnetized planet prevents the penetration of the solar wind magnetic field by the diamagnetic current at the ionopause, which results in the pileup of the IMF outside the ionosphere (e.g., Nagy et al. 2004;Fang et al. 2018). However, the IMF can only be compressed in the direction perpendicular to the solar wind speed, and hence only the IMF components perpendicular to the solar wind velocity, B Y and B Z , participate in the process and determine the size of the pileup region (e.g., Wang et al. 2020a). The pileup process is stronger at a larger IMF intensity. As a result, larger B Y and B Z result in a thicker pileup region and a farther MPB from Mars, which in turn increase the subsolar standoff distance and flaring degree of the MPB. For the IMF B X component, the IMF pileup is weak, so there is little effect of B X on the MPB. The IMF cone angle influence is also mainly due to the changing of the B Y and B Z components. In addition, when the MPB moves away from Mars with increasing intensities of B Y and B Z , the asymmetries of the MPB become smaller.
According to the classic magnetic field draping theory, at Mars the IMF is mainly piled up along the upstream IMF orientation in the plane perpendicular to the solar wind flow direction (e.g., Crider et al. 2004;Brain et al. 2005). Although these pileup and draping processes vary with the altitude and SZA (e.g., Fang et al. 2018), the IMF pileup is the strongest in the IMF clock angle direction (e.g., Luhmann et al. 2004). In this vein, however, the northward and southward B Z should have the same effect on the MPB, which disagrees with our model calculations. This disagreement, however, is explained by the magnetic field topology shown in Fig. 9. Figure 9 shows the magnetic field topology in the X−Z plane for the northward and southward B Z . The intrinsic crustal magnetic field of Mars is complex and the intense crustal sources are mainly concentrated in the area of the southern hemisphere in the Terra Sirenum region, which is located in the longitude range 150 • E to 240 • E and the latitude range 30 • S to 85 • S (e.g., Acuna et al. 1999;Connerney et al. 1999). In this work, for nearly all of the simulations, the subsolar location was fixed to 180 • W and 0 • N, which means that the strongest crustal magnetic field was located on the dayside (Arkani-Hamed 2001; Ma et al. 2004). Panel a of Fig. 9 shows that for the northward B Z , the crustal field and the IMF are basically parallel at SZA ≈ 45 • , which results in a thick pileup region. In contrast, for the southward B Z , shown in panel b, at SZA ≈ 45 • the crustal fields and the IMF are antiparallel, causing magnetic reconnection and leading to the erosion of the crustal field (e.g., Brain et al. 2003). As a result, the MPB for southward IMF is located closer to the planet than for northward IMF. Meanwhile, the northern hemisphere of the MPB is not affected by this as the crustal field is relatively weak there, and its average location is in between the two conditions mentioned above. As a result, for the northward B Z , the MPB in the southern hemisphere is located farther away from the planet than the northern hemisphere and has a larger flaring degree. This also partly agrees with Edberg et al. (2009), who suggested that the southern hemisphere of the Martian MPB reacts differently to the IMF orientation than the northern hemisphere. As the absolute value of B Z increases, this north-south asymmetry decreases as the whole MPB moves away from Mars. Moreover, the effects of the IMF clock angle on the subsolar standoff distance and the flaring degree might also be associated with local magnetic reconnection, which needs further investigation in the future.

Model comparison
In Fig. 10, we compare our parametric model of the Martian MPB to several existing empirical models for different solar wind conditions. It is worth mentioning that Vignes et al. (2000) and Edberg et al. (2008) models do not include any control parameters and only represent the MPB locations under the typical solar wind conditions in panel a. Figure 10 shows that for P d = 1 nPa (panel a), our parametric model is in agreement with the average location of the empirical models in the subsolar region, and is especially close to the Vignes et al. (2000) model. In the mid-latitude region, our model predicts the MPB being farther away from Mars than the other models, and in the tail region, our model also agrees with the empirical models well. For P d = 3 nPa, the results are similar and the two models agree very well on the dayside. Our model is located a bit farther away from Mars in the mid-latitude region when compared to the model of Ramstad et al. (2017). The asymmetries of our model are summarized in Table 2.
The detailed comparison between the subsolar standoff distance r 0 (R M ) and the terminator distance L X=0 (R M ) (measured from the Mars center) for different MPB models is presented in Table 2. For our parametric model, L X=0 values are displayed in the northern (φ = 0 • ), eastern (90 • ), southern (180 • ), and western (270 • ) directions. The IMF conditions used to calculate these values were the same as those used for Fig. 10. For P d = 1 nPa, it appears that the Ramstad et al. (2017) model slightly underestimates the size of the Martian MPB (r 0 = 1.218 R M ), while our parametric model (r 0 = 1.270 R M ) is in closer agreement with the empirical models (1.29 and 1.33 R M ). For P d = 3 nPa, r 0 = 1.190 R M for our parametric model is very close to the 1.187 R M of the Ramstad et al. (2017) model. Our model is the only one so far to include the asymmetries; the obvious asymmetries it shows are the in the north-south (L φ=180 • > L φ=0 • ) and dawn-dusk (L φ=90 • > L φ=270 • ) directions. Asymmetries in the MPB location are expected and have been reported previously (e.g., Crider et al. 2002;Brain et al. 2003;Edberg et al. 2008Edberg et al. , 2009Fang et al. 2017;Matsunaga et al. 2017). In the future, the corresponding crossing data of the Martian MPB from multiple satellites, such as MGS, MEX, MAVEN, and Tianwen-1 can be used to validate and improve our model.

Effects of the crustal magnetic field and the solar cycle
In this section, we make a first attempt at illustrating the effects of the orientation of the crustal magnetic fields and solar activity. We performed four simulations in which the subsolar positions of the crustal magnetic field were set to 180 • W 0 • N and 0 • W A74, page 13 of 15 A&A 664, A74 (2022)  Ramstad et al. (2017) Edberg et al. (2008) 1.33 ± 0.15 1.456 Vignes et al. (2000) 1.29 ± 0.04 1.468  (Ma et al. 2004). Table 3 shows the values of the model parameters r 0 , α 0 , α 1 , and α 2 obtained directly from the simulation data, as well as the average terminator distance of the MPB from the Mars center at X = 0, L x=0 . The other solar wind conditions used for these four cases were P d = 1 nPa (n = 4.0 cm −3 , V = −400.0 km s −1 ), B X = −1.6776 nT, B Y = 2.4871 nT, and B Z = 0 nT. Our results show that the MPB is located farther away from Mars under the solar maximum condition than under the solar minimum conditions, which agrees with Ramstad et al. (2017) and Lentz et al. (2021). Specifically, in the subsolar strong crustal field cases (180 • W 0 • N), L x=0 and r 0 were 3.793 and 0.7087% smaller, respectively, for solar minimum conditions than for solar maximum conditions. In weak crustal field cases (0 • W 0 • N), the differences were 6.246 and 6.544%, respectively. In both the solar maximum and minimum conditions, when the strong crustal field is located in the dayside region, the MPB moves away from Mars. Also, the flaring angle α 0 is enhanced under the solar maximum condition, while it does not change too much under the solar minimum condition. Moreover, in solar maximum conditions, the southern hemisphere asymmetry is more prominent when the strong crustal field is located on the dayside, which agrees well with the results of Fang et al. (2017) and Matsunaga et al. (2017). In solar minimum conditions, the northsouth asymmetry in the simulations is not obvious. In addition, our simulations also show that the effect of the crustal field on the MPB location is larger than that of the solar cycle. In the future, the effects of the crustal fields and the solar cycle need to be studied further and accounted for by Martian MPB models.

Summary and conclusions
Using a 3D multispecies MHD model, we constructed a 3D parametric Martian MPB model that includes the solar wind density and velocity, along with the intensity and orientation of the IMF, as controlling factors. We employed a modified parabola function to describe the shape and size of the Martian MPB. This shape function includes the north-south and dawn-dusk asymmetries, and is controlled by four model parameters. We determined the location of the MPB from the simulation data by combining the pressure balance and the flow speed gradient. A total of 267 individual simulations with different solar wind and IMF conditions were used to produce a database and calculate the four parameters of our MPB model. Conditions corresponding to the solar maximum and the strongest crustal field located on the dayside were used in most cases. Our main results are as follows: (1) the MPB moves close to Mars when the upstream P d increases; the subsolar standoff distance decreases and the flaring degree of the Martian MPB increases with increasing P d according to the power-law relations. Under the same P d , a higher solar wind velocity (and, therefore, a lower density) leads to a farther location of the MPB from Mars and a larger flaring degree. This is explained by the higher solar wind convection electric field and the stronger magnetic pileup process under these conditions. (2) A larger IMF Y or Z component, B Y or B Z , leads to a thicker pileup region and a farther MPB location from Mars, as well as a decrease in the flaring degree. There is little effect of the radial IMF component, B X , on the geometry of the MPB. (3) The conditions corresponding to the strongest crustal field located on the dayside were used in this work. Under these conditions, for the southward IMF, the Martian MPB is located farther away in the northern hemisphere than in the southern hemisphere, and this phenomenon is more prominent for a larger value of the southward IMF. While the well-known north-south asymmetry of the Martian MPB (the southern hemisphere of the MPB being farther away than the northern one) is only observed for other IMF directions. We suggest that the magnetic reconnection of the southward IMF and the crustal field that occurs at the middle latitude of the southern hemisphere results in the different magnetic field topologies and affects the location of the MPB.
The comparison shows a good agreement between our parametric model and the previous empirical models under typical solar wind conditions. There are two aspects of the model that A74, page 14 of 15 should be improved in the future. First, several other factors, such as the location of the crustal magnetic fields (or the seasonal dependency) and the solar cycle (or solar EUV radiation) should be included. Second, the Martian MPB crossing from multiple satellites, such as MGS, MEX, MAVEN, and Tianwen-1 should be used to validate the model.