Influence of an Upstream Dynamic Pressure Pulse on Venus–Solar Wind Interaction by Multispecies Magnetohydrodynamic Simulation

As an unmagnetized planet, Venus does not have a global dipole magnetic field like Earth. The solar wind interacts directly with the Venusian ionosphere. In this paper, the interaction between solar wind and the Venusian upper atmosphere is investigated by conducting a three-dimensional multispecies magnetohydrodynamic (MHD) model, in which the Venusian ionosphere is developed self-consistently by considering chemical reactions among four major ion species (H+, O2 +, O+, and CO2 +). Based on this developed MHD model, we mainly study the effects of the impulsive solar wind dynamic pressure on the interaction process. Simulation results reveal that the locations of the bow shock and magnetic pileup boundary show an abrupt compression, then expansion, and finally shrinkage to reach a new steady state with the influence of the impulsive dynamic pressure. The recovery time of the large-scale magnetic field is a few minutes. In addition, the peak escape fluxes are enhanced by about 1 order of magnitude due to the increased solar wind dynamic pressure pulse.


Introduction
Venus, as a terrestrial planet, possesses neither a global intrinsic magnetic field like a geomagnetic field nor a local crustal field like on Mars (Phillips & Russell 1987).Therefore, different from Earth but similar to Mars, the solar wind can interact with the Venusian atmosphere/ionosphere much more directly, which puts significant impact on its atmosphere and climate changes (Ma et al. 2008).
The supersonic solar wind plasma is thermalized and deflected around Venus when approaching its surroundings, carrying the interplanetary magnetic field (IMF) to pile up around the obstacle.As a result, there are three important plasma boundaries formed in this interaction process: bow shock (BS; Slavin et al. 1979;Xue & Jin 2014;Shan et al. 2015), magnetic pileup boundary (MPB; Bertucci et al. 2003;Kallio et al. 2008), and ionopause (Chong et al. 2017;Han et al. 2020).Previous studies have revealed that the position of these boundaries is dependent on the solar wind dynamic pressure and solar activity (Russell et al. 2006).Analyzing the data observed by Venus Express (VEX), Zhang et al. (2008) suggested that the BS location during solar minimum conditions is 2.14 R V near the terminator and 2.40 R V during solar maximum conditions.The shock location seems not sensitive to the variation of the dynamic pressure, especially during the solar minimum by ASPERRA-4 on VEX (Martinecz et al. 2008).There is a linear dependence of subsolar and terminator position of Venusian BS on sunspot number between different solar cycles (Whittaker et al. 2010).The location of the crossings can be best fitted by the conic section or the parabola (Futaana et al. 2017).The pole-equator asymmetry of the shock can be observed with the solar wind flow nonperpendicular to magnetic field.The north-south asymmetry of the shock is observed for the increase of the obstacle size by the excessive pickup ions accelerated by the convective electric field (Chai et al. 2015).
Early Pioneer Venus Orbiter (PVO) observations revealed that the ionosphere in the terminator region expanded in response to extreme low solar wind dynamic pressure.New observations by VEX completed previous results with the discovery of a teardrop-shaped ionosphere due to the tailward extension of nightside ionosphere (Wei et al. 2012).Statistical studies have found a positive correlation between ion escape rates and upstream solar wind dynamic pressure for Venus (Luhmann et al. 2007;Edberg et al. 2011).The ASPERA-4 Ion Mass Analyzer and Magnetometer on VEX provided observations around Venus encountering the coronal mass ejection (CME) disturbance (Luhmann et al. 2008).The ion escape rate is enhanced by a factor of 1.9 on average during high solar wind dynamic pressure events (Edberg et al. 2011).The outflow of ions is significantly enhanced (10 times higher) in the Martian tail for the response to the impulsive event (Futaana et al. 2008).It was inferred that O + and O 2 + are the dominant heavy ion escape species according to the measurements of the Venusian upper atmosphere/ionosphere and the antisunward ion fluxes in the nightside wake (Brace et al. 1995;Luhmann et al. 2006).The observations made it possible to investigate how the ionized oxygen escape is influenced by the varied solar activity (Persson et al. 2018).The observations and simulations confirm a higher ion-loss rate for the enhanced solar EUV radiation and solar wind dynamic pressure for Mars (Lundin et al. 2008;Ma et al. 2014).Solar energetic particle and corotation interaction region (CIR) events resulted in around a 1 order of magnitude increase of the heavy ion outflow flux from the Martian atmosphere (Dubinin et al. 2009).The Martian magnetosphere and ionosphere were significantly compressed by enhanced solar wind dynamic Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
There are several numerical models studying the interactions of the solar wind with Venus, each of them has their strengths and limitations.The gas-dynamic computational model developed by Spreiter & Stahara (1980) assumed that the ionopause is the obstacle to the solar wind flow.However, Zhang et al. (1991) revealed that the actual obstacle is the combination of the Venus ionosphere and the magnetic barrier, which cannot be simulated by gas-dynamic models because of the lack of a magnetic induction equation (Russell et al. 2006).The hybrid model treats the individual positive ions as particles and the electrons as massless fluid (Kallio & Jarvinen 2012).The position and velocity of each ion are obtained by the integration of the equations of motion.By assuming all the ion species share the same velocity and temperature, the multispecies single-fluid magnetohydrodynamic (MHD) model calculates all the major ion species in the Venus ionosphere, treating the ionosphere as a perfectly conducting sphere in simulation (Terada et al. 2009;Ma et al. 2013Ma et al. , 2020)).The multispecies multifluid MHD model (Dang et al. 2022) solves the continuity, momentum, and energy equations for each ion species considered in addition to the magnetic induction equation, assuming that the plasma is collisional and quasineutral.
Due to the upstream time-dependent solar wind conditions, the Venusian ionosphere is expected to change from one quasisteady state to another due to external drivers, particularly solar wind dynamic pressure.Such a change has been clearly shown by the PVO observations for three consecutive orbits: 175, 176, and 177 (Elphic et al. 1980).Magnetization of the Venus ionosphere occurs under high solar wind dynamic pressure by PVO measurements (Woo & Kliore 1991).The hydrogen ion abundance near solar minimum is more than 1 order of magnitude less than that observed at solar maximum by the Ion Mass Spectrometer on PVO (Donahue & Hartle 1992).Besides, net escape rates tend to be higher in the solar minimum than those in the solar maximum, during which return flows are more frequently observed due to the multiple reversals of the B x field component in the magnetotail (Masunaga et al. 2019).The increased solar wind dynamic pressure tends to decrease the densities of the main planetary ion species (except for CO 2

+
) and enhance the outflow velocity for all the ions with a near-horizontal magnetic field on Mars (Niu et al. 2021).However, the single-spacecraft observation limits a full understanding of how the ionosphere and ion escape change with the solar wind drivers.Ma et al. (2014) estimated the response of the Martian ionosphere to the external solar wind dynamic pressure changes.The evolution of the large-scale magnetic field in the Venus ionosphere using a global multispecies MHD model was investigated by Ma et al. (2020).It should be noted that the identical change in velocity can cause larger enhancement to the solar wind dynamic pressure, and thus leads to more intense influence on Venusian large-scale boundaries and ion escape in the tail region.
This study adopts a global multispecies single-fluid MHD model to investigate the influence of solar wind dynamic pressure pulse on the Venusian plasma environment by increasing only the velocity to 2 times.The study of the changes of boundary positions and tailward ion escape flux is helpful for understanding the impact of dynamic upstream conditions on the solar wind interaction with Venus.In addition, in order to evaluate our MHD model, one more case was conducted to compare with the VEX observations.

