Studying the Spheromak Rotation in Data-constrained Coronal Mass Ejection Modeling with EUHFORIA and Assessing Its Effect on the B_z Prediction

Coronal mass ejections (CMEs) are one of the major sources of space weather disturbances. If the magnetic field inside an Earth-directed CME, or inside its associated sheath region, has a southward-directed north–south magnetic field component (B_z ), then it interacts effectively with Earth's magnetosphere, leading to severe geomagnetic storms, depending on the strength of B_z (Wilson 1987; Tsurutani et al. 1988; Gonzalez et al. 1999; Huttunen et al. 2005; Gopalswamy et al. 2008). Therefore, it is crucial to predict the strength and direction of B_z inside Earth-impacting interplanetary CMEs (ICMEs) in order to forecast their geoeffectiveness. Since the magnetic field of CMEs cannot be reliably measured remotely, and direct in situ measurements of Earth-impacting ICMEs are routinely available only very close to our planet, modeling of CME magnetic properties using near-Sun observational proxies is paramount.

The state-of-the-art global heliospheric MHD models typically implement the axisymmetric spheromak or modified spheromak configurations to characterize the magnetic structure of a CME and simulate its evolution from Sun to Earth (Vandas et al. 1998, 2002; Manchester et al. 2004; Shiota & Kataoka 2016; Scolini et al. 2019; Asvestari et al. 2021; Singh et al. 2022). In a recent study by Asvestari et al. (2022), it was found that a spheromak injected into the inner heliospheric domain gradually rotates, due to the presence of a torque, until it reduces its magnetic potential energy. The rotation stops when its magnetic moment aligns with the ambient magnetic field. This behavior is a reminiscent of the condition known as the "spheromak tilting instability" (Rosenbluth & Bussac 1979; Bellan 2000, 2018; Mehta et al. 2020). According to theory, the torque is larger in the presence of a strong ambient field. Taking into consideration that the magnetic field is indeed stronger near the Sun but drops with heliodistance as 1/r², it is anticipated that the bulk of the spheromak rotation takes place near the Sun and the rotation rate drops with increasing heliodistance. Consistent with theory, in the modeling analysis of Asvestari et al. (2022) it was found that the rotation rate was higher near the inner boundary at 0.1 au and dropped as the spheromak moved away from it. Despite the drop in the rotation rate, the overall rotation angle was still notable at large heliodistances. Observational evidence suggests that some CMEs show clear signatures of rotation in the corona (e.g., Vourlidas et al. 2011; Nieves-Chinchilla et al. 2012; Kay & Opher 2015). CME rotation during its early phases of evolution is also noticed in simulations (Török & Kliem 2003; Lynch et al. 2009). However, such rotation mostly occurs in the low corona (below 0.05 au) and is less likely to happen at larger heliodistances in the heliosphere (Kay & Opher 2015). Therefore, the presence of spheromak rotation even beyond 0.1 au suggests that it may not be a realistic phenomena and thus can affect the performance of space weather forecasting models that use spheromaks as a flux rope.

In the study of Asvestari et al. (2022), a strong or weak unidirectional (inward or outward) magnetic field in the ambient medium was employed to quantify the spheromak rotation. Studying this phenomenon in solar wind conditions representing a specific actual period of time will be an important step forward to assess its effect on space weather forecasting. Event studies of past CMEs within the framework of global MHD models provide an excellent opportunity to compare the model results with in situ observations at different heliodistances. Therefore, using event-based data-driven MHD simulations employing the spheromak model allows us to build a quantitative understanding on how the spheromak rotation could affect the B_z prediction at 1 au.

There are several studies that have previously used spheromaks to study an observed CME event using data-driven MHD models (e.g., Shiota & Kataoka 2016; Scolini et al. 2019; Verbeke et al. 2019; Asvestari et al. 2021). However, a particular focus has not yet been made in any event-based study to address whether the spheromak rotation has significant consequences for B_z forecasts. In this work, we study a CME event on 2013 April 13 using the European Heliospheric Forecasting Information Asset (EUHFORIA) model aiming to understand the consequences of spheromak rotation in physics-based space weather forecasting models.

Previous studies with EUHFORIA use a default uniform mass density (1 × 10⁻¹⁸ kg m⁻³) to specify the mass density inside the spheromak. However, a recent observational study shows that the average CME density at 0.1 au (inner boundary of global heliospheric MHD models) ranges within the order of 1 × 10⁻¹⁸ to 1 × 10⁻¹⁷ kg m⁻³ (Temmer et al. 2021). Therefore, we also explore the effect of spheromak density in data-constrained CME modeling by performing a set of simulations with different spheromak densities that lie within the observational range. With this work, we address the following key scientific questions:

• 1.  
  How does an observationally constrained spheromak evolve/rotate in a data-driven global MHD simulation?
• 2.  
  Until what distance from the Sun does the ambient magnetic field play a major role in the spheromak rotation?
• 3.  
  Does the spheromak density have a considerable effect on its heliospheric evolution?
• 4.  
  Does the spheromak evolution/rotation affect the prediction of B_z at 1 au?

To address these questions, we organize the paper as follows. First, we briefly discuss the CME event on 2013 April 11 in Section 2. The methods to obtain the background solar wind and the observational techniques to constrain the spheromak to mimic the associated CME structure for this event are described in Section 3. Based on our study of the spheromak evolution using EUHFORIA, we present the results in Section 4. Finally, we summarize these results in Section 5 in the context of answering the science questions raised in Section 1.

