Global Simulation of the Solar Wind: A Comparison with Parker Solar Probe Observations during 2018–2022

The computational domain for the 3D MHD simulation is a heliocentric spherical coordinate system (r, θ, ϕ) with the Sun–Earth line in the 180° longitude meridional plane. Note that the latitude of the Earth changes with seasons within ± 72 in this coordinate system, e.g., θ = −72 in early March, θ = 0° in early June and December, and θ = +72 in early September. The simulation domain covers −875 ≤ θ ≤ 875, 0° ≤ ϕ ≤ 360°, and 18 R_⊙ ≤ r ≤ 345 R_⊙. An open (outlet) boundary condition is adopted at θ = ±875 and r = 345 R_⊙, so there are no reflective disturbances. A constant grid size of Δr = 3 R_⊙, Δθ = 5°, and Δϕ = 5° is used, which results in (110 × 36 × 72) grids. Note that a higher resolution is always desirable. The angular resolution used in this study matches the angular resolution of the driving data (e.g., source surface field). There is no significant benefit from having a finer grids in the angular resolution. Although we can reduce the size of radial grids, we have to consider the time it takes to run the simulation. A simulation of this scale (5 yr without gaps) is challenging and has not been performed in the past.

The G3DMHD model is driven by Wilcox Solar Observatory (WSO) solar photospheric maps extrapolated to 2.5 R_⊙ and the expansion factor (Y.-M. Wang & N. R. Sheeley 1990a, 1990b) using the potential field source surface (PFSS) model (M. D. Altschuler & G. Newkirk 1969; K. H. Schatten et al. 1969; J. T. Hoeksema 1984; Y.-M. Wang & N. R. Sheeley 1992). The radial component of the magnetic field (B_r ) at 2.5 R_⊙ is extrapolated to 18 R_⊙ assuming conservation of the magnetic flux (e.g., r² B_r = constant). The radial solar wind velocity (V_r ) is estimated using a modified Wang–Sheeley model (C.-C. Wu et al. 2020a, 2020b): V_r (r = 18 R_⊙) = 150 + 500 fs^−0.4 km s⁻¹. The longitudinal (φ) component of the solar wind velocity is determined with the solar rotation at the equator: V_φ = 33.8 × sin(θ_c ) km s⁻¹, where θ_c is the colatitude and the sine factor approximates the latitudinal distribution. The latitudinal component of the solar wind velocity (V_θ ) is assumed to be negligible (e.g., V_θ = 0). This setup is the same as those used in the ENLIL model (e.g., D. Odstrcil et al. 2020), except the meridional flow speed is set to zero. The solar wind density is estimated assuming mass conservation (e.g., ρ = (ρ_o V_o)/V, where ρ_o = 2.35 × 10⁻⁹ kg m⁻³ and V_o take the average value of V_r at 18 R_⊙). The longitudinal and latitudinal components of the magnetic field are derived assuming a corotating field: B_θ = B_rV_θ /V_r = 0 and B_φ = B_rV_φ /V_r , respectively. Finally, the solar wind temperature is determined assuming conservation of energy, pressure balance, e.g., 2ρRT + ρV²/2 + ∣ B ∣²/2μ₀ + ρg(r–R_⊙) = constant, where R (=8.314 J K⁻¹ mol⁻¹) is the gas constant and g (=274 m s⁻²) is the solar gravitational acceleration constant. At the solar surface, the magnetic intensity (B_o) is a global average value and temperature is (T_o) = 1.5 × 10⁶ K.

To simulate the solar wind and its propagation from the inner boundary (r = 18 R_⊙) to 1 au for one Carrington period (∼27 days) or less, it requires three synoptic maps for the inner boundary input. The extra maps, one before and one after the current map, remove the discontinuity at the start and end times. This ensures that the simulation reaches the steady state of the solar wind. Stitching of synoptic maps is thus required. Quite often, there is discontinuity between them, which may introduce unrealistic disturbances and results in code failures when the disturbances are too large. Smoothing of the maps in longitude in the stitching area is often employed to reduce the discontinuity. For the present simulation, no smoothing is required as the code runs stable across the boundary of the synoptic maps. Note that only two extra maps are needed for the long-term simulation. A total of 68 (e.g., 2267 − 2200 + 1) synoptic maps are used for this study.