Method
In this study, a 3D multispecies MHD model is used to simulate the interaction of the Venusian space environment with the solar wind, in particular the changes caused by a solar wind pressure pulse event.The four main ion species chosen in our model are H + , O 2 + , O + , and CO 2

+
, which account for a large proportion of ion particles in the Venusian ionosphere.In order to solve the mass equation for each ion component separately and the single momentum and energy equation for all ion species, the decoupling method is adopted.The decoupled continuous equation can be expressed as follows: where W s is the mass density of each component and F s is the tensor flux formed by the momentum vector of each component.The expressions of W s and F s can be written as follows: where ρ 1 , ρ 2 , ρ 3 , and ρ 4 represent the mass density of H + , O 2 + , O + , and CO 2 + .u is the average ion velocity.Q s represents mass variations caused by chemical reactions among all the ion species, which is expressed as follows: The governing equation can be expressed as The conservative vector W is composed of eight conserved variables: where ρ is the total mass density of the four ion species , p is the total ion species thermal pressure, and γ is the ratio of specific heats taken to be 5/3.The expressions for the conserved flux F and the source term where G is the Venusian gravity, ν is the ion-neutral collision frequency: ), T o is the temperature of newborn ions assigned the same value as that of the local neutral atmosphere: here, we assume T o = 400 K.The multispecies MHD model in our paper solves the continuity equation for the four species decoupled from the main governing equation, different from solving the total set of equations consisting of all the variables for Venus by Ma et al. (2013).
The neutrals in the Venusian atmosphere considered in this study include O and CO 2 , their densities are based on the model calculations of Fox & Sung (2001).For simplification, spherical symmetric profiles of neutrals are adopted: where H is the altitude from the Venus surface measured in km.
A self-consistent three-dimensional Venusian ionosphere could be generated because our model simulations take consideration of chemical reactions, including photoionization, charge exchange, and dissociative recombination.The associated chemical reactions and corresponding reaction rates refer to Table 1 by Ma et al. (2013), concluding from the work of Schunk & Nagy (2009) and Fox & Sung (2001).
The cosine function was used to obtain the solar zenith angle effect in calculating photoionization rates of the neutrals for the dayside region: = q k CO cos SZA , 11 For the nightside region (SZA > 90°), the photoionization rates are set to zero.Our calculations were performed in the Venus-centered Solar Orbital (VSO) coordinate system, with the x-axis pointing from Venus toward the Sun, the z-axis perpendicular to the x-axis and positive toward the north celestial pole, and the y-axis completing the right-handed coordinate system.The computational domain is set to be −24 R V x 8 R V , −16 R V y, and z 16 R V , where R V is the radius of Venus (R V = 6052 km).
The inner boundary of the calculation grid is set at an altitude of 100 km above the surface of Venus, where O 2 + , O + , and CO 2 + densities are taken to be values in photochemical equilibrium.The density of H + is set to be 30% of the solar wind proton density (Taylor et al. 1980).Assuming that the fluid is inviscid, the normal velocity value at the inner boundary is set to be zero.
In this paper, two simulation cases are conducted.Taking the solar maximum condition as inputs, case 1 is devoted to studying the space environment of Venus, particularly the influence of a solar wind dynamic pressure pulse event.After case 1 achieves a quasi-steady state, we increase the solar wind velocity by 2 times at the simulation time of 01:00:00, as shown in Figure 1.Using VEX observation parameters as another inputs, we run case 2 for validating our simulation model results by comparing with VEX data.VEX is aimed at improving knowledge of the Venus plasma environment with its advantages in instruments, orbit coverage, and lower-altitude observations compared with previous PVO observations.The MAG instrument on board VEX owns a wide dynamic range between ±32.8 and ±8388.6 nT with the corresponding resolution between 1 and 128 pT (Zhang et al. 2006).The magnetometer MAG on board VEX measures the magnetic field with a frequency of 128 Hz and the autonomous mode is the standard mode set at a 1 Hz data rate (Zhang et al. 2006).

