Bow shock is formed when the solar wind plasma flow moving through the interplanetary space is converted from a supersonic flow to a subsonic flow as it meets the obstacle formed by the planet. In the ideal compressible magnetohydrodynamic (MHD), the plasma frozen-in by different magnetic field lines is not mixed with each other. Bow shock, a collisonless shock, is formed as the magnetosphere formed by the earth’s dipolar magnetic field serves as an obstacle that interrupts the supersonic solar wind flow. The thickness of the shock wave is less than the mean free path of plasma. Additionally, by the interaction between the solar wind plasma and the planet that has its own magnetic field, the solar wind diverts around the planet, while the inner space is surrounded by the planetary magnetic field.

On the contrary, in the case of a planet such as Mars that has a very weak magnetic field, the atmosphere of the planet plays the most important role in forming the magnetospheric structure. This phenomenon is the most outstanding in the case of comets, which lack a magnetic field and are surrounded by a thick atmospheric layer. In the case of the planets of which the magnetic field is very weak, the ionosphere and the solar wind wake on the opposite side of the sun are formed by the ions generated by the photoionization of the outer plasma bound to the unstable magnetic field, such as the solar wind, or the initial neutral atmosphere by the solar radiation, or by the new ions generated through the ionization process such as charge exchange and electron impact. The ionosphere serves as a conductor so that the interplanetary magnetic field frozen-in by the solar wind can be modified to form the structures such as bow shock, ionospheric layer interface, and magnetic tail. The Martian bow shock formed in this way is affected by the mass-loading effect, which is different from the planets with a strong magnetic field such as earth (Breus et al. 1989). As the ions of planetary origin generated by the photoionization of the Martian neutral atmosphere by the solar radiation or by charge exchange or electron impact are loaded on the solar wind, the range of the bow shock becomes further from the planet by the mass-loading effect, that is, the increase of the overall solar wind mass and the decrease of the solar wind velocity based on the law of conservation of momentum (Kim 1999, Yi et al. 1999).

In this article, the structure of the Martian magnetosphere was analyzed by applying three-dimensional MHD simulation. Similar to comets, the MHD model including the mass-loading effect has been applied to Mars. To simulate the structure of the Martian magnetosphere, the three-dimensional MHD simulation code of Yi et al. (1996) that includes mass-loading and had been applied to comets, was modified and then applied to the interaction between the solar wind and Mars. Beside the steady solar wind condition, the structure of the Martian magnetosphere was simulated and compared with various cases assuming the initial conditions of the solar wind density and velocity. Based on the previous observation data, we completed the three-dimensional MHD model for the Martian upper atmosphere and magnetosphere and used it to provide the basis for the prediction of the observation results by future Martian probes.

Slavin et al. (1991) introduced the positions of the Martian bow shock found in the observation data from Mariner 4, Mars 2, 3, 5, and Phobos 2. The observed range of the bow shock was 1.3~2.5 R_M (the radius of Mars) in the subsolar direction, and the mean value was 1.56 ± 0.04 R_M. The range of the bow shock in the terminator direction was 2.0~3.5 R_M, and the mean value was 2.66 ± 0.49 RM.

The time-dependent observation data of Phobos 2 shown in Sauer et al. (1992) showed that the x-component of the solar wind velocity, the initial solar wind direction was decreased very much, while the magnitude of the magnetic field was increased after the solar wind passed through the bow shock. Additionally, it was found that, even though the observed number density of protons in the solar wind was reduced after passing through the bow shock, the ratio of O⁺ and H⁺ particles was increased by ten times or even more than one hundred times in the mass-loading region, as shown in the report by Lundin et al. (1990).

In this study, the interactive process between the solar wind and Mars was calculated using the ideal gas MHD equation. In the ideal compressible MHD equation, the mass, momentum, and energy conservation equations are expressed as Eqs. (1-3) (Yi et al. 1996):

Additionally, the Faraday’s law and Ampere’s law for induced current are expressed in Eq. (4):

where ρ denotes the plasma mass density, A the source of mass increase by mass-loading, P the plasma pressure, w the out-flowing velocity, B the solar wind magnetic field, γ (= 5/3) the specific heat ratio when adiabatic process of monoatomic ideal gas is assumed, and v the plasma flow velocity. In this article, the calculation is the single fluid MHD calculation in which the solar wind particle and the planet origin ions are not distinguished, and the source of mass, momentum, and energy is expressed as A.

