A 3D Parametric Venusian Bow Shock Model with the Effects of Mach Number and Interplanetary Magnetic Field

Using global magnetohydrodynamics simulations, we have developed a three-dimensional parametric model for the Venusian bow shock based on a generalized conic section function defined by six parameters, with the effects of the solar wind magnetosonic Mach number (M MS) and the interplanetary magnetic field (IMF) involved. The parametric model’s results reveal the following findings: (1) The size of the Venusian bow shock is primarily determined by M MS. An increase in M MS results in the bow shock moving closer to Venus and a reduction in its flaring angle. (2) Both the subsolar standoff distance and the bow shock’s flaring angle increase with the strength of the IMF components that are perpendicular to the solar wind flow direction (B Y and B Z in the Venus-centered solar orbital coordinate system), whereas the parallel IMF component (B X ) has a limited impact on the subsolar standoff distance but affects the flaring angle. (3) The cross section of the bow shock is elongated in the direction perpendicular to the IMF on the Y–Z plane, and the elongation degree is enhanced with increasing intensities of B Y and B Z . (4) The quasi-parallel bow shock locates closer to the planet as compared to the quasi-perpendicular bow shock. These findings are in alignment with prior empirical and theoretical models. The influences of M MS and IMF on the bow shock’s position and geometry are attributed to the propagation of fast magnetosonic waves, showing the nature of the formation of a collisionless bow shock under the interaction of magnetized flow with an atmospheric object.


Introduction
Unlike our Earth, Venus is an atmospheric body without magnetic field.The solar wind can interact with the Venusian ionosphere directly and produce a bow shock (e.g., Russell et al. 1979;Bertucci et al. 2011;Futaana et al. 2017).The position of the Venusian bow shock is influenced both by upstream solar wind parameters, such as the interplanetary magnetic field (IMF) and Mach number (Zhang et al. 2004;Chai et al. 2014), and by the ionosphere which is contingent on factors like solar extreme ultraviolet (EUV) radiation and solar wind dynamic pressure (Knudsen 1988;Phillips & McComas 1991;Russell et al. 2006).The existence of the Venusian bow shock was initially confirmed by the Mariner 5 spacecraft (Bridge et al. 1967), and over the past five decades, extensive research has been conducted to characterize its size and shape (e.g., Spreiter et al. 1970;Zhang et al. 2007;Shan et al. 2015).
Observations have revealed the independence of Venusian bow shock position from the solar wind dynamic pressure (P d ; Russell et al. 1988;Zhang et al. 2004;Martinecz et al. 2008) because changes in the size of the Venusian ionosphere, resulting from different values of P d , are not substantial enough to influence the standoff location of the bow shock.However, the shock location is sensitive to variations in EUV flux, which controls both the size of ionosphere and the massloading effect (Slavin et al. 1979;Alexander & Russell 1985;Zhang et al. 1990).From the period of low solar activity to high solar activity, the average bow shock location at the terminator increased by approximately 0.35 times the radius of Venus (R V ), marking a significant ∼17% shift (Alexander & Russell 1985).
The primary factor influencing the location of the bow shock is the magnetosonic Mach number.A prior study by Russell et al. (1988) proposed the anticorrelation between the distance to the bow shock and the magnetosonic Mach number.Subsequently, the Venus Express (VEX) spacecraft verified the occurrence of a distant bow shock linked to an extraordinarily low magnetosonic Mach number, close to 1 (Zhang et al. 2008b;Xu et al. 2019).
The influence of the IMF on the Venusian bow shock has been a subject of study.For instance, Zhang et al. (1991) reported an asymmetry in the cross section of the bow shock based on Pioneer Venus Orbiter observations.Specifically, they observed that the bow shock on the quasi-parallel side was positioned closer to the planet than the portion on the quasiperpendicular side.However, the exact cause of this effect, whether it results from kinetic effects as suggested by the hybrid model of Thomas & Winske (1990), or magnetohydrodynamic effects, remains uncertain.Zhang et al. (1991) further quantified this asymmetry as 0.1 R V .In a separate study, Chai et al. (2014) analyzed 1680 bow shock crossings recorded by the VEX spacecraft to investigate the impact of the IMF on the Venusian bow shock.Their findings indicated that the position of the shock depended on the angle θ Bn , which is the angle between the IMF direction and the direction of the shock's normal.Additionally, in the terminator and tail regions, the bow shock was observed to be closer to the planet at the magnetic equator than observations at the magnetic poles (Russell et al. 1988;Chai et al. 2015).This pole-equator asymmetry is attributed to the faster propagation of plasma compressible waves at the magnetic poles where the magnetic field aligns perpendicularly with the wave group velocity (Alexander et al. 1986;Chai et al. 2015).
While investigating the determinants of the Venusian bow shock position, there have been endeavors to model its dimensions and configuration.Nevertheless, the dynamic variation of the bow shock position due to different upstream conditions has not been considered in the previous models.In this study, with the aim of assessing the variability of the Venusian bow shock location from the magnetohydrodynamic point of view, we conducted simulations across 110 cases with different magnetosonic Mach number and IMF.These simulations served as the basis for constructing a comprehensive dynamic Venusian bow shock model that incorporates relevant control parameters.

