CME–CME Interactions as Sources of CME Geoeffectiveness: The Formation of the Complex Ejecta and Intense Geomagnetic Storm in 2017 Early September

The first CME (hereafter CME1) was first observed in the LASCO/C2 coronagraph on September 4 at 19:00 UT as a partial halo with a dominant propagation component toward the southwest. It was associated with an M1.7 class flare (start: 18:46 UT—peak: 19:37 UT—end: 19:52 UT) localized in AR 12673. The CME propagated in LASCO/C2 and LASCO/C3 with an average speed (projected in the plane of sky) of about 600 km s⁻¹, exhibiting a slightly accelerating behavior (from the LASCO CME catalog, https://cdaw.gsfc.nasa.gov/CME_list/). On the day of the first eruption, STEREO-A was separated from Earth by an angle of 128°; therefore Earth-directed CMEs could be well observed by the COR2 instrument on board STEREO-A. In this instrument, the CME was visible starting from 19:39 UT, where it appeared to propagate toward the southwest.

When the leading edge of CME1 was at ∼10 solar radii (R_⊙) as seen by LASCO/C3, it was overtaken by a second, faster CME (hereafter CME2) that was first observed in LASCO/C2 at 20:36 UT. This second CME appeared from Earth as a full halo having an intensified frontal part propagating toward the south, with an average projected speed of about 1420 km s⁻¹. CME2 was associated with an M5.5 class flare (start: 20:28—peak: 20:33—end: 20:37) also localized in AR 12673. In STEREO/COR2-A, the CME was seen to propagate toward the southwest (first appearance: 20:39 UT), catching up with CME1 shortly after 21:00 UT. By 21:42 UT, the leading edges of CME1 and CME2 had completely merged as seen by LASCO/C3 as well, so that the two structures became indistinguishable in both LASCO and COR2 images. The two CMEs erupted from the same AR less than 2 hr apart and exhibited similar coronal signatures, suggesting a sympathetic, (quasi-)homologous nature (Zhang & Wang 2002; Cheng et al. 2006; Wang et al. 2013).

The day after the eruption of CME1 and CME2 (i.e., September 5), little activity was observed in the LASCO and STEREO-A coronagraphs. One faint partial halo CME, most probably also erupting from AR 12673, was visible starting from 17:36 UT in LASCO/C2, but it became too faint to be tracked in the LASCO/C3 field of view. It was at all times barely visible even in running-difference images from STEREO/COR2-A, where its front appeared to propagate predominantly below the ecliptic plane. The faintness of this CME in LASCO and STEREO-A coronagraph images, combined with the limited eruption signatures visible in EUV images of the solar disk, makes the reconstruction of its kinematic, geometric, and magnetic parameters using the techniques presented in Section 3 particularly complicated. Considering also its propagation direction below the ecliptic plane, we have neglected this event in the following analysis (see also Werner et al. 2019, for a similar modeling approach). However, we point out that this CME could have contributed to the complexity of the event observed at Earth by interacting with the preceding and following CMEs.

Finally, on September 6, a full-halo CME, hereafter CME3, was observed entering the LASCO/C2 field of view at 12:24 UT. This CME originated from the same AR as the previous ones, and it was associated with a remarkably intense flare of class X9.3 (start: 11:53 UT—peak: 12:02 UT—end: 12:10 UT). The CME was observed to propagate toward the southwest with a projected speed of about 1570 km s⁻¹, and its leading edge was characterized by a highly elliptical shape tilted by about 45° with respect to the solar equator. In STEREO/COR2-A the CME appeared as a full halo (first appearance: 12:24 UT) characterized by a southwest propagation direction (see Appendix A for additional details).

The source regions of the CMEs discussed above can all be located within AR 12673. This AR presented outstanding levels of CME and flare productivity persisting for more than a full week (Chertok et al. 2018; Redmon et al. 2018; Bruno et al. 2019). The region was first classified as a simple α region (Hale et al. 1919; Künzel 1965) on August 30, when it was rotating toward the solar disk center from the eastern limb. It was then classified as βγ on September 3, that is, the day before the eruption of CME1. The region then developed into a βγδ configuration starting from September 5. Figure 1 shows SDO observations of the AR as observed by the HMI and AIA instruments.

Zoom In Zoom Out Reset image size

From September 4 onward, photospheric magnetograms of the AR show the presence of a complex system of polarity inversion lines (PILs) that evolved and rotated over the days (Figures 1(a)–(c)). Two main PIL systems, one in the southeast part of the AR, characterized by an approximately north–south orientation, and one in the northeast part, exhibiting an approximately east–west direction, are visible.

We use SDO/AIA 171 and 1600 Å images to pinpoint the location of the eruption of the three CMEs within the AR. CME1 erupted in the southeast part of the AR, as indicated by the development of flare ribbons (visible in 1600 Å, Figure 1(d)) and by the southward expansion of coronal loops during the eruption (visible in 171 Å, not shown here). A posteruptive arcade (PEA) (visible in Figure 1(j)) also formed after the first eruption until the onset of the second eruption. Starting at 20:28 UT, associated with the eruption of CME2, more extended flare ribbons developed in the northern part of the AR (Figure 1(e)). These observations suggest that the eruption of CME1 remained confined to the southeast part of the AR, while the eruption of CME2 developed through the whole PIL system up to its northwest end (Figure 1(b)). The formation of a stable PEA is visible in the AIA 193 Å filter (Figure 1(k)), confirming that the magnetic reconnection processes associated with the eruption extended over the whole PIL system elongation. The short waiting time between CME1 and CME2 (less than 3 hr) and their origin from the same AR strongly favor the scenario of quasi-homologous CMEs, where the second eruption is commonly interpreted as a consequence of the flux rope destabilization caused by the rearrangement of coronal magnetic fields following the first eruption (Török et al. 2011; Bemporad et al. 2012; Chatterjee & Fan 2013; Wang et al. 2013; Liu et al. 2017). The eruption of CME3, occurring about 41 hr after CME2, originated from the vertical PIL located around the center of the AR, where changes in the surface magnetic field in 45 s HMI observations are visible starting from 11:54 UT (see also Mitra et al. 2018). Bright flare ribbons visible in the AIA 1600 Å line (Figure 1(f)) and a PEA visible in the AIA 193 Å filter (Figure 1(l)) indicate that magnetic reconnection extended over the whole PIL elongation (Figure 1(c)).

Figure 2 shows 1 minute averaged in situ measurements taken by the Wind and DSCOVR spacecraft during the days following the eruptions, together with the 1 hr Dst index measured on ground (provided by the World Data Center for Geomagnetism, Kyoto; http://wdc.kugi.kyoto-u.ac.jp/dstdir/).

Zoom In Zoom Out Reset image size

Since CME1 and CME2 are observed to merge into a single structure (hereafter CME1+CME2) in the fields of view of LASCO and SECCHI, that is, already in the corona, it is reasonable to expect them to drive a single, common shock (e.g., Odstrcil et al. 2003; Xiong et al. 2007; Lugaz et al. 2013) as they propagate in interplanetary space, most probably the one observed at Wind at 23:13 UT on September 6 (hereafter S1). S1 was followed by a prolonged sheath region (with a duration of ∼21 hr, corresponding to a thickness of ∼0.25 au for a structure moving at ∼500 km s⁻¹) characterized by a fluctuating magnetic field and relatively high density and temperature. Such a structure can imply a spacecraft crossing through a thick sheath region formed by the merging of the CME1 and CME2 sheaths, whose formation is compatible with the early merging of the two CMEs. On September 7 around 20:00 UT a region of low plasma β and smooth magnetic field, compatible with a magnetic ejecta (hereafter E1), was observed. E1, most probably associated with the merged CME1+CME2 ejecta, was also listed in the Richardson & Cane list of interplanetary CMEs (http://www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm; Cane & Richardson 2003; Richardson & Cane 2010) with start time around 20:00 UT and end time around 04:00 UT on September 8, that is, with a duration of 8 hr. In the list, the ejecta was classified with a quality flag equal to 1, indicating that it exhibited only some of the typical characteristics of magnetic clouds, for example, a low plasma β and coherent magnetic field rotation, but lacked some other characteristics such as an enhanced magnetic field magnitude.

Moreover, the ejecta was characterized by the presence of an interplanetary shock (hereafter S2) propagating through it. The shock was observed at 22:38 UT on September 7, that is, 2.5 hr after the start of E1, and it was most likely driven by CME3. S2 compressed the magnetic field of the E1 ejecta, resulting in a significant amplification of the southward B_z from a preexisting value of ∼−10, to ∼−30 nT. The negative z-component of the magnetic field in the shock upstream region triggered the beginning of a geomagnetic disturbance, marked by a decrease in the Dst index to about −50 nT. The further enhancement of the negative B_z in the downstream region led to the development of the first and strongest dip in the Dst profile, reaching −142 nT around 1 UT on September 8. Overall S2 presented several characteristics typical of shocks propagating inside preceding ejecta, including a low β in the upstream and downstream shock regions, and a magnetic field clock angle almost constant across the shock (Lugaz et al. 2015, 2017).

A second period of enhanced magnetic field and low β was observed between 11:00 UT on September 8 and 20:00 UT on September 10. This period was classified in the Richardson & Cane ICME list as a single ejecta with a quality flag equal to 1. However, we note the presence of a region of fluctuating fields and relatively high proton temperature and plasma β between 14:30 UT and 20:30 UT on September 8, which suggests the passage of two separate ejecta regions (hereafter E2 and E3). Both E2 and E3 exhibited low plasma β and enhanced magnetic fields with different levels of rotation. The period marked as E2 is characterized by a rotating magnetic field, a typical characteristic of spacecraft crossings through the flux rope structure. This ejecta period is most probably associated with the passage of CME3 that also drove S2. E3 exhibits typical characteristics of leg encounters, as featured by the large B_x component (not shown), the lack of magnetic field rotations, and a long duration (Marubashi & Lepping 2007; Möstl et al. 2010; Kilpua et al. 2011, 2013; Owens 2016). In view of the coronal and in situ CME/ICME observations, we consider it most probable that E2 and E3 were associated with the same CME (i.e., CME3) at the Sun and simply corresponded to crossings of the spacecraft through different parts of the flux rope. In this picture, E2 would correspond to a passage closer to the apex of CME3 and E3 to the passage through its leg. Additional evidence is provided by the speed profile, which decreases very coherently through and between E2 and E3, as well as by the coherent rotation of the magnetic field vectors. This interpretation would require a bending or deformation of the flux rope global structure as a consequence of its interaction with the ambient solar wind or preceding ejecta (Crooker et al. 1998; Mulligan et al. 1999; Dasso et al. 2007; Marubashi & Lepping 2007). The passage of E2 generated a second dip in the Dst index (minimum of −124 nT on September 8 around 13 UT), that is, the complex ejecta investigated here resulted in a two-step geomagnetic storm (Farrugia et al. 2006; Liu et al. 2014b). The magnetic field in E3 was low (about 5 nT) and pointed primarily to the north. As a consequence, the Dst recovered during E3 without further intensification.

We note that in such high-activity periods, the identification of ICME signatures and their linking with the corresponding CMEs at the Sun become very complex due to interactions among the various structures in the corona and interplanetary space and the elevated number of potential CME candidates. For this reason, we cannot rule out a priori that E2 and E3 were associated with two different CMEs at the Sun, possibly involving the faint partial halo CME discussed in Section 2.1. However, in order to keep the assumptions as simple as possible at the beginning, in our simulations we consider both ejecta to be associated with CME3. This is a reasonable assumption as our primary focus concerns the investigation of the nature and evolution of the main geoeffective structures, that is, the southward field in E1 due to the compression by S2 and the southward field in E2, rather than the origin of E3. We also note that with this interpretation our in situ analysis is consistent with the previous studies by Shen et al. (2018) and Werner et al. (2019).

Between the end of E1 and the start of E2, a region characterized by plasma β ∼ 1, modest and fluctuating magnetic fields, and increasing density suggests the occurrence of magnetic reconnection at the interaction surface between E1 and E2 (Maričić et al. 2014). These in situ signatures exhibit several characteristics of ongoing CME–CME interactions, consistent with a picture in which the interaction of CME1+CME2 and CME3 was still at an early stage at 1 au (Lugaz et al. 2015, 2017).

During the days of the CME eruptions, STEREO-A was located at a longitude of −128° in Stonyhurst coordinates (Thompson 2006), providing a second vantage point to investigate the coronal evolution of the CMEs under study, in addition to the observations made along the Sun–Earth line. To constrain the kinematics and geometry of the CMEs in the corona, we perform a 3D fitting of the events from these two viewpoints (SOHO and STEREO-A) using the graduated cylindrical shell (GCS; Thernisien et al. 2006, 2009) model. The results from the GCS fitting are then used as input for the CME modeling using EUHFORIA (Section 3.3).

The GCS fitting provides as output the following parameters (in Stonyhurst coordinates): CME latitude (θ_CME), longitude (ϕ_CME), front height (h_CME), aspect ratio (κ_CME), half angle (α_CME), and tilt (γ_CME). For ${\alpha }_{\mathrm{CME}}\ne 0^\circ$, the shell of the GCS model corresponds to a croissant-like shape. From the aspect ratio and half angle we derive the edge-on (EO), face-on (FO), and average (AV) CME half-widths (${\omega }_{\mathrm{CME}}/2$) (Thernisien et al. 2009). From a modeling perspective, the determination of a single value for the CME half-width from raw GCS outputs is critical because in EUHFORIA CMEs are initialized with spherical shapes, that is, their cross sections are symmetric in all directions. From the fitting of sequential images, we also derive the deprojected (3D) speed of the CME apex (v_CME) as well as the radial/translational CME speed (${v}_{\mathrm{CME}}^{\mathrm{rad}}$, corresponding to the speed of the center of the croissant tube) and the CME expansion speed (${v}_{\mathrm{CME}}^{\exp }$, corresponding to the rate of increase of the croissant cross-sectional radius) (Scolini et al. 2019, see also Appendix A for the analytical derivation). Each speed was determined by performing a linear fitting of instantaneous results from LASCO/C2–C3 and STEREO/COR2-A images (CME1: between 19:54 UT and 20:39 UT on September 4, for a total of three images; CME2: between 20:39 UT and 21:54 UT on September 4, for a total of five images; CME3: between 12:39 UT and 13:54 UT on September 6, for a total of five images).

Table 1 lists the results of the GCS fitting for the three events under study. Snapshots of the GCS fitting results overplotted on LASCO/C3 and STEREO/COR2-A images are provided in Appendix A. We note that from the extrapolation of the times at which the CME leading edges reached 0.1 au, we ascertain that CME2 reached this distance about 15 minutes earlier than CME1. This indicates that the interaction between CME1 and CME2 occurred at heliocentric distances close to or slightly lower than 0.1 au.

Table 1.  Results from the GCS Fitting of the Three CMEs under Study

Parameter	CME1	CME2	CME3
θCME	0°	−25°	−11°
ϕCME	25°	0°	21°
κCME	0.38	0.50	0.43
αCME	10°	30°	15°
ωCME/2 (EO/AV/FO)	22°/27°/32°	30°/45°/60°	25°/33°/40°
γCME	0°	0°	40°
vCME	960 km s−1	1585 km s−1	1910 km s−1
${v}_{\mathrm{CME}}^{\mathrm{rad}}$	697 km s−1	1057 km s−1	1293 km s−1
${v}_{\mathrm{CME}}^{\exp }$	263 km s−1	528 km s−1	617 km s−1
Time at 0.1 au	2017 Sep 4 23:00 UT	2017 Sep 4 22:44 UT	2017 Sep 6 14:11 UT

Note. EO = edge-on, FO = face-on, AV = average.

Download table as:  ASCIITypeset image

In preparation for the heliospheric CME simulations with the spheromak flux rope model in EUHFORIA (discussed in Section 3.3), we discuss in the following subsections the observational derivation of three key parameters characterizing the magnetic (flux rope) structure of the CMEs under study: their chirality, tilt, and the amount of magnetic flux reconnected during the eruption.

3.2.1. Chirality and Tilt of the Flux Ropes

Chirality of the Flux Ropes—Observationally, the magnetic helicity sign (or chirality) of ARs can be inferred from different morphological features (e.g., Démoulin & Pariat 2009; Palmerio et al. 2017). In the particular case under study, images in the extreme-ultraviolet SDO/AIA filter at 94 Å (Figures 1(g)–(i)) suggest that AR 12673 was characterized by a negative chirality as indicated by the presence of a reverse-S sigmoid in its northern part, which is also consistent with the recent analyses by Mitra et al. (2018), Yan et al. (2018), and Price et al. (2019). Although cases of inconsistency between the chirality of the source region and that of the erupted CME have been observed (e.g., Chandra et al. 2010), for most of the events the two are found to match (Bothmer & Schwenn 1998; Palmerio et al. 2018). For this reason, in the following analysis we start by assuming the erupted structures are characterized by a negative chirality as their source region.

Tilt of the Flux Ropes—In order to estimate the orientation of the flux ropes at the Sun, we use proxies based on the orientation of PEAs and PILs (Möstl et al. 2008; Palmerio et al. 2017, 2018). As shown in Figure 1, the PEA forming after the eruption of CME1 was confined to the southern portion of the AR/PIL and exhibited an approximately north–south orientation. For CME2 and CME3, we observe PEAs developing along the whole PIL structure. Although a global direction from southeast to northwest can be identified, the shape of such PEAs appears to be bent in a reverse-S shape. This reflects the complexity of the underlying PIL system and makes the determination of an unambiguous flux rope tilt based on such observations extremely difficult. Similar conclusions about the initial flux rope tilts can be obtained by considering the locations of coronal dimmings and flare ribbons (Figure 3). Combining these tilts with the information about flux rope chirality and magnetic polarity regions from HMI magnetograms, we recover an ENW flux rope type for CME1 and intermediate ENW–NWS flux rope types for CME2 and CME3 in the lower corona (using the same classification as Bothmer & Schwenn 1998; Mulligan et al. 1998; Palmerio et al. 2018).

Zoom In Zoom Out Reset image size

To constrain the orientation of the flux ropes in the upper corona (∼5 R_⊙ to ∼20 R_⊙), we consider the results from the GCS fittings. The derived tilts suggest that the axial magnetic field of CME1 and CME2 is oriented parallel to the solar equatorial plane (γ_CME = 0°), while that of CME3 has an inclination of ∼40° with respect to the solar equator. According to reconstruction of the heliospheric current sheet (HCS) by potential-field source-surface (Altschuler & Newkirk 1969) models (from the Global Oscillation Network Group (GONG): https://gong.nso.edu/data/magmap/pfss.html), the derived CME tilts are (quasi-) aligned to the HCS (Yurchyshyn 2008; Isavnin et al. 2014). Using the GCS tilt estimates as given in Table 1 and assuming the flux rope chirality at the source region to be preserved in spacetime, two possible flux rope types are possible for each CME. CME1 and CME2 are most probably either SEN flux rope types or NWS flux rope types. CME3 can be associated either with an intermediate ENW–NWS flux rope or with an intermediate WSE–SEN flux rope. We note, however, that the γ_CME parameter is associated with the highest uncertainties (Thernisien et al. 2009), and it is known to be very sensitive to the subjectivity involved in performing the GCS fitting (see, e.g., Figure 5 in Shen et al. 2018, for an alternative fitting of the CMEs using the GCS model, leading to a quite different interpretation of the tilt angles in the corona). For this reason, we consider its reconstruction particularly uncertain.

Overall, the comparison of the CME flux rope types recovered in the lower and upper corona suggests that CME1 and CME2 underwent considerable rotations (rotation ≥90° for CME1, ≥45° for CME2, and ≥0° for CME3). Similar conclusions are reached when considering the GCS fitting performed by Shen et al. (2018). Such estimates are in the upper range of reported values (Lynch et al. 2009; Vourlidas et al. 2011; Isavnin et al. 2014; Kay et al. 2015) but are not surprising given the sympathetic nature of CME1 and CME2, consistent with a scenario of early interactions and strong magnetic forces that may have led to significant CME rotations (Kay et al. 2015).

Association with Interplanetary Structures—In order to further confirm the associations between the CMEs and their interplanetary counterparts described in Section 2.3, we also compare the helicity sign and the flux rope types of the corresponding ejecta at the Sun with the magnetic structures recovered in situ.

The magnetic field rotations observed in association with E1 (most probably associated with the structure resulting from the merging of CME1 and CME2) indicate that the ejecta can be described as a left-handed flux rope characterized by a ∼45° inclination between an SEN- and an ENW-type at 1 au, although the presence of shock S2 propagating through E1 makes the interpretation of the flux rope type particularly problematic because of the disturbed properties in the downstream region. The magnetic field rotations within E2 indicate that the second ejecta can be described as a left-handed low-inclination flux rope of SEN type at 1 au. The lack of magnetic field rotations within E3 makes it difficult to determine the flux rope type of this structure. Overall, the chiralities of the flux ropes recovered from in situ observation at Earth are consistent with the chiralities recovered from source region images. On the other hand, occurrence of significant rotations in the corona and/or interplanetary space is needed to explain the different flux rope orientations recovered at the source region, in the upper corona, and at 1 au. Although deviating from typical scenarios, this aspect can be interpreted as a consequence of the multiple coronal and interplanetary interactions that occurred among the CMEs under study, as well as with other heliospheric structures (e.g., the HCS).

3.2.2. Reconnected Magnetic Fluxes

In order to perform EUHFORIA simulations with the spheromak CME model, next to observational estimations of the flux rope tilts and chiralities, we also need an estimate of the amount of magnetic flux contained within the magnetic structure. As a proxy, we use the flux (φ_r) that reconnected in association with each CME eruption. To have a robust estimate of φ_r, we compare the results obtained from the analysis of a variety of posteruptive signatures (such as flare ribbons, coronal dimmings, and PEAs) in the AR involved in the eruptions.

Statistical Relations—We first estimate the reconnected fluxes using CME–flare statistical relations proposed in previous works. Among all, we focus on the relations between the flare peak intensity in soft X-rays and the reconnected flux derived from flare ribbons and coronal dimmings (e.g., Kazachenko et al. 2017; Dissauer et al. 2018a; Tschernitz et al. 2018) and on the relations between the CME speed and the reconnected flux derived from PEAs (Pal et al. 2018). The major advantages in using statistical relations instead of an in-depth analysis of single events are the applicability to a larger set of events (not restricted to eruptions originating close to the disk center and characterized by specific posteruptive signatures) and the simplicity of use, which makes them potentially suitable for operational forecasting applications, as they can be used to routinely initialize the parameters used by physics-based flux rope CME models running in forecasting mode.

Tschernitz et al. (2018) studied a set of 51 flares ranging between B3 and X17 in GOES class, reporting a very tight correlation (Pearson correlation coefficient r_P = 0.92 in log–log space) between the flare peak intensity I_SXR (in units of W m⁻²) and the reconnected flux φ_r (in units of Mx) estimated from flare ribbons:

$$\begin{eqnarray}&&\mathrm{log}({\varphi }_{r})=24.21+0.58\mathrm{log}({I}_{\mathrm{SXR}}).\end{eqnarray} \tag{ 1 }$$

Considering a larger sample of about 3000 flares ranging from C1 to X5 in GOES class (RibbonDB catalog, http://solarmuri.ssl.berkeley.edu/~kazachenko/RibbonDB/), Kazachenko et al. (2017) reported a correlation of

$$\begin{eqnarray}&&\mathrm{log}({\varphi }_{r})=24.72+0.64\mathrm{log}({I}_{\mathrm{SXR}}),\end{eqnarray} \tag{ 2 }$$

with Spearman's rank correlation coefficient r_S = 0.66. Correcting for the different definitions of φ_r used by Kazachenko et al. (2017) (φ_r = total (unsigned) magnetic flux) compared with Tschernitz et al. (2018) (φ_r = average of the positive and negative fluxes), the relation in Equation (2) becomes

$$\begin{eqnarray}&&\mathrm{log}({\varphi }_{r})=24.42+0.64\mathrm{log}({I}_{\mathrm{SXR}}),\end{eqnarray} \tag{ 3 }$$

where φ_r is now defined to be consistent with the definition used by Tschernitz et al. (2018). Considering the large span of the flare GOES classes associated with the three CMEs under study (M1.7 to X9.3) and the higher correlation coefficients reported by Tschernitz et al. (2018) and Kazachenko et al. (2017) compared with other studies, in the following we use Equations (1) and (3) to identify the most probable range of flare ribbon reconnected flux values φ_r associated with each CME event.

In addition to flare ribbons, we consider coronal dimmings as a secondary signature to estimate the reconnected flux during the eruptions under study based on flare peak intensities. Dissauer et al. (2018a) performed a statistical analysis based on coronal dimming regions observed in association with a set of 62 CME events, reporting a correlation between the flare peak intensity and reconnected flux estimated from coronal dimmings equal to

$$\begin{eqnarray}&&\mathrm{log}({\varphi }_{r})=23.26+0.42\mathrm{log}({I}_{\mathrm{SXR}}),\end{eqnarray} \tag{ 4 }$$

with r_P = 0.62 (in log–log space).

Applying Equations (1), (3), and (4) to the flare peak intensities observed in association with the three CME events under study (obtained from the NOAA SWPC data archive, ftp.swpc.noaa.gov/pub/warehouse), that is, ${I}_{\mathrm{SXR}}=1.7\times {10}^{-5}$ W m⁻² (CME1), ${I}_{\mathrm{SXR}}=5.5\times {10}^{-5}$ W m⁻² (CME2), and ${I}_{\mathrm{SXR}}\,=9.3\times {10}^{-4}$ W m⁻² (CME3), we obtain an estimate of the reconnected fluxes based on statistical relations as listed in Figure 4.

Zoom In Zoom Out Reset image size

In a recent study, Pal et al. (2018) derived statistical relations linking the reconnected fluxes obtained from the Flux Rope from Eruption Data (FRED; Gopalswamy et al. 2017) method, using PEAs as a primary signature to calculate the reconnected flux, and the CME 3D speed v_3D (in units of km s⁻¹) in the corona estimated by applying the GCS fitting technique. Based on 33 CME events, they reported a correlation of (r_P = 0.66)

$$\begin{eqnarray}&&{v}_{3{\rm{D}}}=327\,{\varphi }_{r}^{0.69}.\end{eqnarray} \tag{ 5 }$$

Inverted, the relation above is

$$\begin{eqnarray}&&{\varphi }_{r}={\left(\displaystyle \frac{{v}_{3{\rm{D}}}}{327}\right)}^{1/0.69},\end{eqnarray} \tag{ 6 }$$

which allows us to estimate the reconnected flux once the 3D speed of the CME is known. Using as input the 3D speeds reconstructed from the GCS fitting (v_CME) listed in Table 1, the reconnected fluxes for the three CMEs obtained from Equation (6) are provided in Figure 4.

Reconnected Fluxes from Single-event Analysis—In order to obtain an event-based estimate of the reconnected fluxes and to assess the performance of statistical relations in the specific case of the events considered, we complement the φ_r values recovered from statistical relations with results from the single-event analysis of each of the three eruptions under study. In order to consider a broad spectrum of CME-related signatures, we estimate the reconnected flux involved in each CME eruption using the following methods:

• 1.  
  A method to identify flare ribbon areas based on the work by Kazachenko et al. (2017). Flare ribbons are detected in SDO/AIA 1600 Å images with 24 s cadence by applying a cutoff threshold based on the median image intensity. Images are taken between 30 minutes before and 3 hr after the start of the flare associated with each CME under study.
• 2.  
  A method to identify coronal dimming areas based on the work by Dissauer et al. (2018b). Coronal dimmings are detected based on a thresholding method that is applied to logarithmic base ratio SDO/AIA 211 Å images. Similar to the flare ribbons, dimming pixels are detected between 30 minutes before the flare and up to 3 hr after the flare.
• 3.  
  A method to identify PEA areas based on the FRED method described by Gopalswamy et al. (2017). For each CME under study, we use SDO/AIA 193 Å taken around the moment of maximal extension of the PEAs. We note that instead of performing a manual identification of the PEA areas as in Gopalswamy et al. (2017), we here employ an automatic identification algorithm using a cutoff threshold based on the median image intensity.

To recover φ_r based on the areas identified with the methods above, we use full-disk SDO/HMI line-of-sight magnetograms (hmi.m_720s) taken about 30 minutes before the flare start times (shown in Figure 3).

Results from the single-event analyses are listed in Figure 4, while the location of the various EUV signatures overplotted on HMI magnetograms for the three events is shown in Figure 3. The large spread in the recovered φ_r values from single-event analyses reflects the different areas covered by the signatures considered (ribbons, dimmings, and PEAs), and we note that large uncertainties affect the estimation of φ_r, that is, up to ±50% of the measured value, as reported by various studies (Qiu et al. 2007; Gopalswamy et al. 2017; Pal et al. 2017; Temmer et al. 2017; Dissauer et al. 2018a; Tschernitz et al. 2018). Despite the scatter, the different values recovered can therefore be considered consistent within the (large) error bars. At the same time, these results highlight that when using this methodology one should aim to recover an order of magnitude for the reconnected flux, rather than a precise estimate of it.

When averaging out the variability of the φ_r results obtained from the different methods, we recover very similar results between the single-event analyses and the statistical analyses of CME1 and CME2 (rows in bold in Figure 4). The results from the two approaches for CME3, on the other hand, appear somehow less consistent. We note that a contributing factor to this variability comes from the fact that the CME erupted about 40° away from the solar disk center, which is close to the limit of the applicability of the single-event analysis methods due to the increased projection effects when moving away from the disk center (e.g., Gopalswamy et al. 2017; Dissauer et al. 2018b). This most probably resulted in higher uncertainties affecting the reconstruction of φ_r, that is, the results from single-event and statistical analyses are still consistent among each other due to the larger error bars. Overall, these results indicate that for the specific events considered, using different statistical methods to quantify the reconnected fluxes provides results that are on average consistent with those of more sophisticated single-event analysis methods. Such statistical methods are fast and easy to apply as they only require as input easy to use data products, such as the peak intensity of the CME-associated flares or the 3D speed of the CME in the corona as recovered from the GCS fitting or other reconstruction methods. In the context of operational forecasting operations, these results therefore highlight how statistical methods could represent promising potential alternatives to otherwise time-consuming single-event analysis methods.

Magnetohydrodynamic (MHD) simulations of the evolution of magnetized CMEs in the heliosphere are powerful complementary tools to observations, as they can provide information on the evolution of CME structures in 3D space that is often difficult to infer from remote sensing or in situ observational data alone—particularly in cases limited by single-spacecraft in situ measurements as here. In this work, we investigate the heliospheric propagation and interaction of the three successive CMEs under study using the EUHFORIA model.

EUHFORIA is a recently developed 3D MHD model of the inner heliosphere (Pomoell & Poedts 2018) that allows the modeling of the background solar wind and CME events using a linear force-free spheromak flux rope model (Verbeke et al. 2019b). The model is composed of (1) a semiempirical Wang–Sheeley–Arge-like (Arge et al. 2004) global coronal model that provides the background solar wind parameters at the heliocentric radial distance of 0.1 au, starting from synoptic maps of the photospheric magnetic field and (2) a time-dependent 3D MHD model of the inner heliosphere (between 0.1 and 2 au). In the heliosphere, it is possible to model solely the ambient solar wind or to also include CMEs, which are inserted into the heliospheric domain via time-dependent boundary conditions at 0.1 au. In this work we use a computational domain for the heliospheric model ±180° in longitude (ϕ), ±80° in latitude (θ), and from 0.1 to 2 au in the radial direction (D). The use of a sufficiently high spatial resolution is particularly necessary to better resolve shock structures in the simulation domain, which are extremely important in the context of global CME–CME interactions, as discussed in Section 4. For this reason, our simulations are performed using a homogeneous grid with 1024 cells in the radial direction (corresponding to a radial resolution of ΔD = 0.00186 au = 0.4 R_⊙ per cell) and with 2° resolution in the longitudinal and latitudinal directions. As input for the coronal model we use the magnetogram synoptic map generated by the GONG on 2017 September 4 at 00:04 UT (https://gong.nso.edu/data/magmap/QR/bqs/201709/mrbqs170904/mrbqs170904t0004c2194_055.fits.gz). All simulations are carried out with EUHFORIA, version 1.0.4. In the following, all coordinates are given in the heliocentric earth equatorial coordinate system, unless specified otherwise.

3.3.1. CME Modeling

In order to quantify the resulting geoeffectiveness of the CME events as predicted by EUHFORIA, we calculate the predicted Dst index from the modeled time series at Earth after conversion into geocentric solar magnetospheric (GSM) coordinates, using the AK2 model proposed by O'Brien & McPherron (2000a, 2000b). Such predictions are compared with hourly Dst values from the World Data Center for Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp/dstdir/), and with predictions obtained by applying the same coupling function to 1 minute solar wind time series provided by Wind.

The quality of the CME geoeffectiveness predictions is then quantified by comparing the minimum Dst predicted by EUHFORIA, with predictions obtained from Wind measurements and with actual observations.

An overview of the event in the heliosphere as obtained from EUHFORIA run 01-01-01 is presented in the top panels of Figure 5, showing the radial speed (v_r) in the ecliptic plane and in the meridional plane containing Earth, at three different times in the simulation (the radial speed, number density, and colatitudinal magnetic field plots for all the runs performed are provided in Appendix B). The position of the leading edges of the CMEs as marked by their interplanetary shocks at various times are indicated by the black arrows. The bottom panel of Figure 5 shows the v_r prediction at Earth compared with in situ observations from Wind (a comparison of the EUHFORIA time series at and around Earth for all runs performed and in situ observations from Wind are included in Appendix C). The evolution and interactions of the CMEs in terms of 3D topology of their magnetic field lines, and of the plasma β and B_z (D/1 au)² in the meridional plane containing Earth, are provided in Figure 6. The simulation results show that CME1 and CME2 interacted already during the insertion in the heliospheric domain, as expected given the very close insertion times derived from the GCS reconstruction (listed in Table 2) and from coronagraph images. In the heliosphere, the two CMEs propagate as a merged structure (hereafter CME1+CME2) all the way from 0.1 to 2 au. As further discussed in Section 4.3, its predicted arrival time at Earth supports the interpretation of S1 as the interplanetary shock driven by the CME1+CME2 merged structure. Figure 5 also shows CME3 first propagating through the perturbed solar wind in the wake of CME1+CME2 and then interacting with them. The interaction with CME1+CME2 appears to be still ongoing at the time CME3 reaches 1 au, as indicated by the shocks of CME1+CME2 and CME3 being still distinct as predicted at Earth, also supporting the interpretation of S2 being the interplanetary shock driven by CME3.

Zoom In Zoom Out Reset image size

Because of the intrinsic difficulties in identifying the boundaries of the various ejecta and related structures in 3D, to characterize the phases of the interaction between CME1+CME2 and CME3 as they propagate in interplanetary space, we apply an approach based on 1D cuts taken along the Sun–Earth line (i.e., approximately the direction of propagation of the structures eventually arriving at Earth) at various times in our simulation (run 01-01-01). At each time, we identify the location of the leading and trailing shocks (driven by CME1+CME2 and CME3, respectively) and the boundaries of the two respective ejecta by considering a low-β region to correspond to an ejecta. We also characterize the position of each ejecta in terms of their geometrical center. Figure 7 illustrates the main MHD parameters along the Sun–Earth line, together with the location of the various shock and ejecta structures, at three different times in the simulation, which clearly associate with three of the four typical phases of CME–CME interactions (as defined by Lugaz et al. 2005). In particular, Figure 7(a) associates with phase 1, corresponding to a period before the start of the interaction. Figure 7(b) associates with phase 2, corresponding to the shock–ejecta interaction phase, and Figure 7(c) with phase 3, corresponding to the shock–sheath interaction phase. Phase 4, corresponding to the shock–shock interaction phase, is not present in Figure 7 as this phase occurs after the CMEs have already left the simulation domain. A more detailed investigation of each interaction phase is given in the following paragraphs.

Zoom In Zoom Out Reset image size

We also calculate the speed of the shocks and ejecta centers in the reference frame of the Sun. The propagation of each structure in terms of time–distance and time–speed profiles is shown in Figure 8.

Zoom In Zoom Out Reset image size

At each time available, we characterize the properties of the trailing (CME3) shock by determining its speed in the reference frame moving with the upstream plasma by applying the Rankine–Hugoniot relations (assuming the 1D cut to be parallel to the shock normal) as

$$\begin{eqnarray}&&{v}_{\mathrm{shock}}=\displaystyle \frac{{v}^{\mathrm{down}}\,{\rho }^{\mathrm{down}}-{v}^{\mathrm{up}}\,{\rho }^{\mathrm{up}}}{{\rho }^{\mathrm{down}}-{\rho }^{\mathrm{up}}},\end{eqnarray} \tag{ 7 }$$

where "down" and "up" refer to quantities calculated in the downstream and upstream shock regions, respectively. In addition to the shock speed, we also calculate the Alfvén speed (${v}_{{\rm{A}}}^{\mathrm{up}}$), sound speed (${c}_{{\rm{s}}}^{\mathrm{up}}$), plasma β (${\beta }^{\mathrm{up}}$) in the upstream region, and the shock Alfvén and magnetosonic (fast) Mach numbers (M_A, M_ms), together with the density compression ratio $r={\rho }^{\mathrm{down}}/{\rho }^{\mathrm{up}}$. The evolution in spacetime of all these quantities is shown in Figures 9(a)–(e).

Zoom In Zoom Out Reset image size

To clarify the relation between the evolution of the shock driven by CME3 and the evolution of the two ejecta, we also calculate the average scaled kinetic (${\tilde{\epsilon }}_{\mathrm{kin}}$), magnetic (${\tilde{\epsilon }}_{\mathrm{mag}}$), and thermal (${\tilde{\epsilon }}_{\mathrm{therm}}$) energy densities within the two ejecta as

$$\begin{eqnarray}&&\begin{array}{r}\left\{\begin{array}{l}{\tilde{\epsilon }}_{{\rm{k}}{\rm{i}}{\rm{n}}}= \langle \displaystyle \frac{1}{2}\tilde{\rho }{\tilde{v}}^{2} \rangle \\ {\tilde{\epsilon }}_{{\rm{m}}{\rm{a}}{\rm{g}}}= \langle \displaystyle \frac{{\tilde{B}}^{2}}{2{\mu }_{0}} \rangle \\ {\tilde{\epsilon }}_{{\rm{t}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{m}}}= \langle \displaystyle \frac{2\tilde{n}{k}_{{\rm{B}}}\tilde{T}}{\gamma -1} \rangle \end{array}\right.\end{array},\end{eqnarray} \tag{ 8 }$$

where γ = 1.5 is the adiabatic index (consistent with Pomoell & Poedts 2018), k_B is the Boltzmann constant, $\langle x\rangle$ indicates the average along the radial coordinate taken over the ejecta extension calculated above, and $\tilde{x}$ indicates a scaled quantity. Scaled quantities account for the radial evolution of the various CME parameters and are employed in order to better compare the energy densities at different times in the same run (corresponding to different distances from the Sun) and at the same times in different runs. This is needed because the propagation, that is, the radial distance, of the leading ejecta is greatly affected by the presence (run 01-01-01) or lack (run 01-01-00) of the interaction with CME3, and to compare its properties at the same time in the two runs, a correction for the different distance from the Sun is needed. In particular, we use the following scaling relations: $\tilde{v}=v$, $\tilde{n}=n\,{(D/1\mathrm{au})}^{2.32}$, $\tilde{\rho }=\rho \,{(D/1\mathrm{au})}^{2.32}$, $\tilde{B}=B\,{(D/1\mathrm{au})}^{1.85}$, $\tilde{T}=T\,{(D/1\mathrm{au})}^{0.32}$, $\tilde{P}=P\,{(D/1\mathrm{au})}^{2.64}$, with the exponents derived from statistical studies of CMEs in the inner heliosphere (from Liu et al. 2005; Gulisano et al. 2010). Together with the scaled energy densities, we also calculate the scaled radial size of the ejecta as $\tilde{S}=S\,{(D/1\mathrm{au})}^{0.45}$, where the scaling factor is taken as the lower limit of the range obtained by Gulisano et al. (2010), which proves to scale well the decrease of the radial size of CME1+CME2 in the case without interaction (dashed red line in Figure 9(f)). We also apply the same approach to calculate the average energy densities in the region between the leading edge of the CME1+CME2 ejecta and its shock, corresponding to the sheath ahead of CME1+CME2. The evolution in time of the sizes and scaled energy densities associated with these three structures is shown in Figures 9(f)–(i).

Finally, we calculate the amplification of the helioeffectiveness of CME1+CME2 due to its interaction with CME3 in terms of the magnetic field compression caused by the propagation of the shock driven by CME3 through the CME1+CME2 ejecta. We therefore compute the maximum $\tilde{B}$ (${\tilde{B}}^{\max }$) and minimum $\tilde{{B}_{z}}$ (${\tilde{B}}_{z}^{\min }$) within the boundaries of the ejecta, where ${\tilde{B}}_{z}={B}_{z}\,{(D/1\mathrm{au})}^{1.85}$ assuming that the magnetic field components in a CME scale with the same behavior as that of the magnitude B. The helioeffectiveness amplification factors ($\tilde{A}$) due to the interaction are then calculated as ${\tilde{A}}_{B}={\tilde{B}}_{010101}^{\max }/{\tilde{B}}_{010100}^{\max }$ and ${\tilde{A}}_{{Bz}}={\tilde{B}}_{z,010101}^{\min }/{\tilde{B}}_{z,010100}^{\min }$, that is, taking the ratio of the values from run 01-01-01 and run 01-01-00. Results are shown in Figure 9(j).

Phase 1: Preinteraction—In Figure 7(a), the shock and ejecta associated with CME3 propagate through a high-β (β ∼ 10) solar wind perturbed by the earlier passage of the preceding CME1+CME2 structure. The rear edge of CME1+CME2 is still unaffected by the presence of CME3, indicating no direct interaction has yet occurred. This is also visible from the time–distance profile in Figure 8(a). At this time the shock driven by CME3 is propagating with a speed of ∼2130 km s⁻¹ and its ejecta with a speed of ∼650 km s⁻¹, significantly higher than the shock and ejecta associated with CME1+CME2 (both moving at ∼600 km s⁻¹). CME3 is progressively approaching CME1+CME2, as is clearly visible in Figure 8(a).

During this phase, the scaled energy densities of CME1+CME2 are approximately constant in spacetime (Figures 9(g)–(i)), implying that the (nonscaled) energy densities are decreasing with radial distance due to the known interplay of expansion (which converts magnetic and thermal energy into kinetic energy) and drag (which slows down the CME reducing its kinetic energy) (Cargill 2004; Vršnak et al. 2013). Between 18:00 UT and 21:00 UT on September 6, CME3 is still being inserted in the domain, and it expands (Figures 8(a), 9(f)) due to the propagation through the rarefied perturbed solar wind (Liu et al. 2014b; Temmer & Nitta 2015; Temmer et al. 2017), as indicated by its rapidly growing radial size and the rapid decrease of ${\tilde{\epsilon }}_{\mathrm{kin}}$, ${\tilde{\epsilon }}_{\mathrm{therm}}$, and ${\tilde{\epsilon }}_{\mathrm{mag}}$ (Figures 9(g)–(i)).

Phase 2: Shock–Ejecta Interaction—On September 6 at 22:00 UT, around 0.45 au, the shock driven by CME3 starts interacting with the preceding CME1+CME2 ejecta (see Figures 7(b) and 8(a)). Figures 9(a)–(e) show that as the shock enters the preceding ejecta, the upstream plasma β decreases by a factor of 100, the upstream density decreases by a factor of 8, and the Alfvén speed increases by a factor of 12, while the sound speed remains almost the same. The reduced density and higher magnetic field contribute to increase v_A and to lower M_A; hence the shock quickly accelerates from a speed of ∼1300 km s⁻¹ to a speed of ∼2120 km s⁻¹ on September 6 at 23:00 UT, while the density compression ratio of the shock decreases from r ≃ 5 (with 5 equal to the theoretical maximum for ideal MHD with γ = 1.5) to a value of about 1.9. After reaching the speed of ∼2120 km s⁻¹ , the shock decelerates to a speed of ∼970 km s⁻¹ when reaching the core (densest part) of the CME1+CME2 ejecta. The compression ratio is found to increase again to 2.9 a few hours later. The shock remains at all times a fast shock, with a minimum M_ms of 4.3 in the frame moving with the upstream plasma. Following the passage through the denser ejecta core, the shock propagates through a rarefied region where it accelerates up to ∼2390 km s⁻¹ on September 7 at 10:00 UT, right before exiting the leading edge of the ejecta.

During this phase, the shock–ejecta interaction is most efficient in amplifying the helioeffectiveness of the CME1+CME2 ejecta, as is visible in the ${\tilde{A}}_{B}$ and ${\tilde{A}}_{{Bz}}$ amplification factors in Figure 9(j). The shock driven by CME3 compresses, accelerates, and heats the plasma in the preceding ejecta as it propagates through it, enhancing the density, speed, temperature, and magnetic field in the downstream region. This is visible as an increase in ${\tilde{\epsilon }}_{\mathrm{kin}}$, ${\tilde{\epsilon }}_{\mathrm{mag}}$, and ${\tilde{\epsilon }}_{\mathrm{therm}}$ of CME1+CME2 compared with the simulation without interaction (i.e., red continuous and dashed lines in Figures 9(g)–(i)). The acceleration of CME1+CME2 is also clearly visible in the increased speed of the ejecta center in Figure 8(b). The radial size of CME1+CME2 also rapidly decreases (Figures 8(a) and 9(f)), consistent with a radial compression by the trailing shock and ejecta (Vandas et al. 1997; Schmidt & Cargill 2004; Lugaz et al. 2005; Xiong et al. 2006). During this phase we observe a steady increase of the ${\tilde{A}}_{B}$ and ${\tilde{A}}_{{Bz}}$ amplification factors ("growth phase"), until maximum amplification factors of 3.3 for ${\tilde{A}}_{B}$ and 2.8 for ${\tilde{A}}_{{Bz}}$ are reached around September 7 at 14:00 UT, that is, after the end of this phase ("maximum phase"). The maximum helioeffectiveness amplification occurs when the CME1+CME2 ejecta is around 0.9 au in simulation 01-01-01, that is, close to the location of Earth.

Phase 3: Shock–Sheath Interaction—Eventually, the shock driven by CME3 reaches the front of the CME1+CME2 ejecta and starts interacting with the dense sheath of plasma ahead of it (Figures 7(c) and 8(a)). In our simulations, this occurs around September 7 at 11:00 UT, that is, close to the moment when it passes Earth/1 au. We note that this phase starts at smaller heliocentric distances than in the observed in situ data by the spacecraft at L1, where the shock appears still to be fully embedded in the ejecta. This is most probably due to a combination of two factors: (1) an overestimated drag in our simulations (lower speed and higher density in the solar wind ahead of CME1+CME2 compared with observations, see Appendix C)—this particularly contributed to slowing down CME1+CME2 (postponing its arrival time) and (2) uncertainties in the CME initial speeds based on the GCS and φ_r reconstructions, which may have contributed to predicting an early arrival time of the shock of CME3 compared with in situ observations (bottom panel in Figure 5). During this phase the shock driven by CME3 enters again a high-β (β > 50) environment, characterized by a low Alfvén speed (∼15 km s⁻¹) and a sound speed comparable to that in the ejecta (∼95 km s⁻¹). The density in the sheath is about 2 orders of magnitudes higher than in the ejecta. As the shock enters this new region, its speed with respect to the upstream plasma drops below 850 km s⁻¹, due to the increased upstream density. Its radial speed is still slightly larger than the speed of the leading shock driven by CME1+CME2 (∼800 km s⁻¹ compared with ∼630 km s⁻¹), so the CME3-driven shock continues to get closer to the shock ahead, but at a slower rate. We expect them eventually to merge, but the rate of approach is so low that this would take place only beyond the outer boundary of the simulation domain (i.e., beyond 2 au).

Analyzing the energy of each substructure in this phase, we observe that the CME3 shock starts heating and compressing the CME1+CME2 sheath ahead, as indicated by the decrease in its radial size and by the increase in the associated ${\tilde{\epsilon }}_{\mathrm{kin}}$, ${\tilde{\epsilon }}_{\mathrm{mag}}$, and ${\tilde{\epsilon }}_{\mathrm{therm}}$ compared with the simulation without interaction (green continuous and dashed lines in Figures 9(f)–(i)). At the same time, the scaled radial size of CME1+CME2 (CME3) slowly decreases (increases) until September 7 at 23:00 UT, when the trends invert. Most notably, in this phase we also observe that the ${\tilde{A}}_{B}$ and ${\tilde{A}}_{{Bz}}$ amplification factors follow a decreasing trend starting from September 7 at 14:00 UT. We explain this by noting that while ${\tilde{\epsilon }}_{\mathrm{mag}}$ in CME1+CME2 starts decreasing from September 7 at 14:00 UT, its ${\tilde{\epsilon }}_{\mathrm{kin}}$ and ${\tilde{\epsilon }}_{\mathrm{therm}}$ continue increasing until September 8 at 03:00 UT. The earlier peak of ${\tilde{\epsilon }}_{\mathrm{mag}}$ compared with those of ${\tilde{\epsilon }}_{\mathrm{kin}}$ and ${\tilde{\epsilon }}_{\mathrm{therm}}$ suggests a conversion of the magnetic energy accumulated by CME1+CME2 during the shock–ejecta interaction phase, into kinetic and thermal energy of the same structure. The acceleration of CME1+CME2 due to this energy conversion is particularly visible in Figure 9(b). During this energy-conversion phase, the radial size of CME1+CME2 slightly increases as indicated by the constant scaled radial size (Figure 9(g)), most probably due to the presence of CME3 at the back of CME1+CME2 preventing a further expansion. At the end of this phase CME1+CME2 and CME3 move at almost the same speed (∼670 km s⁻¹, Figure 9(b)). The low expansion of CME1+CME2 is therefore consistent with numerical and observational studies (Xiong et al. 2006, 2007; Gulisano et al. 2010; Lugaz et al. 2012), which highlighted how the expanding behavior of the leading ejecta after the end of the main interaction phase depends on the delicate interplay between the natural tendency for the CME to expand and the compressing action of the trailing ejecta.

The most relevant consequence of the expansion and conversion of magnetic energy into acceleration and heating is that CME1+CME2 progressively "forgets" the amplification of B_z caused by the interaction with CME3, and it slowly returns to B_z conditions similar to the case without interaction ("decay phase"). This suggests that one of the key factors at the origin of the intense geomagnetic storm triggered by the 2017 September 4–6 CMEs was their arrival at Earth during the phase of maximum helioeffectiveness amplification due to the interaction of CME1+CME2 with CME3.

4.3.1. B_z and Arrival Time Predictions

Together with the intensity and orientation of their internal magnetic field at Earth (Savani et al. 2015, 2017; Palmerio et al. 2018), another key CME parameter that the community has long tried to predict is the arrival time (Möstl et al. 2014; Riley et al. 2018; Verbeke et al. 2019a). Current estimates of prediction uncertainties for this parameter are about ±10 hr depending on the exact metric considered and on the model used (Riley et al. 2018; Wold et al. 2018), with numbers increasing for more complex events. In the context of interacting CMEs, the arrival time of individual CMEs can be significantly affected by two different effects impacting the CME kinematics: (1) the preconditioning of the ambient solar wind due to the passage of previous CMEs (Temmer et al. 2017) and (2) direct CME–CME interactions (CME–CME collisions; see, e.g., Liu et al. 2014a; Shen et al. 2017). Given the relevance of CME–CME interactions in affecting the kinematics of individual CMEs and, consequently, the prediction of their arrival times at Earth, we address here how much the B_z and arrival times at Earth of the shocks and ejecta associated with CME1+CME2 and CME3 were affected by their interaction. To do so, we start from the identification of the shocks and ejecta boundaries performed in Section 4.2, and then we extract and compare times at which the shock and leading edge of the ejecta associated with CME3 arrived at 1 au, in simulations with (run 01-01-01) and without (run 00-00-01) the presence of the preceding CMEs. To assess the impact of the interaction on the arrival time of CME1+CME2, we similarly extract and compare the position in time of the shock and leading edge of the ejecta associated with CME1+CME2 obtained from run 01-01-01 and run 01-01-00.

Effect on CME1+CME2—In Figure 10(a) we plot the predicted minimum B_z associated with CME1+CME2 from the EUHFORIA time series at Earth and the surrounding virtual spacecraft, together with the observed minimum B_z associated with E1 (see Figure 2). We observe that the simulation of CME1 alone (run 01-00-00) predicts a modest minimum B_z of −10 nT (−13 to −8 nT considering spacecraft separated by an angular distance Δσ = 5° from Earth and −14 to −4 nT considering spacecraft separated by Δσ = 10°). While the strength of the magnetic field within CME1 at Earth is similar to the unshocked region in E1, due to its relatively low initial speed compared with CME2 (Table 1), the shock driven by CME1 without the inclusion of CME2 in our simulations is predicted to arrive at Earth on September 7 at 09:33 UT, that is, about 10 hr later than the observed S1. When adding CME2 to the simulations (run 01-01-00), the early interaction of CME1 and CME2 results in a merged structure that propagates in the heliosphere with a speed close to that of the fastest CME involved, that is, CME2, arriving at Earth on September 7 at 02:23 UT, that is, only 3 hr later than the actual arrival time of S1. The predicted minimum B_z of this combined structure is −14 nT (−15 to −9 nT considering spacecraft separated by Δσ = 5° from Earth and −16 to −5 nT considering spacecraft separated by Δσ = 10°), less than 5 nT lower than the unshocked (upstream) region observed in E1. The start of the ejecta associated with CME1+CME2 is predicted on September 7 at 23:13 UT, only 3 hr later than the starting time of ejecta E1 based on Wind in situ observations. When including all three CMEs (run 01-01-01), the predicted minimum B_z of CME1+CME2 at Earth drops to −35 nT, that is, very close to the minimum observed value of −32 nT. When accounting for uncertainties related to the initial CME directions of propagation, the predicted value varies between −38 and −19 nT for a Δσ = 5° separation from Earth and between −41 and −4 nT for a Δσ = 10° separation from Earth. Overall, we also observe a good level of agreement between the predicted and the observed B and B_z time series at and around Earth (see Appendix C). While the predicted arrival time of the shock driven by CME1+CME2 is unaffected by the inclusion of CME3 in our simulations, the ejecta associated with CME1+CME2 arrives about 5 hr earlier, that is, on September 7 at 17:43 UT. The ejecta is therefore predicted to arrive about 2 hr earlier than the starting time of ejecta E1 based on Wind in situ observations. Overall, the very close match between the modeled and observed arrival times and minimum B_z in run 01-01-01 provides strong evidence that S1 and E1 were indeed the interplanetary counterparts of CME1+CME2.

Zoom In Zoom Out Reset image size

By comparing the minimum B_z prediction from runs 01-01-00 and 01-01-01 at and around Earth, we conclude that the presence of S2 and E2 contributed to an increase in the minimum B_z associated with E1 by a factor of 2.5 (ranging between 2.1 and 2.5 for a Δσ = 5° separation from Earth and between 0.8 and 2.6 for a Δσ = 10° separation from Earth), a value consistent with the results obtained by Shen et al. (2018) based on observational arguments and with the analysis in Section 4.2.

Effect on CME3—For completeness, in Figure 10(b) we plot the predicted minimum B_z at and around Earth associated with CME3 in a similar way as was already done for CME1+CME2. Simulation results at Earth from run 00-00-01 are consistent with the observed minimum B_z within E2 (−13 nT compared with −16 nT), while inclusion of the preceding CMEs predicts a minimum B_z of −6 nT. By considering virtual spacecraft in the vicinity of Earth, however, we observe that predictions vary between −6 and −3 nT for Δσ = 5° and between −27 and −2 nT for Δσ = 10°. Therefore, although run 00-00-01 gives a slightly better prediction for the minimum B_z within E2 at Earth, results in the vicinity of Earth from both simulations are overall consistent with the value measured in situ. The larger spread in the minimum B_z in the vicinity of Earth predicted in run 01-01-01 compared with run 00-00-01 may reflect the development of finer structures within CME3, or potential deflections in its trajectory, as a consequence of its interaction with CME1+CME2 (or a combination of the two).

Considering its arrival time at Earth, run 00-00-01 predicts the shock of CME3 to arrive at Earth on September 8 at 08:23 UT, followed by the ejecta starting on September 9 at 01:43 UT. As expected, the inclusion of CME1+CME2 ahead of CME3 (run 01-01-01) anticipated the arrival of the shock by about 15 hr to September 7 at 17:03 UT, and the arrival of the ejecta by about 16 hr, to September 8 at 09:43 UT. This was caused by the solar wind preconditioning induced by CME1+CME2 (which also contributed to the early start of the interaction). These final predictions match well the observed arrival time of S2 (September 7 at 22:38 UT, about 6 hr later than predicted) and the starting time of E2 (September 8 at 11:00 UT, about 1 hr later than predicted) at Wind, supporting their association with CME3 at the Sun. At the same time, we note how the predicted long duration of the ejecta associated with CME3 in run 01-01-01 suggests that it was associated with both the E2 and E3 (see Figure 2), supporting the interpretation that they were indeed both associated with the same CME at the Sun.

4.3.2. Dst Predictions and Geoeffectiveness

To quantify the effect of CME–CME interactions on the actual impact of the merged CME1+CME2 structure on Earth (i.e., its geoeffectiveness), we apply the method presented in Section 3.4 to the EUHFORIA time series extracted at and around the location of Earth. To better quantify by how much EUHFORIA-based Dst predictions are off due to mispredictions of the CME impact parameters, and how much of the discrepancy is actually due to limitations of the specific coupling function used, we also apply the coupling function to Wind in situ measurements. As previously done for the minimum B_z, we separately discuss the predictions obtained for CME1+CME2 (Figure 11(a)) and for CME3 (Figure 11(b)).

Zoom In Zoom Out Reset image size

Minimum Dst Associated with CME1+CME2—As already found for the B_z prediction, Dst predictions associated with CME1+CME2 are highly dependent on the simulated CMEs in different runs (Figure 11(a)). While CME1 alone (run 01-00-00) predicts a minimum Dst of −99 nT (−125 to −69 nT considering spacecraft separated by Δσ = 5° from Earth and −135 to −16 nT considering spacecraft separated by Δσ = 10°), the addition of CME2 (run 01-01-00) leads to a predicted minimum Dst of −119 nT (−128 to >0 nT considering spacecraft separated by Δσ = 5° from Earth and −134 to >0 nT considering spacecraft separated by Δσ = 10°) due to the early interaction between CME1 and CME2. The simulation including all three CMEs (run 01-01-01) finally provides us with an estimate of the contribution of CME–CME interaction processes between CME1+CME2 and CME3 to the geoeffectiveness of CME1+CME2. In particular, we note that the minimum Dst predicted at Earth is ∼−215 nT. When accounting for uncertainties related to the initial CME directions of propagation, the predicted value varies between −227 and −123 nT for a Δσ = 5° separation from Earth and between −229 and −33 nT for a Δσ = 10° separation from Earth. By comparing the minimum Dst prediction from runs 01-01-00 and 01-01-01 at and around Earth, we conclude that the presence of S2 and E2 significantly enhanced the geoeffectiveness of E1 by enhancing the minimum Dst by a factor of 1.8 (ranging between 1.7 and 1.8 for a Δσ = 5° separation from Earth and between 1.7 and 1.3 for a Δσ = 10° separation from Earth, consistent with Shen et al. 2018). In terms of absolute values, the minimum Dst recorded on ground was −142 nT, while the predicted minimum based on Wind in situ measurements was −149 nT, that is, the two values are very close. Therefore, we observe that EUHFORIA run 01-01-01 tends to overestimate the Dst at Earth by a factor of 1.5 compared with actual observations (1.4 when compared with Wind predictions). The overprediction is significantly higher than the one observed in B_z (1.1 compared with Wind observations), most probably due to the larger density (entering the calculation of Dst via the solar wind dynamic pressure (P_dyn); see Equation (3) in O'Brien & McPherron 2000b) predicted by EUHFORIA (see Appendix C).

Minimum Dst Associated with CME3—For completeness, in Figure 11(b) we plot the predicted minimum Dst at and around Earth associated with CME3 in a similar way as was already done for CME1+CME2. Dst estimates reflect what is already found for B_z, that is, that the magnetic field strength within CME3 was not significantly altered by the interaction process, as the minimum Dst predicted in runs with (run 00-00-01) and without (run 01-01-01) the preceding CME1+CME2 are consistent with each other and with actual observations and Wind-based predictions. We also observe an anticorrelation between the minimum B_z and the minimum Dst predicted in runs 00-00-01 and 01-01-01: while the minimum B_z increases in the immediate surroundings of Earth, the minimum Dst decreases. This result is most probably due to the combination of the higher density and speed associated with CME3 in run 01-01-01 compared with run 00-00-01 (contributions due to dynamic pressure (P_dyn) and the dawn-to-dusk component of the electric field (${{VB}}_{s}$); see Equations (2) and (3) in O'Brien & McPherron 2000b) and the fact that in run 01-01-01 the Dst dip caused by CME3 started from a condition of highly disturbed Dst already (see also Kamide et al. 1998; Vennerstrom et al. 2016).

From Section 4.2 we concluded that one of the key factors at the origin of the intense storm triggered by the 2017 September 4–6 CMEs was their arrival at Earth during the phase of maximum amplification of the southward B_z (and consequently, of their helio/geoeffectiveness) due to the interaction of CME1+CME2 with CME3. Moreover, Figure 9(j) highlights the existence of a correlation between the evolutionary phase of CME–CME interactions and the amplification of the helioeffectiveness of the leading ejecta involved in the interaction, at least for the specific CME events and the specific (Sun-to-Earth) direction under study.

In general, for two generic CMEs launched in approximately the same direction, the spatial/temporal windows of the various interaction phases depend on three main parameters: (1) the ambient solar wind through which the CMEs (particularly the preceding one) are propagating, (2) the time interval between the eruptions of the individual CMEs involved, and (3) their relative speed. Therefore, we note that different combinations of such parameters will lead the ejecta to reach Earth or any other location in the heliosphere at different evolutionary phases and hence during different phases of helioeffectiveness amplification. An extensive exploration of the parameter space was performed by Xiong et al. (2007) via 2.5D simulations of interacting CMEs varying the time interval between the eruptions (from 10 to 44 hr) and their relative speeds (from 50 to 800 km s⁻¹). By comparing the resulting geoeffectiveness amplification at Earth in the case of interactions that reached 1 au at different evolutionary stages, they similarly suggested that the evolution stage may be a dominant factor in determining the ultimate geoeffectiveness of interacting CMEs, although a comparison of model results with in situ observations of real CME events was not presented. They also suggested that the exact evolution profile in spacetime of the helioeffectiveness amplification may depend on the impact angle between the spacecraft crossing and the CME apex, with spacecraft locations close to the CME nose more likely to maintain the helioeffectiveness amplification due to the persistent push of the trailing ejecta on the leading ejecta (preventing further expansion of the latter) and spacecraft locations close to the CME flanks exhibiting a decay of the helioeffectiveness amplification due to the narrower extension of the trailing magnetic ejecta compared with its driven shock (which induced the compression). An in-depth analysis of such angular dependencies in EUHFORIA simulations is left for future studies.

The existence of a decay phase in the helioeffectiveness amplification starting after the end of the shock–ejecta interaction clearly has strong implications for the impact of CME–CME interactions at various locations in the heliosphere. In particular, each CME–CME interaction event may be associated with a "helioeffectiveness amplification zone," corresponding to the heliocentric distances associated with the maximum amplification phase for a given combination of CME waiting times and relative speeds. The amplification of the helioeffectiveness of individual CMEs will be null for spacecraft locations closer to the Sun than the distance at which the interaction starts, it will progressively increase for spacecraft locations between the start and the end of the shock–ejecta interaction (growth phase), it will be maximal for spacecraft locations at the outer edge of this distance range (maximum phase), and it will progressively decrease for spacecraft locations farther away from the Sun (decay phase). Although a more extensive study of this impact needs to be addressed in a future work, we speculate that the nonuniform probability distribution of the CME waiting times and relative speeds (Wang et al. 2013) may also result in higher probabilities of having helioeffective CME–CME interaction events at specific heliocentric distances than at others.