Application of a Modified Spheromak Model to Simulations of Coronal Mass Ejection in the Inner Heliosphere

The magnetic fields of interplanetary coronal mass ejections (ICMEs), which originate close to the Sun in the form of a flux rope, determine their geoeffectiveness. Therefore, robust flux rope-based models of CMEs are required to perform magnetohydrodynamic (MHD) simulations aimed at space weather predictions. We propose a modified spheromak model and demonstrate its applicability to CME simulations. In this model, such properties of a simulated CME as the poloidal and toroidal magnetic fluxes, and the helicity sign can be controlled with a set of input parameters. We propose a robust technique for introducing CMEs with an appropriate speed into a background, MHD solution describing the solar wind in the inner heliosphere. Through a parametric study, we find that the speed of a CME is much more dependent on its poloidal flux than on the toroidal flux. We also show that the CME speed increases with its total energy, giving us control over its initial speed. We further demonstrate the applicability of this model to simulations of CME-CME collisions. Finally, we use this model to simulate the 12 July 2012 CME and compare the plasma properties at 1 AU with observations. The predicted CME properties agree reasonably with observational data.


Introduction
Coronal Mass Ejections (CMEs) are one of the most explosive events in our solar system, with the kinetic energy release of up to 10 26 J [Forbes, 2000, Vourlidas et al., 2002. During a CME event, the coronal plasma erupts producing a coherent magnetic structure called the flux rope. The standard model of CME eruption states that the pre-eruptive state possesses a core of axial flux that may be either a sheared arcade or a flux rope which exists above the magnetic polarity inversion line [Carmichael, 1964, Sturrock, 1968, Hirayama, 1974, Kopp and Pneuman, 1976. As this flux rope starts to rise due to various possible magnetic force imbalances, the overlying magnetic field lines serve as the dominant source of the poloidal flux in a flux rope [see the review of Chen, 2011]. As a result, a CME carrying mass between 10 14 to 4 × 10 16 gms, starts its journey through the heliosphere with speeds ranging between 10 to >2000 km/s [Hudson et al., 2006]. Most of the time, the CME speed exceeds the fast magnetosonic wave speed in the ambient solar wind, resulting in a shock in front of this CME [Sime andHundhausen, 1987, Raymond et al., 2000].
If a CME impacts Earth, it disturbs the local space environment posing a wide range of risks including harmful effects on space assets, radiation exposure for astronauts and passengers during polar flights, communication losses between satellites and ground receivers, and inducing large currents in long power transmission lines. The ability of a CME to disturb the near-Earth space, also called its geo-effectiveness, is largely due to CME's magnetic field which often reveals itself as a flux rope with a sheath region in front of it. CMEs observed to have flux ropes associated with magnetic clouds at 1 AU form a subset of the total CME variety . If the orientation of the magnetic field in the CME is favorable for magnetic reconnection with the Earth's magnetosphere, such a CME is highly likely to be geoeffective [Burton et al., 1975, Gonzalez et al., 1994. This orientation should involve a negative B z component, since B z is positive in the Earth's day-side magnetosphere.
To mitigate problems caused by geoeffective CMEs with proper precautionary steps, we need to predict accurately the arrival time and magnetic structure of CMEs at 1 AU. Due to the development of supercomputers and highly parallelized codes, magnetohydrodynamic (MHD) simulations of CMEs have become a major tool for CME predictions. A number of authors have shown the application of MHD modeling to CME simulations [e.g. Vandas et al., 1996, Manchester IV et al., 2004, Lugaz et al., 2005, Aulanier et al., 2010, Shen et al., 2014, Titov et al., 2014, Jiang et al., 2016, Shiota and Kataoka, 2016, Pogorelov et al., 2017, Jin et al., 2017, Singh et al., 2018, Scolini, C. et al., 2019]. The cone model of a CME [e.g. Chané et al., 2005, Odstrčil andPizzo, 1999] is already being used for operations by various agencies. This model can predict the CME arrival time, shock and sheath regions while its magnetic structure remains unknown because of the lack of a flux rope treatment. This major disadvantage of the cone model can be addressed by using fluxrope-based models. Many flux rope models have been proposed so far for CME simulations, e.g., 1) the spheromak model [Lites et al., 1995], 2) the Gibson-Low model , 3) the Titov-Demoulin model [Titov and Démoulin, 1999], etc.
For a successful prediction of CMEs, the input parameters in a flux rope model must be constrained by observations as much as possible. Keeping that in mind, Singh et al. [2020] proposed a modified spheromak model based on the observed poloidal flux, toroidal flux, and helicity sign. The parameters of this flux rope can be adjusted so that it erupts with the observed speed, orientation, and direction. In this paper, we show the applicability of this model to simulate CMEs in the inner heliosphere, where the inner boundary of the computation domain is at 0.1 AU. It is important that model CMEs should propagate through a realistic data-driven background solar wind (SW) because the CME-SW interaction plays an important role during the CME propagation. We use the inner-heliospheric background created on the basis of the Wang-Sheeley-Arge (WSA) coronal model as inner boundary conditions [Arge and Pizzo, 2000, Arge et al., 2003, Arge et al., 2005.
Here, we perform a parametric study to analyze how a CME possessing the modified spheromak flux-rope structure evolves through the inner heliosphere. We show the effect of varying poloidal and toroidal fluxes, initial flux rope size, and total energy on the CME dynamics. We also show the applicability of our model to simulations of CME-CME interactions, which are very common and important during the solar activity maxima. We also show the simulation results for the 12 July 2012 CME and compare simulation results with observations at 1 AU.
In Sec. 2, we describe the SW and flux rope models used in this study. We also describe how the flux rope model is inserted into the solar wind background. In Sec. 3, we discuss the results of various CME runs performed in this study. In particular, the CME-CME collision results are presented in this section. We also show the results for the 12 July 2012 CME simulation. Finally, we give our concluding remarks in Sec. 4.

Models
All our simulations are carried out using the Multi-Scale Fluid-Kinetic Simulation Suite (MS-FLUKSS), a highly paralleled collection of codes capable of performing AMR finite-volume MHD simulations of the SW in the presence of neutral atoms and physical processes,e.g. turbulence, induced by charge exchange on adaptive mesh refinement grids . In the following subsections, we describe our SW and flux-rope CME models, and the method of CME insertion into the SW background.

Inner Heliosphere Model
To generate realistic time-varying background solar wind in the inner heliosphere, we couple our heliospheric MHD model with the WSA coronal model at 0.1 AU [e.g., Kim et al., 2019]. As input to the WSA model at 1 R , we use the Air Force Data Assimilative Photospheric Flux Transport (ADAPT) synchronic maps generated from the NSO/GONG magnetograms , Arge et al., 2011, Hickmann et al., 2015. The semiempirical WSA model initially extrapolates the photospheric magnetic field to a source surface, chosen to be at 2.5 R , using the potential field source surface (PFSS) model, and then to the outer boundary at 21.5 R (0.1 AU) using the Schatten current sheet model. In addition to magnetic field strengths, the WSA model provides an estimate of the solar wind speed at the model outer boundary as a function of the flux expansion factor and distance to the nearest coronal hole boundary [Arge et al., 2003, Arge et al., 2005.
While the ADAPT-WSA model provides 12 realizations of magnetic field and solar wind speed maps at 0.1 AU, we select one particular realization that provides us with the best fit to the near-Earth spacecraft data. On interpolating the WSA maps from the original 144 × 72 (2.5 • × 2.5 • ) resolution to a base grid of 256 × 128 (∼ 1.4 • × 1.4 • ) in MS-FLUKSS, we scale the WSA magnetic field strengths by a factor of 3 to compensate for the systematic underestimation of the open magnetic flux at 1 AU [Linker et al., 2016, Wallace et al., 2019. We estimate the radial and azimuthal components of magnetic field at 21.5 r using the local solar wind speed to account for the solar rotation as the solar wind propagates radially outward from 1 to 21.5 r [MacNeice et al., 2011]. Also, we reduce the WSA speeds by 20% to account for the difference in solar wind acceleration between the WSA and MS-FLUKSS models [e.g., MacNeice et al., 2011. To estimate the solar wind density (and temperature) at 0.1 AU, we use empirical correlations between the solar wind speed and density (and temperature) based on OMNI data [Elliott et al., 2016]. With these boundary conditions, we solve the ideal MHD equations on the nonuniform 150 × 256 × 128 (r, φ, θ) fully spherical grid to simulate the solar wind outflow from 0.1 to 1.5 AU. The resulting solar wind solution in the equatorial plane is shown in Fig. 1. Here, we have colored the equatorial plane with radial solar wind speed. We have also shown the magnetic field lines in this plane. Singh et al. [2020] describe the details of the modified spheromak model used in this study. The original spheromak solution was modified so that the poloidal and toroidal fluxes can be set up independently. The analytic solution for the magnetic field inside the flux rope can be given as

Modified spheromak model
where α 0 and r 0 are related as α 0 r 0 = 5.763459, which is the first root of the Bessel function J 5/2 . Here, r 0 is the spheromak radius. The toroidal flux in the spheromak is proportional to the parameter a 1 , whereas the poloidal flux is proportional to a 1 γ. The origin of the spherical coordinate system (r, θ, φ) is placed at the spheromak center. The parameter δ can assume values +1 or -1 for either positive and negative helicity, respectively. The plasma density is distributed uniformly in the spheromak. The thermal pressure is assumed to be proportional to the magnetic pressure through the choice of desired plasma β. Our approach to do this will be discussed in the next subsection. As suggested by , this model can also undergo a stretching operation r → r − a, where a is the stretching parameter to modify the spheromak shape from spherical to a tear-drop. This stretching operation does not violate the solenoidal condition for magnetic field.

Introducing flux rope in solar wind
We propose to introduce a flux rope into the SW in such a way that initially the flux rope is superimposed with the background. This approach is different from the one commonly used, where a flux rope is introduced gradually at the inner boundary [eg. Shiota andKataoka, 2016, Scolini, C. et al., 2019]. We are not introducing the full spheromak into the SW, but rather the top half of it, which resembles the CME flux rope much more accurately, exhibiting a curved front and two legs. We introduce the flux rope parameters as follows: Here, b F R is given by Eq.
(1), ρ is the plasma density, and e is the total energy density. We also note that where p and v are the thermal pressure and bulk velocity, respectively. We have kept the isotropic index γ = 1.5 in this study. We also introduce an energy multiplier factor ξ here, which can be shown to be related to plasma β as β ≈ (γ − 1)(ξ − 1) under a quite realistic assumption that the magnetic pressure created by the flux rope is much greater than the ambient solar wind pressure. Thus we can control the plasma β, as well as the eruption speed of the flux rope, using ξ, as will be shown later. During this insertion procedure, we make sure that total thermal pressure in the domain does not fall below 25% of the background SW thermal pressure [Manchester IV et al., 2004]. Introducing flux rope by modifying the magnetic field in the domain will break the divergence free condition at the flux rope-solar wind interface. This is true even when introducing the flux rope gradually at the inner boundary [Shiota and Kataoka, 2016]. We then rely on divergence cleaning mechanisms like 8-wave approach [Powell et al., 1999], generalized lagrange multiplier approach [Dedner et al., 2002], etc. to convect away and damp the non-zero magnetic divergence.
The result of the spheromak insertion into the SW is shown in Fig. 2. The flux rope is kept in the equatorial plane without any tilt. The introduced poloidal flux is 2 × 10 22 Mx and the toroidal flux is 5 × 10 21 Mx. The helicity sign is set to be positive and the mass of 1.65 × 10 15 g is distributed uniformly inside it. The red sphere shows the inner boundary placed at 21.5 r . The flux rope size parameter r 0 and stretching parameter a are set to be 15 r and 5 r , respectively. The energy multiplier ξ has been set to 2. The flux rope center is kept at the inner boundary. We find that keeping the flux rope center near the inner boundary creates a configuration suitable for eruption. There is a region of high total magnetic pressure inside the flux rope, which can result in its eruption. We notice that for ξ = 2, plasma β is around 0.5, which satisfies our relation β ≈ (γ − 1)(ξ − 1). This value of γ means that magnetic pressure will be the dominant driving force during CME eruption. To visualize the force distribution in this flux rope, we calculated the total force (per unit volume) due to the magnetic pressure gradient, thermal pressure gradient and magnetic tension as shown below: On solving these equations, we get The first term on the RHS is due to the magnetic pressure gradient, similarly to the hoop force described by Kliem and Török [2006] for torus instability. The middle term is due to the magnetic tension force, whereas the last term is due to thermal pressure gradient. The radial component of force F is shown in the bottom right panel of Fig. 2. One can see a large radial force underneath the FR, that makes our model CME to erupt. We will show the erupting stages in the next section.

Results
In this study, we solve the ideal MHD equations in an inertial coordinate system centered at the Sun, whose sidereal rotation period is 25.38 days. We use a TVD, finite-volume Roe scheme with arithmetic averaging to compute the numerical fluxes and the forward Euler scheme for time integration. To remove numerically-induced magnetic field divergence, we use the Powell et al. [1999] approach. We use the WSA solution for the period between 01-June-2012 and 31-Aug-2012 as rotating boundary conditions at the inner sphere to create the time-dependent inner-heliospheric background. We use the spherical domain of size 150 × 128 × 256 in r, θ, and φ directions, respectively. The modified spheromak flux rope, when inserted into the solar wind background as described in Sec. 2.3, erupts immediately due to the outwards forces. The eruption of a flux rope shown in Fig. 2 is demonstrated in Fig. 3. This flux rope was inserted into the domain on 12-July-2012 21:30 UT physical time. We can see both the flux rope expansion and formation of a shock in front of the CME, marked by temperature enhancement. We notice that the shock evolution deviates from the initial nearly spherical shape due to the nonuniform solar wind background. The shock front, as well as the flux rope, protrude along the region of high-speed solar wind seen clearly in this direction (Fig. 1). The shock also changes the direction of the solar wind magnetic field, downstream to it. The speed of the flux rope in the propagation direction was found to be 1172 Km/s. The shock speed, when it reached 50 R height, was 1309 Km/s. We find these values by fitting a quadratic function to the height-time profile up to 70 R .
It should be noted that the initially inserted flux rope detaches from the inner boundary and its legs experience reconnection with the background solar wind magnetic field. Since we introduce the flux rope inside the computational domain, the boundary conditions, which reside in the ghost cells, are not affected by it. Moreover, the solar wind at the inner boundary is already moving faster than the local fast magnetosonic speed. Therefore, no MHD waves can travel back towards the inner boundary.
In Fig. 4, we show the 1 AU simulation data for this CME. The probe, shown by brown circle, is kept at 1 AU in the direction of CME insertion in the equitorial plane, and the solution is probed every hour at this point. In Fig. 4, we can clearly see the passage of shock and the flux rope at the probe. We have marked the shock along with its sheath and the flux rope region with blue and yellow color, respectively. The shock and its sheath reveals themselves as abrupt enhancement in N p and V r and compression of magnetic field. To identify the end of the sheath region and start of the flux rope, we use the probed B N values. In the equatorial plane, our inner heliosphere background model gives B N ≈ 0. Therefore, the smooth change in the B N must be due to the flux rope. This region is also marked by a decrease in plasma density, also commonly observed in in-situ data.

Parametric study
To understand the dependence of CME evolution on its input parameters better, we performed a parametric study, in which we vary some governing parameters keeping the rest of them constant and determining their impact on speed and acceleration of the erupting CME. We simulate 12 cases. The simulation results are given in Table 1. We also probe the CME properties at 1 AU for these cases. We focus on the following parameters in this study: 1. Poloidal flux 2. Toroidal flux 3. Initial size of flux rope 4. The energy multiplier ξ The results are also discussed below in the corresponding subsections.

Relationship between CME speed and acceleration
Observations show that the background solar wind accelerates slow CMEs and decelerates fast ones. For example,  show that the average acceleration of The red line shows the probed data before the CME was launched. The light blue region represents the sheath region, marked by jumps in speed and density. B R and B T are also enhanced in this region due to the shock compression. The yellow region marks the flux rope behind the sheath region, identified by non-zero and rotating B N values. This region is also characterized by very low density, in agreement with observations. The horizontal axis represents time since CME insertion.
CMEs between Solar and Heliospheric Observatory (SOHO) coronagraph field of view and 1 AU is positive for CMEs with initial speed less than 405 km/s and negative for CMEs with initial speed greater than 405 km/s, which is almost equal to the average solar wind speed. To identify this effect in our study, we find the flux rope and shock speeds, and acceleration at 50 R for all cases. This was done by fitting a quadratic function to heighttime data. Figure 5 shows the inverse dependence of the flux-rope and shock acceleration on the CME speed, in agreement with observational data. The speed at which the acceleration changes from positive to negative is around 1000 km/s, which is higher than the ambient solar wind speed. This is stipulated by the method we use for the CME insertion. It results in large acceleration initially, when a CME is just inserted into the solar wind. This is why, acceleration dominates over the solar wind drag at 50 R , where these data are taken.  Figure 5: Dependence of CME acceleration on its speed at 50 R for different runs tabulated in Table 1. Fast CMEs have negative acceleration, whereas the slower ones have positive acceleration, with a clear negative correlation between speed and acceleration.

Effect of Poloidal flux on CME evolution
Keeping r 0 = 15R , r 1 = 26R , a = 5R , ξ = 2, and the toroidal flux = 5 × 10 21 Mx, we vary the poloidal flux within 10-25 ×10 21 Mx. The resulting flux-rope and shock speeds are shown in Fig. 6. We see a strong linear dependence of the CME speed on the poloidal flux.
A similar result was reported by Singh et al. [2019] and Jin et al. [2017] using the Gibson-Low flux rope model. This relationship is observed because by increasing the poloidal flux in a flux rope, we are increasing the magnetic pressure and tension forces responsible for the eruption. We will discuss this point in more detail in the next subsection. A positive correlation between CME speeds and corresponding poloidal fluxes was demonstrated earlier in observations by Gopalswamy et al. [2018]. In Fig. 7, we show the properties of the simulated CME at 1 AU in the equatorial plane and in the direction of CME launch. Here we show the CMEs with poloidal flux of 15 × 10 21 , 20 × 10 21 , and 25 × 10 21 Mx, along with the case when no CME was launched. The magnetic field is given in RTN coordinates, where the B N values are approximately equal to the B z values in the geocentric solar ecliptic (GSE) coordinates. The shock arrival is characterized by a sudden increase in density and radial speed. The flux-rope arrival reveals itself by a large change in B N . The changes in B R and B T , which occur before the change in B N , are due to the shock compression. The results at 1 AU show that CMEs with higher poloidal flux arrive early because they are launched at higher speeds. We also see that the passage of a flux rope is accompanied by a drop in density, as often observed in in situ data [Burlaga et al., 1981].  Figure 6: Dependence of the simulated CME speed at 50 R on the corresponding input poloidal flux. We find a linear trend relating these two properties.

Effect of toroidal flux on the CME evolution
We vary the toroidal flux in the flux rope in the range of 5-12 ×10 21 Mx, while keeping r 0 = 15R , r 1 = 26R , a = 5R , ξ = 2, and poloidal flux = 20 × 10 21 Mx. The effect of this on CME speed is plotted in Fig. 9. Here we see a stark difference in the trend as compared with that of the variation of the poloidal flux. The CME speed is much less dependent on the input toroidal flux. This phenomenon can be understood by looking at the dependence of force distribution inside the flux rope on the varying poloidal and toroidal fluxes. In Fig.  8, we show the radial force acting on the flux rope at the time of its insertion as calculated by Eqn. 4. The three cases shown have the poloidal and toroidal fluxes equal to 10 and 5, 20 and 5, and 20 and 10 respectively. The units used are 10 21 Mx. We show that on doubling the poloidal flux from 1 × 10 22 Mx to 2 × 10 22 Mx while keeping toroidal flux equal to 5 × 10 21 Mx, the force increases considerably. However, when we double the toroidal flux from 5 × 10 21 Mx to 10 × 10 21 Mx while keeping the poloidal flux as 2 × 10 22 Mx, there is very little change in the force distribution. Thus we conclude that in a curved flux rope geometry such used in our model, the poloidal flux is the main contributor of the eruptive force, and hence affects the CME speed stronger than the toroidal flux. Since the observed CMEs also have such curved geometries, this statement applies to them as well [Kliem and Török, 2006]. Observations show that the majority of poloidal flux in a CME is due to the reconnected flux created during the eruption , Longcope et al., 2007. Thus, the major factor that is affecting the CME speed is the magnetic flux which a flux rope gains during its eruption, not the flux it possesses in the pre-eruptive flux rope. Figure 10 shows the probe data at 1 AU for the runs with toroidal flux of 5 × 10 21 , 8 × 10 21 , and 12 × 10 21 Mx. We see that the arrival times of CME shocks do not differ much, but magnetic field magnitudes inside them are very different. This shows that CMEs with the same poloidal flux, speed, orientation, and direction can carry different magnetic field at 1 AU, and hence have different geo-effectiveness. This means that the toroidal magnetic flux in pre-eruptive flux ropes may not affect the CME speed, but it does impact its features at 1 AU.
To demonstrate this point even further, we made two more simulations with the toroidal flux of 5 × 10 21 Mx and 10 × 10 21 Mx, the orientation tilted by 90 degrees, and while the rest of parameters are the same. This should create conditions for the CME toroidal flux to contribute to the B N values at 1 AU, in accordance with the self-similar expansion conditions. Although our simulated CMEs do not expand in a self-similar way, they retain the overall coherent structure, so we should expect the toroidal flux to contribute to B N in any event. Figure 11 shows the comparison of these two CMEs at 1 AU. We see that the toroidal flux makes B N negative. It is worth noticing here is that even though the CMEs arrive at roughly the same time and have the same initial orientation and poloidal flux, the B N values are much more negative for the case of 10 × 10 21 Mx toroidal flux. Thus, this CME is more geoeffective. Therefore, if we wish to use our flux rope-based model for CME predictions, constraining poloidal flux alone is not sufficient. We need to ensure a proper magnitude of toroidal flux as well.

Effect of total energy on CME evolution
As discussed in Sec. 2.3, we use a multiplier ξ while adding the flux rope total energy to the simulation domain. We find that this multiplier can control the eruption speed of the flux rope and therefore can be used to constrain it to desired values. To quantify the effect of this parameter, we vary it between the values of 1 and 6 while keeping r 0 = 15R , r 1 = 26R , a = 5R , and the poloidal and toroidal fluxes equal to as 20 × 10 21 Mx and 5 × 10 21 Mx, respectively. The effect of this parameter on the flux rope and shock speeds at 50 R is shown in Fig. 12. We see a linear dependence of the speed on ξ. This relationship can be used to launch CMEs with desirable speed. However, we note that according to the relation β = (γ − 1)(ξ − 1), the choice of ξ > 3 results in β > 1 in our model, which is an unnatural property for a CME. Therefore, ξ = 3 should be considered as the realistic upper limit for our model.  Figure 12: CME speed dependence on the total energy multiplier ξ. The CME speed is calculated at 50 R using a quadratic fit to height-time data.
3.1.5 Effect of the initial size a flux rope on the CME evolution If we want to change the initial size of a flux rope, while keeping the magnetic flux values unchanged, we need to modify the magnetic field strength in this flux rope. In this way, a  Figure 13: 1 AU probe data for different input energy multipliers ξ. The legend shows the energy multiplier used in simulation. The probe is kept in the equatorial plane in the CME launch direction. The horizontal axis represents time since CME insertion.
bigger initial flux rope will result in a smaller magnetic field strength inside it and decrease the total pressure. This means that larger initial flux rope will erupt at a smaller speed. This effect was also shown in the CME simulations in the solar corona using a modified spheromak model [Singh et al., 2020]. We show that this also takes place in our innerheliosphere simulations by keeping r 1 = 26R , a = 5R , poloidal and toroidal fluxes equal to 20 × 10 21 Mx and 5 × 10 21 Mx, respectively, ξ = 2, while varying r 0 in the range of 10-20 R . The results are shown in Fig. 14. We see a reduction in both shock and flux rope speeds with increasing r 0 . The Earth probe data for these runs are given in Fig. 15.

CME-CME collision
During solar maxima, CME-CME collision scenarios are common due to a large number of CMEs erupting from the Sun. Such collisions are possible if a slower CME is in the path of a faster one. The resulting dynamics of the CMEs can range from purely inelastic to superelastic [Shen et al., 2017, and references therein]. A collision can change the CME direction, as well as its speed, resulting in the CME features observed at 1 AU totally different from such with a single CME [Manchester et al., 2017, and references therein]. This makes  Figure 15: Earth probe data for the cases with different input initial size parameter r 0 . The legend shows the values of r 0 used in our simulations, scaled to R . The horizontal axis represents time since CME insertion.
the collision process very relevant for space weather predictions. A CME model should be able to simulate CME-CME interactions for it to be operationally viable. Flux ropes of each CME play an integral part in the collision process. This is why, we expect flux-rope-based models to be suitable for such simulations. Many authors have shown the ability of different flux rope models to simulate CME-CME collision process [eg. Lugaz et al., 2005, Shen et al., 2016, Shiota and Kataoka, 2016.
To show the applicability of the modified spheromak model to simulate CME-CME interaction, we launch 2 CMEs, with case numbers 11 and 12 in Table 1. The slower CME was launched as in case 12 and had shock speed of 1033 km/s at 50 R . When this CME shock reaches 125 R , 19 hrs after the insertion, a faster CME with parameters of case 11 is launched in the same direction. The flux rope of the slower CME is at 116 R at this time. Owing to its larger speed, this CME catches up the slower CME and the collision process begins. In Fig. 16, we show these two CMEs during the collision phase, 33 hrs after the launch of the slower CME. We can see the axial field lines of two flux ropes and a complex structure made by the poloidal field lines. The shock in front of the slower CME is at 180 R , its flux rope being at 172 R at this time. The shock due to the faster CME has already crossed most of the flux rope of the slower CME by this time.
We show the Earth probe data for this collision in Fig. 17, together with the solution involving only first, slower CME. We find that the shock reaches the probe at the same time in both cases, but the flux rope of the slower CME reaches 1 AU faster when no CME is pushing it from behind. We also find that the density jump, as well as the bulk speed in the case of collision, is stronger. In our example, the orientation of two colliding flux ropes was such that the field lines behind the slow CME and ahead of the faster CME had opposite signs in the z direction. This resulted in magnetic reconnection of these field lines, greatly reducing the magnetic field strength in the z direction. This can be seen as a reduction in B N in the case of collision (see Fig. 17, once again confirming the fact that CME collisions can change the geo-effectiveness of CMEs.

12 July 2012 CME
In this section, we investigate the applicability of our model to simulations of actual CMEs. We have chosen the 12 July 2012 CME for this purpose. The choice was made for several reasons. The speed, direction, and orientation of this CME could be reliably obtained from three different viewpoints of the Solar Terrestrial Relations Observatory (STEREO) and Solar and Heliospheric Observatory (SOHO) spacecraft, since STEREO A and B were making an almost right angle to the Sun-Earth line during this time. This also means that the projection effects can be removed from the mass estimates with higher accuracy. The poloidal and toroidal fluxes were calculated also more accurately, since the source of this CME was near the solar disk center as viewed from the Earth and the solar magnetograms, used to determine these values, are most accurate in this region. Singh et al. [2020] have Shock 1 Shock 2 Axial field of FR 1 Axial field of FR 2 Figure 16: Colliding CMEs discussed in Sec. 3.2, 33 hrs after the launch of the slower CME. We show the equatorial plane colored by plasma temperature along with the magnetic field lines. The axial field lines belonging to both CMEs can be seen. The poloidal field lines become very complicated, owing to the magnetic reconnection between them. The shock of the faster CME, labeled here as "Shock 2" has crossed through the slower CME's flux rope by this time. The image size is 310×310 R . explained the method of determination of these values in detail, using the 12 July 2012 case. They find that this CME had the following properties: 1. The CME speed was found to be 1265 Km/s using a linear fit to the height-time profile.
2. The direction of the CME was found to be at 12 degrees south and 8 degrees west.
3. The CME flux rope orientation, found using the GCS method, was 53 degrees with respect to the solar equator.
4. The poloidal flux of the CME was found to be 1.4 × 10 22 Mx.
5. The toroidal flux of the CME was found to be 2.1 × 10 21 Mx.  Figure 17: Earth probe data from the collision simulations are shown together with the solution for a slow CME only. We find that the collision resulted in the earlier arrival of slow CME. The horizontal axis represents time since CME insertion.
6. The CME flux rope was found to have a positive helicity sign.
7. The mass of the CME was found to be 1.65 × 10 16 g.
As explained in Sec. 2, we can input all these values easily into our model. There is one outstanding problem though. Our model shows a positive acceleration at the initial insertion, which can be at a height of 30-40 R . But the actual CME can be experiencing deceleration at those heights. We somehow need to match the speeds of the simulated and observed CMEs at some height after the initial insertion. In this study, we match the speed of the simulated CME at 50 R with the speed of observed CME at the same height obtained using the drag based model (DBM) [Vršnak, B. and Zic, T., 2007]. Using V CM E = 1265 at 15 R , the drag parameter of 0.1 × 10 −7 and the asymptotic solar wind speed equal to 400 km/s, BDM gives the CME speed at 50 R equal to 1140 km/s. The DBM model also predicts the CME height on 12-July-2012 21:30 UT to be 34 R . The drag parameter and the asymptotic solar wind speed values used here are based on the recommendation by Vršnak et al. [2014]. Using r 0 = 15, r 1 = 26, a = 5 and ξ = 2, the simulated CME acquired the speed of 1160 km/s at 50 R . We insert the CME in the background on 12-July-2012 21:30 UT, which makes the height of the CME at this time to be 36 R . Using this setup, the simulated CME arrives 3 hrs after the observed CME. In a future study, we will try to find the best approach to treat this speed-matching problem using multiple events, so that the arrival time difference can be minimized as much as possible.
In Fig. 18, we show the comparison between the observations and simulations at Earth. We use the 1 hour averaged solar wind data provided by NASA/GSFC's OMNI data set through OMNIWeb [King and Papitashvili, 2005]. The simulation results are probed along the Earth trajectory. Also, to visualize the error that could be due to the error in the initial direction of the CME, we put multiple probes in ± 5 degrees latitude-longitude region and probe the plasma values there as well. They are shown by green lines in Fig. 18. With red lines, we have plotted the simulation data when no CME was launched, to show the modulation of the solar wind by the CME. We note some differences between the model background solar wind and OMNI data ahead of the interplanetary shock arrival on 14 July 2012 ( 2012.5348) that are possible sources of error in the sheath region behind the shock. First, the background magnetic field polarity is different between the model and observations, which could have contributed to the errors in the sheath magnetic fields. Secondly, the model solar wind is slightly denser and faster than observed (i.e., higher dynamic pressure). The propagation of CME is sensitive to the solar wind background [e.g. Temmer et al., 2011, Lee et al., 2013 , so it is important to improve the pre-CME solar wind reconstruction, possibly by further optimizing the WSA parameters.
As for the comparison between CME values, we find that there is a good agreement in the density and B T values, and a fair agreement in B N and speed values. Our model did not capture the negative B R values exhibited by the actual CME. We also notice that the CME properties do not differ much in the ± 5 degrees region. This CME has previously been modeled by Shen et al. [2014] and Scolini, C. et al. [2019] using a magnetised plasma blob model and a spheromak model, respectively. We notice that magnetic field values due to our modified spheromak model agrees much better with observations compared to above mentioned studies. For example, Shen et al. [2014] were not able to match any component of magnetic field whereas Scolini, C. et al. [2019] were able to match only one component of magnetic field in one of their runs. However, the arrival time was missed by about 12 hrs. We believe that constraining both poloidal and toroidal field can be a reason for our model's better results. We have shown in our parametric study that changing these values can have a significant effect on 1 AU properties of CMEs. Another advantage our model has is that we introduce only a part of spheromak in the solar wind, such that it represents a bent tube geometry. Farrugia et al. [1995] have shown that a tube geometry works much better than the full spheromak geometry to reproduce magnetic cloud signatures at 1 AU. While our proposed model and its data constraining technique shows a significant im-provement on existing models, this model still does not capture the magnetic field signature of CMEs perfectly, for example, the negative B R values are not reproduced by this model at the Earth in this case. Several factors, including errors in initial direction, tilt, magnetic fluxes, etc. can be responsible for this, and we plan a detailed study to understand the impact of errors in initial parameters. We also plan to do a detailed study of this model using multiple events to understand its applicability as a operational model for CME predictions. We believe that the application of a flux rope-based model like ours will make a major advantage on existing non-flux-rope models. Figure 18: Comparison of the 1 hr averaged OMNI solar wind data (blue) and simulation results probed at Earth (black). Shock arrival is marked by the vertical dotted lines and the solid lines bound the magnetic cloud. The simulated CME arrived 3 hrs after the observed one. We also show the data probed at ± 5 degrees in latitude and longitude around Earth (green lines), to get an estimate of changes due to the possible error in the initial direction of the CME. Red lines show the solar wind solution when no CME was launched. We would like to emphasize that the B N values in RTN coordinates are approximately equal to the B z values in GSE coordinates

Conclusions
In this study, we demonstrated the application of our modified spheromak model to simulate CMEs in the inner heliosphere driven by the WSA coronal model providing time-dependent boundary conditions. We described a robust technique to introduce the flux rope into the domain. We presented a parametric study to see the effect of model parameters on the simulated CME evolution. They are summarized again below: • Initially, fast CMEs are decelerated and slow CMEs are accelerated. This has been confirmed by the observations as well.
• In our simulations, the CME speed increases with the poloidal flux but does not depend on the toroidal flux. We have shown that this is because the pressure gradient and magnetic tension forces in a bent flux rope geometry, such as ours, are much more dependent on the poloidal flux than on the toroidal flux. Since actual CMEs have similar geometry, we may expect them to show similar behavior.
• The energy multiplier ξ can be used to control the CME speed, since the speed depends on ξ linearly. However, ξ should be kept less than 3 to avoid unnatural thermal pressure in the initial flux rope.; • The simulated CME's speed increases with the decrease in the initial size of the flux rope. This is another parameter constraining the speed of simulated CMEs A CME-CME collision scenario, common during solar maxima can also be simulated using this model. We showed this in Sec. 3.2. Finally, in Sec. 3.3, we simulate the 12 July 2012 CME and compare the 1 AU properties with the observations. We find that our model was able to match the arrival time and plasma parameters with reasonable accuracy. It takes about 30 minutes of computation to simulate a CME propagation to 1 AU, on a 150 × 128 × 256 resolution using 1200 CPUs. With this robustness, this model can be used as an operational tool for space weather prediction.