Simulation Model
The utilization of global magnetohydrodynamics (MHD) simulations enables the systematic examination of the configuration and dimensions of the Venusian bow shock.Here we employ a three-dimensional, multispecies MHD model that was originally developed by Ma et al. (2013) and has been consistently utilized in studies concerning the interaction between the solar wind and Venus, as demonstrated by Ma et al. (2020) and Xu et al. (2022).The input parameters include the following variables within the Venus-centered solar orbital (VSO) coordinate system: solar wind number density (n), the solar wind temperature (T), the components of solar wind velocity (U X , U Y , U Z ), and the components of IMF (B X , B Y , B Z ).In the VSO system, the coordinate origin is situated at the center of Venus, with the X-axis directed toward the Sun, the Zaxis orthogonal to the orbital plane of the planet, and the Y-axis completing the right-handed coordinate system.The numerical simulations are conducted using the Space Weather Modeling Framework (Tóth et al. 2012).The computational domain is defined as km is the radius of the Venus, and the inner boundary is taken to be 100 km above the surface of Venus.The model uses a nonuniform, spherical grid near Venus with the resolution of 2°.5 in longitude and latitude.The radial resolution is set to be 5 km at the lower boundary and increases gradually with altitude.
As the background atmosphere/ionosphere of Venus is associated with the solar EUV radiation or solar cycle (e.g., Fox & Sung 2001;Han et al. 2020), in this model the density profiles of the neutral species of the atmosphere for the two fixed conditions of solar cycle maximum and minimum (Ma et al. 2013) are adopted to study the dependence of bow shock on EUV, while the effects of solar wind parameters are examined for both solar minimum and maximum.
We establish the upstream solar wind temperature at 2.5 × 10 5 K, and assign values of 0 to the solar wind velocities U Y and U Z .Employing these solar wind parameters, we conduct simulations for a total of 110 cases, each varying in terms of solar wind magnetosonic Mach number (M MS ), and the magnitude and orientation of the IMF.These parameters are set to emulate the typical solar wind conditions observed upstream of Venus (Chang et al. 2018;Xu et al. 2021;Rojas Mata et al. 2022).More precisely, within the scope of this study, the solar wind exhibits variations in both density, ranging from 3.0 to 15 cm −3 , and velocity, spanning from 224 to 1000 km s −1 .Simultaneously, the strength of IMF ranges from 2.8 to 20 nT.The IMF orientation is represented by the IMF clock angle, denoted as λ, and the IMF cone angle μ.The IMF clock angle is defined as the angle between the Z-axis and the projection of IMF onto the Y-Z plane in the VSO coordinates,


For B Y = B Z = 0 cases, λ is set to π/2.μ is defined as the angle between the IMF direction and X-axis in the VSO coordinates.
In this work, we separate the IMF cone angle on two planes


Comprehensive information pertaining to all simulations can be found in the supplementary document of Table 1.