Results
Since our global Venus-solar wind interaction model evolves self-consistently, a steady scenario can be obtained through converging in time with specific initial and boundary conditions as described above.The simulated magnetic field configuration and average ion velocity distribution results for case 1 in the x-z and x-y planes are demonstrated in Figure 2.Moreover, the white lines in the magnetic field plots represent the field lines and indicate the streamline of ion flow in the velocity plots; the black dashed line in each panel shows the mean observed BS location from Zhang et al. (1990).As can be seen in Figure 2, the BS boundary is clearly visible in both magnetic field and velocity configurations, which fit quite well with average PVO observed positions (shown by the dashed line in each panel; Zhang et al. 1990).In the contour plot of Figure 2(c), it is obvious that the magnetic field lines piled up to form the magnetic pileup region (MPR) and the current sheet shifts to the duskside (−Y direction) due to the initial 56°P arker spiral IMF.However, the average ion velocity distributions seem to be symmetrical in the x-z and x-y planes.
In order to verify the reliability of the MHD model we conducted in this paper, simulation results were compared with the observation data of VEX before investigating the impact of the solar wind pulse event.Figure 3   simulation results in Figure 4 that there is a sharp increase of B X and B across the BS at around 1:24 UT, which is consistent with the observations.Both observed and modeled B X and B show a gradual increase during 1:24 UT-1:34 UT, indicating that the VEX went into sheath region and crossed the MPB.Around 1:55 UT, the direction of the observed B X suddenly changes from positive to negative, meaning the VEX crossed the Venusian magnetotail current sheet, which was reproduced by our model approximately.From 2:38 UT to 2:50 UT, VEX departed from the southern nightside-induced magnetosphere.Finally, it went into the solar wind at 2:50 UT.The variations of the magnetic field vector were well reproduced by our simulation model results extracted from the whole trajectory, except those between 1:55 UT and 2:38 UT; while the model predicts drops to zero gradually, the observation shows a sharp decrease to zero.This discrepancy can cause a time delay of nearly 20 minutes and also leads to the sharp decrease and increase of B during the time interval from 1:55 UT to 2:17 UT.Meanwhile, simulation results of B y and B z agree well with data observed by VEX.In general, the locations of the inbound/outbound BS and current sheet as well as the peak strength of the magnetic field were all reproduced by our model conducted in this paper.Therefore, it is reasonable and convincing that we use this model to investigate the changes of the Venusian space environment by the solar wind dynamic pressure pulse event.
Since the formation of these boundaries is associated with the solar wind, the solar wind dynamic pressure can play an important role in the position variation of BS and MPB. Figure 5(a) demonstrates pressure distributions along the Sun-Venus line in case 1 simulation results.P D is the solar wind dynamic pressure; P T is the solar wind thermal pressure; P B is the magnetic pressure; P Total is the total pressure (P Total = P D + P T + P B ).The small solid points located at each line with the same color indicate the positions of grid points.Given the gas-dynamic consideration mentioned above, the BS can be observed clearly in the location where the solar wind dynamic pressure balanced with the thermal pressure at 1.4 R V .Moreover, the MPB can form around the region where the solar wind thermal pressure transforms to the magnetic pressure at 1.24 R V .The distance of the BS and MPB represented in Figure 5(a) is basically consistent with the PVO and VEX observations (Zhang et al. 1990) as well as previous simulation results (Ma et al. 2013).
To investigate how the Venusian space environment is influenced by the solar wind dynamic pulse event, we adopt a dynamic pulse as input at 01:00:00 UT, when case 1 is already in a quasi-steady state (shown in Figure 1).In order to examine the variation of boundary locations during the event, Figures 5(b)-(i) show various pressure profiles along the subsolar line for different times, which plots dynamic pressure P D (in green), thermal pressure P T (in red), magnetic pressure P B (in blue), and total pressure P Total (in black) with the grid point (the same color solid circle) overlapped on each line.The equilibrium state in Figure 5(a) is ensured by the basically unchanged positions of the BS and MPB with time.The BS and MPB locations in Figures 5(b)-(i) are obtained during the dynamic process activated by the dynamic pressure pulse.Most panels in Figure 5, except for panel (c), show that the total pressure changes slightly except for abrupt variation in the region around the planet and shock, consistent with previous MHD results (Ma et al. 2004).The total pressure conservation in the BS is broken due to the dissipation caused by magnetic and fluid viscosity (Gombosi 1998).The total pressure near the planet surface is not conserved due to the ionosphere chemical reactions.The influence of chemical reactions may extend between the ionosphere and BS, resulting in a small change of the total pressure.The total pressure within the BS evidently exceeds the solar wind dynamic pressure during the early influence of the pulse as shown in Figure 5(c).This is caused by the transient large shrink of Venusian boundaries due to the abrupt change of external solar wind condition.
Before the arrival of the pressure pulse, the steady locations of the BS and MPB are 1.40 R V and 1.24 R V , respectively (Figure 5(a)).It takes some time (about 1 minute) for the pressure pulse to reach the Venus BS outer boundary (see panel (b)).The sudden increase of dynamic pressure breaks the equilibrium state, forcing the BS and MPB to shrink toward the planet.The solar wind dynamic pressure enhancement causes the locations of the BS and MPB to move inward to 1.25 R V and 1.16 R V in less than 14 s, respectively (see panel (c)).Both the magnetosheath region and MPR respond quickly to the pressure pulse event.Since the solar wind dynamic pressure is almost converted to thermal pressure in the sheath region, the thermal pressure inside the sheath region increases as well as the magnetic pressure in the MPR.The timescale predicted by our model is similar to previous studies (Romanelli et al. 2019;Ma et al. 2020).The location of the BS rebounds a little and settles at 1.35 R V and 1.39 R V , while MPB locates at 1.15 R V and 1.18 R V (see panels (d) and (e)).From panel (a) to panel (i), it can be concluded that when the solar wind dynamic pressure enhances from 4.5 to 18 nPa, the plasma thermal and magnetic pressure shows a sharp increase and then gradually decreases.Figures 5(f)-(i) show that the BS and MPB are finally compressed to reach a new steady state with the BS and MPB locating at 1.35 R V and 1.16 R V , respectively.
Figure 6 presents the distributions of the density of O 2 + and the velocity with θ at altitudes of 200 km and 350 km in the x-y plane, respectively.θ is the angle with a subsolar line with a positive degree for the +y hemisphere and a negative degree for the −y hemisphere of Venus.The simulation results at 1 hr, 1 hr 1 minute 19 s, and 2 hr are represented by red, green, and blue lines, respectively.At the altitude of 200 km, the magnitude of the O 2 + density is about 7 × 10 4 cm −3 around the subsolar point and gradually decreases away from the subsolar point (larger solar zenith angle, SZA).The density is enhanced slightly at large SZA at 1 hr 1 minute 19 s.At the new steady state, the density is increased by a factor of 2 to 3. The velocity at an altitude of 200 km is about 0.2 ∼ 0.5 km s −1 at initial steady state and increases abruptly at 1 hr 1 minute 19 s by a factor of 4 to 6 at small and moderate SZA.The velocity at 2 hr is slightly higher than that at 1 hr.At 350 km, the magnitude of the density is about 2 orders of magnitude less than that at 200 km around the subsolar region.Additionally, the density at 350 km around the moderate SZA (around 50°) is smaller than that at the subsolar point and the terminator by about 1 or 2 orders of magnitude.At 1 hr 1 minute 19 s with the influence of the dynamic pressure pulse, the density is slightly more decreased than the initial state, different from that at 200 km.The velocity is enhanced by about 6 times on average at 300 km, and the variation with time is similar to that at 200 km.
To investigate the variation of the Venusian ionospheric plasma escape, Figure 7 demonstrates the ion fluxes in the logarithmic scale in the x-y plane at 2 hr when a new quasisteady state is reached after the solar wind pulse event.It is obvious that the most significant escape channel is the Venusian magnetotail, accompanied by the dawn-dusk asymmetry.The dawn-dusk asymmetry of the escape flux is caused by the IMF B x component (Kallio et al. 2006).The ion escape fluxes of the dusk-tail flank are much more significant than the counterpart of the dawn-tail.Moreover, O 2 + is the dominant escape ion and CO 2 + is the least significant one.In order to analyze the impact of the solar wind event on the ion escape quantitatively, Figure 8 shows the variation of the ion escape flux along the y-axis at x = −2 R V .It can be seen from each panel that the ion escape flux at the magnetotail is mainly concentrated in −4 R V y 4 R V .Meanwhile, the peak of the escape fluxes is located at y = −0.74R V at 1 hr, i.e., the maximum escape fluxes of O 2 + , O + , and CO 2 + are 1.75 × 10 7 cm −2 s −1 , 1.17 × 10 6 cm −2 s −1 , and 8.72 × 10 4 cm −2 s −1 , respectively.The steady-state magnitude of 10 6 cm −2 s −1 for the escape flux of O + is consistent with the hybrid simulations by Terada et al. (2004).At 01:08:12 of the simulation time, which the new quasi-steady state reached after the solar wind pulse event, the peak O 2 + , O + , and CO 2 + escape fluxes at y = −0.13R V are 4.38 × 10 8 cm −2 s −1 , 2.39 × 10 7 cm −2 s −1 , and 4.45 × 10 6 cm −2 s −1 , respectively.The peak escape ion fluxes are enhanced by about 1 order of magnitude for the three ion species.
Figure 9 presents the density of O 2 + in logarithmic scale and tailward velocity V x along x = −2R V at 1 and 2 hr.It can be seen that more O 2 + ions are concentrated in the central tail region at 1 hr.With the influence of the dynamic pressure pulse, the density in the steady state is increased by 1 order of magnitude around -0.5 R V y 0.3 R V .The velocity is slightly increased in the central tail region by a factor of less than 2. In the region |y| 4 R V , the tailward velocity is close to that of the solar wind.The increase of the ion escape flux by 1 order of magnitude as shown in Figure 8 is mainly contributed by denser ions in the central tail region.

