© The Authors 2024

1. Introduction

The changing physical conditions in the heliosphere induced by the dynamical processes on the Sun, and those in the interplanetary environment, appearing in the solar wind, magnetosphere, ionosphere, and thermosphere, are called space weather. The main drivers of disturbed space weather conditions are coronal mass ejections (CMEs, Howard 2011; Webb & Howard 2012), the giant blobs of magnetised plasma erupting from the Sun into the heliosphere. These large-scale structures propagate through the heliosphere, interacting with the solar wind and often causing geomagnetic disturbances on Earth and other planets and spacecraft.

Interplanetary CMEs (ICMEs) are those propagating beyond the corona into the interplanetary medium. ICMEs manifest lucid magnetic field characteristics that can be explained by a flux-rope structure. Flux ropes are organised bundles of axially twisted magnetic field lines confining plasma within them (Antiochos et al. 1999; Török & Kliem 2005). Flux ropes are ubiquitous, that is, like CMEs, they are also associated with magnetic structures in the heliosphere, streamer blow-outs within the heliospheric current sheet, small-scale structures called plasmoids in the heliosphere, and in solar flares (Nieves-Chinchilla et al. 2023). The CME flux rope is often referred to as a magnetic cloud (MC) or magnetic ejecta (ME), which is characterised by a strong magnetic field, clear rotations in the magnetic field vector, and a low proton temperature (Burlaga et al. 1981; Klein & Burlaga 1982; Kilpua et al. 2017).

The majority of the interplanetary CME reconstruction techniques are built on the principle of a force-free magnetic field configuration of the magnetic cloud. This means, ∇ × B = αB, so that J is parallel to B and J × B vanishes (where B and J are the magnetic field and current density respectively, and α is the force-free constant). This condition implies that the electric current is parallel to the magnetic field only in static conditions. This criterion was incorporated into the first intuitive solution of the cylindrically symmetrical solution with a constant α across the cloud (Lundquist 1951). Another generally assumed characteristic of CME propagation is a self-similar expansion (Poomvises et al. 2010; Davies et al. 2012; Subramanian et al. 2014) in the initial phase of their evolution. However, incorporating force-free and self-similarity characteristics simultaneously was found to be difficult. The cylindrical configuration did not maintain the force-free state once it had begun expanding. Some models changed the geometrical cross-section of the cylinder from circular to elliptical (Hidalgo et al. 2002b; Hidalgo 2003), or kinematically distorted the cylindrical flux-rope geometry during propagation to conserve the force-free assumption during self-similar expansion (Owens et al. 2006; Vandas et al. 2003, 2006). Others proposed non-cylindrical flux-rope models, for example with toroidal geometry (Romashets & Vandas 2003b; Marubashi & Lepping 2007). The orientation of CMEs fitted with torus geometry better matches the in situ observations than the cylindrical geometry and is more consistent with their corresponding source region orientations (Marubashi et al. 2009).

The static CME model must be self-consistently evolved in a realistic solar wind, and its interaction with the large-scale structures in the solar wind and with other transients must be modelled accurately for reliable space weather forecasting. When designing a numerical model of the CME to be propagated in the heliosphere, the model should have the following requirements. It should be consistent with the definition of flux rope emerging from the solar corona and, upon evolving in the heliosphere, it should be able to reproduce the characteristics of the flux ropes observed in the heliosphere. CME evolution and propagation models, such as AWSoM (van der Holst et al. 2014) and CORHEL-CME (Linker et al. 2023), self-consistently erupt the CME from the corona and evolve them to Earth and beyond. However, such modelling is computationally expensive as it covers multiple physical scales and processes. Models like ENLIL (Odstrcil 2003), EUHFORIA (Pomoell & Poedts 2018), SUSANOO-CME (Shiota & Kataoka 2016), MS-FLUKSS (Singh et al. 2018) and so on begin evolving the CMEs from 0.1 au. Although such an initialisation restricts the modelling of the scales and physical processes close to the Sun, it still allows us to understand the heliospheric processes involved in CME propagation. In this work, we use the solar wind and CME evolution model EUropean Heliospheric FORecasting Information Asset (EUHFORIA), which contains functional magnetised FR models. The linear force-free (LFF) spheromak model has shown potential in fitting asymmetric magnetic field observations in situ (Vandas et al. 1991, 1993b,a). Although the spheromak model can fit the central part of the magnetic cloud as well as any flux rope (cylindrical or toroidal models), it is unable to model the ‘edges’ of the cloud (Farrugia et al. 1995). Due to the compact spherical shape and lack of leg-like structures connecting the CME to the surface of the Sun, it is incapable of modelling the flank encounters, that is, CMEs where the Earth (or a satellite) is impacted by one of these legs. Moreover, the spheromak model tilts in the heliosphere to align its symmetry axis with the ambient field to an extent not actually observed in the heliosphere, and therefore leads to erroneous predictions at Earth (Asvestari et al. 2022; Sarkar et al. 2024). To overcome the drawbacks of the spheromak model, the Flux Rope in 3D model (FRi3D; Isavnin 2016) was implemented in EUHFORIA (Maharana et al. 2022). FRi3D has an extended flux-rope geometry with flexible variable cross-sections that incorporate deformations, such as pancaking, flattening, and skewing. This model significantly improved the prediction of CMEs at 1 au compared to spheromak CME predictions (Maharana et al. 2022, 2023; Palmerio et al. 2023). However, the numerical implementation of the deformations and the permanent connection of the CME legs to the Sun gives rise to a higher expenditure of numerical resources. In addition, keeping the legs of one CME attached to the Sun makes it complicated for the numerical injection of the following CME when the legs overlap. Hence, modelling successive CME evolution with the FRi3D model is presently challenging in the framework of EUHFORIA. Therefore, toroidal CME models were developed by Linan et al. (2024) using the analytical magnetic field configurations of the modified Miller-Turner (mMT) and the Soloviev solution of the Grad-Shafranov equation. The aim of this latter work was to simplify the FRi3D geometry while still overcoming the limitations of the spheromak geometry and shortening the wall-clock time of the simulations.

In this study, building on the work by Linan et al. (2024), we focus on validating the modified Miller-Turner (mMT) toroidal CME model and suggest changes to upgrade it into the ‘Horseshoe model’. Due to the multiple non-trivial geometrical parameters of the Soloviev solution, constraining them using observations is not straightforward. Hence, validation of the Soloviev toroidal CME model will be carried out in a future publication. Section 2 introduces the Horseshoe model and its implementation in EUHFORIA. In Sect. 3 we explain the methodologies used to constrain the geometrical and magnetic field parameters of the Horseshoe model from observations. The details of the validation events are provided in Sect. 4. In Sect. 5 we assess the performance of the Horseshoe model and compare it to the other magnetised CME models in EUHFORIA. We summarise our conclusions and discuss our outlook in Sect. 6.

2. Models

In this section we introduce the Horseshoe model, EUHFORIA, and the numerical implementation of the Horseshoe model in EUHFORIA. The comparison of the full torus (Linan et al. 2024) and the Horseshoe implementation of the mMT model is presented.

2.1. Horseshoe CME model

The Horseshoe CME model is a modification of the toroidal CME model introduced by Linan et al. (2024) in EUHFORIA. The geometry of the Horseshoe model is a full torus devoid of the rear part of the torus, and therefore resembles a horseshoe. The rationale behind the new implementation is twofold: (1) It is designed to resemble a realistic flux-rope structure with legs, and (2) to avoid reproducing inaccurate magnetic field configuration in the heliosphere because a real CME does not have the geometry of a full torus. The magnetic field topology, as in Linan et al. (2024), is defined using the modified Miller-Turner model (mMT, Romashets & Vandas 2003b):

$$\begin{array}{cc}\hfill B_{{\rho }_l}& =B_0\frac{R_0-2{\rho }_{l}cos{\theta }_l}{2\alpha R_0(R_0+{\rho }_{l}cos{\theta }_l)}{J}_0(\alpha {\rho }_l)sin{\theta }_l,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill B_{{\varphi }_l}& =B_0(1-\frac{{\rho }_l}{2R_0}cos{\theta }_l)J_0(\alpha {\rho }_l),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill B_{{\theta }_l}& =B_0\frac{R_0-2{\rho }_{l}cos{\theta }_l}{2\alpha R_0(R_0+{\rho }_{l}cos{\theta }_l)}{J}_0(\alpha {\rho }_l)cos{\theta }_l\hfill \\ \hfill & -B_0(1-\frac{{\rho }_l}{2R_0}cos{\theta }_l)J_1(\alpha {\rho }_l),\hfill \end{array}$$(3)

where (ρ_l, ϕ_l, θ_l) are the local toroidally curved cylindrical coordinates (for more information, see Fig. 5 in Linan et al. 2024). Here, R₀ and a are the major and minor radii of the torus, respectively, while B₀ is the axial magnetic field strength. In the local coordinate system, ρ_l extends from [0, a], and ϕ_l and θ_l both cover the [0, 2π] interval. J_m is the Bessel function of m-th order and α is the force-free constant. The advantage of this magnetic field topology is its fully analytical form of the force-free magnetic field (∇ × B = αB) inside the torus. The magnetic field topology is axisymmetric about ϕ_l. To confine the magnetic field inside the torus with a circular cross-section, the magnetic field must be completely poloidal at the torus boundary. B_ρl and B_ϕl are maximal at the centre of the torus cross-section and monotonically drop to zero when αρ_l ≈ 2.405, which gives the first zero of J₀. The parameter α is related to the flux rope minor radius and chirality (C) by:

$$\begin{array}{c}\hfill \alpha =C\frac{2.405}{a}.\end{array}$$(4)

The twist varies from zero at the centre to infinity at the outer boundary of the torus. The mMT configuration is a generalised version of the Miller-Turner topology (MT, Miller & Turner 1981), as it satisfies the exact solenoidality (∇ ⋅ B = 0) for all aspect ratios (R₀/a) of the torus, while MT is solenoidal only for higher aspect ratios (> 10), and coincides asymptotically with the mMT and Lundquist (1951) solutions in that regime. For low aspect ratios (< 3), the magnetic field profiles for the mMT model better match the numerical solutions. Both MT and mMT solutions are approximately force-free for smaller aspect ratios, but the MT solution performs better than the mMT (Vandas & Romashets 2015). The interplanetary flux ropes fitted at 1 au with toroidal CME models point to an aspect ratio of greater than 3 (Marubashi et al. 2015). However, the flux ropes close to the Sun or the locally distorted loop-like flux ropes (Vandas et al. 2002) can have aspect ratios of less than 2 (Romashets & Vandas 2003a). The CMEs initialised at 0.1 au have a low aspect ratio, and therefore we can use both the MT and mMT solutions in our numerical modelling. Confirming our claim, Linan et al. (2024) demonstrated the numerical stability of using tori (with the mMT configuration) of low aspect ratios (between 1 and 3) in the framework of EUHFORIA.

Previous studies (Farrugia et al. 1995; Mulligan & Russell 2001; Hidalgo et al. 2002a; Möstl et al. 2009) suggest that the force-free assumption is mostly invalid, and cannot fully explain the pressure gradients in the core and boundary of the interplanetary CMEs. Similarly, Isavnin (2016) does not consider the non-force-free nature of the FRi3D model as a disadvantage. Although mMT is not force-free for all aspect ratios, it better satisfies the divergence-free condition, and so we consider its adoption in our study. Romashets & Vandas (2003b) suggest the validity of the mMT solution for locally treating a part of an extended flux rope rooted at the Sun as a toroid (a torus is a toroid with a circular cross-section). Hence, even if we deviate from the full torus geometry in the novel horseshoe model, the magnetic field topology of mMT is still applicable.

