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.1943 V 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 time-dependent 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 ₁ = 2/9, a ₂ = 3.8°, a ₃ = 3.6, and a ₄ = 3. The number density and temperature were determined by using an empirical fit to Helios data: N [ c m − 3 ] = 112.64 + 9.49 ⋅ 1 0 7 / ( V r [ k m / s ] ) 2 , T = N ₀ T ₀/N, where T ₀ = 8 ⋅ 10⁵ K and N ₀ = 300cm⁻³. 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 ₁ = 2/9, a ₂ = 2.8°, a ₃ = 1.25, and a ₄ = 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 ₁ = 2/9, a ₂ = 1.9°, a ₃ = 2, and a ₄ = 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. N k T + B r 2 8 π = N slow k T slow , (3) where k is the Boltzmann constant, N _slow = max(N), and T _slow = 5 ⋅ 10⁵ 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: N k T + B r 2 8 π = N ⋆ k T fast + B max 2 8 π , (4) where T _fast = 10⁶ 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. B r = sign B PFSS × 1 2 m e a n B 1  au 1 a u R g b 2 , (5) where R _gb = 0.1 au. The average value of the observed magnetic field, mean (B ^{1 au}), 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 ₁ = 0.22 rad, a ₃ = 1.0, and a ₄ = 1.0 were constants over time, while V _min and a ₂ 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. N = N 1au 1  au R g b 2 V 1au 1 2 V 1au 2 + G M s R s V r 1 2 V r 2 + G M s R s − 1 , (6) 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: T p = 1 2 V r 2 1  au R g b 2 γ − 1 . (7)

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 ₁ = 2/9, a ₂ = 2.8°, a ₃ = 1.25, and a ₄ = 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 NV ², N = N s V s V r 2 . The temperature was set as T = T s N s N . Specifically, V _s = 300 km/s, N _s = 500cm⁻³, 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. (2006, 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 the equation V r f s = V min + V max − V min / 1 + f s f 0 V xp , (8) 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 ₀ 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 ² 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 ₁ = 0.2, a ₂ = 2.6°, a ₃ = 1.25, and a ₄ = 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 B r = sgn B corona B max V r / V max , (9) 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 ₁ = 2/9, a ₂ = 0.02 rad, a ₃ = 1.25, and a ₄ = 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: N = N fsw V max / V r 2 , (10) where N _fsw = 300cm⁻³ 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).

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, 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 ⁻². 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 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) used V r = V min + V del f s V xp + α 1 ⋅ ( θ b f s V xp / 2 − α 2 ) , where V _min = 154 km/s, V _del = 300 km/s, V _xp = 3, α ₁ = 7.4 km/s, and α ₂ = 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 ₁ = 2/9, a ₂ ∈ (2.8°, 3.2°), a ₃ = 2, and a ₄ = 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 data-driven 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 Φ O B = 4 π r 1au 2 m e a n ( 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 ₁ = 2/9, a ₂ = 2°, a ₃ = 1, and a ₄ = 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. ϕ i r a d , N i c m − 3 , V r , i k m / s , B r , i n T , B ϕ , i n T . (11)

The Parker Spiral can be used to obtain their respective cells of origin ( 0.1 a u ; ϕ i ′ ) . ϕ i ′ = ϕ i + Ω s v r , i 1 a u − 0.1 a u . (12)

The radial speed and numerical densities at 0.1 au can be obtained as follows: V r , i ′ = V r , i , V θ , i ′ = 0 , V ϕ , i ′ = − Ω s , r N i ′ = N i 1 a u 0.1 a u 2 .

The treatment of the magnetic field was similar. B r , i ′ = B r , i 1 a u 0.1 a u 2 , B θ , i ′ = 0 , B ϕ , i ′ = B ϕ , i 1 a u 0.1 a u + B r , i Ω s V r , i 1 a u − R s s − 0.1 a u − R s s 0.1 a u 1 a u 2 .

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).

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. V r + 1 , ϕ = V r , ϕ + Δ r Ω s V r , ϕ V r + 1 , ϕ − V r , ϕ Δ ϕ , (13) 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 ₁ = 2/9, a ₂ = median of θ _bE , a ₃ = 1.25, and a ₄ = 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. κ = 5 + 2 π Ω s 21.5 R s V r min . (14)

Ω _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 three-dimensional 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: V = V max − N − N min N max − N min V max − V min , (15) 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.

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 time-dependent 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. E t = ∇ t × ψ r + ∇ t Φ . (16)