Figure 1 shows the time sequence of synoptic maps for B_r input at 2.5 R_⊙ for the years from 2018 to 2023. The x-axis represents the time (in years) from 2018 to 2023. The y-axis represents the solar latitude (in degrees). Every fifth Carrington rotation number is marked on the top of Figure 1. The data resolution is 5° in both longitude and latitude direction for each Carrington map. Each Carrington map lasts for about 27.2753 days. The color bar describes the value of B_r with a range of −0.41–0.43 Gauss (marked on the bottom). Figure 1 clearly shows several solar cycle variations of the magnetic field. First, the Sun was quiet during 2018–2019 and started becoming more active in the middle of 2020 (CR2231). Second, the solar magnetic field polarity flipped signs in late 2022 (between CR2262 and CR2263). Third, the solar minimum occurred in 2020 since the smallest neutral line (B_r = 0, the white solid trace) tilt occurred between 2020 January and April (CR2227–CR2230). After CR2230, the neutral tilt started increasing and reached the south pole at CR2261.

Zoom In Zoom Out Reset image size

Panels from (a) to (e) in Figure 2 shows selected snapshots of G3DMHD-simulated solar wind radial velocity from 2018 to 2022. In each row, there are four maps corresponding to March, June, September, and December at a plane defined by θ = 25N. The exact date and time for each map is marked on the top of each map. The velocity scale is in the range of 250–700 km s⁻¹ (see the color bar on the bottom of each map) for all maps for easy comparisons. The blue curves on each panel represent the possible location of the HCS. In addition, the 1 au distance from the center of Sun is marked as a full circle in the r–φ plane.

Zoom In Zoom Out Reset image size

As shown in Figure 2, there is a clear difference in the flow speed among these maps. The flow speed is faster in 2018 and 2019 and is slower in 2020 and 2021. These velocity maps serve as the gross feature of the solar wind from the late descending phase of solar cycle 24 to the early ascending phase of the solar cycle 25, with 2020 for the solar minimum year. During the descending phase of the solar cycle, the solar wind in the equatorial region is often characterized by corotating interaction regions, where fast flows originating from the coronal holes meet the slow flows originating from the streamers (L. F. Burlaga et al. 1978). During the ascending phase of the solar cycle, the increasing occurrence of CMEs is expected to modify the solar wind significantly. Since CMEs are not simulated in this study, possible CME effects are not present in Figure 2. While including CMEs in our simulation can be beneficial to the present study, it is extremely difficult to implement all CME events in the present study.

As the main objective of the present study, e.g., comparing the simulation result with in situ solar wind measurements, Figure 3 shows the comparison for hourly resolution of G3DMHD simulation results (red dots) with the hourly PSP (black crosses) solar wind plasma data from the Solar Wind Electrons Alphas and Protons (J. C. Kasper et al. 2016) and magnetic field data from the Electromagnetic Fields Investigation (S. D. Bale et al. 2016) during 2018–2022. Panels from top to bottom in Figure 3 show the comparison for the (a) magnetic field intensity (B in nT), (b) solar wind radial velocity (V_r in km s^–1), (c) density (N_p in cm^–3), and (d) temperature (T_p in K). The orbit of PSP (black, blue, and orange lines are the distance from the center of the Sun to the PSP, latitudes, and longitudes of the PSP positions, respectively) is shown in Figure 3(e). Note that, for the comparison, the simulated solar wind parameters are extracted from the grid where the PSP is located for a given time. The total number of data points for the magnetic field and plasma (N_p , V_r , and T_p ) comparison are 28,334 and 9310, respectively. There are large data gaps in N_p , V_r , and T_p and these data gaps occur when PSP is near its aphelion, i.e., close to 1 au.

Zoom In Zoom Out Reset image size