Identification of the Bow Shock
Based on the Rankine-Hugoniot relations and gradients of physical properties, we can pinpoint the locations of the bow shock based on the simulation results, a methodology akin to previous studies conducted for Earth and Mars (Wang et al. 2015(Wang et al. , 2016(Wang et al. , 2020b)).In our research, we implement a modified algorithm that relies on the gradient distribution of solar wind density.To be more specific, within each longitudinal and latitudinal resolution, the bow shock location is determined by identifying the point along the radial direction where the inward gradient of solar wind density is maximized.Figure 1 shows an example of the Venusian bow shock identification from the simulation data on 3D (panel (a)), Y-Z plane (panel (b)), X-Y plane (panel (c)), and X-Z plane (panel (d)) in the VSO coordinate system.Black dots are the identified Note.Column (1) is the number of these cases.Column (2)-( 10) are the input solar wind conditions for each case, including the solar wind number density, n (cm −3 ), the X component of solar wind velocity, U X (km s −1 ), the solar wind dynamic pressure, P d (nPa), the components of IMF, B X , B Y , B Z (nT), the magnitude of IMF, B t (nT), the IMF cone angle, μ (in radians), and the IMF clock angle, λ (in radians), the magnetosonic Mach number, M MS , and solar activity condition, respectively.The upstream solar wind temperature is fixed to 2.505 × 10 5 K.The solar wind velocities U Y and U Z are chosen to be 0. Column (13)-( 18) are the fitting results of model construction parameters, including the bow shock subsolar standoff distance, x 0 (R V ), the flaring degree, e 0 , the elongation degree in the twist direction, e 1 , the north-south asymmetry degree, e 2 , the dawn-dusk asymmetry degree, e 3 , the twist direction, ω (in radians), as well as the coefficient of determination for the case fitting, R 2 , in Column (19).In addition, the eccentricity of the conic section function, ε, is fixed to 1.05.Column (20) and 21 are the coefficient of determination, R 2 model , and the root mean squared error, associated with the identified locations of bow shock from simulations and the final model functions (Equations ( 3)-( 10) and Table 2) for each case.This table shows the information of 10 examples, all the cases are listed in the supplementary document.
(This table is available in its entirety in machine-readable form.)locations in each panel.The identified locations of the Venusian bow shock agree very well with the positions where the radially inward gradient of the solar wind density is maximized.The subsolar standoff distance of the bow shock identified in this instance is 1.348 R V , a value closely aligned with those from previous empirical models.Specifically, it is comparable to the 1.364 R V proposed by Shan et al. (2015), the 1.28 R V reported by Slavin et al. (1984), and the 1.32 R V outlined in Zhang et al. (2008a) at solar minimum.As this simulation corresponds approximately to the solar wind conditions near Venus at solar minimum (Xu et al. 2023), this alignment also serves as an affirmation of the MHD model's capability to faithfully replicate the observed characteristics of solar wind-Venus interactions (Ma et al. 2013).

3D Parametric Model of the Bow Shock
To create the three-dimensional parametric bow shock model, we adopt a versatile conic section function developed by Ramstad et al. (2017), which is represented as follows: where R BS , x, x 0 , e, and α are the radial distance from the symmetry X-axis, the location in X-axis, the subsolar standoff distance, the eccentricity of the hyperbola, and the scale length which is equal to α = 33.5 R V , respectively.In order to account for the three-dimensional irregularity of the bow shock, as described in our previous model for Mars ( Wang et al. 2020bWang et al. , 2022)), we expand upon the parameter governing hyperbola eccentricity, as follows: 2


