Three-Dimensional MHD Modeling of Interplanetary Solar Wind Using Self-Consistent Boundary Condition Obtained from Multiple Observations and Machine Learning

: Three-dimensional (3-d) magnetohydrodynamics (MHD) modeling is a key method for studying the interplanetary solar wind. In this paper, we introduce a new 3-d MHD solar wind model driven by the self-consistent boundary condition obtained from multiple observations and the Artiﬁcial Neural Network (ANN) machine learning technique. At the inner boundary, the magnetic ﬁeld is derived using the magnetogram and potential ﬁeld source surface extrapolation; the electron density is derived from the polarized brightness ( pB ) observations, the velocity can be deduced by an ANN using both the magnetogram and pB observations, and the temperature is derived from the magnetic ﬁeld and electron density by a self-consistent method. Then, the 3-d interplanetary solar wind from CR2057 to CR2062 is modeled by the new model with the self-consistent boundary conditions. The modeling results present various observational characteristics at different latitudes, and are in better agreement with both the OMNI and Ulysses observations compared to our previous MHD model based only on photospheric magnetic ﬁeld observations. F.S.; Y.Y. F.S.; Y.Y. F.S.;


Introduction
Magnetohydrodynamic (MHD) modeling is an important method for space plasma physics research and space weather prediction. In the research field of heliophysics, MHD models have been widely used in the studies of coronal-interplanetary solar wind, the propagation of coronal mass ejections (CMEs), and their interaction process [1]. Obtaining the realistic three-dimensional (3-d) solar wind plasma and magnetic field environment is a significant foundation for studying and predicting the propagation process of catastrophic space weather events such as CMEs and solar energetic particles (SEPs) [2][3][4][5][6][7][8][9]. Therefore, the development of a realistic 3-d MHD model for the interplanetary solar wind has great theoretical and practical value. Different research groups have developed a few reasonable 3-d MHD solar wind modeling systems with relatively good performance at present, for example, the space weather modeling framework (SWMF) developed by the Center for Space Environment Modeling (CSEM) at University of Michigan [10], the Coronal-Heliosphere (CORHEL) model developed by the Center for Integrated Space Weather Modeling (CISM) at Boston University [11,12], the WSA-ENLIL model utilized by the space weather prediction center at National Oceanic and Atmospheric Administration (NOAA) [13], the SIP-CESE model and the COIN-TVD model developed by the National Space Science Center at Chinese Academy of Sciences [14][15][16], the EUHFORIA model developed by the European research community [17], and some models developed by other researchers [18][19][20][21].
Considering solar wind propagates to the interplanetary space through a complex heating and accelerating process at the inner corona, the current 3-d MHD solar wind models usually divide the solar-interplanetary space into two computational regions and deal with them by different methods [22]. One kind of model computes both the coronal and interplanetary regions by the physics-based MHD algorithm, and deals with the heating and accelerating process in the corona by adding mathematical terms to the MHD equations, such as the SWMF model and the SIP-CESE model. The other kinds of model deals with the coronal region by semi-empirical models, and only compute the interplanetary region by the MHD algorism. The second kind of model generally sets the inner boundary of MHD computational region beyond the Alfven critical point (around 0.1 AU from the Sun), and determines the inner boundary conditions with semi-empirical methods, such as the WSA-ENLIL model and the EUFORIA model.
The plasma β is the ratio of the plasma pressure to the magnetic pressure. Due to the plasma β << 1 in the corona, the numerical computation of the first kind of solar wind models requires a large amount of computing resources, while the second kind of model can save the computational cost of the corona and improve the computational efficiency for interplanetary solar wind greatly. Since the heating and accelerating mechanisms of the solar wind have not been totally revealed, the models based on empirical boundary conditions are no less effective than the models fully based on physical equations in reproducing large-scale structures of interplanetary solar wind, as a study that compared the two kinds of models in detail presented [23]. Therefore, the interplanetary solar wind model based on empirical boundary conditions has important value in solar wind prediction in terms of both efficiency and accuracy. For example, the WSA-ENLIL model based on the WSA empirical function for solar wind velocity has been utilized as an operational model by the SWPC/NOAA.
The performance of the hybrid empirical-MHD solar wind model is highly dependent on the initial boundary conditions used to solve the MHD equations, thus the realistic initial boundary conditions are significant to the simulation accuracy of the interplanetary solar wind [24][25][26][27][28][29]. Since the magnetic field is the main factor controlling the motion of solar wind plasma, most advanced solar wind MHD models utilize observational data of the photospheric magnetic field and some theoretical extrapolation models to obtain the initial boundary conditions of the magnetic field. For example, the potential field source surface (PFSS) model can extrapolate the coronal magnetic field from the photospheric magnetograms, and is widely used in current solar wind models. The boundary condition of velocity is usually given by the WSA model [30][31][32] based on an empirical relation between solar wind velocity and magnetic field expansion factor. In addition, the initial boundary conditions of density and temperature are generally derived by the assumption of constant momentum flux and the total pressure (sum of thermal and magnetic pressures) balance [19,28,33,34]. Thus, these two parameters are closely related to the derived magnetic field and velocity.
However, the coronal magnetic field extrapolated through theoretical models with some ideal assumptions, such as the potential field assumption, may be biased from the reality. Therefore, some additional types of observational data can be considered to improve the boundary conditions of solar wind MHD models. For example, some recent studies have tried to correct the coronal magnetic field extrapolated from the PFSS model by the white light brightness observation, according to the magnetic freezing effect of the solar wind plasma [35,36]. In addition, the empirical functions for solar wind velocity in the WSA model are derived only from ecliptic solar wind observations, so it may have uncertainties for providing the global distribution. Since the global structure of solar wind velocity can be obtained from the interplanetary scintillation (IPS) observations, some researchers have utilized the IPS data to give the boundary conditions of solar wind MHD models [25,[37][38][39]. However, the IPS observation stations cannot operate every winter due to thick snow, and the data also could not cover all the global regions for other seasons due to the lack of proper radio sources [40].
With the development of computational technology and the accumulation of observational data, machine learning techniques have been introduced to predict the solar wind condition in recent years. For example, different machine learning algorisms have been utilized to predict the solar wind speed near the Earth [41,42]. However, these models can only give the prediction of solar wind speed near the Earth, and cannot give the 3-d distribution of all the solar wind parameters. In our recent work [43], we have developed a new method using multiple observations and the artificial neural network (ANN) machine learning technique to obtain the self-consistent global distribution of the magnetic field, electron density, flow velocity and plasma temperature on the coronal source surface. We established an empirical model for solar wind velocity based on the magnetic field properties and the electron density derived from the observations of photospheric magnetogram and polarized brightness. The new empirical correlation is determined from global observations at all latitudes as it utilized the IPS solar wind velocity data. The plasma temperature distribution is deduced by solving a 1D MHD system at the source surface. Therefore, the global distribution of different solar wind parameters is deduced self-consistently based on multiple observational data.
In this paper, we introduce a new 3-d MHD interplanetary solar wind model that coupled the method of obtaining self-consistent global distribution of solar wind parameters at the coronal source surface [43] and the interplanetary total variation diminishing (IN-TVD) MHD solar wind model [29]. By analyzing the modeled global structure of solar wind plasma and magnetic field, the new model is validated. Then, the modeling results are compared with in situ observations both at the ecliptic plane and at higher latitudes along the Ulysses spacecraft trajectory. In addition, the comparison with our previous MHD model shows the advantages of the multiple observation driven model.
The outline of this paper is as follows: in Section 2, we describe the data and methods that we used to establish the new MHD model. In Section 3, we first make analyses for the modeled global structure of different solar wind parameters during one Carrington rotation (CR), then we validate the model by comparing the modeling results with the Ulysses observations for six CRs from CR 2057 to CR 2062, i.e., the period from May 24th to Nov 4th of 2007. In Section 4, we discuss the results and draw some conclusions.

The MHD Numerical Model
Our interplanetary solar wind model is based on the 3-d ideal MHD equations. Because the solar wind and interplanetary magnetic field are rigidly co-rotating with the Sun, a co-rotating coordinate system can facilitate the processing of the bottom boundary conditions. In the co-rotating coordinates, the ideal MHD equations in the conservation form can be written as: ∂ρV ∂t where the physical parameters to be solved are the plasma density ρ, pressure p, velocity vector V and magnetic field vector B. Equation (1) is the continuity equation. Equation (2) is the momentum equation, where µ 0 is the vacuum permeability, I denotes the unit tensor, G is the universal gravitational constant, M s is the mass of the Sun, ω is the angular velocity of the solar rotation, and f = −ω × [ω × r + 2ω × V] is the centrifugal force introduced because of the adoption of the co-rotating coordinate system. Equation (3) is the simplified electromagnetic equation considering the ideal MHD condition. Equation (4) is the pressure equation used instead of the energy equation, where γ is the polytrophic index. The use of pressure equation is beneficial to ensure the stability of the numerical scheme [13].
In our model, the Lax-Friedrich finite difference scheme is used with the total variation diminishing (TVD) feature to solve the above 3-d time-dependent MHD equations in all the computational regions. In addition, since the numerical calculation cannot guarantee the physical condition that the divergence of the magnetic field is zero everywhere, so the divergence of the magnetic field may accumulate at some places and lead to the collapse of the calculation, and some methods have to be adopted to limit the divergence of the magnetic field. In this paper, we have adopted the diffusion method [44,45], which added a diffusion term to the equation to make the magnetic field divergence generated by numerical calculation spread rapidly to the whole calculation domain, so that the stability of the calculation can be ensured. Other details of the numerical scheme can be referred to the previous papers [15,46].
Our calculation domain is set to be from 0.1 AU to about 2 AU from the Sun's center, exactly [21.5R s , 450R s ] in units of solar radius (R s ). The meshes in the radial direction are divided by equal ratio, so the meshes become larger with the increase in heliocentric distance. The radial mesh size increases gradually from 0.37R s to 3.6R s . The allowable time step will be larger as the grid is larger, according to the Courant-Friedrichs-Lewy (CFL) conditions of numerical calculation. Therefore, we divide the whole computing space into six regions in the radial direction, and each region uses its own allowable time step, so that asynchronous parallel method can be adopted in the calculation to save the calculation time [15]. The division of radial grids is shown in Figure 1a. Different colors represent six different regions divided radially. There are 48 grid point in each region, and 4 of them overlap with the nearby region. The number of grids in the radial direction can be adjusted as needed to obtain the appropriate calculation area. On the spherical shell of our computation domain, a six-component grid system is adopted to avoid the singularity in the spherical coordinate system and facilitate parallel computing with high performance. The six-component grid system consists of six identical regions, four of which are arranged in order along the longitude in the middle and low latitudes, and the other two are used to cover the North and South poles. Each component intersects with the nearby ones to cover an entire sphere. Figure 1b shows a diagram of the six-component grid systems in 3-d with each component represented by a different color.
In the six-component grid system, each component has an independent coordinate, and the coordinate range of each region is as follows: where ∆ represents the size of the overlap between the components on a spherical shell. Each component is divided according to the following rules: where n is the number of grids in the direction of θ or φ. In this work, we adopted a mesh with the resolution of 1 • × 1 • on a sphere, and the overlapping parts of different components are set to 7 grids, so the value of n is 97. During the calculation process, the values of physical parameters for the inner grids of each component are solved by the numerical scheme, and the values on the boundary of each components are obtained by interpolation of the adjacent overlapping grids.

Boundary Conditions Obtained from Multiple Observations and Machine Learning
As mentioned above, the inner boundary of the calculation domain in this paper is located at 21.5R s , which belongs to the supersonic and super-Alfven-velocity region. Thus, all the characteristic waves on the inner boundary are pointed to the calculation region, so all the physical values in the MHD equations need to be set empirically. The outer boundary adopts no reflection condition, which can be obtained by linear extrapolation of the points in the calculation region.
If realistic boundary conditions are given to the MHD model, the simulated interplanetary solar wind can be well consistent with actual observations. In order to provide a more realistic inner boundary, here we try to obtain the global distribution of various physical parameters using multiple observational data. Since there is no directly available interplanetary observation data at the inner boundary location, we firstly determine the global distribution of the physical parameters at the coronal source surface (2.5R s ) using the observational data, and then extrapolate the global distributions to the inner boundary of the MHD model according to the Parker spiral [47]. The following describes how the boundary conditions for magnetic field, density, velocity, and temperature are given based on multiple observations. Our boundary condition of radial magnetic field is given by the PFSS model [48,49], which can extrapolate the coronal magnetic field structure using the observational data of photospheric magnetograms. Although the PFSS model is relatively simple, it can still predict the polarity of the interplanetary magnetic field well [50], so it is commonly used in a lot of solar wind models. Currently, there are many data sources for the photospheric magnetograms data. According to a comparative study [51,52], we used the data provided by NSO/GONG (http://gong.nso.edu/, accessed on 30 September 2021). Previous studies have demonstrated that the nondipole harmonics of the solar magnetic field can be significant [53]. In this paper, we use a truncation number of 11 for the spherical harmonic expansion in the PFSS model. After the calculation of the magnetic field intensity, we can also obtain some other magnetic characteristic parameters by tracing field lines. The expansion factor f s on the source surface can be derived by the function: where B 0 and B ss are the magnetic field strength at the photosphere and at the source surface, respectively. The parameter θ b is defined as the minimum angular distance between an open field foot point and its nearest coronal hole boundary. The f s and θ b are good indicators for solar wind speed as they are both used in the WSA model. In addition, we define L m as the latitude difference between a point at the source surface and the heliospheric current sheet (HCS), i.e., the magnetic latitude. The L m is also well related to the solar wind structure, as has been studied by many researchers [54][55][56]. Another feature we derived is the latitude of the magnetic foot-point (defined as L 0 ), since it may also have a relation with the solar wind velocity [57]. The six magnetic field features named B ss , B 0 , f s , θ b , L m and L 0 which are obtained from the PFSS model were utilized to predict the solar wind speed at 1AU in our previous paper, and achieved a good performance [58]. Thus, we intend to use these features to derive the global distribution of solar wind velocity in this study.
During the propagation of solar wind in interplanetary space, the dynamic pressure ( 1 2 ρV 2 ) dominates the plasma flow in the region with a weak magnetic field. Thus, it is very necessary to give a reasonable density boundary condition for the 3-d MHD solar wind model. According to the Thomson scattering theory, the coronal electron density can be determined from the white light polarization brightness (pB) observation. Therefore, we obtain the global distribution of electron density at the source surface using the pB data from LASCO/C2 coronagraph and a coronal electron density inversion methods developed by Van de Hulst et al. [59]. Although the LASCO/C3 coronagraph can provide observations at the inner boundary location of our model directly, it cannot be used to invert the electron density distribution, because only the pB observation at low corona can avoid the noise caused by the interplanetary stray light [60]. The boundary condition of solar wind radial velocity is very important in solar wind MHD modeling, because the interaction between high and low velocity flows during the propagation of interplanetary solar wind has a great influence on the evolution of all other solar wind parameters. The current MHD solar wind model generally gives the radial velocity boundary value by the empirical WSA model, which can be expressed in the following formula: where f s is the expansion factor, θ b is the angular distance between the foot point of the magnetic field line and the coronal hole boundary, V s and V f are free parameters used to define the minimum and maximum velocities, and a 1 to a 4 is the free parameters that control the influence of f s and θ b . There are variants of similar formula used in different studies [23,30,61], and we chose the above function because it has reduced free parameters and was validated in our previous study [29]. As mentioned above, the empirical formula of WSA velocity is obtained from the in situ observational data at 1AU in the ecliptic, so the initial boundary conditions of the high latitude region given by the WSA formula may have uncertainties. Moreover, the WSA model only relies on the photospheric magnetic observation. Therefore, it may be better if some other types of observations can be added to constrain the value of the velocity. For these reasons, we have used an artificial neural network (ANN) machine learning technique to establish an empirical relationship among the source surface velocity, six magnetic field parameters based on photospheric magnetic observations and the electron density inverted from pB observations. The new empirical relationship was established using IPS observations of global velocity, so it could be more reasonable than the WSA formula at high latitudes. The ANN model we used is a feed forward network with three layers, including the input layer, the hidden layer and the output layer. The architecture of our ANN model is shown in Figure 2. The calculation of the ANN can be expressed as follows: where x i is the input features; w ji is the weights between the neuron i of input layer and node j of hidden layer, v j are the weights between the hidden layer and the output; b j and b 0 represent the biases; and y is the target output. Besides, f is a bipolar sigmoid nonlinear activation function used in the hidden layer, and g is a linear function applied for the output layer. The function of f is: As shown in Figure 2, we use seven features derived from multiple observations as input and 150 neurons in the hidden layer. Using the above method, we can obtain a velocity boundary condition that are related to both photospheric magnetic observations and coronal pB observations. Figure 3 shows the comparison of (a) the IPS observed solar wind velocity distribution, (b) the ANN model derived velocity distribution, and (c) the WSA model derived velocity distribution during CR2057 to CR2062. We use this period to validate our model because there are few CMEs during this period, so it is proper to examine the model performance for background solar wind.
Our temperature boundary condition is derived by a self-consistent approach based on the already obtained magnetic field intensity and density distributions. It assumes that solar wind plasma satisfies the ideal gas conditions and the statistical results for the constant thermal pressure and plasma beta β at different regions, which can be presented as the following: where T is the temperature, N is the density, p is the dynamical pressure, B is the magnetic field intensity, and the value of plasma β is given by: and the value of β c and C are given empirically according to the previous study [62]. Since the magnetic field and density are derived from photospheric magnetograms and pB brightness observations, respectively, the temperature computed from them is also related to the two kinds of observations self-consistently. In addition, the longitudinal and latitudinal components of the magnetic field vector and velocity vector are given according to the corotation as follows: Above all, we have obtained the boundary conditions of all the solar wind physical parameters according to three kinds of observational data, including the photospheric magnetogram, the pB brightness and the IPS velocity data. Then, we can drive the 3-d MHD model introduced in Section 2.1 and achieve the steady-state solution of the 3-d interplanetary solar wind after sufficient time relaxation iteration.
The left column of Figure 4 shows the global distributions of (a) radial velocity, (b) electron density, (c) temperature and (d) radial magnetic field at the inner boundary of our 3-d MHD model obtained from multiple observations and the ANN technique during CR2062. The right column is a comparison with the inner boundary condition based solely on the the WSA model, which we previously used. As the figure showed, the heliospheric current sheet during CR2062 is mainly distributed in low latitude near the equator and is well presented by both methods. However, the global distributions of density and temperature presented in the right column have exactly the same distribution structure with the velocity, while the global distribution of velocity, density and temperature in the left column have many more longitudinal and latitudinal variations, especially for the regions at the high latitudes.

Results
Utilizing the new 3-d MHD interplanetary model driven by multiple observations and ANN, we simulated the 3D interplanetary background solar wind during six CRs, which are CR2057 to CR2062. The Ulysses spacecraft had experienced its third fast passage from the solar south to the north, so we could examine the modeling performance at high latitude during this period. Besides, this period was in the year 2007, which is close to the minimum of the 23rd solar cycle, so there were few coronal mass ejections, meaning this time periods would be very suitable for the verification of a background solar wind model.

Analyses of the Modeled Solar Wind Parameters at 1AU
Firstly, we analyzed the modeling results for the global distribution of different solar wind parameters at 1AU during CR2062 as an example.
The left column of Figure 5 shows the global distribution of the interplanetary solar wind parameters at 1AU during CR2062 simulated by the MHD model with our new boundary condition settings. From a comprehensive view of various physical parameters, we can see that in the polar region the velocity is higher, the density is lower, the temperature is higher and the magnetic field is also higher than that in the equatorial region. These characteristics are consistent with the actual observation of large-scale solar wind structures near the minimum activity years in a solar cycle.
By comparison with Figure 4, it can be found that all the solar wind parameters at 1AU had obvious changes compared with the inner boundary. The magnetic field and density at low latitudes showed some inhomogeneous structures due to the interaction between high and low velocity flows, and a few sudden increases emerged due to the compression effects. Besides, the average solar wind speed increased by tens of km/s while the magnetic field and temperature decreased exponentially relative to the inner boundary. This means that when the solar wind travels through interplanetary space, the magnetic energy and internal energy are converted into kinetic energy. In addition, we can intuitively see that the large-scale configuration of all the solar wind parameters at 1AU is shifted to the left about 50 degree in the longitudinal direction relative to the inner boundary, which is due to the co-rotation effect of the solar wind [47].
The right column of Figure 5 shows the global distribution of solar wind parameters at 1AU during CR2062 simulated by the MHD model with our previous boundary condition treatment [29]. The comparison between the left and right columns in Figure 5 shows the modeling results from both models are similar on the large scale, but the new model driven by multiple observations presents more small variations for all the solar wind parameters. In the right column, the values of solar wind plasma parameters along the heliospheric current sheet are equal, while in the left column the values can be different. For example, the densities along the heliospheric current sheet are all above 10 cm −3 according to the results of our previous model, while the new model indicates the density can be less than 10 cm −3 . Besides, there are more complex distributions of the parameters at high latitude in the left column. For example, the velocity value is almost the same above 60 degrees according to the result of our previous model, while the new model shows some variations.
In order to examine if the small scale variations simulated by our new MHD model are closer to the actual solar wind environment than the previous model, we compared the modeling results with the in situ observations at 1AU provided by OMNIweb (http: //omniweb.gsfc.nasa.gov, accessed on 30 September 2021), and the results are shown in Figure 6. and Table 1.   Figure 6 shows the comparison of our MHD simulation results of solar wind (a) velocity, (b) density, (c) temperature and (d) radial interplanetary magnetic field (IMF) at 1AU with OMNI hourly averaged in situ observational data. The blue lines in Figure 6 show the results from the new model driven by multiple observations and ANN (MHD1), and the cyan lines present the results from our previous MHD model driven by the WSA-based boundary condition (MHD2). The overall trend of all simulated solar wind parameters at 1AU is generally in agreement with the observed value for both the MHD models. The two high-speed flows larger than 400 km/s are relatively accurate. The simulated values of density are generally consistent with the observation, but some peak values are missed. The two periods of temperature simulation values greater than 10 5 K are generally in agreement with the OMNI observations. The radial magnetic field also agrees with the observations at large scale, but the small fluctuations cannot be simulated by both the MHD models since they may come from the turbulence of the solar wind. We can see intuitively that the lines of the new model are closer to the observations.  Table 1 presents the correlation coefficient (CC) and the root mean squared error (RMSE) between the MHD models and the OMNI observation. We can see that both the MHD models performed well, especially for the solar wind velocity. Additionally, all the parameters simulated by the new model were more accurate than the previous model quantitatively for both CC and RMSE.

Comparison with the Ulysses Observation
Comparison between the modeling results and the OMNI in situ observational data at 1AU can only verify the performance of our MHD model on the equatorial plane. In order to further verify the performance of the model at high latitudes, we also compared the results with the observation data of Ulysses. The trajectory of Ulysses spacecraft during CR2057 to CR2062 is shown in Figure 7 as the black curve. The background is the modeled 3-d distribution of velocity during CR2062. We can see that the Ulysses passed a large latitude range from the Southern Hemisphere to the Northern Hemisphere during CR2057 to CR2062, which is very suitable for comparison with our simulation results at high latitudes. The cone-shaped zones with a reduced speed that can be seen in the polar regions may reflect the presence of conic polar current sheets, which are introduced by a previous research [63]. Figure 8 presents the solar wind parameters simulated by the MHD models and the six-hour averaged in situ observational data from Ulysses. The x axis at the bottom of each panel is the date of the year 2007, and the x axis at the top of each panel is the latitude of the Ulysses spacecraft. The green dashed lines are used to separate the six CRs. We can see that Ulysses traveled from about -50 degrees to above 55 degrees during this period.
Firstly, we examined the general trend over the whole period, and found that both the MHD models indicate that the solar wind velocity and temperature increase with the increase in latitude, while the density decreases with the increase in latitude, which is generally in agreement with the observed characteristics.
Then, if we focus on the latitudes above ±15 degree, the modeling results of our new MHD model coincide with the Ulysses observations better than the previous MHD model for all the solar wind parameters, intuitively. The most obvious improvement is in the temperature at high latitudes, since both the observation and the new MHD model presented the temperature as around 2 ×10 5 below −30 latitude degrees and above 35 degrees, but the previous model presented much lower values. However, the two large dips in speed and the corresponding impulses of density observed by Ulysses during late September to the beginning of October in Figure 8a,b were not well reproduced by both our MHD models. According to the theoretical studies on the quasi-stationary current sheets [64,65], these differences might be caused by the large-scale current sheets at midheliolatitudes that are not reproduced by our models, as we can see the inversions of the radial magnetic field direction are only presented by the observational data in Figure 8d.
In order to validate the modeling results more quantitatively, we also calculated the CC and RMSE between the MHD models and the OMNI observation. Overall, both the MHD models perform well, but the solar wind parameters simulated by the new model are more accurate than the previous model.
From the comparison of Tables 1 and 2, we can find that the metrics of the new MHD model using multiple observations and ANN technique stay nearly the same as at 1AU on the ecliptic plane as at high latitudes, indicating that the model can properly reproduce the solar wind plasma and magnetic field structures globally.

Discussion and Conclusions
In this paper, we introduced a new 3-d MHD interplanetary solar wind model using boundary conditions obtained from multiple observations and ANN technique. The model utilizes the TVD Lax-Friedrich scheme to solve ideal MHD equations in a parallelhelpful grid system with six components on a spherical shell and six radial partition. In the new boundary conditions, the magnetic field is extrapolated from the PFSS model using photosphere magnetic map observations, density was inverted from pB observation, velocity is obtained by an ANN, and temperature follows a self-consistent method that can be derived from the magnetic field and density. With the help of the ANN machine learning technology, the new velocity empirical relationship established from IPS data should be more reliable at high latitudes.
The new model is used to simulate the 3-d interplanetary solar wind during CR2057 to CR2062. Various analyses of the simulation results show that the new model using boundary conditions derived from multiple observations and ANN can reproduce the large-scale observational characteristics of the 3-d interplanetary solar wind better than our previous MHD model, especially at high latitudes. The modeled parameters are generally in agreement with the in situ observation data at 1AU and at the Ulysses orbit.
However, although the IPS data can provide the global velocity distribution including high latitudes, the error of IPS velocity is a problem we have not considered in this work. Besides, the ANN velocity function in this paper is trained from limited IPS data during the year 2007, so it may have larger errors during other periods. Nevertheless, utilizing more IPS data from different solar cycle phases to train the ANN in the future may improve the long-term performance of the ANN velocity function over the whole solar cycle. Since the Solar Orbiter has launched and will provide observations of solar wind plasma and the magnetic field at high latitudes in the coming years, we may extend this work to use more reliable data.
In addition, the numerical scheme can be improved to further improve the computational efficiency, so that the computational region can start directly from the near corona to the Earth orbit, which is more conducive to the simulation research and prediction of CME events. Data Availability Statement: The data used for this work are available at the official websites of NSO/GONG, SOHO/LASCO, NASA/GSFC and IPS group at Nagoya University. The Global Oscillation Network Group (GONG) magnetograms were downloaded from the website of the National Solar Observatory (NSO, http://gong.nso.edu/, accessed on 30 September 2021). The SOHO/LASCO data used here were produced by a consortium of the Naval Research Laboratory (USA), Max-Planck-Institut fuer Aeronomie (Germany), La-boratoire d'Astronomie (France), and the University of Birmingham (UK). The polarized bright-ness (pB) data can be downloaded from the LASCO instrument website at NRL (http://lasco-www.nrl.navy.mil/content/retrieve/polarize/, accessed on 30 September 2021). The IPS data were provided by the ISEE solar wind group at Nagoya University, and can be downloaded from the official website (http://stsw1.stelab.nagoya-u.ac.jp/ ips_data-e.html, accessed on 30 September 2021). We acknowledge the use of the OMNIWeb service provided by the NASA/GSFC Space Physics Data Facility. The Ulysses and near-Earth in situ data in this work were both obtained from the web interface at http://omniweb.gsfc.nasa.gov (accessed on 30 September 2021).

Conflicts of Interest:
The authors declare no conflict of interest.