Qualitatively, the modeled solar wind parameters compare reasonably well with the PSP observations. The values of the Pearson correlation coefficient (cc = $\sum \left({x}_{i}-\bar{x}\right)\left({y}_{i}-\bar{y}\right)/\sqrt{\sum {\left({x}_{i}-\bar{x}\right)}^{2}}\sqrt{\sum {\left({y}_{i}-\bar{y}\right)}^{2}}$, where x_i and y_i are the individual sample points and $\bar{x}$ and $\bar{y}$ are the sample mean of variables x and y, respectively) are calculated and listed in Table 1. The Pearson cc, which ranges between −1 and 1, measures the direction and strength of the tendency of two linearly dependent variables. Therefore, a value closer to 1 means better correlation between the modeled and observed variables. The correlation in B and N_p is especially good. The values of the cc for the entire data is 0.75 for B and 0.73 for N_p . However, the correlation for V_r and T_p is poor, cc = 0.15 and 0.28, respectively. Since correlation is not a measure of the accuracy of forecasts but a measure of the trend between two variables, we also calculate the mean absolute scaled error (MASE = $\sum \left|{y}_{i}-\widehat{{y}_{i}}\right|/\sum \left|{y}_{i}-\bar{y}\right|$, where y is the actual value, is the predicted value, and $\bar{y}$ is the mean of the data; R. J. Hyndman & A. B. Koehler 2006) to represent the forecase accuracy, and the results are listed in Table 1. MASE measures the mean of the absolute error of model prediction against the mean of the absolute error from the mean of all data. Therefore, a value of MASE smaller than 1 indicates a better prediction against the mean of the data. Based on this metric, the model performs best in predicting the density (MASE = 0.5), followed by the magnetic field intensity (MASE = 0.6), temperature (MASE = 1.1), and radial velocity (MASE = 1.3). Note that although the correlation for the magnetic field intensity is the highest among the four parameters, there is a clear negative offset or underestimate in the predicted values. In addition, the offset is not random but has a clear trend. It is the largest in 2018 and decreases toward 2022 (see Table 1). For the solar wind density, the correlation is the best in 2018 (cc = 0.86) and becomes slightly weaker toward 2022 (cc = 0.67). There is no systematic trend in the values of MASE, which varies in the range of 0.4 and 0.6. For solar wind radial velocity and temperature, there is no systematic trend in the yearly values cc and MASE. There are years (2019–2021) in which correlation for V_r is poor.

It is well known that data smoothing can affect (improve) the values of cc. In this study, we use the hourly solar wind data obtained from the Space Physics Data Facility at the NASA/Goddard Space Flight Center to calculate cc and MASE without smoothing the data. In addition, there is a concern about the data gaps in the calculated cc and MASE. The presence of data gaps in the PSP data reduces the "within-series dependence" and thus improves the quality of the calculated cc. For calculating MASE, there are two approaches for estimating the naïve forecast. For continuous time series data, the mean absolute error from the (one-step) prior period are used for the naïve forecast. On the other hand, for non-time-series data, the average data value is used as the base predictor (R. J. Hyndman & G. Athanasopoulos 2021). We applied the latter approach as it is more appropriate to the gappy PSP data.

We have performed MHD simulations of the solar wind in the inner heliosphere (18 R_⊙–345 R_⊙ ) for the period of 2018–2022. The simulated results (B, V_r , N_p , and T_p ) are compared with in situ measurements of the solar wind from the PSP to evaluate/validate the G3DMHD model. While the comparison is limited to PSP's orbit, the result shows some qualitative agreement. The major finding is that the model performs better in reconstructing solar wind density and magnetic field intensity than in reconstructing the other two parameters (V_r and T_p ). Solar wind B, N_p , T_p , and V_r are in the ranges of 1–1000 (nT), 5–3000 (cm^–3), 104–106 (K), and 200–800 (km s^–1), respectively. The factor between minimum and maximum are ∼ 1000, 600, 102, and 4 for B, N_p , T_p , and V_r , respectively. While our simulation results of these four solar wind parameters generally match the observations, the large dynamic range of B may be the major cause for better cc for B.

The generally good agreement between modeled and observed B is expected as it is the only measurement available for the boundary conditions. On the other hand, the good agreement in N_p and the poor agreement in V_r are totally not unexpected. The boundary value for N_p is not measured but derived based on the conservation of mass, and an empirical value for the initial mass flux is assigned at the solar surface. Therefore, the value of N_p is determined by the value of V_r at the inner boundary. Thus, an inaccurate V_r can lead to an inaccurate N_p . In other words, an accurate N_p implies an accurate V_r . This is obviously not consistent with the present finding, and we do not have an explanation for this result.