Here, f is the azimuthal angle measured from the positive Zaxis of the VSO coordinate system in polar coordinates.As such, six parameters, denoted as x 0 , e 0 , e 1 , e 2 , e 3 , and ω, are employed to ascertain the bow shock locations.Specifically, x 0 stands for the bow shock subsolar standoff distance which is the distance from the point where the bow shock intersects the Venus-Sun line to the center of Venus.The parameter e 0 represents the degree of flaring of the bow shock tail, functioning in a manner akin to its role in the model presented Ramstad et al. (2017).The parameter e 2 characterizes the north-south asymmetry of the bow shock.If e 2 > 0, the flaring degree in the northern hemisphere surpasses that in the southern hemisphere.Conversely, when e 2 < 0, the situation is reversed.Similarly, e 3 characterizes the dawn-dusk asymmetry.Parameters ω and e 1 incorporate quadratic terms into the asymmetry.In the case of a bow shock exhibiting rotational asymmetry, we suggest that its cross section is not circular but rather takes on the shape of an ellipse, with the elongation along the major axis of the ellipse.The parameter ω represents the angle between the major axis of the elliptical cross section and the X-Z plane within the VSO coordinate (or equivalently, the Z-axis within the Y-Z plane).Here e 1 quantifies the extent of elongation.Note that a positive e 1 , with a corresponding ω, represents the same elliptical cross-section shape as a negative e 1 of the same magnitude paired with ω ± π/2.In this study, we chose to set e 1 > 0 to better elucidate the influence of solar wind parameters on ω.We use a set of the parameters of x 0 , e 0 , e 1 , e 2 , e 3 , and ω to describe one of the shape and size of the Venusian bow shocks with the north-south, dawn-dusk, and rotation asymmetries.At last, according to Wang et al. (2020b), we suggest that the effect of the IMF cone angle on the XY plane μ XY on the dawn-dusk asymmetry e 3 corresponds to the μ XZ effect on e 2 .That is, for the same value of μ XY and μ XZ , e 3 equals to e 2 .Moreover, the influences of B Y and B Z on x 0 , e 0 , and e 1 should also be the same.Hence, in this work, we mainly focus on μ XY (or B Y ) effect on the geometry of the bow shock, and the result of the μ XZ (or B Z ) effect is inferred from μ XY (or B Y ).
Subsequently, we apply the model function to correlate it with the bow shock positions under different solar wind conditions.It is worth mentioning that, to avoid any the deviation at the subsolar region, the value of x 0 is extracted directly from the determined bow shock location.The remaining parameters, e 0 , e 1 , e 2 , e 3 , and ω, are obtained by fitting the bow shock locations using this x 0 .
The determination coefficient, denoted as R 2 , quantifies the extent to which the variance in the dependent variable can be explained by the independent variable(s).A value of R 2 = 1 represents a perfect fit.For the entire set of 110 simulations in this study, R 2 quals 0.9969, indicating a high level of agreement between the chosen bow shock shape and the MHD model results.In the following section, we delve into the precise influence of solar wind parameters on the dimensions and shape of the Venusian bow shock.