The three-dimensional simulation code that had been applied to comets (Yi et al. 1996) was modified and applied to the interaction between Mars and the solar wind. Fig. 1 shows the three-dimensional MHD simulation box. The number of grids in the simulation box was (x, y, z) = (126, 62, 62), and the grid width was 120 km in the x direction, and 200 km in the y and z directions. The size of the simulation box was -2.1 R_M ~ 2.38 R_M in the x direc-

krtion, and -3.56 R_M ~ 3.56 R_M in the y and z directions. For the calculation, the volume of 1/4 of the entire box was simulated and the result was symmetrically reflected to the planes, z = 0 and y = 0.

Initially, the x direction of the simulation box was set to the initial solar wind direction, and the interplanetary magnetic field direction was set to the +y direction. The simulation was continued until the steady state where the whole structure of the magnetosphere is stabilized. The assumed initial conditions of the solar wind were the solar wind density ρ_sw of 5 m_H/cm³, the solar wind velocity v_sw of 350 km/sec, the interplanetary magnetic field strength Bsw of 5 nT, and the solar wind temperature T_sw of 105 K. The radius of Mars, R_M, was assumed to be 3,400 km, and the inner boundary which is the bottom of the exosphere was set as R_s = R_M + 300 km. The generation rate of the ions that are loaded on the solar wind by the mass-loading effect was set as A = σΣ_i n_i, considering only the photoionization of the Martian exospheric neutral atmosphere by the solar radiation. The photoionization rate σ was assumed to be 10^-6/sec. The simulation was performed varying the initial conditions, the solar wind density and velocity, to analyze how the structure of the Martian magnetosphere is affected by the dynamic pressure of the solar wind.

The ions that are the source of mass loaded on the solar wind are generated by photoionization of the Martian neutral atmosphere. The Martian neutral atmospheric model of Nagy et al. (1990) was applied to the three-dimensional simulation calculation. Fig. 2 shows the altitude distribution of the thermal hydrogen and oxygen as well as non-thermal hydrogen and oxygen.

Among these, the oxygen components whose mass is larger than that of hydrogen mainly serve as the source of mass in the mass-loading effect. Since all the thermal oxygen components are distributed at lower altitude, the ions that become the main source of mass density are formed by the non-thermal oxygen that is uniformly distributed over the entire altitude. The respective density and scale height of the neutral hydrogen and oxygen are shown in Table 1. The altitudinal distribution of the atmospheric components was calculated by Eq. (5):

where N_i (r) denotes the number density of each component, R_S the distance from the center of Mars to the height posiof exosphere R_s = R_M + 300 km, hs_i the scale height of each component, N_0i the initial neutral number density at the bottom of the exosphere.

Simulation was performed using the three-dimensional MHD model to view the structure of the Martian magnetosphere formed by the interaction between Mars and the solar wind. Under the solar wind conditions of the density ρ_sw of 5 m_H/cm³ and the velocity v_sw of 350 km/s, the time required to the steady state where there was no change in the plasma distribution in the subsolar and terminator directions was 108.0 seconds. Fig. 3 is the simulation result of the interaction between Mars and the solar wind in the simulation box filled with the solar wind of which initial mass density ρ_sw and velocity vsw were 5 m_H/cm³ and 350 km/s, respectively, and the plasma density distribution around Mars is shown in the interval of 12.0 sec from t₁ = 12.0 sec to t₉ = 108.0 sec, on the xy plane. The individual contour lines represent the density distribution from 15 m_H/cm³ to 150 m_H/cm³ or higher in the

interval of 15 m_H/cm³. Fig. 4 is the simulation result of the plasma pressure distribution in the range between 0.05 to 3.5 nPa on the xy plane under the same conditions and in the same time of Fig. 3. The plasma density and pressure were higher in the subsolar direction as the position was nearer to the planet, as shown in Figs. 3 and 4.

Fig. 5 shows the plasma density and pressure on the xy plane under the normal solar wind conditions of the density ρ_sw of 5 m_H/cm³ and the velocity v_sw of 350 km/s, as the result of the three-dimensional MHD model simulation of the interaction between the solar wind and Mars. The figure on the left shows the density distribution from 15 m_H/cm³ represented by white color to 120 m_H/cm³ represented by black color in the interval of 15 m_H/cm³. The figure on the right shows the contour lines for the pressure distribution from 0.1 nPa to 0.8 nPa in the 0.1 nPa interval. The region where the contour lines that were in the uniform interval started to be dense is the position of the shock wave. It was found that solar wind density around Mars was increased by three times from the initial solar wind density of 5 m_H/cm³ in the region before the shock wave to 15 m_H/cm³ after passing through the shock wave, and it was increased as the posi-