The input for V_r is based on the modified Wang–Sheeley (WS) model (C.-C. Wu et al. 2020a, 2020b), which is tuned to match in situ observations at 1 au using 3D MHD simulations. The poor agreement in V_r may suggest a defect in the empirical formula for the solar wind speed. This defect may be due to the oversimplified formula. Other more sophisticated empirical solar wind speed models such as the Wang–Sheeley–Arge (WSA; C. N. Arge et al. 2003) and the coronal hole boundary (DHCB; P. Riley et al. 2001) models take into account the shortest distance from the open–closed boundary. Some previous works seem to suggest that both WSA and DHCB perform better than the WS model (e.g., P. Riley et al. 2015). However, C.-C. Wu et al. (2020a) demonstrated that the modified WS model can outperform the WSA model for solar wind at 1 au during solar minimum. In addition, the modified WS model is derived based on data acquired in one Carrington rotation (C.-C. Wu et al. 2020b) and may not be suitable for other solar activity phases. Other models have considered solar activity phases by using different coefficients that are used in their empirical solar wind speed model (e.g., P. Riley et al. 2015). If the phase of solar cycle is important, one should see a worse result for nonsolar minimum years. Indeed, Table 1 clearly shows such a result but only moderately. Obviously, one cannot rule out that the effect of CMEs, which are not modeled in the present study. However, we believe that such an effect should be small as we are able to model B and N reasonably well without including CME events. The solar wind temperature for the boundary is also not measured but derived from the assumption of conservation of energy. This formula depends on the solar wind speed. Thus a poor modeled T_p may be attributed to a poor specification of the solar wind speed at the boundary. However, we cannot rule out the oversimplified solar wind temperature at the solar surface as the cause. A single temperature value is clearly not realistic. Works on improving the inner boundary conditions for the solar wind parameters (e.g., V_r and T_p ) mentioned above will certainly improve G3DMHD and benefit the community.

Another finding worthy of mention is that the performance of our simulations as indicated by the cc and MASE metrics does not show clear systematic changes over the phase of a solar cycle. This is contrary to our intuition. Our MHD simulations are designed to reconstruct the background solar wind without transit events such as CMEs and shocks. Therefore, it is expected that the best performance would occur in the minimum year of the solar cycle (i.e., 2019 December ³ ) because the occurrence of CME events is least frequent. While both cc and MASE values for N_p and V_r in 2018 suggest a better agreement than other periods, we do not see such an effect in B and T_p . It is reasonable to argue that a better agreement could be achieved if CME events were excluded in the PSP data and not compared. Unfortunately, classification of the PSP data into CME and non-CME categories is not a trivial task and will not be addressed here. In addition, there are two possible factors that affect the comparison result. For example, the ideal MHD model, thus our simulation, is not perfect. PSP often observed large Alfvénic waves and turbulence closer to the Sun. This finding has led to some workers to incorporate MHD turbulence to their simulation models (e.g., F. Fraternale et al. 2022). Another possible factor is that the assigned boundary values are not correct. After all, only one parameter (B in the photosphere) is measured, and others are based on some conservation laws. We believe this is a major factor that will be explored in our future work.

In this section, we intend to address the issue of the underestimate of the total magnetic field. According to Figure 3(a), the underestimate decreases with time, such that it is largest in 2018 and smallest in 2022. This is in good agreement with some previous results (e.g., Y.-M. Wang & N. R. Sheeley 1988; Riley et al. 2015; L. K. Jian et al. 2015; J. A. Linker et al. 2017; S. Wallace et al. 2019; S. T. Badman et al. 2021) and is referred to as the open flux problem by J. A. Linker et al. (2017) although these studies are based on 1 au data. To further demonstrate the systematic yearly offset in B, additional simulation runs are performed by increasing the inputs of the B_r derived from the PFSS model by 2, 3, and 4 times at the inner boundary of the simulation domain. Quantitatively, there is little change in the cc for B, N_p , and T_p . Although the correlation for V_r becomes slightly weaker for increasing B, all coefficient values are too small to be correlated. The MASE values technically remain the same for the three MHD parameters except for the magnetic field intensity.

We may determine the scaling factor for B that best fits the data on a year-by-year basis, and since the magnetic field offset is scalable, one can improve reconstruction of B by scaling up the PFSS magnetic field intensity input at the inner boundary. The smallest MASE value for B are 0.54, 0.67, 0.57, 0.56, and 0.43 for the year 2018 (4B), 2019 (3B and 4B), 2020 (3B), 2021 (2B), and 2022 (1B), respectively. To avoid the discontinuity between years, a continuous scaling factor, e.g., time-dependent scaling factor, is obtained by fitting a straight line to these yearly scaling factors. The time-dependent scaling factor (or scaling function) is expressed as SB = 4.63–0.767 × (f_yr − 2018), where f_yr is the fractional year. A simulation run using the scaling factor SB is performed, and the result is shown in Figure 4. The results of cc and MASE are listed in Table 2. It is shown that the agreement between the observed and simulated B improves significantly. While the rescaling of B does not improve the correlation, the overall agreement improves significantly, with the value of MASE reducing to 0.54. While this rescaling of B improves the performance of our simulation result, one has to be cautious about the use of the scaling function as it has not been tested with more data. For example, the time-dependent scaling factor is derived from data acquired during the ascending phase of current solar cycle. It is clearly not useful for the descending phase as the expected trend in B might require a switch sign for the fitting slope. In conclusion, this type of B rescaling has no physical ground other than a demonstration that B is scalable.

