Numerical modeling of solar wind and coronal mass ejection in the inner heliosphere: A review

The predictions of plasma parameters in the interplanetary medium are the core of space weather forecasts, and the magnetohydrodynamics (MHD) numerical simulation is an important tool in the prediction of plasma parameters. Operational space weather forecasts are commonly produced by a heliosphere model whose inner boundary is set at 18 R s or beyond. Such predictions typically use empirical/physics-based inner boundary conditions to solve the MHD equations for numerical simulation. In recent years, significant progress has been made in the numerical modeling of the inner heliosphere. In this paper, the numerical modeling of solar wind and coronal mass ejection in the inner heliosphere is reviewed. In particular, different inner boundary conditions used in the simulation are investigated since the MHD solutions are predetermined by the treatment of the inner boundary conditions to a large extent. Discussion is made on further development of the heliosphere model.


Introduction
Over the past few decades, the effects of space weather on technology have become an important field of research. The most serious space weather impact in near-Earth space is usually related to interplanetary coronal mass ejections (ICMEs), which are affected by the surrounding solar wind when propagated to Earth orbit (Wu et al., 2006;Kilpua et al., 2019;Li et al., 2020). For this reason, solar wind condition is a necessary prerequisite for the propagation of CMEs, which governs the propagation of ICMEs in the interplanetary space. The prediction of solar wind and ICME plasma properties plays a crucial role in the space weather forecast. However, the in situ observations of plasma parameters are only applicable to several points where spacecraft are located, for example, the Advanced Composition Explorer (ACE) (Stone et al., 1998), Wind (Bochsler et al., 1995), Parker Solar Probe (PSP) (Fox et al., 2016), Solar Orbiter (SO) (Müller et al., 2020), and Solar Terrestrial Relations Observatory (STEREO) spacecraft (Howard et al., 2008;Kaiser et al., 2008). Therefore, we have to rely on numerical technology to predict the plasma parameters.
Currently, a wide variety of techniques have been developed to generate plasma parameters in the interplanetary space. In particular, magnetohydrodynamics (MHD) simulations are important tools in this endeavor (Pizzo, 1982;Usmanov, 1993;Odstrcil, 2003;Detman et al., 2006;Nakamizo et al., 2009;Hayashi, 2012;Riley et al., 2012;Feng et al., 2014;Zhang and Feng, 2015;Shiota and Kataoka, 2016;Zhang and Xueshang, 2016;Pomoell and Poedts, 2018;Shen et al., 2018;Scolini et al., 2019;Feng, 2020;Scolini et al., 2020;Singh et al., 2020). The MHD simulations can be classified into two types: the full physics-based models and the hybrid empirical/physics-based models. The full physics-based models solve the MHD simulation from the solar surface to 1 au or beyond, and on the other hand, the inner boundary of the hybrid empirical/physics-based models is usually beyond the corona region, with a set of empirical formulas to specify the solar wind distribution at the lower boundary. Since the plasma β in the lower corona is low, the time step determined by the Courant-Friedrichs-Levy condition is small in the corona region. In this sense, the full physics-based models are computationally expensive and the hybrid empirical/physics-based models have their merits from the perspective of forecasting application because of the inexpensive computational cost and relative simplicity of implementation. Also, Owens et al. (2008) showed that the hybrid empirical/physics-based MHD model can reproduce large-scale solar wind structures comparable with the full physics-based models. The purpose of this study is to review the development status of numerical simulation for solar wind and coronal mass ejection in the inner heliosphere, and in particular, we investigate the different inner boundary conditions used in the models.

Numerical simulation of solar wind in the inner heliosphere
The lower boundary of the heliospheric models is usually located around 0.1 au, where the solar wind becomes supersonic and super-Alfvénic. Thus, the MHD solutions are predetermined by the treatment of the inner boundary conditions to a large extent. Different kinds of lower-boundary conditions existed in the inner heliosphere models.
In general, the radial component of the magnetic field B r at the lower boundary of the heliospheric models is usually obtained by the potential field source surface (PFSS) model (Altschuler and Newkirk, 1969;Schatten et al., 1969;Mackay and Yeates, 1999). In the PFSS model, a source surface (around 2.5 R s ) is defined, and the magnetic field between the solar surface and the source surface is assumed to be current-free and becomes open and purely radial above the source surface. With the magnetograph observations of the photosphere as the input, the PFSS model can rapidly extrapolate the large-scale coronal magnetic structures with reasonable accuracy. However, magnetic fields derived from the PFSS model have no thin heliospheric current sheet (HCS) or Parker spiral in the interplanetary space. To solve this problem, some authors combined the PFSS model with the Schatten current sheet (SCS) model (Schatten, 1971). The SCS model takes the absolute value of the (radial) field on the source surface obtained from the PFSS solution as a new lower-boundary condition, the new potential field is solved between the source surface and infinity, and then, the proper sign is restored. Another method for reconstructing the coronal magnetic field is the current sheet source surface (CSSS) model (Zhao and Hoeksema, 1995) that explicitly takes into account additional sheet currents. The CSSS model includes a cusp surface height R CS and a source surface height R SS , where the magnetic field becomes open above R CS and becomes radial above R SS . The Ulysses observations show that the strength of the interplanetary magnetic field B r has no dependence on latitude. However, the largescale structures of the coronal magnetic field derived from the PFSS model have a systematic gradient of B r in the latitudinal direction. The CSSS model can better reproduce the latitude-invariant nature of B r .
Based on the coronal magnetic topology parameters obtained from the previous models, the Wang-Sheeley (WS) or Wang-Sheeley-Arge (WSA)-type models (Arge and Pizzo, 2000;Arge et al., 2004;Mcgregor et al., 2011) are the most widely used empirical models that predict the solar wind velocity at the lower boundary. The specific form of V r (km/s) in the WS model can be written as where f s is the magnetic expansion factor which reads where B R s and B R are the magnetic field strength at the solar surface and at the source surface, respectively. The distribution of f s at the lower boundary of the heliospheric models is the same as its distribution at the source surface. V min sets the minimum possible solar wind speed, while V max sets the maximum possible solar wind speed. α determines the effect of f s on V r . By including an additional variable θ b , the specific form of V r (km/s) in the WSA model can be written as where θ b denotes the minimum angular separation between an open-field foot point and its nearest coronal hole boundary. The angle a 2 and exponent a 3 determine the angular extent and influence of the open-flux boundary layer, respectively. Any of the parameters in Eqs 1, 2 can be modified. The other components of the magnetic field and velocity, such as the meridional magnetic field B θ and the meridional flow velocity V θ , are always set to 0. The azimuthal flow velocity V ϕ is always set to zero in the inertial frame and set as Ω s R gb sin θ in the rotating frame, with Ω s denoting the solar angular speed and R gb standing for the lower boundary of the heliospheric models. The azimuthal magnetic field B ϕ is determined by Other solar wind parameters, including the density and temperature, are always prescribed by an assumption of the constant momentum flux and total pressure (sum of thermal and magnetic pressures). In the following section, the lower-boundary conditions in the existing inner heliospheric models are introduced in detail.