tion was nearer to the planet surface in the subsolar direction. The solar wind pressure was also increased drastically from the initial solar wind pressure, 0.069 nPa, to 0.1 nPa after passing through the shock wave, and it was increased as the position was nearer to the planet surface in the subsolar direction.

Fig. 6 shows the distribution of the plasma velocity and the magnitude of the magnetic field on the xy plane around Mars, as the result of the simulation under the same solar wind conditions in Fig. 5. The figure on the left shows the velocity distribution from 300 km/s to 20 km/s in 40 km/s interval. It shows that the plasma velocity was drastically decreased from the initial solar wind velocity 350 km/s to 100 km/s or less after it passed through the shock wave in the subsolar direction. The figure on the right shows the distribution of the magnetic field strength around Mars, indicating that it was increased from the initial solar wind magnetic field strength 5 nT to 24 nT or higher after it passed through the shock wave in the interval of 2 nT.

Fig. 7 shows the plasma velocity and magnetic field vectors around Mars on the xy plane under the same solar wind conditions in Fig. 5. It shows that the velocity vector of the inflowing solar wind in the subsolar direction (the +x direction) steered by the planet as the +x direction component was decreased and the ±y direction component was increased around the leading edge of the shock wave. The magnitude of the solar wind plasma velocity that seemed to stop in the subsolar direction was increased again as the ±y direction velocity component was increased as the solar wind moved to the terminator direction of the ±y direction. The figure on the right side shows that the strength of the interplanetary magnetic field by the solar wind which was initially in the +y direction was increased in the inner side near to the shock wave and the ±x component was increased except in the subsolar direction so that the solar wind was steered by the planet.

Fig. 8 shows the plasma velocity and magnetic field vectors around Mars on the xz plane under the same solar wind conditions in Fig. 5. The plasma velocity on the left side shows the same trend found in the xy plane (the left figure on Fig. 7). The magnetic field shows that the ±x and ±z direction components of the interplanetary magnetic field that had only the +y component were increased in the inner side near to the shock wave.

Fig. 9 schematically illustrates the structure in which the supersonic flow of the fluid became subsonic in the part perpendicular to the shock front that was formed as the fluid flow meets a round-shaped obstacle in the

front but the velocity is accelerated to supersonic velocity again as the fluid approached the terminator region. Fig. 10 is the simulation result of the interaction between Mars and the solar wind under the solar wind conditions that are 5 m_H/cm³ of density and 350 km/s of velocity, showing the Mach number distribution around Mars on the xy plane. The Mach number, M, is defined as Eq. (6):

where v, ρ and P denote the velocity, density and pressure of the plasma, respectively. As the shock wave strength is stronger, the Mach number in the shock wave position is increased in the upstream but decreased in the downstream. In this article, the change of the shock wave strength was shown using the Mach number in the upstream. Fig. 10 well shows that the plasma traveled between the shock wave and the planet in the subsolar direction, as in Fig. 9, in a subsonic velocity of which Mach number was less than 1, while the Mach number was increased again to the supersonic velocity as it approached to the terminator.

Fig. 11 shows the change of the bow shock position as the distance was varied from -2.1 R_M to -1.1 R_M in the subsolar direction and the solar wind density was varied to 2, 5, 10, 15 and 20 m_H/cm³. The figure shows the change of the density (a), velocity (b), pressure (c), and magnetic field (d), depending on the distance expressed by the radius of Mars. As the solar wind density was higher, the position of the shock front was closer to the center of Mars. However, the change of the shock wave position was relatively small when the solar wind density was high as 10, 15 and 20 m_H/cm³ than when it was low as 2 and 5 mH/cm3. The figure also shows that the plasma velocity (b) was reduced in the position where the plasma density (a), pressure (c) and magnetic field strength (d) were increased.