Model Results
Figures 2(a1)-(d1) show the effects of the magnetosonic Mach number on the subsolar standoff distance x 0 (squares in panel (a)), the flaring degree e 0 (circles in panel (b)), the elongation degree in the twist direction e 1 (upward triangles), the degree of north-south asymmetry e 2 (downward triangles), the degree of dawn-dusk asymmetry e 3 (right triangles in panel (c)), and the twist direction ω (diamonds in panel (d)), respectively.In Case 1-40, the magnetosonic Mach number varies with solar wind density and velocity.The IMF conditions are B X = 3.6 nT, B Y = −5.4nT, and B Z = 0 nT for the cases shown in Figures 2( a1)-(d1).Additional information regarding the simulations and fitting procedures can be found in the supplementary document.The solid lines illustrate the expected relationships between the solar wind parameters and the bow shock construction parameters of x 0 , e 0 , e 1 , e 2 , e 3 , and ω.In these figures, it is evident that x 0 , e 0 , e 1 , and e 3 decrease with increasing M MS , following hyperbolic tangent functions.Meanwhile, the twist direction ω displays minimal sensitivity to changes in M MS .It is important to note that, in panels ((a1)-( d1)), where B Z = 0 nT for all cases, the degree of north-south asymmetry e 2 keeps close to 0. As mentioned above, if the B Z component has a nonzero value, we speculate that the corresponding μ XZ would induce a nonzero value of e 2 , which should also follow hyperbolic tangent relationships with M MS .In addition to analyzing the effects of M MS for this specific IMF condition, we have also conducted a study for other IMF conditions (refer to the supplementary document).The final relationship functions will be provided in the following Equations (3)-(8).
Figures 2(a2)-(d2) show the influences of the IMF components on x 0 (squares in panel (a)), e 0 (circles in panel (b)), e 1 (up triangles), e 2 (down triangles), e 3 (right triangles in panel (c)), and ω (diamonds in panel (d)).All these simulations maintain a constant solar wind density of n = 9.0 cm −3 and a velocity of U X = −448.0km s −1 spanning from Case 41 to 64 (refer to the supplementary document for details).The different colors of red, green, and blue correspond to IMF B t (Case 41-48), B X (Case 49-56), and B Y (Case 57-64), respectively.The figures reveal that x 0 , e 0 , and e 1 increase with the growing intensity of B t and B Y , while the effect of B X is negligible.The parameter e 3 also exhibits an increase with the intensification of B Y when B X has a nonzero value.Similar to panel (1), e 2 remains close to 0 as B Z = 0 nT for the cases.Meanwhile, ω seems to be more closely related to the IMF direction rather than the IMF magnitude.
Figures 2( a3)-(d3) show the effects of the IMF cone angle μ XY (in radians) on the bow shock parameters.Parameters x 0 , e 0 , and e 1 can be well approximated by quadratic dependence on IMF cone angles μ XY .That is, x 0 , e 0 , and e 1 are all enhanced with the increase of μ XY ä ( −π, −π/2), and decreased with the increase of μ XY ä ( −π/2, 0).There seems to be no obvious dependence of parameters of ω on μ XY .In addition, μ determined by B Y and B Z has different influences on the asymmetry of the bow shock.That is, the dawn-dusk asymmetric degree e 3 has an obvious sine and cosine correlation with μ XY , and a similar relationship is also expected between the north-south asymmetric degree e 2 and μ XZ .
Figure 3 displays the impact of variations in solar activity on the configuration of the bow shock.The parameter x 0 -solarmax linearly increases with x 0 -solarmin, while a similar relationship is observed between e 0 -solarmax and e 0 -solarmin.However, the other parameters e 1 , e 2 , e 3 , and ω do not display a distinct correlation in the context of solar minimum and maximum conditions; hence, we do not present them here.
Finally, we compile the outcomes of our analysis based on the entire set of 110 simulations, revealing the correlations between the bow shock parameters and various solar wind and IMF conditions.The results include the best-fitting functions for these parameters: The numerical values of the 27 coefficients in these equations are given in Table 2.Here B X , B Y , and B Z in nT, λ and μ in radians, x 0 in R V .The fitting determination coefficients R 2 for Equations (3)-( 10) are 0.9648, 0.9165, 0.8042, 0.8712, 0.8712, 0.9975, 0.8822, and 0.6981, respectively.In addition, we also calculate the R 2 model and the root mean squared error (RMSE) associated with the identified locations of bow shock from simulations and the final model functions (Equations ( 3)-( 10) and Table 2) for all the 110 cases which Figure 3. Solar activity effect on the geometry of the bow shock.Parameter x 0 -solarmax is the subsolar standoff distance under the solar maximum condition in Case 87-110, and x 0 -solarmin is the subsolar standoff distance under the same solar wind condition as that of Case 87-110 but for solar minimum condition.Parameters e 0solarmin and e 0 -solarmax are the flaring degree under the corresponding conditions.3)-( 10) and Table 2) for all the 110 cases.The black dashed lines represent the average values.
have also been listed in the supplementary document as shown in Figure 4.The RMSE represents the square root of the second sample moment of the differences between predicted values and the observed values or the quadratic mean of these differences.A smaller RMSE generally indicates a better fit and a larger R 2 .In our analysis, it is noticeable that although R 2 model is correspondingly less than R 2 which is directly fitted from simulation in most cases, the average value of R 2 model reaches 0.9961, as indicated by the black dashed line of Figure 4(a).This value indicates a good fit and the overall accuracy of the model results.Moreover, the average value of RMSE is 0.07986, as denoted by the black dashed line in Figure 4(b).This signifies that the average quadratic mean of the differences between the bow shock locations and the model's predictions is 0.07986 R V .It's worth noting that, in very few cases, both R 2 model and RMSE deviate significantly from the average, which may be attributed to errors in the simulation or the fitting process.In general, both R 2 model and RMSE demonstrate a good fitting of the model result.
The dawn-dusk asymmetry in the bow shock position is mainly controlled by the orientation of the IMF when the IMF lies on the X-Y plane, as shown by Figure 5(a1) and (a2).Specifically, the distance of the parallel bow shock is closer to the planet compared to the perpendicular bow shock.Furthermore, Figure 5(c3) and (c4) along with Equations (8) demonstrate that the cross section of the bow shock tail is elongated in the direction perpendicular to IMF.From Figure 5((a2)-( c2)), we deduce that the IMF component B X has a relatively small influence on the geometry of the bow shock.In these cases, the flaring degree e 0 slightly increases with the magnitude of B X .However, the degree of north-south asymmetry (e 2 ) and dawn-dusk asymmetry (e 3 ) are barely affected by B X .In Figure 5((a3)-( c3)), it is evident that the bow shock moves farther away from Venus, and the flaring degree increases as B Y or B Z increases.For B Y dominating cases the cross section of the bow shock is elongated in the Z direction.