We consider the CME event that erupted on 2013 April 11. This event has been studied extensively in previous works, for example, Sarkar et al. (2020) and Vemareddy & Mishra (2015). The CME appeared as a halo and was first observed at 07:24 UT in the field of view of the LASCO/C2 coronagraph. The eruption was associated with an M6.5 class flare that occurred in NOAA active region 11719 at around 06:50 UT at the position of N07E13. Figure 1 depicts the multiwavelength observation of the CME source region in AIA passbands and the CME morphology in white light as observed by LASCO. Based on the multispacecraft observations and reconstruction techniques, it was derived that the CME approximately propagates along the Sun–Earth line and was related to a well-defined magnetic cloud (MC) observed at 1 au (Vemareddy & Mishra 2015). In situ observations from WIND reveal that the leading edge of the MC arrived at L1 on 2013 April 14 at around 17:00 UT (Sarkar et al. 2020). Using the observations from heliospheric imagers on board STEREO, the interplanetary propagation of the CME was studied in Vemareddy & Mishra (2015), and it was found that the CME did not interact with any other CMEs during its propagation from Sun to Earth. The availability of cradle-to-grave observations of this event provides an excellent opportunity to model this event using near-Sun observational inputs and assess the model results at 1 au.

Zoom In Zoom Out Reset image size

The underlying methodology of the modeling performed in this work is to first generate a steady-state ambient solar wind condition in the corotating frame, representing the large-scale heliospheric plasma environment prior to the observed eruption, and then inject an observationally constrained spheromak into the medium to describe the propagating disturbance. The detailed methodology for the abovementioned steps are discussed below.

3.1. Constructing the Background Solar Wind

3.2. Constraining the Spheromak from Observations

Information on several characteristics of the CME is needed to constrain the spheromak model that is injected at the inner boundary of the heliospheric model at 0.1 au. All the required parameters and the methods of how they are determined are described in the following subsections. We first describe the geometric parameters (tilt and radius) followed by the kinematic (speed and propagation direction) and magnetic (helicity sign and magnetic flux) parameters. Finally, we discuss the density and temperature to be used for the spheromak. The majority of the parameters that we use to constrain the model are based on the multiwavelength and multispacecraft observational analysis reported in Sarkar et al. (2020) and Vemareddy & Mishra (2015). Notably, these studies determine the CME parameters assuming that the underlying flux rope structure is rooted in the Sun, which has a different geometry than that of the magnetically isolated spheromak. We emphasize below those differences and present the methods on how those CME parameters can be utilized to constrain a spheromak. The values used in the simulation runs are collected in Table 1.

Table 1. CME Input Parameters for Simulation Run 1

Parameter	Value
Insertion time	2013-04-11T11:12:00
Speed (vradial)	400 km s−1
Longitude ϕ	0°
Co-latitude θ	0°
Radius	7.2 RS
Density	1 × 10−18 kg m−3
Temperature	0.8 × 106 K
Helicity	−1 (left handed)
Tilt τ	−70°
Toroidal flux	2.4 × 1013 Wb

Download table as:  ASCIITypeset image

3.2.1. Tilt Angle