Fig. 12 shows the Mach number of the plasma in the subsolar direction for each case of Fig. 11. Since it is subsonic in the inner side of the shock wave in the subsolar direction, as explained in Fig. 9, the position of the bow shock in each case of Fig. 12 is the place where the Mach number was declined to less than 1. Table 2 shows the position of the bow shock in each case when the solar wind velocity was constant as 350 km/s and the density was varied to 2, 5, 10, 15 and 20 m_H/cm³. Here, R_M, which is 3,400 km, denotes the radius of Mars. In each case, the position of the shock wave fell within the observed Martian bow shock range, 1.3~2.5 R_M (Slavin et al. 1991).

Fig. 13 shows the distribution of the density (a), veloc-

ity (b), pressure (c), and magnetic field (d) of the plasma around Mars in the subsolar direction for each case when the solar wind density was 5 m_H/cm³, and the solar wind velocity was 300, 350, 450, 500, 600 and 700 km/s. As the solar wind velocity was higher, the position of the shock front was closer to the center of Mars. However, the change of the shock wave position was relatively small when the solar wind velocity was high as 450, 500, 600, 700 km/s than when it was low as 300 and 350 km/s. The figure also shows that the plasma velocity (b) was reduced in the position where the plasma density (a), pressure (c) and magnetic field strength (d) were increased.

Fig. 14 shows the Mach number of the Martian plasma in each case when the solar wind density was 5 m_H/cm³, and the solar wind velocity was 300, 350, 450, 500, 600 and 700 km/s. Table 3 shows the position of the bow shock in each case depending on the solar wind velocity change. In each case, the position of the bow shock fell within the range of the observed shock wave positions. In comparison to Fig. 12, the change of the bow shock position by the solar wind velocity variation was relatively smaller than that by the density variation.

Fig. 15 shows the dynamic pressure of the solar wind calculated by Eq. (7) and position of the Martian bow shock at the corresponding dynamic pressure when so-

lar wind density and velocity were varied for each case. The dynamic pressure and the bow shock position were marked as ‘triangles’ in the cases where the solar wind density was varied to 2, 5, 10, 15 and 20 m_H/cm³ under the assumption that the solar wind velocity is constant as 350 km/s, while they were marked as ‘x’ in the cases where the solar wind velocity was varied to 300, 350, 450, 500, 600 and 700 km/s under the assumption that the solar wind density is constant as 5 m_H/cm³. The result that the gradient of the bow shock position change by the solar wind density variation was larger than that by the solar wind velocity variation even when the dynamic pressure of the solar wind was the same indicates that the change of the solar wind density has a greater effect on the change of the bow shock position.

If it is taken into consideration that the momentum, ρv, is increased as the density, ρ, is increased at the same dynamic pressure, ρv², the result that the shock front position is pushed into the center of Mars as the solar wind density is increased can be explained.

Different from the earth of which own magnetic field is strong, the Martian magnetic field is formed by the direct interaction between the solar wind and the Martian atmosphere. The most important role in the formation of the Martian bow shock is played by the mass-loading effect by which the planet origin ions are loaded on the solar wind flow so that the overall solar wind density is increased and the solar wind velocity is decreased and thus the position of the bow shock becomes farther from the planet. We simulated the Martian magnetic field depending on the changes of the solar wind density and velocity using the three-dimensional MHD model prepared by modifying the comet code (Yi et al. 1996) that includes such mass-loading effect. The Martian exosphere model of Nagy et al. (1990) was employed as the Martian atmosphere model, and the photoionization by the solar radiation was assumed to be the only ionization process of the neutral atmosphere.

In the three-dimensional simulation result under the normal solar wind condition, the position of the Martian bow shock in the subsolar direction was well consistent with the result of Liu et al. (1999). In the result of the three-dimensional simulation performed by varying the solar wind density and velocity, all the Martian bow shock positions fell within the range of the observed positions in the subsolar direction, 1.3~2.5 R_M (Slavin et al. 1991). The simulation result also showed that the change of the bow shock position by the solar wind density variation was greater than that by the velocity variation.

In order to accurately predict the Martian bow shock position based on our results, not only the photoionization but also other ionization processes should be included in the mass-loading procedure, and the model needs to be improved to describe in detail the conditions that were simplified in this study as the conditions that are more similar to the actual conditions, including the consideration of the asymmetric neutral atmospheric model caused by the rotation of Mars. We expect that this article, as the initial result of the model, to be useful in analyzing the data observed by the future Martian probes and understanding the interaction of the solar wind and the Martian atmosphere.