2.2. EUHFORIA

EUHFORIA is a physics-based, data-driven magnetohydrodynamic (MHD) model of solar wind and CME evolution used for both space weather research and forecasting. It has a modular design, enabling the combination of different coronal and heliospheric models. The coronal model extends from the photosphere to the low corona up to 0.1 au and provides the boundary condition for the heliospheric domain. The heliospheric domain normally covers the region between 0.1 au and 2 au (expandable) and solves the MHD equations on the 3D grid. We use the semi-empirical, modified Wang–Sheeley–Arge model (WSA, Arge et al. 2004; Pomoell & Poedts 2018) as the coronal model. The WSA model incorporates the potential field source surface (PFSS, Altschuler & Newkirk 1969) extrapolation initialised by photospheric magnetograms from sources; for example, the Global Oscillation Network Group (GONG, Harvey et al. 1996) and SDO’s Helioseismic and Magnetic Imager (HMI, Schou et al. 2012) in the low corona up to 2.6 R_⊙. This is followed by the Schatten Current Sheet (SCS, Schatten et al. 1969) model, which then radially extends the magnetic field lines up to 21.5 R_⊙ (0.1 au). Finally, the plasma properties (speed, density, temperature) are computed at 0.1 au using empirical functions of the flux tube areal expansion factor and the distance of the foot point of the open field lines to the coronal hole boundary (more details are provided in Pomoell & Poedts 2018; Asvestari et al. 2019). After the solar wind relaxation phase (i.e. filling the heliospheric domain with steady co-rotating solar wind plasma and magnetic field), CMEs are injected by means of time-dependent boundary conditions at 0.1 au and are then self-consistently evolved in the heliosphere by solving the MHD equations. In this work, the simulations with the Horseshoe model in EUHFORIA are performed using a spherical grid in the Heliocentric Earth EQuatorial (HEEQ) coordinate system, with a spatial resolution of 1.6 R_⊙ in the radial direction (0.1–2 au), and 4° angular resolution in the latitudinal (−70° to 70°) and longitudinal (0°–360°) directions. In this work, the resolution of the EUHFORIA time series output is 10 minutes.

2.3. Numerical implementation of the Horseshoe model in EUHFORIA

The Horseshoe CME is injected into the heliospheric domain of EUHFORIA as a time-dependent boundary condition starting from a predetermined initiation time (based on observations of the CME speed). There are two layers of ghost cells before the actual computational domain of the heliosphere. The CME values are first assigned at the interface between the ghost cells and the in-domain cells of the heliosphere and are then interpolated to the in-domain cells (see details in Pomoell & Poedts 2018). The implementation follows the same methodology as the torus models as detailed in Sect. 3.2 of Linan et al. (2024). A mask function is defined to identify the grid cells on the inner boundary (a spherical surface at 0.1 au), which intersect with the CME upon its insertion. Then, ρ_l is computed, and all points on the boundary – where this radius is less than the minor radius of the torus – are considered part of the CME. As the magnetic field for the mMT configuration within the torus is defined in the local curved cylindrical coordinates (ρ_l, ϕ_l, θ_l), the magnetic field is derived in the local spherical system (ρ, θ, ϕ), where ρ is the radius, θ is the co-latitude, and ϕ is the longitude. The two sets of coordinates are related to each other by the relations:

$$\begin{array}{cc}\hfill {\rho }_l& =\sqrt{{\rho }^2+R_0^2-2\rho R_{0}sin\theta }\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill {\theta }_l& =\{\begin{array}{cc}{sin}^{-1}\frac{\rho cos\theta }{{\rho }_l},\hfill & \text{if}\rho sin\theta -R_0\ge 0\hfill \\ \pi -{sin}^{-1}\frac{\rho cos\theta }{{\rho }_l},\hfill & \text{otherwise}.\hfill \end{array}\hfill \end{array}$$(6)

More details on the coordinate systems and the transformations can be found in the Appendix of Linan et al. (2024).

Instead of pushing the full torus of radial size 2(a + R₀), we cut the torus at 1.98 R₀, which is slightly less than 2 R₀ (a + R₀ (front half of the torus) + R₀ − a (from centre until the last half of the torus)) in order to exclude the injection of the trailing part of the torus. The modified torus of size 1.98 R₀, defined as the Horseshoe model, is illustrated schematically in Fig. 1. At each time step, the centre of the torus is advanced with a uniform radial speed and the 3D mask region is updated. In the mask region, the solar wind magnetic field values are substituted with the mMT magnetic field. We incorporate a uniform mass density and temperature in the Horseshoe model, similar to the full torus model. We performed a comparison of the EUHFORIA profiles at 1 au for the full torus and the Horseshoe models. We launched the CME along the Sun–Earth line with zero tilt to get a simplified magnetic field profile at Earth. The time series plot in Fig. 1 shows the contrasting features in the magnetic field components. The double polarity feature resulting from the front and the back part of the full torus is reflected in the B_y component. The 3D visualisation in Panels 1–3 of Fig. 2 shows the evolution of the full torus in the equatorial plane, 4 and 11 hours after the start of the CME injection into the heliospheric domain, respectively. The magnetic field lines are colour-coded with the total magnetic field strength, hence distinguishing the background solar wind from the strong field inside the CME. The divergence of the velocity (∇ ⋅ v) is plotted in the background to show the flow of plasma – the blue spectrum corresponds to the accumulation of matter in the sheath ahead of the CME, and the red region inside the CME depicts the outflow of mass. Due to the uniform injection speed, the torus is deformed during the injection. More precisely, the torus gets flattened, and the rear part becomes bean shaped. Moreover, the rear part propagates much faster than the front part as the front part cleans up the background wind so that the density behind it is much lower. Consequently, the Alfvén speed increases and the rear part evolves faster than the front part and catches up with it. These features suggest the possibility of erroneous space weather predictions if a full torus model is used. Panels 2–4 of Fig. 2 depict the evolution of the Horseshoe CME and the mitigation of the double B_y polarity signature. The B_z profile at 1 au in the full torus has a shorter duration as a result of the trailing part of the torus compressing the front part and merging with it, whereas the CME expands more freely in case of the Horseshoe.

Fig. 1. Horseshoe model implementation. Left: schematic representation (not-to-scale) of the Horseshoe model propagating into the EUHFORIA heliosphere domain in the equatorial plane (x − y); Right: comparison of the EUHFORIA predictions at 1 au between the Horseshoe and the full torus geometry with the mMT magnetic field configuration. The red arrow in the By profile shows the double polarity signatures of the full torus. The Bz profile of the full torus is compressed and merged with the trailing part of the torus.

Fig. 2. Equatorial snapshot of EUHFORIA simulations of the full torus and the Horseshoe implementation of the mMT magnetic field configuration. The evolution of the internal magnetic field lines of the CMEs (colour-coded with the magnetic field strength), 4 hours (panels 1–2), and 11 hours (panels 3–4), respectively, after the start of the injection at 0.1 au, is illustrated. In the background, the divergence of the speed is plotted. The sheath ahead of the magnetic cloud is depicted by the negative divergence (accumulation of plasma) region, showing a clear envelope around the CME. In panels 1 and 3, the rear part of the full torus is seen injected, whereas panels 2 and 4 show the horseshoe-like geometry creating CME leg-like structures connected to the inner boundary.

In the current version of the Horseshoe model in EUHFORIA, we also implemented the MT magnetic field topology. MT is as stable as mMT, and qualitatively, it gives similar results for low aspect ratios. We also noticed that the aspect ratios obtained from the 3D reconstruction of CMEs from white-light images within 0.1 au are usually low (< 2). In that case, both mMT and MT will perform similarly. However, due to the more general divergence-free nature of the mMT solution, we further develop the methodologies to constrain the magnetic field parameters and optimise our simulations for real events for the mMT implementation in the Horseshoe geometry.

3. Constraining the CME parameters

3.1. Magnetic field parameters

To provide an early warning of CME impact, we constrain the CME parameters from remote sensing observations of the photosphere and corona before, during, and after the time of eruption. Estimating these parameters is not straightforward for CMEs without low-coronal signatures, such as the streamer blow-out CMEs (Robbrecht et al. 2009). However, for the classical CMEs associated with flares, the source region features allow estimation of the magnetic configuration of the flux rope (Hudson & Cliver 2001; Palmerio et al. 2017). The magnetic field input parameters required for the initialisation of the CME at 0.1 au are the axial magnetic field strength (B₀), the chirality (C), and the orientation of the magnetic axis (tilt).

Axial magnetic field strength. The field strength is quantified through the amount of the magnetic flux released during the eruption. Hence, we determine the toroidal and poloidal flux as a function of B₀. The toroidal flux, ϕ_t, is the magnetic flux passing across the cross-section of the torus. Considering flux conservation and substituting with the mMT magnetic field equations from Eqs. (1)–(3),

$$\begin{array}{cc}\hfill {\phi }_t& ={\displaystyle {\int }_S}\mathit{B}·\mathrm{d}\mathit{S}={\displaystyle {\int }_0^a}{\displaystyle {\int }_0^{2\pi }}{B}_{{\varphi }_l}{\rho }_l\mathrm{d}{\rho }_l\mathrm{d}{\theta }_l\hfill \end{array}$$(7)

$$\begin{array}{cc}& ={\displaystyle {\int }_0^a{\int }_0^{2\pi }}{B}_0(1-\frac{{\rho }_{l}cos{\theta }_l}{2R_0})J_0(\alpha {\rho }_l){\rho }_l\mathrm{d}{\rho }_l\mathrm{d}{\theta }_l,\hfill \end{array}$$(8)

where S is the cross-section of the torus and magnetic field B is integrated over the crossing area. The radius of the surface S is the minor radius of the torus (a). Here, $\alpha =\pm \frac{a_0}{a}$, where a₀ is the first root of J₀, that is, J₀(a₀) = 0 ⟹ a₀ = 2.41 a. The sign of α represents the chirality of the flux rope. Inverting the above integral yields

$$\begin{array}{c}\hfill B_0=\frac{{\phi }_t\alpha }{2\pi a{J}_1(\alpha a)}.\end{array}$$(9)

The poloidal flux, φ_p, is computed by integrating over the magnetic field lines crossing the θ_l = 0 plane of the torus. Hence, we integrate B_θl (θ_l = 0) over r and ϕ:

$$\begin{array}{cc}\hfill {\phi }_{\mathrm{p}}& =\frac{1}{2\pi }{\displaystyle {\int }_0^a{\int }_0^{2\pi }}{B}_{{\theta }_l}{\rho }_l\mathrm{d}{\rho }_l\mathrm{d}{\varphi }_l\hfill \end{array}$$(10)

$$\begin{array}{cc}& =B_0{\displaystyle {\int }_0^a}(\frac{R_0-2{\rho }_l}{2\alpha R_0(R_0+{\rho }_l)})J_0(\alpha {\rho }_l)\hfill \\ \hfill & -(1-\frac{{\rho }_l}{2R_0})J_1(\alpha {\rho }_l){\rho }_l\mathrm{d}{\rho }_l.\hfill \end{array}$$(11)

This integral can be computed numerically. The value of φ_p is obtained from observations, and hence the value of B₀ can be obtained by inverting Eq. (10). The values of the geometrical parameters a and R₀ – obtained when the torus is self-similarly expanded up to 0.1 au – are substituted in Eq. (10).

One of the fastest and most accessible methods to obtain ϕ_p from observations is using the empirical relation, that is, using the peak X-ray intensity of the flare. Kazachenko et al. (2017) from a catalogue of 3000 flares (RibbonDB catalogue¹) found a correlation between the reconnection flux (ϕ_r in Mx) and the flare peak intensity (I_SXR, in W m⁻²). The relation obtained by these authors was log₁₀(ϕ_r) = 24.42 + 0.64log₁₀(I_SXR), with a Spearman’s rank correlation coefficient of 0.66. We set the value of ϕ_p to half of the total unsigned ribbon reconnection flux provided by the catalogue. The catalogue of Kazachenko (2023) provides the observed ϕ_p (hereafter, ϕ_p, o) and the error associated with each event (Δϕ_p, e) for 479 flares. The relative error, RE = Δϕ_p, e/ϕ_p, o × 100, of the sample set ranges between 7 and 60%. The average RE (RE_avg) from the catalogue is around 23%. We derive the absolute average error, Δϕ_p, avg, from RE_avg = Δϕ_p, avg/ϕ_p, o × 100. While performing EUHFORIA simulations, the error estimates will be considered for making ensemble runs. The flares associated with the validation events (Sect. 4) considered in this study are available in the RibbonDB catalogue. Hence, we use the observed poloidal flux to constrain B₀.

Chirality. With this parameter, we incorporate the sign of the helicity. The orientation of the magnetic flux rope can be interpreted from the hemisphere rule (Pevtsov & Balasubramaniam 2003) or EUV/soft-X-ray sigmoids (Titov & Démoulin 1999). A comprehensive list of different observational proxies of chirality can be found in Palmerio et al. (2017). A right-handed sigmoid is associated with a positive chirality and a left-handed one with a negative chirality.

Tilt. The flux-rope orientation is generally inferred to be parallel with the polarity inversion line (PIL). We define the east–west pointing flux rope – placed parallel to the solar equator – as the reference zero tilt for the Horseshoe model. The acute angle is assigned a negative (positive) sign when rotated anticlockwise (clockwise).

3.2. Geometry

The geometrical parameters for the Horseshoe model can be constrained equivalently to the full torus, as its geometry is not altered during the numerical implementation; rather, only a part of the full torus is injected. The simplest CME geometry is generally defined using the angular half-width (φ_hw) and angular half-height (φ_hh), which define the maximum extent of the CME (relative to the plane containing its magnetic axis) in the azimuthal and polar directions, respectively. This approach has been adopted in 3D reconstruction techniques of the Graduated Cylindrical Shell (GCS, Thernisien et al. 2006) and Flux Rope in 3D (FRi3D, Isavnin 2016) models. For constraining the geometrical parameters to initialise CMEs in heliospheric models like EUHFORIA, the 3D reconstruction of a CME is performed as long as the CME structures can be identified distinctly in white-light coronagraph images. Multi-point observations provide the view of the CME in the upper corona until somewhere below 21.5 R_⊙. Such data are currently available from the C2 and C3 instruments of the Large Angle and Spectrometric COronagraph (LASCO, Brueckner et al. 1995) on board the Solar and Heliospheric Observatory (SOHO), and the COR-2 instrument on board the Sun–Earth Connection Coronal and Heliospheric Investigation (SECCHI) package of the twin-spacecraft Solar Terrestrial Relations Observatory (STEREO, Kaiser et al. 2008), the Metis instrument (Antonucci et al. 2020) on the Solar Orbiter (Müller et al. 2020), and others. Therefore, we determine the geometry of the torus assuming a self-similar expansion of the CME beyond the last observed time frame until 21.5 R_⊙. We inject the already expanded torus directly at 21.5 R_⊙. The time of arrival of the CME at 21.5 R_⊙ is calculated using the speed obtained in the last frame and assuming a uniform propagation of the CME up to 21.5 R_⊙. In this section, we define two geometries (Case 1 and Case 2) for establishing the association between the torus parameters a, R₀, and the centre of the torus (T_c) with the observed angular CME parameters, φ_hw and φ_hh.

• Case 1: This geometry is represented by the red curve in the equatorial and meridional planes as shown in Figs. 3a and b, respectively, when T_c = 21.5 R_⊙. The φ_hw and φ_hh are defined as

  $$\begin{array}{c}\hfill sin({\phi }_{\mathrm{hw}})=\frac{R_0+a}{T_{\mathrm{c}}},\end{array}$$(12)

  and

  $$\begin{array}{c}\hfill tan({\phi }_{\mathrm{hh}})=\frac{a}{T_{\mathrm{c}}+R_0}.\end{array}$$(13)

  Fig. 3. Torus geometry in the equatorial and meridional planes for Case 1 (a, b) and Case 2 (c, d), where geometries are represented with red curves. The horizontal black solid line is the Sun–Earth line. The blue lines on the equatorial and meridional plots have angular extents of φhw and φhh, respectively, about the Sun–Earth line. In this example, φhw = 40° and φhh = 12°. In Case 1, R0 = 7.6 and a = 6.2, making the aspect ratio 1.23. In Case 2, R0 = 11.1 and a = 6.9, making the aspect ratio 1.6.

  Using Eqs. (12) and (13), we derive R₀ and a:

  $$\begin{array}{c}\hfill R_0=T_{\mathrm{c}}\frac{sin({\phi }_{\mathrm{hw}})-tan({\phi }_{\mathrm{hh}})}{1+tan({\phi }_{\mathrm{hh}})},\end{array}$$(14)

  and

  $$\begin{array}{c}\hfill a=(T_{\mathrm{c}}+R_0)tan({\phi }_{\mathrm{hh}}).\end{array}$$(15)

• Case 2: This geometry is represented by the red curve in the equatorial and the meridional planes as shown in Figs. 3c and d, respectively, when T_c = 21.5 R_⊙. The φ_hw and φ_hh are defined as

  $$\begin{array}{c}\hfill tan({\phi }_{\mathrm{hw}})=\frac{R_0+a}{T_{\mathrm{c}}},\end{array}$$(16)

  and

  $$\begin{array}{c}\hfill tan({\phi }_{\mathrm{hh}})=\frac{a}{T_{\mathrm{c}}+R_0}.\end{array}$$(17)

  Using Eqs. (16) and (17), we derive R₀ and a:

  $$\begin{array}{c}\hfill R_0=T_{\mathrm{c}}\frac{tan({\phi }_{\mathrm{hw}})-tan({\phi }_{\mathrm{hh}})}{1+tan({\phi }_{\mathrm{hh}})},\end{array}$$(18)

  and

  $$\begin{array}{c}\hfill a=(T_{\mathrm{c}}+R_0)tan({\phi }_{\mathrm{hh}}).\end{array}$$(19)

Figure 4 shows the distribution of R₀ and a of the torus as a function of φ_hw and φ_hh for the two geometries, and the allowed parameter space where the torus exists. The distribution of the R₀ in the φ_hw − φ_hh space is shown in Figs. 4a,b for the two geometries. R₀ becomes negative when sin(φ_hw) < tan(φ_hh) for the Case 1 geometry, and when φ_hw < φ_hh for the Case 2 geometry. The Case 2 geometry results in a comparatively larger equatorial cross-section, resulting in a higher aspect ratio. With Case 1, we get R₀ < 21.5 R_⊙ even for wide CMEs, whereas Case 2 estimates a very large R₀ that cannot be injected at the EUHFORIA inner boundary. The distribution of a (Figs. 4c,d) is similar for both geometries for the low φ_hw − φ_hh space. However, for the values in the higher extreme of the parameter space, Case 2 estimates higher values of a as compared to Case 1. The next step is to ensure the allowed pairs of (a, R₀) such that the torus is physical, that is, R₀ > a. This criterion is depicted in the boolean plots as shown in Figs. 4e and 4f for Case 1 and Case 2, respectively. The white space (True) in the plots corresponds to the allowed (φ_hw, φ_hh) values for both geometries. The boundary separating the allowed and non-allowed values can be fit with a univariate cubic spline (piecewise polynomial function of degree 3) plotted in yellow. This diagnostic shows that the maximum φ_hh with Case 1 geometry is less than 20° even for a CME of as wide as 70°. Hence, with this geometry, we can construct a thin cross-section of CMEs, which could imply a potential underestimation of CME mass. On the other hand, with Case 2 geometry, the range of obtained maximum φ_hh is higher as compared to Case 1. We can say that the Case 2 geometry performs better at estimating the CME cross-section (a), although it overestimates R₀. Therefore, depending on the observed event, either Case 1 or Case 2 can be chosen to constrain (a, R₀) from (φ_hw, φ_hh). In general, we require the strong condition 2(a + R₀) < 21.5 R_⊙ to hold, which in turn obliges us to constrain a small a for high-aspect-ratio cases.

Fig. 4. Contours of major radius (R0) as a function of φhw and φhh for (a) Case 1 and (b) Case 2 geometries. Contours of the minor radius (a) as a function of φhw and φhw for (c) Case 1 and (d). The parameter space where R0 > a is identified for (e) Case 1 and (e) Case 2, geometries, respectively.

The distribution of the aspect ratio in the φ_hw − φ_hh parameter space in Fig. 4 is plotted in Figs. 5a,b. The black dashed and solid curves are the cubic spline defining the boundary for the physical torus parameters for Case 1 and Case 2 geometries, respectively (same as in Figs. 4e,f). These figures suggest the range of aspect ratios that can be allowed in both geometries. Moreover, it can be observed that Case 2 favours higher aspect ratios for φ_hw = 20° as compared to Case 1. The volume of the Horseshoe geometry is computed as 3/4 ⋅ 2π²R₀a². Figure 5 illustrates the extent of the torus volume for the two geometries in φ_hw − φ_hh space. The colour map covers the full range of positive volume. The order of magnitude of the volume is in the range 10²⁷ − 10³¹ m³ for Case 1 and is up to 10³² m³ for Case 2, corresponding to the range modelled by FRi3D and spheromak models in EUHFORIA as discussed in Maharana et al. (2022). The volume of the FRi3D model (with a flexible and extended flux-rope geometry) belongs to the range 10²⁹ − 10³¹ m³. The spheromak model can also be modelled using a sine geometry and a tan geometry (Scolini et al. 2019), similar to the Case 1 and Case 2 geometries defined for the torus model in our case. The range of volume is 10²⁹ − 10³¹ m³ and 10²⁹ − 10³² m³ for the sine and tan geometries, respectively. The Horseshoe volume profiles of Case 1 and Case 2 are similar up to φ_hw = 20°, beyond which the Case 2 geometry gives a higher estimate. Depending on the volume estimated from observations of a CME, the choice of geometry and therefore the volume can be constrained accurately. Using the average CME mass density of 10⁻¹⁷ kg m⁻³ at 21.5 R_⊙ as suggested by Temmer et al. (2021), the torus mass for Case 1 geometry lies in the range 10¹⁰ − 10¹⁴ kg and that for Case 2 lies within 10¹⁰ − 10¹⁵ kg.

Fig. 5. Distribution of aspect ratio (a) Case 1 and (b) Case 2 geometries and the volume (c) Case 1 and (d) Case 2 geometries. The colour bars show the magnitude of the aspect ratio and the volume in the respective plots. The black dashed and the solid curves demarcate the allowed range of φhw − φhh space for Case 1 and Case 2 geometries, respectively.

Statistical studies (Marubashi & Lepping 2007; Marubashi et al. 2015) that reconstructed in situ observations of magnetic clouds at 1 au with toroidal magnetic field configurations suggest CMEs with high aspect ratios greater than 5 (the sample average aspect ratio is ∼9) as compared to the estimates from cylindrical configuration fitting. It must be noted that these reconstruction techniques estimate the global flux-rope properties based on the localised fitting of a part of it. Hence, errors associated with the aspect-ratio estimates are possible. We note that due to the expansion of the flux rope during propagation, the aspect ratio can increase. Vandas et al. (2002, 2003) used an aspect ratio of 3 − 4 for their toroidal MT CME launched into the solar wind in MHD simulations, which is not as high as the estimates obtained from the in situ reconstructions, and obtain a reasonable magnetic field at 1 au. In EUHFORIA simulations, the CME injection happens at 21.5 R_⊙, and considering the mentioned studies, it is reasonable to choose lower aspect ratios to constrain the Horseshoe model close to the Sun.

3.3. Speed

Applying a 3D reconstruction model to the white light images, we fit the φ_hw, φ_hh and the leading edge (LE; the bright front):

$$\begin{array}{c}\hfill \mathrm{LE}=T_{\mathrm{c}}+R_0+a,\end{array}$$(20)

where T_c denotes the centre of the torus. To obtain the total speed (v_3D) at the CME apex, which is the sum of the radial speed, v_rad (rate of change of the T_c), and the expansion speed, v_exp (rate of change of the torus cross-section), we take the time derivative of the above equation and substitute Eqs. (14) and (15) for R₀ and a,

$$\begin{array}{cc}\hfill v_{3\mathrm{D}}=\frac{\mathrm{dLE}}{\mathrm{d}t}& =\frac{\mathrm{d}{T}_{\mathrm{c}}}{\mathrm{d}t}+\frac{\mathrm{d}{R}_0}{\mathrm{d}t}+\frac{\mathrm{d}a}{\mathrm{d}t}\hfill \end{array}$$(21)

$$\begin{array}{cc}& =(\frac{\mathrm{d}{T}_{\mathrm{c}}}{\mathrm{d}t}+\frac{\mathrm{d}{R}_0}{\mathrm{d}t})(1+tan({\phi }_{\mathrm{hh}}))\hfill \end{array}$$(22)

$$\begin{array}{cc}& =\frac{\mathrm{d}{T}_{\mathrm{c}}}{\mathrm{d}t}(1+sin({\phi }_{\mathrm{hw}})).\hfill \end{array}$$(23)

The radial speed is given by

$$\begin{array}{c}\hfill v_{\mathrm{rad}}=\frac{\mathrm{d}{T}_{\mathrm{c}}}{\mathrm{d}t}=\frac{v_{3\mathrm{D}}}{1+sin({\phi }_{\mathrm{hw}})}.\end{array}$$(24)

4. Validation events

In this work, we validate the Horseshoe model with two events. For each event, we present a brief introduction, the details of the parameters constrained by the observations, and the results of the EUHFORIA simulations. The disturbance storm time (Dst) index is computed using the empirical AK2 model of O’Brien & McPherron (2000a,b) on the plasma and magnetic field properties at 1 au. The Dst model itself has some associated errors; hence, we first apply it to the observed data for each event. This serves as the reference modelled Dst. We then apply the Dst model to the results of the EUHFORIA simulations and assess the performance of EUHFORIA results in predicting Dst as compared to the reference modelled Dst.

4.1. Event 1: 12 July 2012

The aim is to validate the ability of the Horseshoe model to predict the arrival time and magnetic field profile caused by the impact of a CME at 1 au. The first studied event is an isolated, non-interacting Earth-directed CME. This textbook event of 12 July 2012, which was triggered by an X1.4 flare from NOAA active region (AR) 11520 located at S17E06, is one of the most popular studied geo-effective events (Hu et al. 2016; Webb & Nitta 2017; Gopalswamy et al. 2018). The interplanetary propagation of this CME was previously modelled in the framework of EUHFORIA using the magnetised CME models – spheromak (Scolini et al. 2019) and the FRi3D (Maharana et al. 2022). A shock driven by the fast CME (with an average projected speed of 890 km s⁻¹) arrived at Earth on 14 July 2012 at 17:39 UT. The magnetic cloud signatures were recorded between 15 July 2012 at 06:00 UT and 17 July 2012 at 05:00 UT, with a temporally long negative B_z signature (minimum B_z = −18 nT). The prolonged southward B_z resulted in an intense storm with the main phase starting around 10.00 UT on 15 July 2012 and maintained a Dst < −100 nT until around 9.00 UT on 16 July 2012 (minimum Dst = −122 nT on 15 July). The observed in situ signatures are plotted in Figs. 6 and 7 along with the simulation results.

Fig. 6. Results of Horseshoe ensemble simulations for Event 1 with varying ϕp obtained using EUHFORIA. Panels from top to bottom: speed (v), proton number density (np), Bx, By, Bz, magnetic field strength (|B|), plasma beta (β), and Dst. The observations from the WIND spacecraft are plotted in black in all the plasma and magnetic field properties panels. In the Dst panel, the observed Dst from WDC, Kyoto, is plotted in red, the reference modelled Dst using the observed data in black, and the Dst using EUHFORIA simulations in other colours. The solid black vertical line marks the shock arrival time, and the dashed magenta lines show the start and end of the magnetic cloud as reported in the WIND ICME catalogue (https://wind.nasa.gov/ICME_catalog/ICME_catalog_viewer.php).

Fig. 7. Results of the best Horseshoe ensemble simulation of Event 1 obtained using EUHFORIA at Earth (solid blue line) and at virtual satellites in the 5 − 10° latitudinal and longitudinal offset around Earth (shaded region). All other plot details are similar to Fig. 6. The grey bar at the top of the Bz panel depicts the time span of the dynamic time warping analysis in Sect. 5.

4.1.1. Observationally constrained parameters for the EUHFORIA simulation

The 3D reconstruction of the CME, performed using the FRi3D reconstruction tool (Isavnin 2016), is provided in Appendix B.1. Applying the geometrical transformations given by Eqs. (14) and (15) to the fitted φ_hw = 48° and φ_hh = 11°, the a and R₀ are obtained as 6 R_⊙ and 10 R_⊙, respectively. Temmer et al. (2021) estimate the mass and density of the CME in Event 1 as 1.84 ⋅ 10¹³ kg and 1.75 ⋅ 10⁻¹⁷ kg m⁻³, respectively, resulting in a volume of 3.2 ⋅ 10³⁰ m³. Using the constrained R₀, a, and the uniform density 10⁻¹⁷ kg m⁻³, the volume of the CME modelled by the Horseshoe model is 1.8 ⋅ 10¹³ kg, which is in agreement with the observed estimate. The reconstruction provides us with the v_3D, which is converted to v_rad using Eq. (24). From the fitting, the geometrical inclination (unsigned tilt) is 51°, which is consistent with the PIL orientation at the source region. The magnetic field properties of the CME could be extracted given the clear association between the flare and the CME. The chirality is right-handed and the polarity is west to east. The direction of the axial magnetic field by incorporating the polarity and chirality is 129° for the Horseshoe model. The signed φ_p is 4.32 ⋅ 10²¹ Mx with an error of 0.66 ⋅ 10²¹ Mx (15%). The B₀ is computed using the observed ϕ_p value in the inverted Eq. (10). All the parameters used for the EUHFORIA simulation using the Horseshoe model are reported in Table 1.

4.1.2. EUHFORIA simulations

We created an ensemble of simulations based on the varying B₀ obtained from ϕ_p, o, ϕ_p, o ± Δϕ_p, e and ϕ_p, o ± Δϕ_p, avg. For Event 1, ϕ_p, o ± Δϕ_p, e and ϕ_p, o ± Δϕ_p, avg are associated with RE = 15% and RE_avg = 23%, respectively. The results of the ensemble run, in comparison with observations (plotted with the black curve), are presented in Fig. 6. The first two panels show the speed (v) and proton number density (n_p), which are qualitatively modelled in a similar way for all the ensembles with a difference in the arrival time of the shock, characterised by the sharp peak in v. Upon interaction with the solar wind, a non-zero Lorentz force (j × B) can develop due to local misalignment between the currents and magnetic fields within the evolving flux rope. The higher the magnetic flux, the higher the Lorentz force, leading to enhanced CME acceleration and expansion and earlier arrival time at Earth (Subramanian & Vourlidas 2009). Panels 3, 4, 5, and 6 of Fig. 6 show B_x, B_y, B_z, and the total magnetic field (|B|). The magnetic field components show a similar behaviour to the observations, that is, they match qualitatively. The magnetic field profiles show not only the difference in arrival time of the ensemble runs but also the expansion of the magnetic cloud. The distinction between the sheath and the magnetic cloud regions can be made using the plasma beta (β, shown in panel 7), that is, the ratio between the plasma pressure and the magnetic pressure. The sheath region is characterised by β ≫ 1, and the magnetic cloud boundary begins when the β starts falling towards 1, and it stretches through the β < 1 region. Finally, panel 8 shows the predicted Dst computed with the empirical model (O’Brien & McPherron 2000a; O’Brien & McPherron 2000b) using the EUHFORIA output. The observed Dst provided by the World Data Center for Geomagnetism, Kyoto², is plotted in red, and the reference modelled Dst (computed using observations) is plotted in black. The coloured lines represent the results of the ensemble runs, as mentioned in the legend.

The simulation with ϕ_p, o + Δϕ_p, avg showed the best performance in predicting the arrival time 14 July 2012 at about 18:13 UT, which is only 34 minutes later than the observed shock arrival, and with an accurately predicted magnetic field profile as well. The most accurately predicted simulation time profiles of the Horseshoe model are plotted in Fig. 7. Virtual spacecraft are placed at an offset of latitude and longitude around Earth (σ_θ, ϕ) with a step of 5° in order to capture the variability of the plasma and magnetic field properties in the vicinity of Earth. The profiles at the location of these virtual spacecraft are plotted in shades of dark and light blue for σ_θ, ϕ = ±5° and σ_θ, ϕ = ±10°, respectively. The percentage difference in minimum B_z is computed as

$$\begin{array}{c}\hfill \mathrm{\Delta }\min (B_z)=\frac{\min (B_{z,\mathrm{euh}})-\min (B_{z,\mathrm{obs}})}{\min (B_{z,\mathrm{obs}})}\times 100,\end{array}$$(25)

where B_z, euh and B_z, obs are the EUHFORIA-simulated and in situ observed time profiles, respectively.

The best Horseshoe simulation results in Δmin(B_z) = − 35%, which means the simulated minimum B_z is underestimated by 35% compared to observations. In the vicinity of Earth (at all the virtual spacecraft), Δmin(B_z) lies in the range [−60%, − 15%]. For the comparison of the results from different ensemble runs of Event 1, the Δmin(B_z) at Earth and the virtual spacecraft are plotted for each ensemble run in Fig. 8a. The red dot and the orange line in boxes represent the Δmin(B_z) at the location of Earth and the median Δmin(B_z) of the distribution, respectively. These values lie within the interquartile range of the data for all the ensemble runs. The distribution of the absolute median Δmin(B_z) reduces with increasing B₀, which is consistent with the result that increasing ϕ_p leads to more negative B_z values (closer to the observations) in the case of Event 1. The absolute median Δmin(B_z) obtained with the poloidal fluxes ϕ_p − Δϕ_p, avg and ϕ_p − Δϕ_p, e are lower than that of ϕ_p, o. This is an exception to the trend whereby the absolute Δmin(B_z) decreases with increasing ϕ_p. However, the whiskers (extent of the farthest data point from the box) of those two cases still lie within the range of whiskers of the ϕ_p, o case. This implies that the predicted B_z values at different virtual spacecraft were more localised around the predicted B_z at Earth for those two cases and within the uncertainty range of ϕ_p, o, and are therefore consistent with the trend mentioned above.

Fig. 8. Boxplot of the spread of the relative minimum Bz (Δmin(Bz)) values at Earth and the virtual spacecraft at an offset of latitude and longitude of 5–10° – (a) Event 1, (b) Event 2. The red dot and the orange line in the boxes represent the Δmin(Bz) at the location of Earth and the median Δmin(Bz) of the distribution, respectively. The dashed blue line indicates Δmin(Bz) = 0 as a reference.

To assess the performance of the Horseshoe model, we refer to the results of Event 1 obtained using the spheromak model (Scolini et al. 2019) and the FRi3D model (Maharana et al. 2022), as illustrated in Fig. 9. The quantitative analysis of the performance of the CME models in predicting the full B_z profile at Earth is discussed in Sect. 5. We introduce Fig. 9 here to set the context for the minimum B_z analysis. All the Horseshoe ensemble runs perform better than the spheromak model as shown in Fig. 8a in predicting the minimum B_z. Although the FRi3D model has an absolute median Δmin(B_z) closer to zero, the positive median Δmin(B_z) implies that the minimum B_z is overestimated as compared to observations in that case. The best simulation with the Horseshoe model, that is, the ϕ_p + Δϕ_p, avg case, has overlaps between its first quartile and the third quartile of the spheromak model and between its third quartile and the first quartile of the FRi3D model. This implies that the minimum B_z prediction capability of the Horseshoe model is intermediate between those of the spheromak and FRi3D models. The relevant information regarding the Δmin(B_z) from the boxplot analysis for Event 1 is provided in Table 2.

Fig. 9. Comparison of the predicted CME profiles modelled with spheromak, FRi3D, and Horseshoe models for Event 1. The description of plot panels is the same as in Fig. 6. The grey bar at the top of the Bz panel depicts the DTW analysis period in Sect. 5.

The last panel of Fig. 6 shows the results of Dst modelled for the ensemble runs compared to the measured Dst. The reference modelled min(Dst) = − 176 nT) is overestimated (Dst_obs = −139 nT). The min(Dst) modelled by the ϕ_p + Δϕ_p, avg ensemble run (−133 nT) is the best prediction and is underestimated by 24% with respect to the reference prediction. As a comparison, the FRi3D model predicts the min(B_z) = − 159 nT and the spheromak model predicts −75 nT, which are underestimated by 9% and 57%, respectively, with respect to the reference prediction (Fig. 9). It can be inferred that the Horseshoe model is an improvement over the spheromak model in predicting the minimum B_z and Dst and that, with further optimisation, it has the potential to match the FRi3D model.

4.2. Event 2: 10 September 2014

This event involves two successive CMEs erupting between 8 and 10 September 2014. On 8 September 2014, around 23:12 UT, the first CME (hereafter, CME1) was associated with the M4.6 flare originating from the NOAA AR 12158. The second CME (hereafter, CME2) erupted on 10 September 2014 at 17:21 UT from the same AR 12158, and was associated with the X1.6 flare. Both CMEs were detected by C2 and C3 instruments on board LASCO and COR2B instruments on board STEREO-B. Previous studies have addressed the details of CME eruption, propagation, and geo-effectiveness of the CME2 in detail (Cho et al. 2017; Webb & Nitta 2017; An et al. 2019; Maharana et al. 2023). The arrival of CME1 was not considered in these studies as it was just a flank hit at Earth and was not reported in ICME catalogues (e.g. Richardson & Cane 2010; Nieves-Chinchilla et al. 2018). However, CME2 was a candidate of the ISEST VarSITI campaign³ and was used to perform the exercise of real-time forecasting. Due to its head-on impact with a strong negative B_z, the CME2 was forecasted to cause severely disturbed geomagnetic conditions. However, it arrived at Earth with a positive B_z in its magnetic cloud and a negative B_z in the sheath ahead, causing only moderate storm (Dst ∼ −88 nT) conditions. Later, a thorough analysis of the remote observations of the source region (Vemareddy et al. 2016; Dudík et al. 2016; Zhao et al. 2016) and the in situ measurements of the interplanetary CME (Marubashi et al. 2017; Kilpua et al. 2021; An et al. 2019) suggested a possible northward deflection and significant rotation of the CME in the low corona, which were not captured in observations. A close-to-flank passage of the flux rope at Earth therefore led to erroneous space weather predictions. The comprehensive story of the erroneous space weather prediction associated with this event (CME1+CME2) is provided in Maharana et al. (2023). These authors analysed the orientation of CME2 at the source region close to 0.1 au and at 1 au and hypothesised that the rotation must have happened in the low corona (see Fig. 8 and associated text). The authors corroborated the claim with simulations and reproduced the negative B_z in the sheath created because of the interaction between CME1 and CME2.

The shock driven by CME1 (S1) was observed on 11 September at 22:50 UT based on the IPShock catalogue⁴ (Kilpua et al. 2015). The start of the magnetic ejecta associated with CME1 (ME1) was identified on 12 September at 8:45 UT. It continued until the shock associated with CME2 (SE2) arrived at L1 on 12 September at 15:17 UT according to the WIND ICME catalogue⁵. The magnetic ejecta associated with CME2 (ME2) passes through L1 between 12 September at 21:36 UT (corrected by Maharana et al. 2023) and 14 September at 11:38 UT, as reported in the WIND ICME catalogue.

In Maharana et al. (2023), this event was modelled using the combination of the spheromak and FRi3D models for CME1 and CME2, respectively. This is because combining two spheromaks resulted in erroneous B_z profiles at Earth, and using the FRi3D model for two successive CMEs was numerically challenging. For the consistent analysis of the performance of different flux-rope models in EUHFORIA, we performed a simulation where CME1 was modelled with spheromak and CME2 was modelled with the horseshoe model. We keep the parameters of CME1 the same as in Maharana et al. (2023) and focus on optimising the Horseshoe model parameters of CME2 through a parametric study by varying the ϕ_p in the error range. In addition, we provide results of the first attempt to model two successive CMEs with the Horseshoe model. The Horseshoe parameters for CME1 will be further optimised in future studies as there is significant scope for improvement.

4.2.1. Observationally constrained parameters for EUHFORIA simulation

Similar to Event 1, we performed the 3D reconstruction of CME2 with the FRi3D reconstruction tool to fit a toroidal geometry. The details can be found in Appendix B.2. Our fitting of CME2 for constraining the Horseshoe model is consistent with the parameters obtained by Maharana et al. (2023), who derived the parameters for constraining the FRi3D model. Overall, the CME is fitted to be thinner (low φ_hh), keeping in mind the toroidal geometry for the Horseshoe model, as compared to the reconstructed φ_hh for the spheromak and FRi3D models. We obtain R₀ = 10 R_⊙ and a = 6 R_⊙ applying the geometrical transformations given by Eqs. (14) and (15) on φ_hw = 52° and φ_hh = 12°. The v_rad is obtained from the reconstructed v_3D using Eq. (24). The ϕ_p, o values are taken from the Kazachenko (2023) catalogue. The unsigned φ_p for CME2 is 1.22 ⋅ 10²² Mx with an error of 1.21 ⋅ 10²¹ Mx (10%). Considering the left-handed chirality and the E–W polarity of the CME2 inferred by Maharana et al. (2023) and the geometrical orientation obtained from the reconstruction, we used a tilt of −35° for the Horseshoe model. The parameters of CME1 are the same as when modelled by spheromak in Maharana et al. (2023). We note that the CME1 was shifted in longitude towards the Sun–Earth line to model the interaction with CME2 in EUHFORIA. As the spheromak has no legs, the shift was necessary to simulate the flank encounter (Maharana et al. 2023). Further, we also constrained CME1 for the Horseshoe model. Applying the same methodology as that applied for CME2, we obtain R₀ = 7 R_⊙ and a = 7 R_⊙. CME1 has an unsigned poloidal flux of φ_p = 8.4 ⋅ 10²¹ Mx with an error of 1.58 ⋅ 10²¹ Mx, and a left-handed chirality. The parameters of the Horseshoe model used for CME2 are listed in Table 3 along with the spheromak and FRi3D parameters for comparison.

4.2.2. EUHFORIA simulations

In this section, we discuss the results of the ensemble modelling of CME2 based on changing the B₀ parameter and keeping the properties of CME1 constant, as shown in Fig. A.1. The arrival time and n_p of the ensembles qualitatively match the observations. CME arrives earlier for the higher poloidal flux, ϕ_p, values because of the higher Lorentz force, a trend similar to that seen for Event 1. All ensemble modelling results show an underestimation in the B_x and B_y profiles of ME2. The B_z inside ME2 is positive and increases with increasing ϕ_p. There is an overestimation of B_y in the sheath (β ≫ 1 region before the magnetic cloud). On the other hand, the coherent positive and negative features of the B_z component in the sheath region are well reproduced by the model. The negative value of the B_z component during this event is formed in the sheath of CME2 because of the interaction between CME1 and CME2 when the CME2 compresses the negative B_z part of ME1 in the trailing part of CME1. We note that the higher the difference in the speed of the CMEs, the more negative the minimum B_z. When the CME2 is slower (low ϕ_p), it creates a continuous compression of CME1 over a longer period before reaching Earth. Hence, if we want to model the most negative minimum B_z, then we predict a delayed arrival of CME2. There is therefore a trade-off between obtaining the best arrival time or the minimum B_z for this case based on only varying ϕ_p. In future work, a more comprehensive parameter study must be performed by varying the φ_hh and φ_hw to check the contribution of the extent and mass of the CME2 in the formation of the compressed features in the sheath ahead of it. For the current study, we chose, as the best case, the base run with a poloidal flux ϕ_p, 0 as it models the most accurate arrival time and also the negative B_z in the sheath reasonably well. The time series plot of the best case is shown in Fig. A.2 with the spatial variability of the plasma and magnetic field at σ = ±5° and σ = ±10° around Earth.

The boxplots of Δmin(B_z) (see Eq. (25)) for the Horseshoe ensemble runs, along with the simulations of CME2 performed with the spheromak and the FRi3D models (Maharana et al. 2023) are presented in Fig. 8b. A comparison of the performance of the three magnetised CME models in modelling Event 2 is provided in Sect. 5. The trends are opposite to those obtained for Event 1, as the lowest ϕ_p results in the lowest Δmin(B_z), as explained above. The range of Δmin(B_z) of the spheromak model is extended over [−65%, 25%], but this does not mean that the model captures the observed minimum negative B_z; it is rather the result of the wrongly predicted B_z component using the spheromak model, which is possibly due to spheromak rotation in the heliosphere (Asvestari et al. 2022). All the ensembles of the Horseshoe model must rather be compared with the FRi3D simulation as the B_z profile is more correctly modelled by the latter. For the best run, Δmin(B_z) lies in the range [−68%, −46%] at the position of all the virtual spacecraft and is −60% at Earth, meaning that min(B_z) is underestimated everywhere. The boxplots of the Horseshoe ensembles overlap with that of the spheromak and the FRi3D models and imply that the Horseshoe model performs better than both models in predicting the min(B_z) in the case of Event 2 (Fig. 8b). The reference modelled min(Dst) = − 47 nT is underestimated (min(Dst_obs) = − 88 nT). The min(Dst) modelled with the best Horseshoe ensemble run (ϕ_p, 0 case) is −17 nT, which is an underestimation by 63% as compared to the reference prediction. The Horseshoe result is close to that of the FRi3D prediction of −18 nT (underestimated by 61% compared to the reference model). The relevant information about the Δmin(B_z) analysis using boxplots for Event 2 is summarised in Table 2.

Horseshoe+Horseshoe simulation. With Event 2, we validated the use of the Horseshoe model in combination with other magnetised CME models. We go a step further, and model both CMEs with the Horseshoe model. The location of the CME source region (inferred from the 3D reconstruction) is kept the same as in the spheromak model. The geometrical and the B₀ values are constrained for the Horseshoe model, following the methodologies described in Sect. 3. The results are plotted in Fig. A.3 (orange dashed line). The CME1 shock arrival predicted by the Horseshoe+Horseshoe simulation is 2 hours before the observed arrival time. On the other hand, in the case of the spheromak+Horseshoe simulation, the arrival time was delayed by ∼2 hours. The Horseshoe+Horseshoe simulation predicts the Dst (see the last panel of Fig. A.3) better than the spheromak+Horseshoe case. These results highlight the possibility of using the Horseshoe model to model successive CMEs. This mitigates the limitation of the FRi3D model, namely that it cannot be used to numerically inject two successive CMEs, and improves the prediction of magnetic field components as compared to the spheromak model. The Horseshoe simulations will be further optimised in future studies, and more work will be carried out to model successive CMEs.

5. Analysis and discussion

In this section, we present an evaluation of the performance of the three magnetised CME models (spheromak, FRi3D, and Horseshoe) based on three criteria. First, we assess the ability of the models to predict the arrival time of the CME shock. Further, we evaluate the accuracy of the models based on their B_z prediction capabilities, and then also based on the speed of their numerical computations. Metrics provide a way to quantitatively assess which model provides the best results for improving space weather forecasting. We begin with identifying metrics for evaluation based on criterion 2. The metrics of mean absolute error (MAE) and root mean square error (RMSE) provide the averaged uncertainties for the whole time series, whereas advanced metrics such as dynamic time warping (DTW; Keogh & Pazzani 2001; Górecki & Łuczak 2013; Laperre et al. 2020; Samara et al. 2022) determine how similar two time series are by performing optimal temporal alignments between the common features (Muller 2007).

To accurately predict the geomagnetic indices, it is crucial to obtain not only the minimum negative B_z caused by a CME impact but also the commencement and the duration of the negative B_z values. Hence, we employ the DTW technique to assess the accuracy of the EUHFORIA-modelled B_z time profiles compared to the observations. The DTW technique is a distance measure similar to the Euclidean distance. It estimates the similarities between two time series containing similar patterns but differing in time. This technique has been previously applied to assess the EUHFORIA-modelled solar wind speed profiles with observations by Samara et al. (2022), and the Dst computed using EUHFORIA-modelled data by Maharana et al. (2024). We employ the open-source code developed by Samara et al. (2022)⁶ for our analysis.

We first applied DTW between the observed B_z data from WIND and the B_z modelled by employing different CME models in EUHFORIA. We ensure the following for the correct application of the algorithm: (a) the first and the last points of one sequence are matched with the first and last points of the other; (b) the mapping is monotonically increasing in time; and (c) there is no data gap, that is, every point in the two sequences is matched with at least one point in the other. To limit the ‘pathological alignment problem’, which creates singularities, that is, when one point of a sequence is matched with multiple points in the other, we apply a time ‘window’ for the alignment. Windowing restricts the mapping of the points to a certain time window and restricts the number of singularities. We also optimally smooth the observed one-minute cadence data at L1 (containing high-frequency fluctuations) in order to better match the trends in the smooth-modelled data. The fluctuations serve as local minima and maxima and influence the DTW results by creating more singularities and increasing the DTW cost. However, the limit of smoothing should be carefully determined in order to avoid missing the important features in the data. To apply the algorithm, first the DTW cost matrix is computed based on the following equation:

$$\begin{array}{c}\hfill D(i,j)=\delta (s_i,q_i)+\text{min}\{D(i-1,j-1),D(i-1,j),D(i,j-1)\},\end{array}$$(26)

where D(i, j) is the cumulative DTW cost or distance, and δ(s_i, q_i) = |s_i − q_i| corresponds to the Euclidean distance between the point s_i from one time series and the point q_i from the other time series. The first element of the array D(0, 0) is equal to δ(s₀, q₀). The last element of this cost matrix, the DTW score, is presumed to be a quantification of alignment between the two sequences. To evaluate the performance of different sequences with respect to a single time series, we calculate the sequence similarity factor (SSF) for each DTW analysis. The SSF quantifies how close the modelled result is to an ideal (observations) and a non-ideal (B_z = 0, reference) prediction scenario; it is the ratio between the DTW score of the observed and modelled B_z time series and the DTW score between the observed and reference scenario time series, and is defined as

$$\begin{array}{c}\hfill \mathrm{SSF}=\frac{{\mathrm{DTW}}_{\mathrm{score}}(O,M)}{{\mathrm{DTW}}_{\mathrm{score}}(O,N)},\mathrm{SSF}\in [0,\infty ),\end{array}$$(27)

where O, M, and N represent the observed, modelled, and non-ideal cases, respectively. We evaluated the CME models based on their SSF values for each validation event. In addition, we present, for each event, the plots of the DTW alignment maps, the histograms representing the distribution of the time differences (Δt), and amplitude differences (ΔB_z) between the observed and the modelled sequences (Δ = Observed − Modelled).

To assess the performance of the CME models based on criterion 3, we compared the computational time of the simulations. As the goal is to improve the models for operational space weather purposes, it is necessary to reduce the required computational time as much as possible. For consistency, all the simulations were performed using EUHFORIA (ver 2.0) on the wICE cluster of the Vlaams Supercomputer Centrum⁷ using two nodes with 72 cores per node (144 parallel processes).

5.1. Event 1

The in situ plasma characteristics, magnetic field, and the Dst predictions at Earth are presented in Fig. 9. Of the three considered models, the Horseshoe model predicts the most accurate shock arrival time with 34 minutes of delay compared to observed shock arrival. Unlike in the case of the FRi3D model, the n_p profile obtained with the Horseshoe model decreases to the ambient solar wind density values after the CME passage as in the observations. The DTW alignments between the observed B_z and the corresponding modelled profiles using the Horseshoe, FRi3D, and spheromak models for Event 1 are illustrated in Fig. 10 (panels a, d, and g), respectively. DTW is applied for the time period between 06:00 UT on 14 July 2012 and 12:00 UT on 18 July 2012, which covers the B_z profile of interest (the shaded grey region in the B_z panel of Fig. 9). We applied a window of 600 minutes for this event by visual inspection as most of the features could be matched within that time frame. Smoothing of the observed data is done over 500 minutes. Panels b, e, and h of Fig. 10 show the histograms of the Δt for Horseshoe, FRi3D, and spheromak models, respectively. The Δt is mostly negative for the FRi3D and the spheromak models, which means that they predict the features to occur later than they do in observations. The Horseshoe model has quite a flat distribution of the Δt, with most of the alignments in the positive spectrum. This implies that the features observed by the model occur earlier. The ΔB_z histograms for the Horseshoe, FRi3D, and spheromak models are shown in Fig. 10 (panels c, f, and i), respectively. Most of the ΔB_z alignments lie between ±2.5 nT for all three models. However, there are multiple alignments with ΔB_z in the range [ − 10, −5] nT in the case of the spheromak model, implying the underestimation of the modelled B_z.

Fig. 10. DTW analysis of Event 1 for all CME models. Rows 1 (a, b and c), 2 (d, e and f), and 3 (g, h and i) show the results for the Horseshoe model, the FRi3D model, and the spheromak model, respectively. Columns 1, 2, and 3 depict the DTW alignment between the observed (blue) and modelled (red) time series, histograms of the time differences between the aligned points, and histograms of the Bz differences between the aligned points, respectively.

The SSF (see Table 4) for the Horseshoe model is lower than in the case of the Spheromak model, meaning that the Horseshoe model better reconstructs the overall B_z profile. At the same time, the SSF for the Horseshoe model is only slightly higher than it is in FRi3D, with the temporal alignments quite homogeneously distributed for the Horseshoe model. The min(B_z) might not be as negative as predicted by FRi3D. However, the time alignment of different features in the B_z is better for the Horseshoe CME.

The Horseshoe model performs the best according to criterion 1 as it gives the best estimate of the arrival time. Considering criterion 2, the SSF of Horseshoe is quite close to the FRi3D model, making its predictions more reliable than the spheromak model. Using the same number of processing cores and resolution, the computational time for the simulations with the spheromak model is 20 minutes, for the Horseshoe model it is 3 hours 24 minutes, and for the FRi3D model it is more than three times longer, amounting to 9 hours and 2 minutes. Evaluating the models based on criterion 3, the Horseshoe can be placed as intermediate between the spheromak and the FRi3D models. Therefore, by combining evaluations based on both criteria, the Horseshoe CME model presents promising capabilities for reliable B_z predictions.

5.2. Event 2

The in situ plasma characteristics, magnetic field, and the Dst predictions at Earth are provided in Fig. A.3. To compare the performance of Horseshoe, FRi3D, and spheromak models in modelling CME2, we consider the simulations where CME1 is modelled with spheromak for all the cases. With the Horseshoe model, the predicted arrival time of S2 (shock associated with CME2) matches the observations. FRi3D predicts the arrival time ∼2 hours in advance with respect to the observed arrival times, and spheromak estimates the arrival time to be delayed by ∼2 hours.

DTW is applied for the time period between 18:00 UT on 11 September 2014 and 12:00 UT on 14 September 2014 (the shaded grey region in the B_z panel of Fig. A.3). This is a complicated case for applying DTW because of the observed temporally fast fluctuations in the sheath that are close to the features we want to capture. Smoothing of those features in observed time series, that is, to the extent required to be mapped with the modelled sequences, results in reduction of the magnitudes of the peaks of the positive and negative B_z in the sheath region. Choosing a smaller smoothing window increases the number of singularities in the alignment. Hence, we optimally smooth over 150 minutes to preserve the magnitude of the negative B_z in the sheath region. This allows us to compare the minimum negative B_z predicted by the CME models in EUHFORIA. Panels (b, e, and h) of Fig. A.4 show the distribution of Δt in the alignments of the Horseshoe, FRi3D, and spheromak models, respectively, for Event 2. The maximum number of alignments for both Horseshoe and Spheromak models lies within [300, 400] minutes (Δt), which implies that the modelled profiles predict the majority of features earlier than observed. For the FRi3D model, although the maximum alignments fall within [300, 400] minutes, the distribution of Δt has an increasing trend towards 400 minutes (i.e. the number of early predictions are more distributed over various Δt as compared to the other CME models). The histograms of ΔB_z for the Horseshoe, FRi3D, and the spheromak models are shown in Fig. A.4 (panels c, f, and i), respectively. Most of the alignments lie in the ΔB_z range [−2.5, 2.5] nT in the case of both the Horseshoe and the FRi3D models. For spheromak, the alignments are between [−5, 5] nT of ΔB_z, and also a large number of alignments are close to 10 nT. This is a consequence of the erroneous prediction of negative B_z in ME2 by the spheromak model as opposed to the observed positive B_z.

The Horseshoe model performs the best according to criterion 1, that is, it gives the best estimate of the arrival time. The SSF for the Horseshoe-modelled profile is better than in the case of both FRi3D and spheromak models and therefore the Horseshoe model also shows the best performance according to criterion 2. The Horseshoe simulation best reproduces the positive-to-negative switch in the sheath B_z. However, FRi3D matches the B_z strength in the magnetic cloud better, except for the initial overestimation of the positive B_z after the passage of the sheath. The computational times of the simulations with the spheromak, Horseshoe, and FRi3D models are 32 minutes, 8 hours 8 minutes, and 15 hours 22 minutes, respectively. For the FRi3D model, the majority of the total computational time is spent on the calculation of the magnetic field at the EUHFORIA inner boundary (i.e. mask computation). In the case of the horseshoe model, however, due to its complete analytical form, the time spent on the mask computation is reduced. The spheromak model still requires the least amount of time to compute the mask because of its purely force-free nature during injection into the heliospheric domain. The force-free nature is not completely satisfied when changing the geometry of the full torus into Horseshoe. This results in some unrealistic high-speed parcels near the inner boundary, which emerge from the cells where the legs of the Horseshoe CME are connected. Hence, pronounced gradients are formed in the speed profile in the computational domain close to the inner boundary. To maintain stability according to the Courant-Friedrichs-Lewy (CFL) condition, the allowed numerical time step, dt, becomes much lower for the Horseshoe model. This is why a greater number of iterations are necessary to reach the end of the simulations, which extends the overall computational duration. Such problems are reported in previous works, such as Regnault et al. (2023) and Linan et al. (2024). Therefore, despite the analytic formulation of the Horseshoe magnetic field, the computational time is still high because of how we disconnect the CME. As the high-speed artefacts occur closer to the boundary, they do not affect the predictions at Earth. However, we acknowledge that this is a crucial problem in the case of successive CME injection or when predicting the impact of space weather at locations closer to the inner boundary. Detailed future work will therefore be dedicated to mitigating this issue.

The Horseshoe model shows an intermediate performance between the other two models based on criterion 3. It meets both our criteria – speed of running the simulations and accuracy of prediction – reasonably as compared to the other CME models in the framework of EUHFORIA. This highlights the potential of the Horseshoe model in providing reliable and efficient operational space weather forecasting.

5.2.1. General remarks

Isavnin (2016) points out that the FRi3D model overestimates the magnetic flux budget of CMEs due to underestimation of the magnetic field line twist near the flux-rope boundary. The fact that FRi3D performs better in matching the total magnetic field is because of the overestimation of B_y and B_z components (B_x is not properly modelled for the events in this study). The Horseshoe model, with the current method of constraining magnetic field properties, provides reasonable magnitudes of B_y and B_z but still fails to model the B_x component. We note that B_x is not modelled accurately by any of the models we studied. The reasons for not capturing B_x could be the erroneous reconstruction of the CME geometry and the location of launch, which means that the impact of the CME is not through its flank. Other reasons could be because of the circular cross-section of the models, or that the CME was deflected during its propagation between 0.1 and 1.0 au, which our simulations do not capture self-consistently.

6. Summary and outlook

We summarise the key findings of this study as follows:

1. We implemented the Horseshoe CME model in EUHFORIA. Our Horseshoe model has a modified Miller-Turner magnetic configuration and a modified torus-like geometry that better mimics the CME leg structures. We point out the differences in the magnetic field profiles predicted by the full torus and the Horseshoe models at 1 au. Based on our results, we are able to recommend the use of the Horseshoe model in space weather forecasting for its better computational performance and modelling accuracy.

2. We designed methodologies to constrain the geometrical, plasma, and magnetic field parameters of the Horseshoe model from observations, and demonstrate their efficiency with two examples.

3. We validate the Horseshoe model with two events: first, a single non-interacting CME event of 12 July 2012 (Event 1), and second, a CME–CME interaction event that took place between 8 and 10 September 2014 (Event 2). We use ensemble modelling to address the variability in the predicted plasma and magnetic field time profiles due to errors in observational constraint of the axial magnetic field strength. In addition, we analyse the changes in the predicted min(B_z) compared to the observations at spatial locations around Earth. The CME shock arrival time to Earth and minimum negative B_z component predicted by the Horseshoe model are comparable to those predicted by the FRi3D model (a realistic flux-rope model) in both events. The computational time required by the Horseshoe model is intermediate between those required by the spheromak and the FRi3D models.

4. We make the first ever attempt to test the capability of the Horseshoe model to inject two successive CMEs into EUHFORIA using Event 2. This mitigates the limitation of the FRi3D model and provides reasonable results, meaning a significant advancement in comparison to the spheromak model. Further validation of the Horseshoe model in modelling multiple CMEs must be carried out.

5. The B_z component of the magnetic field obtained with the Horseshoe model is compared with those obtained with the FRi3D and spheromak models in the framework of EUHFORIA using the dynamic time warping (DTW) technique.

Overall, the Horseshoe model shows the great potential for predicting CME arrival time and magnetic field properties in an operational setup. We note that this model is numerically more stable when launching initially force-free magnetic field configurations (spheromak) than the non-force-free models. The full toroidal mMT topology is divergence-free and approximately force-free initially. Although not the most accurate in reproducing the magnetic field profiles in the heliosphere, it is numerically more stable and computationally less expensive. On the other hand, incomplete injection of the force-free magnetic field through the inner boundary gives rise to residual j × B forces that create excess ∇ ⋅ B. As the Constrained Transport methodology in EUHFORIA is not efficient in cleaning the ∇ ⋅ B, it is difficult to sustain the stability of CME models with legs, which continuously inject additional magnetic field. One of the next steps is to implement the Horseshoe model in Icarus (Verbeke et al. 2022; Baratashvili et al. 2022), which is an advanced version of EUHFORIA implemented in the framework of MPI-AMRVAC (Keppens et al. 2003; Xia et al. 2018) and uses a parabolic divergence cleaning method, which is more efficient than the Constrained Transport in EUHFORIA. Additional future work will be to contract the back part of the torus and inject it fully through the inner boundary in order to conserve the force-free fields. The commonly used constant-alpha (force-free) magnetic field models used to fit interplanetary CMEs exhibit an increasing twist outward towards the flux-rope boundary (e.g. the mMT topology as mentioned in Sect. 2.1). However, most of the observed interplanetary CMEs have a uniform twist distribution (Hu et al. 2015). Hence, Vandas & Romashets (2019) prescribe a flat twist profile to a cylindrical flux rope, which renders the flux rope non-force-free. On large scales, the assumption of a force-free magnetic field is generally invalid (Kilpua et al. 2011). Hence, non-force-free magnetic field topology with a uniform twist profile can be used for interplanetary CME models, although their numerical implementation might pose challenges. In future works, such a uniform twist distribution must be developed for toroidal flux ropes in order to obtain a more accurate modelling of CME evolution.

Acknowledgments

The authors acknowledge support from the projects C14/19/089 (C1 project Internal Funds KU Leuven), G0B5823N and G002523N (WEAVE) (FWO-Vlaanderen), 4000134474 (SIDC Data Exploitation, ESA Prodex), Belspo project B2/191/P1/SWiM, the FEDtWin project PERIHELION and the financial support received to attend the third BINA workshop under the Belgo-Indian Network for Astronomy and astrophysics (BINA) project approved by the International Division, Department of Science and Technology (DST, Govt. of India; DST/INT/BELG/P-09/2017) and the Belgian Federal Science Policy Office (BELSPO, Govt. of Belgium; BL/33/IN12). We thank Dr. Christine Verbeke for insightful discussions about CME model implementation in EUHFORIA, and Dr. Brigitte Schmieder for their constant support and suggestions for the improvement of this model.

References

Appendix A: Additional figures

All the figures associated with Event 2, discussed in Sections 4.2 and 5 are depicted included in this section as Figs. A.1,A.2, A.3 and A.4.

Fig. A.1. Results of Horseshoe ensemble simulations for Event 2 with varying ϕp obtained using EUHFORIA. The detailed plot description is the same as for Fig. 6.

Fig. A.2. Results of the best Horseshoe ensemble simulation of Event 1 obtained using EUHFORIA. The detailed plot description is the same as for Fig. 7. The grey bar at the top of the Bz panel depicts the timespan of the dynamic time warping analysis in Sect. 5.

Fig. A.3. Comparison of the predicted CME profiles modelled with spheromak, FRi3D, and Horseshoe models for Event 2 (CME1 modelled with spheromak). An additional simulation with both CME1 and CME2 modelled with the Horseshoe model is provided (dashed orange line). The description of the plot panels is the same as for Fig. 6. The grey bar at the top of the Bz panel depicts the time span of the dynamic time warping analysis in Sect. 5.

Fig. A.4. DTW analysis of Event 2 for all CME models. Rows 1, 2, and 3 show the results for the Horseshoe model, the FRi3D model, and the spheromak model, respectively. Columns 1, 2, and 3 depict the DTW alignment between the observed (blue) and modelled (red) time series, histograms of time differences between the aligned points, and the histograms of the Bz differences between the aligned points, respectively.

Appendix B: Geometrical reconstruction of the events

In this section, we provide the details of the geometrical reconstruction of the CMEs from the white light images. As this procedure can be subjective and, hence, introduces uncertainties in CME arrival time predictions (Verbeke et al. 2023), we provide our rationale behind the features we fit. We use the white-light image tool of Flux Rope in 3D model (FRi3D, Isavnin 2016).

B.1. Event 1: July 12, 2012

Although this event has been reconstructed before in previous studies (Scolini et al. 2019; Maharana et al. 2022), we repeat the procedure to match the constraints of a toroidal geometry to validate the Horseshoe model. The spacecraft STEREO-A (STA) and STEREO-B (STB) observed the CME as a limb event in their field of view (FOV). They were separated by 125° in longitude and 10° in latitude (Fig. B.1). COR2A and CORB are the coronagraphs onboard STA and STB, respectively. LASCO spacecraft observed the CME as a halo event. However, the coronagraph data was unavailable during the event. We temporally fitted the flux-rope mesh between 17:24 and 18:24 on 12 July 2012, when the full CME was visible clearly in the coronagraph FOV. The 3D reconstruction of the CME was done considering the following aspects noticed in the observations:

• The streamer is observed to be bent from the COR2A perspective, and there is movement through the streamer. However, there are no such signatures from the COR2B perspective, which means that the CME interacted with the streamer on its western limb during the evolution. Hence, the southern part of the CME is fitted in such a way that it looks like it fits part of the streamer, but it is the flux rope passing through it.

• From COR2A FOV, we fit the angular half-width of the CME in the northern limb to not fit the northern streamer completely, which is on the way of the CME and is then deflected.

• For fitting the leading edge, we fit the COR2A observation of the magnetic cloud (darker region) as it is brighter than the COR2B counterpart. Physically, the CME front is closer to the COR2A’s perspective, so it is better to rely on it. The leading edge looks underestimated from the COR2B perspective.

• The stand-off distance between the shock and the CME is insignificant, making it difficult to distinguish the CME parts for reconstruction. These observations point to a very twisted flux rope with local deformations that are not straightforward to capture with the smooth global mesh structure of the existing reconstruction models.

• Flux-rope models like FRi3D or the Torus model represent a more realistic geometry of the CME represented by the reconstruction shown in Fig. B.2. Whereas to constrain parameters for a spheromak model, the CME volume is overestimated and does not represent an extended flux-rope geometry (see, e.g. Fig. 4 of Scolini et al. (2019).)

B.2. Event 2: September 10, 2014

The CME2 was tracked between 18:54 and 19:54 on 10 September 2014, when all the features of the CME were visible in the white light images recorded by LASCO C3 and STB COR2B. The 3D reconstruction of this event is based on the following distinct features:

• The CME interacts with two streamers during its evolution in the coronagraph FOV. From the perspective of LASCO C3, the CME pushes through the streamer visible on the north-western limb (Streamer 1) and the other close to the north pole (Streamer 2). The CME shock and the magnetic cloud interact with the streamers as distinct brightening of the streamers, which in turn fades out when the CME has completely passed. We fit the parts of the Streamer 2 on the northern flank, and part of the Streamer 1 lies within the southern extent of the CME on its western limb.

• It is the front part of the CME which interacts with Streamer 1. Although the shock is observed to disturb the streamer from the COR 2A perspective, the bulk of the CME moves away from the streamer root. As it is a back-sided halo event from COR 2A, we fit the features as per C3 features. Hence, it might seem that the features are ill-fitted in COR 2A FOV.

• Streamer 2 is brightened throughout the complete propagation of the CME, which means that the streamer is present on the front side of the Sun.

• Both the shock front and the twisted magnetic cloud structure have a wavy morphology, which we are trying to fit with a smooth croissant-like structure. The brighter structure seen close to the CME nose in COR 2 FOV could be the twisted extension at the CME back, and hence, it is less bright from the C3 perspective. However, it isn’t very easy to do justice to fit all the features of the CME with simplified reconstruction models. Hence, we fit the features that are best visible from the Earth-directed FOV.

We also performed the reconstruction of CME1, which erupted more than a day before CME2, with the FRi3D model and the results are presented in Fig. B.5.

Fig. B.1. Position of the STEREO-A and STEREO-B spacecraft during the CME eruption on 12 July 2012. The image was created using the freely available ‘Where is STEREO’ tool (https://stereo-ssc.nascom.nasa.gov/cgi-bin/make_where_gif).

Fig. B.2. White light images as observed from COR2-A and COR2-B and the same images overlaid with FRi3D reconstructed wire frame in the bottom panels. The arrows show different projected white light features used as proxies to reconstruct the 3D structure as discussed in the text.

Fig. B.3. Position of the STEREO-A and STEREO-B spacecraft during the CME eruption on September 10, 2014. The image was created using the ‘Where is STEREO’ tool as in Fig. B.1.

Fig. B.4. White light images observed from COR2-B and C3 and the same images overlaid with FRi3D reconstructed wire frame in the bottom panels. The arrows show different projected white light features used as proxies to reconstruct the 3D structure, as discussed in the text.

Fig. B.5. White light images observed from COR2-B and C3 and the same images overlaid with FRi3D reconstructed wire frame in the bottom panels. The arrows show different projected white light features used as proxies to reconstruct the 3D structure, as discussed in the text.

All Tables

All Figures

Fig. 1. Horseshoe model implementation. Left: schematic representation (not-to-scale) of the Horseshoe model propagating into the EUHFORIA heliosphere domain in the equatorial plane (x − y); Right: comparison of the EUHFORIA predictions at 1 au between the Horseshoe and the full torus geometry with the mMT magnetic field configuration. The red arrow in the By profile shows the double polarity signatures of the full torus. The Bz profile of the full torus is compressed and merged with the trailing part of the torus.

In the text

Fig. 2. Equatorial snapshot of EUHFORIA simulations of the full torus and the Horseshoe implementation of the mMT magnetic field configuration. The evolution of the internal magnetic field lines of the CMEs (colour-coded with the magnetic field strength), 4 hours (panels 1–2), and 11 hours (panels 3–4), respectively, after the start of the injection at 0.1 au, is illustrated. In the background, the divergence of the speed is plotted. The sheath ahead of the magnetic cloud is depicted by the negative divergence (accumulation of plasma) region, showing a clear envelope around the CME. In panels 1 and 3, the rear part of the full torus is seen injected, whereas panels 2 and 4 show the horseshoe-like geometry creating CME leg-like structures connected to the inner boundary.

In the text

Fig. 3. Torus geometry in the equatorial and meridional planes for Case 1 (a, b) and Case 2 (c, d), where geometries are represented with red curves. The horizontal black solid line is the Sun–Earth line. The blue lines on the equatorial and meridional plots have angular extents of φhw and φhh, respectively, about the Sun–Earth line. In this example, φhw = 40° and φhh = 12°. In Case 1, R0 = 7.6 and a = 6.2, making the aspect ratio 1.23. In Case 2, R0 = 11.1 and a = 6.9, making the aspect ratio 1.6.

In the text

	Fig. 4. Contours of major radius (R0) as a function of φhw and φhh for (a) Case 1 and (b) Case 2 geometries. Contours of the minor radius (a) as a function of φhw and φhw for (c) Case 1 and (d). The parameter space where R0 > a is identified for (e) Case 1 and (e) Case 2, geometries, respectively.
In the text

	Fig. 5. Distribution of aspect ratio (a) Case 1 and (b) Case 2 geometries and the volume (c) Case 1 and (d) Case 2 geometries. The colour bars show the magnitude of the aspect ratio and the volume in the respective plots. The black dashed and the solid curves demarcate the allowed range of φhw − φhh space for Case 1 and Case 2 geometries, respectively.
In the text

Fig. 6. Results of Horseshoe ensemble simulations for Event 1 with varying ϕp obtained using EUHFORIA. Panels from top to bottom: speed (v), proton number density (np), Bx, By, Bz, magnetic field strength (|B|), plasma beta (β), and Dst. The observations from the WIND spacecraft are plotted in black in all the plasma and magnetic field properties panels. In the Dst panel, the observed Dst from WDC, Kyoto, is plotted in red, the reference modelled Dst using the observed data in black, and the Dst using EUHFORIA simulations in other colours. The solid black vertical line marks the shock arrival time, and the dashed magenta lines show the start and end of the magnetic cloud as reported in the WIND ICME catalogue (https://wind.nasa.gov/ICME_catalog/ICME_catalog_viewer.php).

In the text

	Fig. 7. Results of the best Horseshoe ensemble simulation of Event 1 obtained using EUHFORIA at Earth (solid blue line) and at virtual satellites in the 5 − 10° latitudinal and longitudinal offset around Earth (shaded region). All other plot details are similar to Fig. 6. The grey bar at the top of the Bz panel depicts the time span of the dynamic time warping analysis in Sect. 5.
In the text

	Fig. 8. Boxplot of the spread of the relative minimum Bz (Δmin(Bz)) values at Earth and the virtual spacecraft at an offset of latitude and longitude of 5–10° – (a) Event 1, (b) Event 2. The red dot and the orange line in the boxes represent the Δmin(Bz) at the location of Earth and the median Δmin(Bz) of the distribution, respectively. The dashed blue line indicates Δmin(Bz) = 0 as a reference.
In the text

	Fig. 9. Comparison of the predicted CME profiles modelled with spheromak, FRi3D, and Horseshoe models for Event 1. The description of plot panels is the same as in Fig. 6. The grey bar at the top of the Bz panel depicts the DTW analysis period in Sect. 5.
In the text

Fig. 10. DTW analysis of Event 1 for all CME models. Rows 1 (a, b and c), 2 (d, e and f), and 3 (g, h and i) show the results for the Horseshoe model, the FRi3D model, and the spheromak model, respectively. Columns 1, 2, and 3 depict the DTW alignment between the observed (blue) and modelled (red) time series, histograms of the time differences between the aligned points, and histograms of the Bz differences between the aligned points, respectively.

In the text

	Fig. A.1. Results of Horseshoe ensemble simulations for Event 2 with varying ϕp obtained using EUHFORIA. The detailed plot description is the same as for Fig. 6.
In the text

	Fig. A.2. Results of the best Horseshoe ensemble simulation of Event 1 obtained using EUHFORIA. The detailed plot description is the same as for Fig. 7. The grey bar at the top of the Bz panel depicts the timespan of the dynamic time warping analysis in Sect. 5.
In the text

Fig. A.3. Comparison of the predicted CME profiles modelled with spheromak, FRi3D, and Horseshoe models for Event 2 (CME1 modelled with spheromak). An additional simulation with both CME1 and CME2 modelled with the Horseshoe model is provided (dashed orange line). The description of the plot panels is the same as for Fig. 6. The grey bar at the top of the Bz panel depicts the time span of the dynamic time warping analysis in Sect. 5.

In the text

Fig. A.4. DTW analysis of Event 2 for all CME models. Rows 1, 2, and 3 show the results for the Horseshoe model, the FRi3D model, and the spheromak model, respectively. Columns 1, 2, and 3 depict the DTW alignment between the observed (blue) and modelled (red) time series, histograms of time differences between the aligned points, and the histograms of the Bz differences between the aligned points, respectively.

In the text

	Fig. B.1. Position of the STEREO-A and STEREO-B spacecraft during the CME eruption on 12 July 2012. The image was created using the freely available ‘Where is STEREO’ tool (https://stereo-ssc.nascom.nasa.gov/cgi-bin/make_where_gif).
In the text

	Fig. B.2. White light images as observed from COR2-A and COR2-B and the same images overlaid with FRi3D reconstructed wire frame in the bottom panels. The arrows show different projected white light features used as proxies to reconstruct the 3D structure as discussed in the text.
In the text

	Fig. B.3. Position of the STEREO-A and STEREO-B spacecraft during the CME eruption on September 10, 2014. The image was created using the ‘Where is STEREO’ tool as in Fig. B.1.
In the text

	Fig. B.4. White light images observed from COR2-B and C3 and the same images overlaid with FRi3D reconstructed wire frame in the bottom panels. The arrows show different projected white light features used as proxies to reconstruct the 3D structure, as discussed in the text.
In the text

	Fig. B.5. White light images observed from COR2-B and C3 and the same images overlaid with FRi3D reconstructed wire frame in the bottom panels. The arrows show different projected white light features used as proxies to reconstruct the 3D structure, as discussed in the text.
In the text