Zoom In Zoom Out Reset image size

The systematic IMF B offset suggests the source of properties varying with the solar cycle. An underestimate of B predicted by the PFSS model has been reported previously. Y.-M. Wang & N. R. Sheeley (1988) extrapolated magnetic field strength (27 day average) from the WSO potential field to 1 au and compared it with in situ observations for the period of 1976–1986. They found a factor of ∼ 1–2 difference in B when radius of source surface (R_SS) is either set to 1.6 or 2.5 R_⊙. This work is extended to include data from 1968 to 2022 and eight different synoptic maps recently (C. C. Wang et al. 2022). They concluded that the total open fluxes, extrapolated with R_SS = 2.5 R_⊙, underestimate the observed radial component of B by factors of ∼2–5 for most synoptic maps.

According to the magnetic field observations from the Ulysses spacecraft, the total open flux of B_r is distributed isotropically at 1 au (A. Balogh et al. 1995; E. J. Smith & A. Balogh 2008). The isotropization of the flux could have occurred within r ∼ 10–15 R_⊙ due to the zonal flows (e.g., Y.-M. Wang 1996; X. P. Zhao & J. T. Hoeksema 2010; O. Cohen 2015). Therefore, for each map, we average the magnitude of B_r along a fixed longitude and apply the mean value to that longitude at the inner boundary (18 R_⊙). This procedure is carried out for all synoptic maps for the entire studied period (e.g., 2018–2022) with no change in the V_r and N_p (note that T_p changes with B). Figure 5 shows the simulation result using the latitudinally averaged magnetic maps, and the ccs and MASE values for B, N_p , V_r , and T_p are listed in Table 3. The correlation for B is improved significantly (from cc = 0.75 to cc = 0.83), and the MASE value is reduced by 20% (from 0.60 to 0.49) comparing with simulation results without any modification for B. Note that the extremely good correlation (cc = 0.93) for B for the year of 2021 and the value of MASE drops from 0.68 to 0.43 for the year of 2021. Therefore, using constant ∣B_r ∣ at the inner boundary is a good approach for realistic solar wind simulation.

Zoom In Zoom Out Reset image size

Since PSP orbits cover a wide range of the radial distance from the Sun, it would be important to explore if the performance of G3DMHD is radial dependent. We divided PSP's radial component orbit into three regions: r < 0.3 au, 0.3 au < r < 0.6 au, and 0.6 au < r < 0.9 au. The calculated cc and MASE for these three regions are listed in Table 4. In general, there is no significant change in the performance G3DMHD in reconstructing the solar wind parameters, e.g., the modeled B and N_p match the observed B and N_p than V_r and T_p do. With regard to the radial dependence, the table indicates that when there is a good correlation between the modeled and observed parameters (B and N_p ), the performance of the G3DMHD model reduces with increasing radial distances from the Sun. This result is not surprising because a better agreement between the modeled and observed results closer to the inner boundary is expected if the boundary conditions are correctly prescribed. The values of cc for B and N_p are moderate (0.76 and 0.62, respectively), whereas the values of cc for V_r and T_p are very small (0.09 and −0.05, respectively). Similarly, the values of MASE for B and N_p are only moderate (0.67 and 0.76), whereas the values of cc for V_r and T_p are greater than 1 (1.77 and 1.89). This may suggest that none of the boundary conditions are correctly specified, especially the V_r and T_p parameters. It is also possible that the G3DMHD is not perfect for solar wind modeling. We tend to believe that the former is the major cause because the G3DMHD has been shown to be able to reconstruct the solar wind in some cases (e.g., C.-C. Wu et al. 2016, 2019, 2022).

Note that Table 4 also indicates that the performance of G3DMHD reduces significantly after r > 0.6 (cc = 0.18 and MASE = 1.18 for B and cc = 0.17 and MASE = 0.98 for N_p ). We will show later that the performance of the G3DMHD model also depends on the longitude of observations from the Sun–Earth line.