Conclusion and Discussion
In this paper, the multispecies MHD model is established to investigate the Venusian space environment, especially its changes due to the solar wind pulse event.The model is under consideration of production and loss of the main ion species in the Venusian ionosphere, associated with chemical reactions among all species.Our simulation results effectively reproduce the well-known large-scale structures, such as BS, MPB, and MPR.The comparison between the magnetic field result from the MHD model and VEX observation along a single orbit reaches a good agreement.
BS and MPB experience compression, slight expansion, and a decay to the steady state.The sudden increase of dynamic pressure breaks the equilibrium state, forcing the BS and MPB to shrink toward the planet.At the same time, the kinetic energy transits to thermal and magnetic energy.The increased magnetic and thermal pressure pushes BS and MPB to bounce back after the compressibility reaches its limit.Meanwhile, previously stored magnetic and thermal energy are released to the solar wind.In the final shrinking stage, there is a slight increase of the thermal pressure.When the dynamic pressure pulse event arrives at the Venusian surroundings, the influences on the magnetosheath and MPR are very rapid.Moreover, the magnetosheath and MPR need 3 minutes 51 s to achieve a new steady state.This vibration of the Venusian plasma obstacle has  + in logarithmic scale and the velocity with θ at altitudes of 200 km and 350 km in the x-y plane, respectively.θ is the angle with the subsolar line with a positive degree for the +y hemisphere and a negative degree for the −y hemisphere.The red, green, and blue lines represent the simulation results at 1 hr, 1 hr 1 minute 19 s, and 2 hr, respectively.also been reported on Mars, stimulated by interplanetary shock with a quasiperiod of 1 minute and recovery time of several minutes (Futaana et al. 2006).In particular, Romanelli et al. (2019) found that the recovery timescales of the dayside Martian magnetosphere range between 8 s and 11 minutes for about 90°IMF rotation that lasted 50 s.Observations and simulations revealed the shrunken and less flared BS with a high Mach number for Venus and Mars (Ma et al. 2013;Edberg et al. 2010).Mach number also has an effect on the physical process related to the transport and dissipation of energy in the Mars solar wind interaction (Halekas et al. 2017).More investigations are needed to examine whether the recovery time of the BS and MPB under dynamic pressure pulse is related to Mach number.
Day to night transport is an important source for the nightside ionosphere and ion escape.The density of O 2 + and velocity in the dayside region are increased a few times by the dynamic pressure pulse.The ion escape flux through the Venusian magnetotail presents obvious dawn-dusk asymmetry.The ion escape fluxes in the dusk-tail flank are much more significant than those in the dawn-tail.Previous hybrid simulations obtained that the dawn-dusk asymmetry of the ion species is stronger at Venus than that at Mars (Jarvinen et al. 2016).O 2 + is the dominant escape ion and CO 2 + is the least significant one.Meanwhile, the increase of the solar wind dynamic pressure can lead to about 1 order of magnitude enhancement of the peak ion escape fluxes through the Venusian magnetotail at the steady state.The outflow of ions is significantly enhanced (10 times higher) in the Martian tail in the response to the impulsive event (Futaana et al. 2008).Our simulation results indicate that the increase of the peak ion escape flux by 1 order of magnitude in the Venusian tail is mainly contributed by 10 times denser ions in the central tail region.
The relations between ion escape rate and solar wind condition (e.g., the energy flux) can be used to evaluate the escape rate in the early time of Venus evolution (Persson et al. 2020).Astronomical observations indicated that the early Sun was generally more active than the present.The solar EUV condition was 10-1000 times higher (Ribas et al. 2005) and the solar wind was about 300 times thicker than today (Terada et al. 2009).The ionosphere might have been greatly heated and expanded by the extreme high solar EUV flux, causing a much more severe neutral escape of hydrogen and oxygen than the case for the cooled atmosphere today for Venus (Gillmann & Tackley 2014).Thus, the MHD simulation model can provide a more effective method to solve these problems, which merits further investigation.