Here, ψ is determined by solving the Poisson equation following Faraday’s law. ∇ t 2 ψ = ∂ B r ∂ t . (17)

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: B ϕ = B r V ϕ V r + E θ V r , (18) B θ = − E ϕ V r . (19)

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. θ − θ CME 2 + ϕ − ϕ CME 2 < α t 2 , where the width as a function of time is α t = ω CME 2 sin π 2 t − t onset / t half , where t _half = R _gb  tan (ω _CME/2)/v _CME, θ _CME, and ϕ _CME were the co-latitude 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⁻¹⁸ kgm ⁻³ and 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.

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 ̃ , θ ̃ , ϕ ̃ = 2 B 0 / α r ̃ j 1 α r ̃ cos θ ̃ , B ̃ θ ̃ r ̃ , θ ̃ , ϕ ̃ = − B 0 j 1 α r ̃ α r ̃ + j 1 ′ α r ̃ sin θ ̃ , B ̃ ϕ ̃ r ̃ , θ ̃ , ϕ ̃ = H ⋅ B 0 j 1 α r ̃ sin θ ̃ , where j ₁ and j ₀ 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 R cut ′ t = R cut ′ t 0 − v CME t − t 0 , where R cut ′ ( t 0 ) = 30 R s defines the nose of the deformed spheromak. t 0 = t onset + ( 30 R s − r mc ) v CME is the start time when the CME appeared at the inner boundary. Finally, at the inner boundary, the magnetic field is given as B ′ 30 , θ ̃ , ϕ ̃ = B ′ R cut ′ t , θ ̃ , ϕ ̃ = B ̃ R cut ′ t − a , θ ̃ , ϕ ̃ , where a = 30R _s − r ₀ − r _mc is the radial shift. The CME radius was r _mc = 1R _s ⋅ sin ω/2, and r ₀ = 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. B r = 0 , B ϕ = H B 0 j 1 α r ̃ , B z = B 0 j 1 α r ̃ .

The axial component of the magnetic field was assumed to be zero at the edge of the flux rope; thus, they could get α R = 2.405 , where R stands for the radius of the cylindrical shell, which was obtained from the GCS model. The maximum of magnetic field B ₀ is defined as B 0 = 2.405 H m 4 π L R j , with j = ∫ 0 R j 1 2 ( α r ̃ ) d r , 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: ρ CME = ρ max 1 − d 2 2 r 2 , v CME = v max 1 − d 2 2 r 2 , T CME = T max 1 − d 2 2 r 2 .

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.

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 x CME − x bound 2 + y CME − y bound 2 + z CME − z bound 2 ≤ r m c , (20) 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⁻¹⁸ kgm ⁻³ and temperature T _CME = 0.8 MK. The initial speed required for the spheromak CME model should reflect only the radial speed. v rad = v cme ⋅ 1 − 1 1 + 1 2 1 + cot ⁡ w . (21)

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 ₀. Φ = ∬ B ̃ ϕ ̃ r ̃ d r ̃ d θ ̃ = 2 B 0 α 2 ⋅ − sin α r 0 + ∫ 0 α r 0 sin ⁡ t t d t . (22)

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¹⁴ 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).

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⁵ 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⁻¹⁸ kgm ⁻³ 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. (2019, 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 = d h front d t .

Then, the radial and expansion contribution were obtained as follows: v rad = 1 1 + κ d h front d t , v exp = κ 1 + κ d h front d t .

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: ϕ RC = 1 2 ∫ PEA | B los | d A = 1 2 | B los | tot A PEA .

The axial field strength B ₀ was estimated with the following formula: B 0 = α 3 2 π ϕ p r ⋆ r ⋆ sin α r ⋆ − α r ⋆ ⁡ cos α r ⋆ , 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 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⁻¹⁷ kgm ⁻³ 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. b ⃗ = 1 r ⁡ sin ⁡ θ ζ 1 r ∂ A ∂ θ r ̂ − ζ ∂ A ∂ r θ ̂ + H α 0 A ϕ ̂ , A = 4 π a 1 α 0 2 r 0 2 g α 0 r − r 2 sin 2 ⁡ θ , g α 0 r = sin α 0 r α 0 r − cos α 0 r , where α ₀ 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: b ⃗ final = b ⃗ F R , ρ final = ρ F R + ρ S W , e final = ξ e F R + e S W , 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¹⁵ 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²² Mx and 5 × 10²¹ 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. (2019, 2020). With the launch of the PSP and SO into the inner heliosphere, the ICME will be constrained by more observations.