Boundary conditions based on the PFSS model
Using the PFSS + WS model, Odstrcil (2003) modeled the ambient solar wind from 0.1 au to 1 au using the WS-ENLIL mode. The solar wind velocity at 0.1 au was calculated according to Equation 1 with V min = 285 km/s, V max = 575 km/s, and α = 1/1.7. The mass density ρ was determined by an assumption of the constant momentum flux, and the temperature T was chosen to assure the total pressure is uniform on the source surface. Hayashi (2003) reconstructed three-dimensional solar wind structure from 50 R s to 1 au. The radial flow velocity V r at the inner boundary was determined by solar wind data obtained from interplanetary scintillation (IPS) observations, and the number density and temperature were determined by the empirical relation obtained from the Helios data through least-squares fitting, N(V r ) = 62.98 + 866.4( V r −154.9 100 ) −3.402 and T(V r ) = −0.455 + 0.1943V r 100 . They found that the solar wind research and space weather prediction algorithms can be enhanced by the combination of MHD simulation and IPS tomography. Based on that work, Hayashi (2012) modeled timedependent solar wind from 50 R s up to 1,250 R s . The solar wind was simulated with the time-dependent boundary condition. They applied temporal linear interpolation from one Carrington rotation to the next Carrington rotation. Shiota et al. (2014) developed an inner heliospheric MHD solar wind model for practical use in real-time space weather forecasting. The numerical domain was set to 25R s ≤ r ≤ 425R s . The radial flow velocity was determined using Equation 1 with V min = 267.5 km/s, V max = 410 km/s, and α = 0.4. The number density and temperature were obtained using empirical relationships that have been used by Hayashi (2003) and Hayashi (2012), respectively. Since Hayashi (2003) and Hayashi (2012) showed the relation at 50 R s , their model specified the boundary condition at 25 R s . For modification, the number density was multiplied by 4 assuming constant velocity propagation from 25 R s to 50 R s , and the temperature was multiplied by 4 γ−1 with the assumption of polytropic expansion.
Based on the PFSS + WSA model, Merkin et al. (2011);Pahud et al. (2012); and Merkin et al. (2016) simulated the solar wind from 0.1 au to 1 au or beyond using the Lyon-Fedder-Mobarry (LFM) MHD simulation code. This new version of the code is referred to as LFM-helio. The values of V r at the inner boundary were all calculated according to Equation 2 but with different parameters. In Merkin et al. (2011), V min = 200 km/s, V max = 750 km/s, a 1 = 2/9, a 2 = 3.8°, a 3 = 3.6, and a 4 = 3. The number density and temperature were determined by using an empirical fit to Helios data: N[cm −3 ] = 112.64 + 9.49 ⋅ 10 7 /(V r [km/s]) 2 , T = N 0 T 0 /N, where T 0 = 8 ⋅ 10 5 K and N 0 = 300cm −3 . Figure 1 shows the radial velocity and magnetic field at the inner boundary of the simulation (0.1 au) (Merkin et al., 2011). The magnetic field at 0.1 au was scaled down by a factor 2/3 × 100.
The inner boundary conditions used by Pahud et al. (2012) were similar to those used by Merkin et al. (2011), except for the following parameters: V min = 240 km/s, V max = 675 km/s, a 1 = 2/9, a 2 = 2.8°, a 3 = 1.25, and a 4 = 3. By comparing the MHD results with MESSENGER and ACE spacecraft observations, Pahud et al. (2012) found that the uncertainty in the specification of the boundary conditions, rather than a poor performance of the solar wind model, led to the discrepancies between in situ measurements and simulations. Merkin et al. (2016) used the LFM-helio MHD model to explore the heliospheric consequences of time-dependent changes at the Sun. The Air Force Data Assimilate Photospheric flux Transport (ADAPT) model was used to obtain daily updated photospheric magnetograms. These time-dependent magnetograms were then used to obtain the solar wind parameters at 21.5 R s with the WSA model. In particular, the changes of the longitudinal and latitudinal components of the magnetic field at the inner boundary were induced by the tangential electric fields, which were derived by a Poisson equation. The V r at the inner boundary was calculated according to Equation 2 with V min = 240 km/s, V max = 675 km/s, a 1 = 2/9, a 2 = 1.9°, a 3 = 2, and a 4 = 3. The number density was determined in the same way as in Merkin et al. (2011). The temperature was defined from the approximate full pressure balance in the angular directions.
where k is the Boltzmann constant, N slow = max(N), and T slow = 5 ⋅ 10 5 K is the nominal slow wind plasma temperature. In some cases, one can specify the fast wind temperature parameter T fast instead to avoid a negative temperature: where T fast = 10 6 K, B max is the maximum magnetic field located somewhere in the fast solar wind, and N ⋆ represents the plasma density at that location. Shen et al. (2018) and Wang et al. (2020) used a new boundary treatment to model the solar wind in interplanetary space with MHD equations. The magnetic field was derived by the PFSS model combining with the OMNI database. where R gb = 0.1 au. The average value of the observed magnetic field, mean(B 1au ), was determined by the past three Carrington rotations (CRs) at 1 au from the OMNI database. The solar wind velocity at 0.1 au was obtained from Equation 2. V max = 675 km/s, a 1 = 0.22 rad, a 3 = 1.0, and a 4 = 1.0 were constants over time, while V min and a 2 can be adjusted with different solar cycles. The number density was calculated by assuming that the solar wind energy flux is constant over the solar cycle.
where V 1au = 750 km/s is velocity at 1 au and N 1au is deduced from the average solar wind energy flux based on the OMNI observations. Various studies showed that there are linear or quadratic relationships between T and V at 1 au (Lopez and Freeman, 1986;Richardson and Cane, 1995;Elliott et al., 2005;Verbanac et al., 2011;Chat et al., 2012). The relationship of T p and V r at the lower boundary could be obtained as follows: Then, the pressure was obtained by p = 2NKT. Figure 2 presents the maps of solar wind parameters for CR2053 at the lower boundary of 0.1 au (Shen et al., 2018). Gonzi et al. (2021) investigated the impact of inner heliospheric boundary conditions on solar wind predictions at Earth. The radial component B r was computed using two methods, the PFSS model or the magnetofrictional (MF) model. The solar wind velocity was obtained from three different methods. The first was the WSA model in Equation 2 using V min = 240 km/s, V max = 675 km/s, a 1 = 2/9, a 2 = 2.8°, a 3 = 1.25, and a 4 = 3. The second was a modified DCHB model V r (θ b ) = V min + (V max − V min )(θ b ⋅ ω) δ , and the parameters were chosen by comparing model-predicted solar wind speed histograms at 1 au with observations V min = 200 km/s, V max = 700 km/s, ω = 7 radian, and δ = 1.25. The third was determined by the interplanetary scintillation solar wind data. The number density was determined from the interplanetary scintillation solar wind data or the assumptions of the constant momentum flux The temperature was set as T = T s N s N . Specifically, V s = 300 km/s, N s = 500cm −3 , and T s = 500 kK in the slow stream. The results showed that the selection of the coronal magnetic field model had a significant impact on the model results, and IPS data had a high success rate in reproducing high-speed solar wind. Detman et al. (2006Detman et al. ( , 2011 introduced a Hybrid Heliospheric Modeling System (HHMS) for background solar wind to aid in the operational forecasting of geomagnetic activity, where the lower boundary was also set at 0.1 au. The radial magnetic field was provided by the PFSS + SCS model. The principle in determining the solar wind speed at 0.1 au was similar to the WS model given by FIGURE 2 Maps of solar wind parameters at the lower boundary of 0.1 au for CR 2053. (A-D) The panels show the radial velocity Vr (km/s), number density N (cm −3 ), temperature T (K), and radial magnetic field Br (nT) (Shen et al., 2018).

Boundary conditions based on the PFSS + SCS model
where V max and V min were set at 700 and 200 km/s, respectively. The parameter V xp = 1.5 determined the sharpness of the fast-slow transition in the relationship. The parameter f 0 scaled the parameter f s , which was adjustable to match the observed solar wind speeds. The mass density ρ = F mass V r was determined by the assumption of an constant mass flux, F mass was obtained by taking observed values of ρV r at Earth and applying an r 2 scaling to 0.1 au. The temperature T = p total − B 2 8π 2ρR gas was determined by an assumption of the uniform total pressure p total on the source surface, and p total was a free parameter that required tuning. Figure 3 shows the block diagram of the HHMS in Detman et al. (2006). Narechania et al. (2021) described a 3D MHD-based heliospheric model based on a semi-empirical data-driven approach. The inner boundary of the heliospheric model was set at 25 R s , the magnetic field was derived by the PFSS + SCS model as in Detman et al. (2006), and the velocity was provided by Equation 2 with V min = 250 km/s, V max = 850 km/s, a 1 = 0.2, a 2 = 2.6°, a 3 = 1.25, and a 4 = 2.5. The remaining parameters were determined as in Shiota et al. (2014). Pomoell and Poedts (2018) presented a new space weather forecasting-targeted inner EUropean Heliospheric FORecasting Information Asset (EUHFORIA) model. The plasma parameters at 0.1 au were obtained by using the empirical solar wind model. The radial component B r was computed as where sgn(B corona ) is the sign of the magnetic field as given by the PFSS + SCS model. B max = 300 nT. The boundary value for the solar wind velocity was calculated according to Equation 2 with V min = 240 km/s, V max = 675 km/s, a 1 = 2/9, a 2 = 0.02 rad, a 3 = 1.25, and a 4 = 3. A constant value of 50 km/s is subtracted from the solar wind velocity to avoid systematically overestimating the wind speed. The number density is given as follows: where N fsw = 300cm −3 is the number density of the fast solar wind. The temperature is given as T = T fsw (ρ fsw /ρ) 2 with T fsw = 0.8 MK. Figure 4 shows the sequence of steps that provides the boundary conditions for the heliospheric MHD model (Pomoell and Poedts, 2018).

Boundary conditions based on the CSSS model
The inner heliospheric solar wind was simulated from 21.5 R s to 1 au with the MHD code CRONOS in Wiengarten et al. (2013) and Wiengarten et al. (2014). In Wiengarten et al. (2013), the CSSS model was used to obtain B r at 21.5 R s . The magnetic flux at the solar surface was obtained by the solar surface flux transport (SFT) model, which was then extrapolated to 10 solar radii by the CSSS model with R CS = 1.55R s and R SS = 10R s . Then,

FIGURE 3
Block diagram of the HHMS (Detman et al., 2006). the radial magnetic field at the source surface was scaled to the lower boundary at 21.5 R s by using the scale factor r −2 . In Wiengarten et al. (2014), the magnetic field at 0.1 au was derived by the PFSS + SCS model as in Detman et al. (2006). A Finite Difference Iterative Poisson Solver (FDIPS) was used to obtain the PFSS solution to avoid the numerical artifacts generated by the usual spherical harmonics expansion of the coronal potential field. At 0.1 au, the radial magnetic field was given a value of 300 nT, while keeping the

FIGURE 4
Sequence of steps that provide the boundary conditions for the heliospheric MHD model (Pomoell and Poedts, 2018).
orientation derived by the PFSS + SCS model. The radial velocity at 0.1 au was obtained with different forms for Wiengarten et al. (2013) and Wiengarten et al. (2014). Wiengarten et al. (2013) .4 km/s, and α 2 = 3.5. The velocity V r in Wiengarten et al. (2014) was calculated according to Equation 2 with V min = 200 km/s, V max = 675 km/s, a 1 = 2/9, a 2 ∈ (2.8°, 3.2°), a 3 = 2, and a 4 = 3. The remaining MHD quantities at the inner boundary were all similar to that described in Detman et al. (2006). Li et al. (2020) simulated the interplanetary B z using a datadriven heliospheric solar wind model. Following Merkin et al. (2016), the time-dependent boundary conditions were also used in their study. The radial component B r at the inner boundary was derived from the CSSS model with R CS = 2.5R s and R SS = 15R s . The value of B r given by the CSSS model was scaled by a factor b up = Φ OB /Φ CSSS to provide enough magnetic flux at the inner boundary, where Φ CSSS was the magnetic flux from the original CSSS solution and Φ OB = 4πr 2 1au mean(B 1 au )/ √ 2 was the observed average magnetic flux at 1 au derived from OMNI data. The WSA model (Equation 2) was used to determine the solar wind velocity with V min = 240 km/s, V max = 560 km/s, a 1 = 2/9, a 2 = 2°, a 3 = 1, and a 4 = 3. As in Detman et al. (2006), the density was obtained by an assumption of a constant mass flux. The distribution of temperature at the heliobase was similar to that used in Shen et al. (2018) and Wang et al. (2020). Biondo et al. (2021) described a new approach to determine the plasma parameters at the inner boundary of 0.1 au. This scheme applied a back reconstruction technique to remap the in situ measurements acquired at 1 au into the inner heliosphere. The observation data at 1 au were thought as a ring of L cells.

Boundary conditions based on in situ measurements
The Parker Spiral can be used to obtain their respective cells of origin (0.1au; ϕ ′ i ).
The radial speed and numerical densities at 0.1 au can be obtained as follows: The treatment of the magnetic field was similar.
This method could generate an interplanetary spiral reconstruction well with observations, and additionally, many small-scale features were also generated. Figure 5 presents the steps involved in RIMAP for the reconstruction of heliospheric plasma conditions (Biondo et al., 2021).

Numerical simulation of solar wind with the Heliospheric Upwind eXtrapolation model
To provide a computationally efficient solution, the Heliospheric Upwind eXtrapolation (HUX) model (Riley and Lionello, 2011) was widely used to extrapolate velocity from near the Sun to 1 au without providing physical insight. The HUX model was essentially a 1D extrapolation with velocity in a radial direction.
with Δr and Δϕ representing the grid spacing in radial and longitude directions, respectively. A comparative study of extrapolation models and empirical relations in forecasting solar wind was made in Kumar et al. (2020). The PFSS + SCS model was used to extrapolate magnetic fields. The velocity was obtained using the WS model or the WSA model. Then, the HUX model and a physics-based model PLUTO (Mignone et al., 2007) were used to extrapolate these velocities into the inner heliosphere zone, and the different magnetic field extrapolation models combined with velocity formulations were compared to predict the solar wind properties at L1. Figure 6 shows the sequence of steps used for each of these models. They showed that the WSA model captures the contrast between the slow and fast solar wind better than the WS model. Additionally, the PFSS + SCS magnetic field extrapolation model combined with the WSA model had the best performance in all the cases. Mayank et al. (2022) used a physics-based inner heliospheric model to forecast the ambient solar wind from 0.1 au to 2.1 au. The three-dimensional MHD equations were solved using the PLUTO code. The properties in the inner heliospheric model were derived from the PFSS + SCS model in combination with the WSA model (Equation 2). Additionally, the HUX algorithm was used to find the optimal set of values of independent parameters in the WSA model. The following parameters were chosen: V min = 250 km/s, V max = 650 km/s, a 1 = 2/9, a 2 = median of θ bE , a 3 = 1.25, and a 4 = 3. θ bE was calculated on the field lines that reached the location of Earth. Then, the WSA speed profile was rotated in the longitudinal direction by angle κ to account for solar rotation.
Ω s ranged from 27.21 to 27.34 days, depending on the location of Earth in its orbit. The radial component B r and plasma number density were similar to those used by Pomoell and Poedts (2018). The plasma thermal pressure was chosen to be constant on the boundary and equal to p = 6.6 nPa. Figure 7 presents the process flow diagram of the proposed solar wind model showing the range of numerical models involved in the inner heliospheric model.
The time-dependent form of the HUX model, which was referred to the HUXt model, was used to give solar wind velocity at Earth in Owens et al. (2020) and Bunting and Morgan (2022). Owens et al. (2020) set the inner boundary at 30 R s , and the computed V r at 30 R s from the Magnetohydrodynamics Algorithm outside a Sphere (MAS) coronal model serves as the input to HUXt. They suggested that HUXt can act as a surrogate for full threedimensional MHD models when very large ensembles are required. In Bunting and Morgan (2022), the inner boundary condition of the HUXt model was based on the coronal plasma density gained from coronagraph observations. The solar wind velocity was obtained using a simple linear relationship as follows: where N max and N min are the maximum and minimum densities in the equatorial plane, respectively. The ambient solar wind velocity at Earth was obtained by using the HUXt model. Compared to an HUXt model with a traditional MAS as the input, the results showed that the tomography/HUXt model can predict solar wind velocity much better.

FIGURE 6
Flowchart of the models utilized in the work (Kumar et al., 2020).

FIGURE 7
Process flow diagram of the proposed solar wind model showing the range of numerical models involved in the inner heliospheric model (Mayank et al., 2022).
Overall, various inner boundary conditions have been developed in recent years, which provide important input conditions for space weather prediction. Table 1 summarizes the various methods to set the solar wind parameters at a lower boundary. The MHD model and the HUX model have been widely used for numerical simulation of solar wind. The HUX model can act as a surrogate for full three-dimensional MHD models when very large ensembles are required. However, it can only be used to obtain the velocity in the interplanetary space. Thus, the MHD model is necessary when we want to get other parameters, such as magnetic field, temperature, and density. For MHD simulation, these variables are all needed to be determined at the lower boundary initially.The magnetic field and velocity distribution are important because temperature and density are always determined from them. To calculate the magnetic field distribution, PFSS + SCS may be a better choice because the magnetic fields derived from the PFSS model have no thin HCS or Parker spiral in the interplanetary space. To better reflect the latitude-invariant nature of B r , the CSSS model can be used. For calculating velocity distribution, the WSA model is better than the WS model in capturing the contrast between the slow and fast solar wind. The IPS observation has also been used to obtain velocity distribution, and the results show that IPS data has high success rate in detecting high-speed solar wind and the space weather prediction can be enhanced by the combination of MHD simulation and the IPS data.
It can be seen that the distribution of the inner boundary is fixed in many studies, the numerical processing on this fixed boundary is simple, and the steady-state solar wind can well represent the background conditions of the short-term transient phenomena using this fixed boundary. Since some studies used time-dependent boundary (Hayashi, 2012;Merkin et al., 2016;Li et al., 2020), it is expected that time-dependent boundary conditions can simulate the inner heliosphere more accurately and in real time. Fixed boundary conditions can accommodate observational data made in only one CR period. However, time-dependent boundary conditions can determine the MHD state of the solar wind at distant regions from the Sun, which are determined in a longer time period than one CR. Also, the heliospheric consequences of time-dependent conditions can be explored, that is, what kind of plasma and magnetic field structures are created due to the changing boundaries of coronal holes. Merkin et al. (2016) showed that time-dependent boundary conditions can reproduce the gross-scale structure of the heliosphere with higher fidelity, and provide important insights to interpret data on smaller spatial and faster time scales (e.g., 1 day). Thus, the timedependent boundary condition is a promising direction of research both for space weather applications and fundamental physics of the heliosphere. For the time-dependent boundary, how to maintain the divergence-free condition is one difficulty. Since the time-dependent radial magnetic field ∂B r ∂t is always obtained by observations, we need to derive an electric field tangential to the inner boundary spherical surface that would evolve the radial magnetic field. To do that, the Helmholtz theorem is applied on a sphere and the tangential electric field is expressed as the sum of the curl and the gradient of scalar potentials.
Here, ψ is determined by solving the Poisson equation following Faraday's law.
Because there is no available information to determine Φ, Φ is always set to 0. Then, the transverse components of the magnetic field in the corotational frame can be calculated as follows: 3 Numerical simulation of the coronal mass ejection in the inner heliosphere CMEs are giant clouds of magnetized plasma erupting from the Sun. When CMEs propagate through the corona and the interplanetary space, they interact with the solar wind plasma or other ICMEs and undergo changes. At larger heliocentric distances, CMEs are known to exhibit rotations and deflections, which can play a crucial role in their geomagnetic impact (Maharana et al., 2022). MHD simulations are often used for CME modeling purposes. There are two categories of MHD simulation models, one is anchored to the low corona and the other is fully embedded into the middle or upper corona. The first category is built through a large number of observations and enables the incorporation of more detailed physics, but the total computation time is much larger. On the other hand, the second category is suitable for space-weather forecasts because of the inexpensive computational cost and relative simplicity of implementation. Also, with the launch PSP and SO, the ICME will be constrained by more observations. In the past decades, many different models have been developed to simulate CME propagation in the heliosphere, such as ENLIL (Odstrcil, 2003), EUHFORIA (Pomoell and Poedts, 2018), SUSANOO-CME (Shiota and Kataoka, 2016), and MS-FLUKSS (Singh et al., 2020). In these models, CMEs are inserted at 0.1 au and evolved self-consistently by solving the MHD equations and a suitable CME initialization model is necessary. In this section, we focus on the initial CME parameters that need to be determined and extrapolated to the lower boundary.
The initial CME parameters rely heavily on CME reconstruction methods. The cone model (Zhao et al., 2002;Xie et al., 2004;Michaek, 2006) and the graduated cylindrical shell (GCS) model (Thernisien et al., 2006(Thernisien et al., , 2009Thernisien, 2011) are widely used to estimate the CME kinematic and geometric parameters in the past years. The cone model represents CME as a hydrodynamic pulse with a constant angular width, propagation direction, and speed. The angular width and source position of the CME are obtained by matching the modeled halos with the observed halos. The radial velocity and acceleration of the CME are determined by a series of times and radial distances. Figure 8 shows the cone model topology and relationship between the heliocentric coordinate system and the cone coordinate system (Michaek, 2006). The heliocentric coordinate system is defined such that x h points to Earth, z h points north, and y h , z h defines the sky plane. The cone coordinate system has the origin at the apex of the cone, z c is the cone axis, and x c , y c defines the plane parallel to the base of the cone. The geometry of the cone can be determined by heliographic longitude ϕ, latitude θ, and angular width ω.For the GCS model, the CME is described as a flux-rope-like structure that expands in a self-similar fashion. The three-dimensional shape of the CME consists of two legs and a curved front. Figure 9 shows a face-on schematic of the GCS model (Thernisien et al., 2009) and schematic of the detailed geometric parameter (Thernisien, 2011). The dashed-dotted line shows the axis of the model, and the solid line the represents of the shell where the density is placed. The cross section of the model is a group of circular annuli with a gradually varying radii r 1 = κr, where r is the distance from a point on the shell to the center of the Sun and κ is a constant depending on the studied event. The angle between the axes of two conical legs is 2ω, and the height of the cone is h. The geometry of the shell can be completely defined with these three parameters. In the  Shen et al., 2018;Wang et al., 2020 In situ measurements

In situ measurements
In situ measurements Biondo et al., 2021 Frontiers in Astronomy and Space Sciences

FIGURE 8
Cone model topology and relationship between the heliocentric coordinate system and the cone coordinate system (Michaek, 2006).
following section, lower boundary conditions in the existing inner heliospheric models are introduced in detail.

Numerical simulation of coronal mass ejection with hydrodynamic plasma cloud
Based on the cone model, the ENLIL model was used to simulate transient heliosphere disturbances (Odstrcil, 2003;Odstrcil et al., 2004;Odstrcil, 2005). A CME was simulated by launching a time-dependent hydrodynamic plasma cloud through the inner boundary of the heliospheric. The cone model was used to describe the location, speed, and angular width of the CME. The CME disturbance had uniform velocity, density, and temperature. The CME's density (temperature) was taken to be four times larger than (equal to) the mean values in the fast stream. With the simple geometry, the ENLIL model had been successfully used to study the global evolution of CMEs in the heliosphere and predicted the CME arrival times at Earth. Figure 10 presents the distribution of solar wind parameters with the introduction of the input pulse on 12 May 1997 at 1800 UT (Odstrcil, 2005).
With the development of the EUHFORIA heliospheric model (Pomoell and Poedts, 2018), a new CME propagation model simulation was performed. The cone CME model was used to determine the CME angular width, propagation direction, and speed. Then, the CME was introduced to the inner radial boundary at 0.1 au as a time-dependent boundary condition. Specifically, the velocity, density, and temperature for solar wind values were replaced by constant values at R gb = 0.1 au.
where the width as a function of time is where t half = R gb tan (ω CME /2)/v CME , θ CME , and ϕ CME were the colatitude and longitude of the propagation direction of the CME center, ω CME was the angular width of the CME, and v CME was the velocity of the CME. These parameters were provided by the cone model. For simplicity, the density and temperature of the modeled CMEs were taken to be density ρ CME = 10 −18 kgm −3 and FIGURE 9 (A) Face-on schematic of the graduated cylindrical shell model (Thernisien et al., 2009), (B) Schematic of the detailed geometric parameter (Thernisien, 2011).

FIGURE 10
Distribution of solar wind parameters with the introduction of the input pulse on 12 May 1997 at 1800 UT (Odstrcil, 2005).
temperature T CME = 0.8 MK, respectively. Scolini et al. (2018) used the EUHFORIA heliospheric model to test the effect of different CME shapes on simulation outputs, such as CME radius, the CME cross-section, and velocity initialization at the inner boundary. They showed that these geometrical parameters particularly affect predictions at locations hit by the CME, and the definition of the CME radius was the parameter having the greatest impact on simulation outputs. Verbeke et al. (2022) presented a new inner heliospheric model, ICARUS, for the simulation of the propagation and evolution of superposed CMEs. The cone CME model described in Scolini et al. (2018) was employed. For both scientific and forecasting purposes, the authors showed that the properties of radial grid stretching and solution adaptive mesh refinement were useful to save computational time.

Numerical simulation of coronal mass ejection with the flux-rope model
The geoeffectiveness of a CME is largely due to its magnetic field. If the CME contains a positive Bz component, it will be favorable for magnetic reconnection with Earth's magnetosphere and results in a strong geomagnetic storms. The model described so far has not yet predicted the magnetic field structure of CMEs arriving at Earth. In order to overcome the limitations of the cone model, recent work has focused on modeling the CME using more realistic flux-rope models, such as spheromaks or toroidal-like structures.
The SUSANOO-CME model, developed by Shiota and Kataoka (2016), was used to simulate the interplanetary propagation of multiple CMEs with an internal magnetic flux-rope structure. The CME passed through the inner boundary at 30 R s as the time evolution of the MHD variables. They adopted a spheromak-type magnetic flux rope, and radial compression was adopted for the flux rope to make the CME pass quickly. First, the CME model was defined as a simple linear force-free spheromak model near the Sun.
B̃r (r ,θ ,φ ) = (2B 0 /αr ) j 1 (αr ) cosθ , where j 1 and j 0 are the first-order and zero-order spherical Bessel function, respectively. The CME was defined in local spherical coordinates (r ,θ ,φ ), the center of the CME was chosen to be the origin. Then, the CME model was shifted outward in the radial direction from the Sun to the inner boundary. The spherical cut R ′ cut of the shifted magnetic configuration moved across the computational boundary as where R ′ cut (t 0 ) = 30R s defines the nose of the deformed spheromak.
is the start time when the CME appeared at the Frontiers in Astronomy and Space Sciences  inner boundary. Finally, at the inner boundary, the magnetic field is given as where a = 30R s − r 0 − r mc is the radial shift. The CME radius was r mc = 1R s ⋅ sin ω/2, and r 0 = 1R s + r mc was the distance between the central position of the spheromak and the heliospheric center. The model had 10 free parameters. t onset was the onset time and v CME was the propagation speed, which can be found from the LASCO CME catalog; λ s and ϕ s were heliographic latitude and longitude, respectively, which came from the flare list in National Geophysical Data Center (NGDC, https://www.ngdc.noaa.gov/stp/ spaceweather.html); τ was the tilt angle of the spheromak (±90°with Hale-Nicholson law), χ was the inclination angle of the spheromak, H was the chirality of helicity in the spheromak (1, set −1 if opposite to the Bothmer-Schwenn rule), Φ mag was the magnetic flux within CME, and w = 60°was the angular width of CME. The density inside the CME was in proportion to the magnetic field strength. The CME pressure was assumed to be a constant, which was defined as ρ(v CME − v SW ) 2 /2. Because of the constant pressure and high density, the temperature was the lowest at the center of the flux rope. Liu et al. (2019) used a three-dimensional flux-rope CME initialization model based on the GCS model to simulate the propagation and deflection of a CME from 21.5 R s to 244 R s . The Lundquist force-free CME model in a cylindrical geometry (r, ϕ, z) was used.
The axial component of the magnetic field was assumed to be zero at the edge of the flux rope; thus, they could get where R stands for the radius of the cylindrical shell, which was obtained from the GCS model. The maximum of magnetic field B 0 is defined as with j = ∫ R 0 j 2 1 (αr )dr, H m is the magnetic helicity of the Lundquist flux rope and L = 2ω(H front − R front ) is the flux-rope length, and H front and R front are the maximum height and maximum radius of the front, respectively, that were obtained from the GCS model. The initial density, velocity, and temperature at point P in the CME are given as follows:

FIGURE 11
Magnetic field lines depicting the structure of the LFFS CME model with a 90.0°tilt (Verbeke et al., 2019).
Within the GCS model, r represents the radius of the cross section and d stands for distance from point P to the center point at the front part, while defined as the distance from point P to the inner side of the cone in the cross-section plane. Table 2 shows the initial CME parameters (Liu et al., 2019). Shen et al. (2021) made a numerical research on how the initial CME parameters can affect simulation results. They found that the initial density and geometric size of the CME both had an effect on its arrival time at Earth. The initial magnetic field had a large effect on the CME's geomagnetic effect. All of these results confirmed that the initial geometric and physical parameters had an important role on space-weather research and forecasts.

FIGURE 12
Three-dimensional view of the magnetic field for CASE-1 and CASE-2 at different times (Zhang et al., 2021). The EUHFORIA CME model in Pomoell and Poedts (2018) and Scolini et al. (2018) did not contain an intrinsic magnetic field and was injected at the inner radial boundary with the enhancement of velocity, density, and pressure. Verbeke et al. (2019) introduced a linear force-free spheromak (LFFS) solution to the EUHFORIA CME model. The CME was considered to be a sphere with a radius r mc = 0.1au ⋅ sin ω/2. The CME latitude and longitude (θ CME , ϕ CME ), angular width ω, and speed v CME were the same as in Pomoell and Poedts (2018). Located at (0.1au − r mc , θ CME , ϕ CME ), the CME moved through the 0.1 au boundary with speed v CME . A point at (0.1au, θ, ϕ) satisfied the following formula was inside the CME where (x CME ; y CME ; z CME ) and (x bound ; y bound ; z bound ) were the coordinates of the CME center and the considered point on the boundary in Cartesian HEEQ coordinates, respectively. The magnetic field was defined in a local spherical coordinate system (r ,θ ,φ ) as in Shiota and Kataoka (2016), with the center of the spheromak as the origin. As in Pomoell and Poedts (2018), density ρ CME = 10 −18 kgm −3 and temperature T CME = 0.8 MK. The initial speed required for the spheromak CME model should reflect only the radial speed.
The remaining magnetic input parameters required in the spheromak CME model were the flux-rope handedness, tilt, and the toroidal magnetic flux at 0.1 au. The total toroidal flux Φ was used to define the magnetic field strength B 0 .
In that paper, a +1 handedness with a 90.0°tilt angle were selected to obtain a negative Bz magnetic field when the CME was arriving at Earth. Lastly, the total toroidal flux was set to be 10 14 Wb to produce the best results they performed so far. Figure 11 shows the magnetic field lines depicting the structure of the LFFS CME model with a 90.0°tilt (Verbeke et al., 2019).

FIGURE 13
FRi3D flux rope emerging out of the EUHFORIA's inner heliospheric boundary at 0.1 AU (Maharana et al., 2022). Zhang et al. (2021) performed a numerical study of two injection methods for CME in the inner heliosphere, where CASE-1 CME was introduced into the inner boundary with a radial compression (Shiota and Kataoka, 2016) and CASE-2 CME was introduced into the inner boundary without a radial compression (Scolini et al., 2019;Scolini et al., 2020;Verbeke et al., 2019). The CME angular width was obtained by using the sky-plane-projected speed range (Gopalswamy et al., 2010). The CME mass density and temperature were set to be homogeneous as in Pomoell and Poedts (2018); Verbeke et al. (2019). The CME position, velocity, and magnetic input parameters were estimated from image modeling and geometric triangulation analysis following Liu et al. (2010). They found that CASE-2 CME overestimated the radial extension at 1 au and the modeled magnetic fields at 1 au were lower compared to CASE-1 CME. Figure 12 shows the three-dimensional view of the magnetic field for CASE-1 and CASE-2 at different times. The top row shows the three-dimensional view of the magnetic field at 0.4 h (left), 2.0 h (right) after the addition of the CME into the ambient solar wind for CASE-1, while the bottom row shows the magnetic field at 5 h (left), 15 h (right) after the addition of the CME into the ambient solar wind for CASE-2. Solid black lines display the magnetic field lines and the color code represents the magnetic field B in units of 10 5 T. The blue sphere represents the 20 R s (Zhang et al., 2021).
Cone CMEs in Pomoell and Poedts (2018) were initialized using a set of seven input parameters during the CME insertion at the inner boundary of the heliospheric model. For example, the CME insertion time, speed v CME , source (θ CME , ϕ CME ), and half angular width ω/2 at 0.1 au were usually derived from cone model. In addition, the CME mass density and temperature were set to be homogeneous with density ρ CME = 10 −18 kgm −3 and temperature T CME = 0.8 MK. When simulated CMEs use the spheromak model as in Verbeke et al. (2019), three additional input parameters were needed: the helicity sign (chirality), the tilt, and the toroidal magnetic flux at 0.1 au. Since the spheromak CME in Verbeke et al. (2019) was initialized from observations partially, Scolini et al. (2019Scolini et al. ( , 2020 aimed to constrain all the CME input parameters from remote-sensing observations. The CME geometric and kinematic parameters were derived using the GCS model: the CME direction, the height of the CME apex h front , the tilt angle around the axis of symmetry γ (with respect to the solar equator), the half angle between the legs ω, and the half angle of the cone δ related to the "aspect ratio" κ by the relation κ = sin δ. The 3D speed at the CME apex was derived from the GCS model with a sequence of images.
v CME = dh front dt .
Then, the radial and expansion contribution were obtained as follows: The flux-rope handedness was determined from pre-eruption EUV sigmoids observations. The flux-rope tilt angle/orientation was inferred from the orientation of the source region polarity inversion line (PIL) and/or from the orientation of the post-eruption arcades (PEAs). The flux-rope toroidal magnetic flux was based on the reconnected flux, which was computed by the FRED method described by Gopalswamy et al. (2017). Using the PEA area, the FRED method computed the total (unsigned) magnetic flux from line-of-sight magnetic field data. Then, the reconnected flux ϕ RC could be obtained as follows: The axial field strength B 0 was estimated with the following formula: where r ⋆ is the distance from the center of the spheromak, on the plane θ = π 2 , where the magnetic field becomes completely axial. ϕ p (r ⋆ ) was the poloidal magnetic flux that is transferred from the reconnected flux. Maharana et al. (2022) introduced a fully analytical "Flux Rope in 3D" (FRi3D) CME model to EUHFORIA (Pomoell and Poedts, 2018) to improve the modeling of CME flank encounters. The geometrical parameters for the FRi3D flux rope contained the half-width φ hw in the azimuthal direction, the angular half-height φ hh in the polar direction, the toroidal height R t defined the heliocentric distance to the apex of the CME axis, a coefficient n defined the deformation at the CME front, pancaking φ p meant the deformation due to radial expansion, and skew φ s defined the deformation due to solar rotation. The geometrical parameters were obtained by using the forward modeling tool of FRi3D to remote multi-viewpoint observations. The Lundquist model in cylindrical geometry was used to get the magnetic field configuration of the CME. The magnetic field parameters, such as the tilt, the magnetic flux, the twist, the chirality, and the polarity, were obtained

FIGURE 14
Magnetic field line configuration of the modified spheromak inserted at the inner boundary R = 0.1 AU, shown with the red sphere (Singh et al., 2020).
by fitting the FRi3D CME model to WIND spacecraft at 1 au (Isavnin, 2016). The FRi3D CME model was also filled with uniform density plasma with density ρ CME = 10 −17 kgm −3 and temperature T CME = 0.8 MK. Figure 13 shows a FRi3D flux rope emerging out of EUHFORIA's inner heliospheric boundary at 0.1 au. The color of the field lines is based on the magnetic field strength. The field lines are twisted and deformed as per a particular CME geometry and have maximum strength close to the axis which reduces outward (typical Lundquist model characteristics) (Maharana et al., 2022). Singh et al. (2020) applied a modified spheromak model for simulations of the CME in the inner heliosphere. The simulations were carried out using a MultiScale Fluid-Kinetic Simulation Suite (MS-FLUKSS) code. They introduced a modified spheromak flux rope into the solar wind, where the poloidal and toroidal fluxes were set up independently.
where α 0 r mc = 5.763, 459. The origin of the spherical coordinate system was placed at the spheromak center. The flux-rope parameters were introduced into the solar wind as follows: where ξ was related to plasma β as β ∼ (γ − 1) (ξ − 1) indicated that the magnetic pressure in the flux rope was larger than the pressure in the solar wind. The mass inside the spheromak was a constant with 1.65 × 10 15 g. The thermal pressure was assumed to be proportional to the magnetic pressure. The poloidal flux and the toroidal flux were set to be 2 × 10 22 Mx and 5 × 10 21 Mx, respectively. A positive helicity sign was used. Figure 14 shows the magnetic field line configuration of the modified spheromak inserted at the inner boundary (Singh et al., 2020). Table 3 summarizes the various methods to set the CME parameters at lower boundary. Since the ability of a CME to disturb the near-Earth space is largely due to its magnetic fields, the flux-rope model is widely used in recent years because of its internal magnetic fields structure. How to constrain the initial CME parameters by observations is one difficulty. The cone model and GCS model are widely used to estimate the CME kinematic and geometric parameters. The magnetic field parameters, such as handedness, tilt, and toroidal magnetic flux, are deduced from remote-sensing observations in Scolini et al. (2019Scolini et al. ( , 2020. With the launch of the PSP and SO into the inner heliosphere, the ICME will be constrained by more observations.

Conclusion and discussion
Overall, the three-dimensional MHD numerical model has become an important tool to report interplanetary solar wind and the CME. The MHD heliospheric model is relatively cheap in terms of calculation, and it takes tens of minutes to hours to simulate the CME propagation from the Sun to Earth. The operational spaceweather forecast is usually generated by the heliospheric model. The resulting forecast errors are largely dependent on the uncertainty of the inner boundary conditions, rather than the heliospheric model.
For the solar wind simulation in the inner heliosphere, the solar wind density and temperature at the inner boundary are always derived as a function of the solar wind speed and magnetic field; thus, the speed and magnetic field can significantly affect the accuracy of the simulated results. The speed specified by the WSA formula given in Eq. 2 seems to be the most efficient formula. However, there are several free parameters in the WSA model, whose uncertainties may result in the deviations between the simulated results and observations. Thus, how to make parameter optimization is a major problem. The dynamic time warping (DTW) algorithm described in Bunting and Morgan (2022) can be used to derive the optimized values for the speed in the WSA formula. DTW algorithm is effective to quantify the agreement between two time series and has recently been used in a variety of fields such as economics, biology, and space weather. Here, DTW distance is used to measure the difference between the model and in situ solar wind velocities, and a smaller DTW distance means that the modeled data is in good agreement with the in situ data. The initial parameters of the WSA formula were adjusted to minimize the DTW distance, and then, we can obtain the optimized parameters in WSA formula. Also, the IPS observation is a good choice to obtain the velocity distribution, and the results show that IPS data have a high success rate in detecting high-speed solar wind and the space-weather prediction can be enhanced by the combination of MHD simulation and the IPS data. For the magnetic field at the inner boundary, the PFSS model is always used. However, the PFSS model underestimated the magnetic flux in the heliosphere. Shen et al. (2018) only kept the polarity of the magnetic field from the PFSS model and used the observational data at 1 au in the immediate past to limit the value of the magnetic field, and the results showed that the modeled magnetic field at 1 au agrees with the in situ observations much better. Additionally, the magnetic fields derived from the PFSS model have no thin HCS or Parker spiral in the interplanetary space; thus, the PFSS + SCS model may be a good choice.
For the CME simulation in the inner heliosphere, CMEs are always initialized by flux-rope models, such as spheromaks or toroidal-like structures, and the CMEs are introduced to the inner radial boundary as a time-dependent boundary. Scolini et al. (2019Scolini et al. ( , 2020 suggested that the prediction of the CME arrival time at Earth was found to be highly dependent on the CME model and CME input parameters used. Thus, how to obtain the CME parameters at the inner boundary is important. At present, CME geometric and kinematic parameters, such as CME position, radius, and velocity, are always derived from the GCS model or the cone model. The density and temperature of the CMEs are assumed to be constant. To improve the prediction of the CME internal magnetic field at Earth, it is necessary to constrain the flux-rope magnetic parameters from the remote-sensing observations of the corona. Scolini et al. (2019Scolini et al. ( , 2020 constrained flux-rope magnetic parameters from remotesensing observations, and the results showed that the current model can improve the predictions of B z up to 22-40 percentage points compared to a cone model.
At present, the numerical simulation in the inner heliosphere has a large dependence on the empirical model. How to use the observations to constrain the uncertainties in the empirical models is something we must consider in future. Maybe we can use the in situ data in the past CR to correct the uncertainties in the empirical models at the current CR, but this can only be applicable to the solar wind. With the multi-spacecraft observations in the future, such as SO and PSP, the detailed measurements of the inner heliosphere will be provided. All of these will provide us with more realistic inner boundary conditions, and we will have a deeper understanding of the triggering and propagation evolution of solar wind and disturbance.

Author contributions
XF designed the research. MZ and HL participated in data analyses and the writing of the manuscript. MX, FS, LY, XZ, YZ, and XL participated in data interpretation and discussion. All authors read and agreed to the published version of the manuscript.