Discussions
In hydromagnetic theory, the position of a planetary bow shock is determined by the ratio of the upstream density to the downstream density, as described by the Rankine-Hugoniot relations.This relationship is associated with Mach numbers.Additionally, the limiting Mach cone angle of a bow shock is expected to decrease as the Mach number increases, as suggested by earlier studies (e.g., Spreiter et al. 1966;Slavin & Holzer 1981;Verigin et al. 2001).Therefore, it is widely accepted that when the Mach numbers increase, the location of the bow shock tends to move closer to the planet.This phenomenon has been observed in various studies (e.g., Edberg et al. 2010).
The influence of IMF on the bow shock is connected to both Mach numbers and the pileup of IMF around Venus.First, a weaker IMF corresponds to a higher Mach number, and this relationship can lead to the compression of the bow shock (Zhang et al. 2004;Chai et al. 2014;Signoles et al. 2023).This influence is evident on the change of B Y and B Z .However, the B X is independent of the MHD fast mode speed at the subsolar region, as discussed by Wang et al. (2020b).Therefore, x 0 remains relatively constant with variations in B X .In the tail region, the response to B X is different, as the bow shock moves outward when B X increases.Second, the phenomenon of IMF pileup around the ionosphere primarily depends on the values of B Y and B Z since the IMF can only be compressed in the perpendicular direction.This, in turn, affects the size of the obstacle, namely the induced magnetosphere (Xu et al. 2022).However, an increase in B X does not enhance the IMF pileup around Venus to the same extent, resulting in a relatively minor effect on the bow shock.
The influence of the IMF orientation can also be explained by the fast mode speed v MS .This speed is highest when the IMF direction is perpendicular to the local shock normal but lowest for the parallel situation (Russell et al. 1988;Wang et al. 2020a).As a consequence of this, on the Y-Z plane, the cross section of the Venusian bow shock tail becomes elongated in the direction perpendicular to the IMF projection vector.Furthermore, the degree of this elongation (e 1 ) becomes more pronounced with increasing strength of the IMF components B Y and B Z , but it decreases with increasing M MS .
The impact of these effects is prominently demonstrated in Figure 5. Notably, the asymmetry in the bow shock is particularly evident on the terminator plane, as depicted in Figure 5(c1)-(c4).On this plane, the cross section of the bow shock is notably elongated in the direction perpendicular to the IMF direction.This phenomenon is analogous to what is observed in Earth's bow shock (e.g., Sibeck & Lin 2014;Wang et al. 2016Wang et al. , 2018)).In our parametric model, the degree of asymmetry is approximately 0.5 R V , which is somewhat larger than the roughly 0.1 R V reported by Zhang et al. (1991).
In Figure 6, our parametric model is compared to the empirical models for both solar maximum and minimum conditions.Our parametric model demonstrates a reasonably good agreement with the average location of the empirical models, with the exception of a slightly larger flaring angle at solar maximum.
Compared with solar minimum, the bow shock moves farther away from Venus for solar maximum condition, aligning with prior research findings (Shan et al. 2015;Signoles et al. 2023).Specifically, L x=0 and x 0 exhibit a decrease of 16% and 7%, respectively, when transitioning from maximum to minimum.The difference of L x=0 in our study is similar to the approximately 14% change reported by Alexander & Russell (1985), but greater than those found by Shan et al. (2015) and Signoles et al. (2023).