Figure 1 .
Figure 1.Solar wind velocity used in the simulation.
shows the 3D trajectory of VEX on 2006 June 1 also chosen by Ma et al. (2013).It can be seen that VEX flies from point A to B and point C marks the time of closest approach at 125 km above the Venusian north pole.The comparison of the simulation results with VEX magnetometer data is demonstrated in Figure 4.The red lines represent the observation data of B x , B y , B z , and |B| as a function of time.The blue lines show simulation results extracted from the same satellite trajectory.It can be seen from

Figure 2 .
Figure 2. Contour plots of average ion velocity and magnetic field strength in the x-z and x-y planes.The dashed black curves show the average bow-shock location from Zhang et al. (1990).

Figure 3 .
Figure 3. Orbit of VEX on 2006 June 1 with the spacecraft position in the x-yz coordinate.

Figure 4 .
Figure 4. Comparison of simulation results with the VEX observations along the trajectory.

Figure 5 .
Figure5.Snapshots of pressure profiles along the subsolar line during the solar wind pressure pulse event.Pressure profiles of thermal pressure (red), dynamic pressure (green), magnetic pressure (blue), and total pressure (black) are plotted in each panel.

Figure 6 .
Figure 6.The distributions of the density of O 2+ in logarithmic scale and the velocity with θ at altitudes of 200 km and 350 km in the x-y plane, respectively.θ is the angle with the subsolar line with a positive degree for the +y hemisphere and a negative degree for the −y hemisphere.The red, green, and blue lines represent the simulation results at 1 hr, 1 hr 1 minute 19 s, and 2 hr, respectively.

Figure 7 .
Figure 7. Contour plots of escape flux of different ion species in the x-y plane.

Figure 8 .
Figure 8. Variation of escape flux of different ion species along x = −2R V during the solar wind pressure pulse event.