While the entire photospheric magnetic field is measured "instantaneously," the making of synoptic maps is not. The synoptic map combines the latitudinal stripe of daily magnetic field data near the Sun–Earth line to construct the synoptic maps by assuming the solar magnetic field does not change significantly from one solar rotation to another. Therefore, the reliability of the synoptic maps as the boundary condition for B may be in question because the solar magnetic field is not static. As a result, the performance of G3DMHD should reveal a longitudinal dependence. Here, we divide the longitudinal domain into 12 sectors with 30° for each sector. The values of cc and MASE are calculated for data falling within each sector and are listed in Table 5. The best performance occurs, as expected, when PSP is within 15° (e.g., 345°–15°) from the Sun–Earth line. The values of cc are in the range of 0.64–0.93 for B and 0.45–0.86 for N_p . There is no systematic reduction in cc with the longitudinal separation between the PSP and Sun–Earth line for these two parameters. A modest improvement in the correlation (cc = 0.46) in the modeled solar wind speed occurs near the Sun–Earth line. The result for T_p is consistently poor. We noticed that the cc value is 0.62 at 75°–105°. However, this improvement comes with a large MASE value (3.04). Therefore, we do not think this is a good match.

The present study demonstrated that G3DMHD is better suited for reconstruction of the solar wind density and magnetic field intensity. However, it performs poorly in simulating the V_r and T_p . It is reasonable to attribute the poor agreement to lack of a realistic input from measurements. All current data-driven MHD simulation models are subject to the same vexing drawback. In addition, the poor agreement may be attributed partially to the lack of CME simulation. However, the good agreement in the modeled B and N _p with observations and lack of clear solar cycle variation may suggest that effects of CME might be small. Therefore, we believe that the boundary condition must be the major cause. This suggestion is also reinforced by the result of radial and longitudinal analysis as shown in Section 4.3. As part of our future efforts, we will point out a few possible directions to improve the G3DMHD model. First, the simple model of Y.-M. Wang & N. R. Sheeley (1990a) is used for specifying the solar wind speed at the inner boundary. Although the model has been modified using MHD simulations (C.-C. Wu et al. 2020a, 2020b), it is for solar quiet times only. This solar wind speed model may be improved by considering other solar conditions. Second, C. N. Arge (2003) proposed an improved model (known as the WSA model) of the solar wind that incorporates the angular distance of the field line from the nearest coronal-hole boundary. The WSA-ENLIL model has been adopted by CCMC and has been used in many studies. We will incorporate the WSA model into G3DMHD to test its performance. Third, the solar wind speed and temperature are usually well correlated at 1 au (M. Neugebauer & C. W. Snyder 1966; I. B. Strong et al. 1966; L. F. Burlaga & K. W. Ogilvie 1970; R. E. Lopez 1987; H. A. Elliott et al. 2010). However, no study was performed sunward of 1 au. We anticipate that the agreement in T will improve once the agreement in V is improved. On the other hand, the PSP data can be used to test if the T–V relationship still holds close to the Sun. This will allow us to build an empirical T–V relationship for MHD modeling of the solar wind.

As a final note, the intent of this study is to evaluate the performance of the G3DMHD model for the background solar wind. No solar disturbance (e.g., CME) has been applied to the inner boundary. This could be one of the reasons that contribute to some of the poor results. Introducing CME events into solar wind simulation has been performed in some single-event studies (e.g., C.-C. Wu et al. 2016, 2020a, 2022). Modeling multiple CMEs is still rare (e.g., C.-C. Wu et al. 2022). Some studies model the CME with a pressure pulse (e.g., B. E. Wood et al. 2011, 2012; C.-C. Wu et al. 2016, 2020a) and some with a flux rope (e.g., W. B. Manchester et al. 2008; F. Shen et al. 2012) or two flux ropes (G. J. Koehn et al. 2022). Either approach is still difficult currently without mentioning that it is impractical to include all CME events because identifying CME events also constitutes another difficulty. While the present study suggests that some improvements are required, it encouraged us to implement flux ropes into the G3DMHD model (e.g., C.-C. Wu et al. 2024b). Once the flux-rope model + G3DMHD model is able to simulate flux-roped CMEs and its driven shock correctly, we will try to add all available CMEs into the inner boundary and rerun the yearly simulation in the future.