Conclusions
In this work, we have developed a three-dimensional parametric model of the Venusian bow shock by employing a 3D single-fluid MHD model of Ma et al. (2013).This model incorporates the Mach number and the IMF as controlling The solar wind conditions is set to be n = 17.0 cm −3 , U X = −400.0km s −1 and B = 12.1 nT at solar maximum (Xu et al. 2022), and n = 6.9 cm −3 , U X = −384.0km s −1 and B = 6.5 nT at solar minimum (Xu et al. 2023).The IMF Parker spiral angle is considered as 42°near Venus (Chang et al. 2019).Black solid lines and red dashed lines represent our parametric model results on the X-Y and X-Z planes.The results of Slavin et al. (1980), Russell et al. (1988), Zhang et al. (1990) factors.To derive this parametric model, we conducted a series of 110 simulations, each under different solar wind and IMF conditions.The solar wind parameters were explored within the ranges of M MS ä [2.8, 10.1] and IMF strength ä[2.5, 20] nT.
Our model has provided several key insights into the size and geometry of the Venusian bow shock: (1) The size of the Venusian bow shock is primarily influenced by the M MS .An increase in M MS results in the bow shock moving closer to Venus and a reduced flaring angle of the shock.
(2) Enhancement of B Y or B Z lead to the bow shock moving farther away from Venus with an increased flaring degree.Although the influence of B X is relatively small, it does impact the flaring degree.
(3) The cross section of the bow shock tail is elongated in the direction perpendicular to the IMF.This elongation degree increases with higher intensities of B Y and B Z .In cases where B Y (B Z ) dominates, the cross section of the bow shock is elongated in the Z (Y) direction, contributing to the magnetic pole-equator asymmetry as previously reported by Chai et al. (2015).The dawn-dusk asymmetry is linked to the IMF cone angle when the IMF lies in the X-Y plane.These effects are attributed to the properties of the fast magnetosonic wave, which align with previous discussions by Zhang et al. (1991) and Russell et al. (1988).
In summary, our parametric model demonstrates strong alignment with prior empirical models under average solar wind conditions.However, our study only focused on modeling the Venusian bow shock under solar maximum and minimum conditions.To create a more comprehensive model, the continuous influence of EUV radiation should be incorporated into the model, as it significantly affects the Venusian ionosphere and thus the bow shock position.The model currently does not address asymmetry along the convection electric field, a factor that plays a crucial role in the Venus-Sun interaction.To account for this, future studies could involve multifluid models (Li et al. 2023), which are more adept at capturing such complex interactions.The parametric model is to be further refined in the future.

Figure 1 .
Figure 1.Example of identified locations of the Venusian bow shock from simulation data on 3D (panel (a)), Y-Z plane (panel (b)), X-Y plane (panel (c)), and X-Z plane (panel (d)) in the VSO coordinate system.The solar wind conditions for this example are n = 6.0 cm −3 , U X = −400.0km s −1 , B X = 3.6 nT, B Y = −5.4nT, and B Z = 0 nT.The color code is the solar wind density (cm −3 ) in Figure (b), (c), and (d).

Figure 2 .
Figure 2. The dependence of the bow shock parameters x 0 (squares in panel (a)), e 0 (circles in panel (b)), e 1 (up triangles, panel (c)), e 2 (down triangles, panel (c)), e 3 (right triangles in panel (c)), and ω (diamonds in panel (d)) on the magnetosonic Mach number M MS (panel (1)), IMF intensity (panel (2)), IMF cone angle μ XY (panel (3)), and IMF clock angle λ (panel (4)).The magnetosonic Mach number effects are plotted from Case 1 to 40 in the supplementary document, the IMF intensity effect shows Case 41-64 (In Case 41-48 B t varies with B X and B Y when B Z = 0, in Case 49-56 B X changes when B Y = B Z = 0, and B Y varies in Case 57-64 when B X = B Z = 0), the IMF cone angle effect displays Case 65-74, the IMF clock angle effect is from Case 75-86.Symbols show the results of MHD simulations and the lines are the potential fits.

Figure 4 .
Figure 4. R 2 model (a) and RMSE (b) associated with the identified locations of bow shock from simulations and the final model functions (Equations (3)-(10) and Table2) for all the 110 cases.The black dashed lines represent the average values.

Table 1
Supplementary Information of the Simulation Cases

Table 2
The Coefficients of the Parametric Venusian Bow Shock Model Function in