The insertion tilt of the spheromak is defined as the orientation angle (τ) of its symmetry axis measured clockwise from the meridional direction in the tangent plane to the inner spherical boundary (r = 0.1 au) at the point (θ, ϕ), which corresponds to the colatitude and longitude of the insertion direction of the spheromak (Verbeke et al. 2019; Asvestari et al. 2022). For example, for an insertion direction of θ = 0°, ϕ = 0° in the HEEQ coordinate system, τ = 0° corresponds to a scenario where the symmetry axis of the spheromak is oriented along the z-axis (along the solar rotation axis) so that its magnetic axis aligns along the y-direction. Under the abovementioned circumstance, the direction of the spheromak magnetic axis at the apex point (distant point from the Sun's center) on its magnetic axis curve points toward the positive y-axis.

On the other hand, the axis direction of a croissant-like flux rope obtained from a morphological fit to the white light observations is different from the tilt angle defined above. Observationally, the axial field direction of a croissant-like flux rope is determined from the near-Sun magnetic proxies and geometrical tilt angle of the flux rope obtained using the graduated cylindrical shell (GCS) model (Thernisien et al. 2009; Thernisien 2011). However, this axis orientation obtained from GCS refers to the orientation of the magnetic axis (relative to the solar equator) of a spheromak that is perpendicular to its symmetry axis (Asvestari et al. 2021). Therefore, one will need to translate the flux rope tilt extracted from observations with the GCS model to the insertion tilt of a spheromak in EUHFORIA.

Based on the near-Sun magnetic proxies (orientation of the magnetic connectivity between the two foot points of the preeruptive sigmoid) and the GCS fitting, the axis orientation at the apex of the CME flux rope under study is reported as northward with an angle of 70° at approximately 10 solar radii (R_S) (Sarkar et al. 2020). This angle (θ) is measured counterclockwise from the solar equator as shown in the inset of Figure 2. As the most significant rotation of a CME mostly occurs close (below 10 R_S) to the Sun (see Section 1), we use the axial direction of the CME obtained at ≈10 R_S as the input tilt angle (θ) at 0.1 au, that is, we assume that no further change of direction of the magnetic field structure occurs between 10 and 21.5 R_S.

Zoom In Zoom Out Reset image size

In order to estimate the associated angle τ, we first transform the orientation angle (θ = 70°) of the magnetic axis measured counterclockwise from the solar equator to that (20°; 90° −θ) measured clockwise from the meridional direction. Notably, the symmetry axis (shown by the yellow dashed arrow in Figure 2) of a spheromak is 90° apart (counterclockwise) from the direction of its magnetic axis (shown by the white dashed arrow in Figure 2) as mentioned above. Therefore, applying a counterclockwise rotation of 90° to the orientation angle (20°) of the magnetic axis of the spheromak, we determine the orientation angle of its symmetry axis as −70°. We use this value (−70°) as the input tilt angle (τ) for the spheromak. Figure 2 shows the orientation of the magnetic axis of the spheromak during its insertion phase, which is consistent with the orientation (see the inset of Figure 2) as reported in Sarkar et al. (2020).

3.2.2. Radius

The geometrical reconstruction of the 3D morphology of a CME can be approximated using the GCS model. The geometry of the GCS model provides both the edge-on and face-on radii of a CME flux rope that resembles a croissant-like structure. To transfer this shape to the spherical case, we average the edge-on and face-on radii of the GCS structure to determine the spheromak radius. Notably, we assume self-similar expansion of the CME up to 21.5 R_S from the last location where the 3D reconstruction was performed. Following the aforementioned assumption, we consider that the edge-on and face-on angular width of the CME as obtained from GCS fitting remains constant up to 21.5 R_S . Using the fitted edge-on and face-on angular width as reported in Sarkar et al. (2020), we obtain the radius of the spheromak as 7.2 R_S if its leading edge locates at 21.5 R_S .

3.2.3. Speed and Propagation Direction

The LASCO CME catalog provides a linear speed of the CME as 861 km s⁻¹, projected in the plane of sky as viewed from Earth. Applying the GCS fitting to the multi–vantage point white light observations of the CME at different instances of time, we find the deprojected 3D speed of the CME to be ≈930 km s⁻¹. Notably, the speed (V_3D) obtained from the GCS fitting can be decomposed as the translational speed (V_radial) and the expansion speed (V_exp) of the structure. The expansion of the linear force-free spheromak due to the pressure imbalance with the ambient solar wind is modeled self-consistently with EUHFORIA. Therefore, it is required to only use the V_radial as the input speed to run the simulation. We utilize the empirical relationship, V_radial = 0.43 × V_3D, as used in Scolini et al. (2019) to obtain the input translational speed of the linear force-free spheromak as 400 km s⁻¹.

Following the GCS results on this event studied by Sarkar et al. (2020) and Vemareddy & Mishra (2015), we notice that the propagation direction of the CME ranges between 9°N and 15°S in latitudinal and 0° and 13°E in longitudinal direction in HEEQ. Throughout the paper, we guide the reader about the spheromak evolution by presenting its cross-sectional map on the HEEQ equatorial plane and take the Sun–Earth line as the reference line with respect to which we show any change in orientation of the spheromak magnetic axis curve projected on the equatorial plane. Therefore, we select 0° longitude and 0° latitude in HEEQ (within the abovementioned observational range) as the insertion direction of the spheromak, as in such a scenario the equatorial plane approximately cuts through the spheromak centroid and propagates nearly along the Sun–Earth line. We further assess the effects of uncertainties related to the direction of propagation by comparing the model results from multiple virtual spacecraft located within 15°E–15°W and 15°N–15°S at 1 au (see Section 4.3).

3.2.4. Helicity Sign

Observations of the inverse S-shaped morphology of the preflare sigmoid structure as reported in Vemareddy & Mishra (2015) Sarkar et al. (2020) suggest a left-handed chirality of the associated flux rope. This helicity sign is also consistent with the hemispheric helicity rule as the source active region of the CME was located in the northern hemisphere. We therefore set the chirality of the spheromak to be negative.

3.2.5. Magnetic Flux

Observational studies show that the poloidal flux of a CME can be obtained from the estimation of reconnection flux as computed from the cumulative flare ribbon area (Kazachenko et al. 2017) or posteruption arcades (Gopalswamy et al. 2018a). In order to constrain the magnetic flux of a spheromak in EUHFORIA, Scolini et al. (2019) equates the observationally estimated poloidal flux to that of a spheromak. However, we notice that the magnetic axis of a spheromak resembles that of a full torus-like structure without having any anchoring points on the Sun. On the other hand, the observed poloidal flux is believed to be distributed over a nearly half-torus-like structure with an angular extent similar to the side-on angular width of a CME. Therefore, equating the observed poloidal flux to that of a spheromak would underestimate the field strength at its magnetic axis due to the distribution of observed poloidal flux over a comparatively large volume of a spheromak.

In order to avoid the abovementioned underestimation, we follow a different approach to constrain the magnetic flux of the spheromak. Instead of directly equating the observed poloidal flux with that of the spheromak, we equate the field strength (B_spheromak) at the magnetic axis of the spheromak to the axial field strength (B₀) obtained from the flux rope eruption data (FRED) technique (Gopalswamy et al. 2018b). Conceptually, FRED uses the information on poloidal flux and a Lundquist flux rope model (Lundquist 1950) constrained with realistic CME geometry to estimate the axial field strength of a CME. Further, using the value of B_spheromak, we estimate its toroidal flux content, which is used as input for the EUHFORIA simulation. Therefore, the axial field strength (B_spheromak) of a spheromak remains identical to that obtained from FRED and does not depend on the spheromak volume.

Sarkar et al. (2020) reported the average observed poloidal flux of the associated CME under study as 2.1 × 10¹³ Wb and the axial field strength (B₀) at 10 R_S as 5200 nT. Applying the FRED technique further out to the inner boundary (21.5 R_S) of EUHFORIA, we obtain B₀ at 14.3 R_S (21.5 R_S − r_sph) as 2500 nT. Notably, each point on the magnetic axis of a spheromak is not equidistant to the Sun center due to the circular shape of the axis. Therefore, we use the radial distance (14.3 R_S) of the spheromak center (21.5 R_S − r_sph) for deriving B₀, which is the average distance of all the points on its magnetic axis. Equating the estimated B₀ value (2500 nT) to B_spheromak, we determine the toroidal flux of the spheromak as 2.4 × 10¹³ Wb and use this value as input for EUHFORIA simulation.

3.2.6. Temperature and Density

The default values for the CME mass density (1 × 10⁻¹⁸ kg m⁻³) and temperature (0.8 × 10⁶ K) are widely used in studies employing EUHFORIA and are set to be uniform inside the CMEs (Pomoell & Poedts 2018; Scolini et al. 2019; Asvestari et al. 2021). For the initial simulation run (Run1), we use these default values of the CME mass density and temperature as input for the spheromak.

We further study the dependence of the CME evolution on the choice of the initial mass density by conducting a set of EUHFORIA simulations with density values higher than the default one. Five simulations (Run2, 3, 4, 5, and 6) are performed by selecting the density values within the observational range (1 × 10⁻¹⁸–1×10⁻¹⁷ kg m⁻³) at 0.1 au (Temmer et al. 2021) in steps of 2 × 10⁻¹⁸ kg m⁻³ (see Table 2). We discuss the results based on these runs in Section 4.2. We notice that the spheromak centroid height locates at ≈14 R_S when we start to insert it at the lower boundary (0.1 au) of the simulation. The observational upper limit of the CME density at ≈14 R_S turns out to be ≈5 × 10⁻¹⁷ kg m⁻³ as reported in Temmer et al. (2021). Considering this higher-density value, we perform a further simulation (Run7) and refer to it as the "high-density (HD) run." In Section 4.2, we compare the results of the HD run with Run6 and Run1, which are hereinafter referred to as the "moderate-density (MD) run" and "low-density (LD) run," respectively. For all runs, we keep the rest of the input parameters unchanged.

Table 2. Model Runs with Different Density Values for the CME

Run Number	Spheromak Density
Run1 (LD)	1 × 10−18 kg m−3
Run2	2 × 10−18 kg m−3
Run3	4 × 10−18 kg m−3
Run4	6 × 10−18 kg m−3
Run5	8 × 10−18 kg m−3
Run6 (MD)	10 × 10−18 kg m−3
Run7 (HD)	50 × 10−18 kg m−3

Download table as:  ASCIITypeset image

For all the runs in this work, a uniform mesh grid is used with a 2° angular resolution and a radial grid spacing of Δr ≈ 0.0037 au ≈0.8 solar radii (512 cells in radial direction).

4.1. Spheromak Rotation

The evolution of the spheromak in the heliospheric domain as obtained from the simulation results for model Run1 is illustrated in Figure 3. In order to visualize the 3D orientation of its magnetic axis, we plot a selection of field lines that wrap around the magnetic axis of the spheromak. We do not plot the less-twisted magnetic field lines passing close to the symmetry axis of the spheromak. Therefore, in the visualization of Figure 3, a void structure appears at the central part of the spheromak, shaping it to a torus-like structure (e.g., see panel (a) or (l)). Tracking this structure at different instances of time during its propagation from Sun to Earth, the image shows that the axis orientation of the spheromak changes considerably from that of its insertion at 0.1 au.

Zoom In Zoom Out Reset image size

In the course of the insertion, as depicted in the right-most column (Figures 3(d), (h), and (l)), the spheromak orients in such a way that the central void part of its analogical torus-like shape is fully observable from the side view on the meridional plane (see panel (l)) and its magnetic axis curve as projected on the equatorial plane nearly aligns to the x-direction, that is, the Sun–Earth line (see panel (h)). As the spheromak propagates further out in the heliosphere, it rotates in such a way that the central void part of the torus-like shape starts to disappear from the side view of the meridional plane (see the evolution from right to left in the lower row of Figure 3, i.e., panels (l), (k), (j), and (i), respectively) and becomes partial to fully visible along the Sun–Earth line (upper row of Figure 3, i.e., panels (d), (c), (b), and (a), respectively). The evolution of its projected axis on the equatorial plane as shown in the middle row of Figure 3 clearly depicts that it exhibits a clockwise rotation (as viewed from the top of the equatorial plane) of ≈90° in the interplanetary space and becomes almost perpendicular to the Sun–Earth line (x-direction) when it arrives at 1 au. We present a more quantitative analysis of the changes in the orientation angle of the spheromak with increasing heliocentric distances from 0.1 to 1 au in Section 4.2.

The significant rotation of the spheromak as obtained from Run1 is consistent with the finding of spheromak rotation as reported in Asvestari et al. (2022). Indeed, the symmetry axis of the spheromak rotates in such a way that it aligns approximately along the Sun–Earth line. As a result, the virtual spacecraft at Earth encounters the least twisted part of the spheromak along its symmetry axis and misses the part of the magnetic structure that is constrained from observations (see Figure 2). This can be a significant issue for space weather modeling when using spheromak as CMEs.

Notably, we performed Run1 with the default density value for spheromak that has been used in previous studies with EUHFORIA. We further present the results of a set of simulations (see Table 2) performed with higher-density values of spheromak in Section 4.2.

4.2. Spheromak Evolution with Higher Densities

We find that the spheromaks with higher densities result in higher magnetic field strengths at the position of Earth. We present this result in Figure 4, which shows the temporal profiles of spheromak field strength (top panel) at 1 au from a selection of runs performed with different densities. Notably, a spheromak with higher density is expected to undergo higher expansion due to the increased internal thermal pressure. Indeed, the comparison of the model results at the position of Earth as obtained from Run1 and Run6 (Figure 4) depicts that the spheromaks with initial higher density become larger in size (enclosed by the two vertical dotted lines in red in the figure) compared with the ones with initially lower density (enclosed by the two vertical dotted lines in magenta) based on their appearance in the in situ observations by a virtual spacecraft at Earth. Therefore, due to the higher internal expansion, one would expect that the HD spheromaks will result in a correspondingly lower magnetic field strength upon its arrival at 1 au, which is in contrast to our results as mentioned at the beginning of this section (top panel of Figure 4). Notably, different trajectories of the spacecraft through the spheromak can significantly contribute to the changes in its magnetic strength profile and the inferred size. If the amount by which the spheromak rotation varies due to varying the initial mass density, then the virtual spacecraft at Earth would probe different parts of the spheromak for different cases. Therefore, we further compare the results of Run1 and Run6 for several virtual spacecraft placed at different HEEQ longitudes and latitudes within ±20° from the nominal position of Earth at 1 au (see Figure 5). The plasma beta of the spheromak for Run1 and Run6 in both the equatorial (Figures 5(a) and (b)) and meridional planes (Figures 5(d) and (e)) clearly shows that the size of the spheromak is larger for Run6 (indicated by the red dashed contour) compared with Run1 (indicated by the yellow dashed contour), which is in agreement with our results based on Figure 4. Figures 5(c) and (f) shows that the magnetic field strength profiles obtained at the different virtual spacecraft for Run6 (orange shaded profile) are predominantly higher compared with that for Run1 (cyan shaded profile). This again indicates that the HD spheromaks lead to higher field strengths at 1 au as also seen in Figure 4. However, the considerable variability in magnetic field strength observed as virtual spacecraft traverse distinct regions within the spheromak for both Run1 and Run6 (Figures 5(c) and (f)) highlights the notion that if the spheromak undergoes varying degrees of rotation due to different densities, it would result in diverse magnetic profiles observed at Earth. On the other hand, the predominantly higher field strength obtained for the higher-density run (Run6) incorporating the different spacecraft crossing through the spheromak suggests that the variance in spheromak rotation alone cannot account for these observations, implying the existence of an additional compressional influence acting upon the spheromak. Therefore, we suggest the following two possible scenarios that contribute to the diverse spheromak profiles observed at Earth under varying density conditions:

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

(1) The spheromaks rotate differently when their density increases. Therefore, the virtual spacecraft located at the position of Earth encounter different parts of the spheromak for the runs performed with higher-density values, which results in different field strength at 1 au.

(2) The spheromaks with higher density undergo compression, which leads to higher field strength at 1 au when compared with the LD cases.

We explore the abovementioned two possibilities based on our results in the following Sections 4.2.1 and 4.2.2, respectively.

4.2.1. Effect of Density on Spheromak Rotation

The 3D configuration of the magnetic field of the CME around the magnetic axis at the time when the leading edge of the CME reaches 1 au for the LD (panel (a)) and HD (panel (b)) runs is visualized in Figure 6. The white dashed lines overplotted on each panel of Figure 6 depict the approximate projection of the magnetic axis curve on the equatorial plane. The image clearly illustrates that at 1 au the orientation of the magnetic axis for the HD run is significantly different from that for the LD run. Interestingly, we find that the magnetic axis associated with the HD run undergoes less rotation, as the projected orientation of the magnetic axis curve at 1 au is nearly aligned with its initial orientation (approximately along the Sun–Earth line; see Figure 3(h) and Figure 6(b) for comparison) at the time of insertion.

Zoom In Zoom Out Reset image size

We further present a quantitative analysis of the change in the axis orientation of the modeled CME as it evolves from 0.1 to 1.0 au launched with different densities. A detailed description of the methodology to estimate its orientation angle is given in the Appendix. Figure 7 depicts the evolution of the magnetic axis orientation (projected on the equatorial plane) with increasing heliocentric distances for the HD, MD, and LD runs. Note that the orientation angle is measured with respect to the Sun–Earth line, which is the approximate direction of propagation for this event. It is clearly delineated in Figure 7 that while the magnetic axis undergoes a large rotation (≈90°) for the LD run, the amount of rotation is significantly lower for the MD (≈70°) run and in particular so for the HD (≈35°) run before it arrives at 1 au. Importantly, a significant part of the rotation for all the runs occurs below 0.3 au. This indicates the role of the strong ambient magnetic field that exerts higher torque on the ejected structure when it is close to the Sun (Asvestari et al. 2022). Our finding of the density dependence on the change of orientation suggests that the ambient magnetic field becomes less effective in rotating the more massive spheromaks due to the higher linear momentum of the spheromak toward its direction of propagation. Therefore, the high-density spheromaks exhibit less rotation than the low-density ones.

Zoom In Zoom Out Reset image size

4.2.2. Effect of Density on the Internal Expansion of the Spheromak

The choice of the initial plasma density in the spheromak not only changes the amount of rotation it experiences during its interplanetary propagation but also affects its internal expansion as well. The magnetic field magnitude associated with the spheromak and the north–south component (B_z in HEEQ) of the magnetic field in the HEEQ equatorial plane are shown in Figure 8. The B_z component for the LD, MD, and HD runs as shown in the bottom panels of Figure 8 shows two separate regions of oppositely directed B_z that are made more evident by the enclosed red and blue dashed contours drawn at levels ±2.5 nT, respectively. For the sake of further explanation, we refer to these two regions of dominant positive and negative B_z as red and blue regions following the color map as depicted in the figure. In particular, these red and blue regions as shown on the cross-sectional slices of the spheromak contain the twisted field lines (see Figure 7), centering its magnetic axis, which is pointing in and out of the equatorial plane, respectively.

Zoom In Zoom Out Reset image size

Notably, the magnetic axis orientation of the spheromak in the red region bearing the positive B_z in Figure 8 points toward the northward direction, which has been constrained from the observed axial direction (see Figure 2) of the CME. Being a closed curve (see Figures 3(a) and 3(b)), the magnetic axis of the spheromak turns its direction from the red region to the blue region and becomes southward in the blue region. During the insertion phase of the spheromak, the red region enters into the simulation first, followed by the blue region (see Figures 3(h) and (g) sequentially) as the spheromak itself is an isolated magnetic structure and does not have any legs attached to the Sun. However, instead of propagating as one followed by another, the two aforementioned regions of the spheromak for the LD run approximately move parallel to each other (see Figure 8(d)) along the direction of propagation due to its significant rotation below 0.3 au. Therefore, during the majority of the propagation path (from approximately 0.3 au onward), both regions remain in contact with the solar wind ahead. This causes similar interplanetary evolution of both regions, which results in their similar shapes at 1 au (see Figure 8(d)).

However, for the HD run in the absence of any significant rotation, the regions of opposite B_z polarity move approximately one after another along the propagation direction (see Figure 8(f)). Therefore, in this case, the northward and southward field portions in the spheromak interact differently with the background solar wind as the negative B_z portion no longer remains in contact with the solar wind ahead. The drag force due to the upstream solar wind ahead acts almost entirely on the positive B_z portion (enclosed by the red dashed contour) at the front, which results in a curved and compressed structure of the frontal positive B_z part (see Figure 8(f)). Due to this compression, the field strength at the frontal part of the spheromak enhances as depicted by the color map in Figures 8(c) and (f). More precisely, the maximum field strength of the spheromak as it passes through the virtual spacecraft at 1 au for this case (Run7) reaches 19.3 nT, which is almost 2 times larger than that (11.3 nT) as obtained in Run6 (see the top panel of Figure 4). On the other hand, the rearward negative B_z portion of the spheromak gets shielded from the upstream solar wind by the frontal portion and therefore does not undergo any significant compression or change in shape.

An intermediate situation arises for the MD run (Figures 8(b) and (e)) where the magnetic field structure undergoes less rotation compared with the LD run but experiences a larger rotation compared with the HD run. In this scenario, the region bearing the positive B_z value partially shields the other region located at a smaller heliocentric distance and experiences less compressive interaction with the upstream solar wind compared with the HD run. Therefore, the results presented in this section clearly explain why the HD spheromaks result in higher magnetic field strength at 1 au as depicted in Figure 4.

4.3. In Situ Comparison

The comparison of the in situ virtual spacecraft observations at Earth with the simulation results in Figure 4 indicates that the observed magnetic field strength of the ICME most closely resembles the results of the MD run (Run6). Notably, the input density (1 × 10⁻¹⁷ kg m⁻³) used in Run6 is in the same order of magnitude as that (2.2 × 10⁻¹⁷ kg m⁻³) estimated for this event by employing the observational techniques in Temmer et al. (2021). To further assess this particular run, the in situ results obtained from Run6 are compared with the corresponding observations from the WIND spacecraft for all components of the magnetic field vector as well as the plasma density and speed (right panel of Figure 9). Apart from the magnetic field strength, the modeled speed profile (right panel, last row) also shows a remarkably good agreement with that of the observed ICME and the preceding solar wind stream. Notably, the consistency between the observed and modeled arrival time of the high-speed stream ahead of the ICME was achieved by using the optimized rotation angle (ψ) for the background solar wind solution as discussed in Section 3.1.

Zoom In Zoom Out Reset image size

It is important to note that the magnetic field structure of the spheromak for this particular run undergoes a rotation during its propagation such that both the northward and southward field (in HEEQ) portions (indicated by the red and blue domains) are intersected by the virtual spacecraft at Earth (see the left panel of Figure 9). However, only the magnetic axis orientation in the northward field portion of the spheromak is consistent with the axis orientation of the associated near-Sun flux rope as discussed in Section 4.2.2. The two regions enclosed by the blue and red dashed lines in the left panel of Figure 9 clearly depict that the rotation of the horizontal field (${B}_{x}{\hat{x}}_{\mathrm{HEEQ}}+{B}_{y}{\hat{y}}_{\mathrm{HEEQ}}$) component (indicated by the white arrows in the left panel of the figure), as well as the B_z component in those two regions, is in a different direction. Indeed, these two regions in the cross-sectional plane of the spheromak resemble the cross section of two high-inclination cylindrical-like flux ropes with oppositely directed axial magnetic fields. This indicates that in such a scenario the in situ spacecraft encounters a double flux rope signature when the whole spheromak passes through it. Therefore, we only take into account the northward field portion of the spheromak that is constrained from the observations while comparing the model results with in situ observations. In this case, the domain enclosed by the red dashed line (as depicted in the left panel of Figure 9) is part of the spheromak that bears the flux rope signature constrained from the observations. The gray shaded region in the right panel of Figure 9 indicates the temporal passage of the abovementioned part of the spheromak where the B_z component is throughout positive. Therefore, we compare the magnetic vector profiles (red solid lines) of the spheromak inside the gray shaded region with those (black solid lines) observed inside the MC boundary as indicated by the two vertical blue dashed lines (selected based on the ICME plasma parameters as reported in Vemareddy & Mishra (2015) for this event). This approach helps us to avoid any misinterpretation in the in situ assessment of the modeled magnetic vectors due to the presence of double flux rope signatures inside a spheromak.

A smooth rotation from positive to negative as observed in the B_y component and a predominant positive B_z profile of the observed MC (see the black solid curves in the B_y and B_z plot within the vertical blue dashed lines in the right panel of Figure 9) are well captured by the model results within the gray shaded region. Notably, the leading front of the spheromak arrives earlier than that of the observed MC. Therefore, for the sake of comparison, we shift the modeled magnetic vectors by 9 hr 40 minutes so that the front edge of both the modeled and observed flux rope temporally coincide. The time-shifted modeled magnetic profiles of B_y and B_z as indicated by the cyan solid lines depict a good agreement between the observation and model output. In particular, the modeled B_y component shows excellent correspondence with the observed profile, while a decent match is achieved for the B_z component. However, the modeled B_x component turns out to be opposite to the observed one as the virtual spacecraft at Earth intersects through a different part of the spheromak. In the absence of any significant rotation, we show in Figures 10 and 11 that all three components (B_x , B_y , and B_z ) of the flux rope magnetic field can be well captured by the model within the uncertainty limit in obtaining the CME's direction of propagation. Notably, the amount of rotation exhibited by the spheromak is the least for the HD case (Run7) and maximum for the LD case (Run1) (see Figure 7) compared with the other runs performed in this work. Therefore to assess how the spheromak rotation affects the uncertainty in predicting the magnetic vectors of ICMEs, we present a comparative uncertainty analysis for these two extreme cases (Run1 and Run7) in Figures 10 and 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The in situ assessment made for the modeled results obtained from different virtual spacecraft placed within ±15° longitude (see Figure 10) and latitude (see Figure 11) gives insight into how the spheromak rotation affects the B_z prediction at 1 au. In the presence of significant rotation of the spheromak, the modeled B_z component becomes highly sensitive to the spacecraft location as shown in Figures 10(b) and 11(b). Indeed, the figures show that completely opposite profiles of B_z can be obtained at nearby spacecraft as the set of modeled B_z in such a scenario can start with positive as well as negative values as depicted in panel (b). This indicates that the spheromak rotation leads to a large uncertainty in B_z prediction for this event.

On the other hand, in the absence of significant spheromak rotation, the uncertainties in B_z prediction significantly reduce as shown in Figures 10(d) and 11(d). All the possible B_z profiles at different nearby virtual spacecraft show predominant positive values of B_z , which is in agreement with the observed B_z profile at 1 au. These results suggest that the prediction efficacy improves when a spheromak undergoes minimal rotation effect.

We also note that the simulation outputs of EUHFORIA capture well the formation of the shock as well as the sheath structure ahead of the spheromak as illustrated in Figure 12. However, a detailed study of the sheath region requires a higher spatial resolution in the simulation, which is outside the scope of the current study. In future work, we plan to explore the formation and evolution of sheath regions in detail.

Zoom In Zoom Out Reset image size

In this work, we have assessed the applicability of spheromaks to be used as a model for the magnetic field of CMEs in global MHD simulations for space weather forecasting. In order to test the performance of a spheromak model in real-event forecasting, we carry out a set of data-constrained MHD simulations with EUHFORIA for an Earth-impacting CME event on 2013 April 11. As the spheromaks are prone to the "tilting instability" that causes inherent rotation of its magnetic axis (Asvestari et al. 2022), we perform a quantitative analysis of the resultant rotation experienced by the spheromak to understand how that affects the B_z prediction at 1 au. Further, we study the role of spheromak density to mitigate the rotation of the spheromak that may reduce the uncertainty in forecasting B_z . The important findings of this study are discussed below in order to answer the key scientific questions as mentioned in Section 1.

• 1.  
  Our simulation results confirm the presence of spheromak rotation, reminiscent of the tilting instability in data-constrained MHD simulations. Indeed, we find that a significant rotation up to 90° may occur during the heliospheric propagation of a spheromak. The 3D visualization of the spheromak rotation reveals that this rotational motion is different from the CME rotation that may occur in the lower corona due to the unwrithing motion (Lynch et al. 2009; Zhou et al. 2022). During the lower coronal evolution, the magnetic axis of a CME may rotate about its direction of rise or propagation (Zhou et al. 2022). In contrast, we find that the rotation of the spheromak magnetic axis as observed for this event occurs approximately about the line perpendicular to the direction of propagation. Moreover, the sense of CME rotation in the lower corona follows its chirality, that is, CMEs with negative/positive chirality rotate counterclockwise/clockwise (Green et al. 2007; Lynch et al. 2009; Zhou et al. 2022). However, the sense of rotation of the magnetic axis of the spheromak is not solely dependent on its chirality but also depends on the direction of the ambient magnetic field (Asvestari et al. 2022).
• 2.  
  The most significant part of the spheromak rotation is observed to take place close to the Sun below 0.3 au. The underlying reason for this is that the ambient magnetic field is stronger close to the Sun. Therefore, the torque force exerted on the spheromak becomes less effective above 0.3 au.
• 3.  
  In the presence of significant spheromak rotation, the predicted magnetic field topology of the ICME at 1 au is expected to show poor results when the simulation output is compared with the observational signatures. Interestingly, we find that the spheromak density has a major role in mitigating its rotation effect. Running a set of simulations by using different density values within the observed range, we find that the spheromak rotates less when the density is higher. This can be explained as a consequence of the increased moment of inertia (I) of the more dense spheromaks. Assuming the torque (T) exerted by the background wind remains the same, the angular acceleration (α) of the spheromak needs to decrease since T = I × α. Notably, the insertion speed of the spheromak also has a role (Asvestari et al. 2022) in the amount of rotation experienced by the spheromak. Therefore, a similar effect could probably be achieved by not changing the density but by instead changing the speed. However, density is the focus in the paper as in general the speed of the CME is much more well constrained from remote observations than the mass density.
• 4.  
  As the degrees of spheromak rotation significantly vary between high-density and low-density spheromaks, different portions of the spheromak are being probed by a virtual spacecraft at Earth under varying density conditions. For the low-density spheromak that exhibits the maximum rotation (≈90°), the virtual spacecraft at Earth approximately crosses through the symmetry axis of the spheromak, thereby missing the twisted flux rope part. On the other hand, in the case of high-density spheromaks with minimal rotation, the virtual spacecraft at Earth predominantly passes through the twisted part within the spheromak. Therefore, spheromaks with different densities lead to different in situ magnetic profiles at Earth. In addition to this, the part of the spheromak that accounted for modeling the observed flux rope undergoes significant compression in the high-density case compared with the low-density scenario, leading to further changes in the magnetic strength of the spheromak under varying density conditions. As a combined result of the two aforementioned phenomena, the high-density spheromaks exhibit higher field strength as probed by the virtual spacecraft at Earth.
• 5.  
  Our assessment of the simulation results with the in situ observations for the event under study shows that the spheromak rotation leads to large uncertainties in B_z predictions, whereas the prediction efficacy significantly improves in the absence of any significant rotation of the spheromak.

Our results imply that a different part of the spheromak than that constrained from the observations may arrive at 1 au due to its rotation in the interplanetary domain. As a consequence of this rotation, the prediction of CME magnetic vectors at 1 au can be largely affected. Therefore, care must be taken when using the spheromak in global MHD models for space weather forecasting. In particular, the magnetic configuration of the simulated flux rope should be visualized in 3D to check the effect that the changes in tilt during the propagation have on its final orientation at 1 au.

Intuitively, inserting only half of the spheromak in the heliospheric MHD domain may result in a different scenario compared with that reported in this study. In such a scenario, when the peripheral magnetic field of the partially inserted spheromak gets peeled off due to erosion, it may mimic a flux rope structure with two legs attached to the inner boundary of the simulation. A similar situation would be expected to arise in the case of inserting a half-torus-like flux rope. Under those circumstances, the rotation of the flux rope magnetic axis may not be as significant as that for a fully inserted spheromak or a torus. However, as claimed in Asvestari et al. (2022), a flux rope with two legs attached to the inner boundary may still have a magnetic moment that will try to align with that of the ambient magnetic field. As the legs of the flux rope remain fixed on the inner boundary, the manifestation of the force acting on the flux rope under that circumstance may mostly result in a deflection, in contrast to the large unrealistic rotation as expected for the cases of fully inserted spheromaks. It is possible that such deflections could be realistic and therefore would be important to incorporate into global MHD models from the perspective of space weather forecasting. As the ambient magnetic field is largely effective below 0.3 au (see Section 4.2.1), the MHD modeling approaches that superimpose a flux rope in background wind when the leading-edge height is already at approximately 0.3 au (e.g., see Singh et al. 2022) are not capable of capturing such an effect happening within 0.1–0.3 au. Therefore, we emphasize that in contrast to the method of superposing a flux rope, the method of inserting that from the inner boundary of the MHD model is capable of capturing its early evolution close to the inner boundary (0.1 au) of the heliospheric domain. In future work, we plan to implement the insertion of a half-torus-like flux rope with two legs attached to the inner boundary of the heliospheric model in EUHFORIA. This would allow us to explore whether such an implementation technique still leads to some amount of deflection and rotation of the flux rope, which can then be confirmed as the realistic scenario present in the CME evolution.

This research has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 870405 (EUHFORIA 2.0). R.S. acknowledges support from the project EFESIS (Exploring the Formation, Evolution, and Space Weather Impact of Sheath Regions), under the Academy of Finland grant No. 350015. J.P. acknowledges Academy of Finland Project 343581. E.A. acknowledges support from the Academy of Finland (Postdoctoral Researcher grant No. 322455). E.K., E.A., and J.P. acknowledge the ERC under the European Union's Horizon 2020 Research and Innovation Programme Project 724391 (SolMAG). The work at the University of Helsinki was performed under the umbrella of the Finnish Centre of Excellence in Research of Sustainable Space (Academy of Finland grant No. 312390, 336807). N.W. acknowledges support from NASA program NNH17ZDA001N-LWS and from the Research Foundation—Flanders (FWO-Vlaanderen, fellowship No. 1184319N). A.M. and S.P. acknowledge support from the projects C14/19/089 (C1 project Internal Funds KU Leuven), G.0025.23N (WEAVE FWO-Vlaanderen), SIDC Data Exploitation (ESA Prodex-12), and Belspo project B2/191/P1/SWiM. Open access is funded by Helsinki University Library.

We estimate the trend of spheromak rotation for different density runs as shown in Figure 7 by identifying the orientation of its magnetic axis curve during different phases of its evolution in the heliospheric domain. For this purpose, we developed a method to estimate the projected orientation of the magnetic axis curve of a spheromak on any plane that passes close (well within the radius of the magnetic axis) to its centroid and is perpendicular to its magnetic axis.

In order to describe the method, we present an example of a cross-sectional plane of the analytical linear force-free spheromak used in EUHFORIA as shown in Figures 13(a) and (d). The spheromak solution on this cross-sectional plane is expressed in Cartesian coordinates (B_x , B_y , and B_z ), where B_x and B_y lie in the plane of the figure and B_z is perpendicular to the depicted plane. The strength of B_h ($=\sqrt{{B}_{x}^{2}+{B}_{y}^{2}}$) as shown in Figure 13(a) depicts the strength of the poloidal magnetic field, whereas the B_z map as shown in Figure 13(d) represents the toroidal field. The overplotted black streamlines delineate the direction of the poloidal magnetic field, which swirls around a common center (marked by the white dots in Figure 13(d)) indicating the location of the magnetic axis as it crosses through the cross-sectional plane. Therefore, identifying the locations of those two swirling centers of the poloidal field would allow us to estimate the projected orientation of the magnetic axis curve by connecting those two points on the aforementioned plane. The B_h map as shown in Figure 13(a) shows that the value of B_h attains a minima at the two swirling centers of the poloidal field as indicated by the color map. Notably, it can be further observed that the B_h value also decreases toward the two poles of the symmetry axis of the spheromak. Therefore, we apply a mask on the B_h map as indicated by the regions enclosed by the red dashed contours (see Figure 13(a)) within which the strength of B_z is greater than 0.6 times the maximum strength of B_z within that plane. This selection of mask helps us to discard the regions of lower B_h value toward the poles of the symmetry axis. Therefore, within the masking area, identifying the two points where the B_h value becomes minimum gives us the locations of the two points where the magnetic axis of the spheromak passes through the plane. Finally, connecting those two points with a straight line gives us the orientation of the projected magnetic axis curve of the spheromak.

Zoom In Zoom Out Reset image size

We apply this method to track the change in rotation angle of the spheromak magnetic axis in the simulation outputs obtained from EUHFORIA. As a first step, we identify the plane in which the magnetic axis curve of the spheromak is clearly seen to rotate during its heliospheric evolution. Based on the 3D visualization of the spheromak magnetic axis as shown in Figures 3 and 6, we identify that the rotation of the spheromak for the runs performed in this work is well seen in the HEEQ equatorial plane. Therefore, we use the 2D map of B_x , B_y , and B_z in HEEQ on the equatorial plane to track the orientation of the magnetic axis curve of the spheromak. Before applying the method, we multiply the components of the magnetic field at each grid point of the 2D map with r², where r is the distance of each grid point from the Sun center. This helps us to remove the gradient of magnetic field strength that decreases with larger helio distances (r) following the 1/r² relation.

We show the results of our tracking method during two different time frames for Run7 as shown in Figures 13(e) and (f). Figures 13(b) and (c) show the HEEQ B_h ($=\sqrt{{B}_{x}^{2}+{B}_{y}^{2}}$) component of the spheromak magnetic field at the equatorial plane. The overplotted red contours on the aforementioned B_h maps enclose the masking area within which we identify the local minima of B_h . We further connect the identified conjugate locations of the B_h minima by a white dashed line as plotted on top of the B_z maps shown in Figures 13(e) and (f). This white dashed line represents the projected magnetic axis curve as identified from the tracking method. Comparing this with the background streamline plot of the equatorial component of the magnetic field shows that the locations of two swirling points of the poloidal magnetic field inside the spheromak cross section are correctly identified by the two endpoints of the white dashed line. We consider the centroid (as indicated by the red dot) of the two endpoints of the projected magnetic axis curve as the reference point to measure the distance of the spheromak from the Sun center. The rotation angle (α) of the projected magnetic axis curve as plotted in Figure 7 is measured with respect to the reference line (parallel to the X_HEEQ or Sun–Earth line) as indicated by the yellow dashed line in Figure 13 (f).

Notably, the tracking method applied in this work is different from the one employed in Asvestari et al. (2022), which tracks the symmetry axis instead of the magnetic axis of the spheromak to estimate its orientation angle. However, tracking the magnetic axis as demonstrated in this work turns out to be a useful alternative tool incorporating the cases of HDs (e.g., Run7), where it becomes difficult to define the symmetry axis due to the high compression at the frontal part of the spheromak (see Figures 8(f) and 13(f)).