Study of Evolution and Geo-effectiveness of Coronal Mass Ejection–Coronal Mass Ejection Interactions Using Magnetohydrodynamic Simulations with SWASTi Framework

Most of the previous studies (N. Lugaz et al. 2017 and references therein) have focused on defining a single stage for the entire CME–CME interaction structure. However, given the relatively small scale of Earth's magnetosphere compared to the CME structure, examining the evolution of these stages at different longitudes is critical. For instance, the structures passing through the 1 au locations marked by the three dots in Figure 2 at 0° and ±10° longitudes are at different stages. The trailing shock does not extend to the sheath of CME1 along the −10° longitude in any of the cases. The radial evolution of this shock is greater along 0° and even more so along +10°, where it interacts with the sheath and merges with the first shock in some instances. The uneven evolution of these stages is primarily caused by the deformation of CME1 due to the ambient SW preceding it, particularly due to the SIR in front of it in these specific situations. This SIR positions the top flank of the CME predominantly within the fast-SW stream, while the bottom flank experiences a significantly stronger antiradial drag force. This dynamic results in an overexpansion of the upper flank and an underexpansion of the lower flank, as depicted in Figure 1.

This four-stage evolution concept can be applied both globally across the entire structure and locally along different longitudes to study the progression of CME–CME interactions comprehensively. To facilitate a robust comparison of different cases with varying initial properties, it is crucial to identify both global and local stages. For this purpose, we utilize the concept of virtual spacecraft to observe the in situ plasma conditions at the front and rear of CME1 along three longitudes (0° and ±10°). Six comoving virtual spacecraft were positioned, with three at the rear and three at the front of CME1. Figure 3 displays the in situ speed profiles collected by these spacecraft for the HSHDLF0 case. In Figure 3(a), the solid lines represent data from the rear of CME1, while the dashed lines are from the front. At the rear of CME1, the trailing shock arrives with a time difference of 1 hr between the ±10° longitudes, resulting in a similar speed jump along these longitudes due to the equivalent shock propagation time. However, at the front of CME1, the arrival time difference of the trailing shock between these longitudes is much larger. This is because the radial width of CME1 varies significantly between +10° and −10° longitudes, leading to an enhanced temporal gap in shock arrival. Furthermore, since the radial extent of CME1 is greater along −10°, the trailing shock travels a longer distance and loses more energy, resulting in a smaller speed jump compared to +10° longitude.

Zoom In Zoom Out Reset image size

For comparison across all ensemble cases, we focused on +10° longitude to study the properties of shock propagation. This approach examines the time it takes for the trailing shock to travel from the rear to the front of CME1, a duration primarily influenced by two factors: the radial width of CME1 and the strength of the trailing shock. Given that the CME1 properties remain consistent across all cases, the strength of the trailing shock emerges as the sole determinant of this duration. The bottom left and right subplots of Figure 3 illustrate the arrival and departure times of the trailing shock from the CME1 structure, revealing two distinct behaviors influenced by the initial density of CME2. The speed jump due to the shock is consistently higher in all HD cases compared to all LD cases. Although the HD cases have similar arrival times, the sharpness of the shock is greatest in the HSHDHF case and least in the LSHDLF case. As the trailing shock reaches the front of CME1, the temporal differences between the cases become more pronounced, with the shock propagation duration for the strongest-shock (HSHDHF0) and weakest-shock (LSLDLF1) cases widening from 13.3 to 28.26 hr.

The trailing shock propagation time (Δ) for all ensemble cases is presented in Table 2. The duration is shortest for the HSHDHF case, in both the with and without relative tilt scenarios, and it progressively increases as the initial flux, speed, and density decrease, with the longest duration observed in the LSLDLF case. Cases with a higher initial density consistently exhibit shorter Δ values. Since this duration is directly correlated with the shock strength, it suggests that higher-density CMEs generate stronger shocks. Elevated density leads to an increase in internal pressure within the CME, which amplifies the pressure differential with the ambient SW. This, in turn, results in an accelerated rate of CME expansion. Moreover, greater momentum allows the CME to plow more effectively through the SW plasma, especially in the LD environment of the preconditioned SW. Therefore, a combination of a faster expansion rate and more efficient plowing contributes to the observed trend. Furthermore, cases with 0° tilt consistently show smaller Δ values compared to those with 45° tilt. This could potentially be due to the preconditioned SW: CME2, having no relative tilt, follows the same trajectory as CME1 and encounters less drag force from the LD SW swept by CME1. However, with nonzero relative tilt, the structure of CME2 does not completely align with the cleared path and faces greater resistance from the ambient SW, potentially leading to a reduction in shock strength.

As the trailing shock progresses inside the first magnetic cloud, the rear CME begins colliding with the first. This collision starts in stage 2, where the trailing shock accelerates the rear of CME1. Depending on the advancement of the shock, the difference in speed between the rear of CME1 and the front of CME2 can exceed or fall below the local fast magnetosonic speed. When this speed is exceeded, the trailing CME pushes the plasma faster than the leading CME can smoothly adjust, leading to the formation of a shock wave directed toward CME2.

Fast reverse shocks associated with SIRs and their corresponding in situ properties have been well studied (E. Kilpua et al. 2015; D. M. Oliveira 2016). However, their formation due to CME–CME interaction has not been reported until very recently, by D. Trotta et al. (2024) using Solar Orbiter observations. Similar in situ features have also been observed through virtual spacecraft in our simulation. The bottom row of Figure 4 (subplots (a2), (b2), and (c2)) demonstrates the existence of fast reverse shocks. The characteristic feature of anticorrelation between plasma speed and magnetic field and density can also be seen. As this reverse fast-mode shock wave propagates antiradially, it compresses the downstream plasma inside CME2, resulting in higher plasma density in the downstream region compared to the upstream region. Conversely, the effective plasma speed will be lower in the downstream region due to the antiradial direction of the shock, leading to an anticorrelation between speed and density profiles.

Zoom In Zoom Out Reset image size

Given the nonuniform interaction demonstrated in earlier sections, the deformation of CME1 causes the collision to initiate in one or more localized areas rather than across the entire interface simultaneously. This nonuniformity leads to the generation of independent shock fronts in each of these areas. The interactions among these shocks can lead to complex dynamics, including the merging of shocks, the amplification or attenuation of wave fronts, and the formation of complex wave patterns. Figure 4 illustrates this complexity, showing the reverse shock's propagation in the rear magnetic cloud for the HSHDHF0 case. Inside CME2, a complex pattern of alternating compressed and rarefied zones is visible, manifesting as ripple-like structures in the in situ plots. Such variations hold the potential to influence the duration of geomagnetic storms.

The subplots (a1), (a2), (d1), and (d2) of Figure 5 illustrate the percentage change in the total radial momentum (TRM) of CME1, which is calculated using following equations:

$$\begin{eqnarray} &&\% \,\mathrm{Change}=100\times \left(\displaystyle \frac{{\mathrm{TRM}}_{i}-{\mathrm{TRM}}_{0}}{{\mathrm{TRM}}_{0}}\right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\mathrm{TRM}=\displaystyle \sum _{j}{\rho }_{j}\,{v}_{{r}_{j}}\,{{dV}}_{j}.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

Here, the subscript "i" represents cases 1–16, while "0" denotes case 0. ρ_j , ${\rm{}}{v}_{{\rm{}}{r}_{j}}$, and dV_j are the mass density, radial velocity, and volume of the grid cell "j," respectively. The TRM is the summation of the radial momentum at each grid cell. Based on the slope of the momentum profile, we interpret the temporal evolution of momentum change in CME1 in two distinct phases: the Rising Phase and the Diminishing Phase.

The Rising Phase commences when the trailing shock from CME2 impacts the rear of CME1. During this phase, the effective total momentum of CME1 starts to increase, due to the local velocity increase of the CME1 plasma. As the shock propagates deeper into CME1, the total momentum continues to rise. This momentum increase of CME1 solely due to the trailing shock occurs very briefly (1–2 hr), followed by the collision of the trailing sheath and magnetic cloud with CME1. This collision initiates a substantial momentum exchange between the two CMEs. The high-velocity trailing CME2 exerts significant force on the slower-moving CME1, resulting in a notable increase in CME1's total momentum. During this phase, the rate of momentum increase (indicated by the slope of the momentum curve) achieves its highest value, reflecting the extent of the direct and robust transfer of momentum from CME2 to CME1.

In the final Diminishing Phase, CME2 has transferred a considerable amount of its momentum to CME1, significantly reducing its own momentum. This momentum exchange decreases the relative speed between the two CMEs. Despite the reduction, CME2 still maintains a higher speed than CME1 and continues to push it forward, albeit with diminishing force. Consequently, the rate of change in momentum for CME1 gradually decreases, indicating that the interaction force is waning and the system is approaching a new equilibrium state.

The transition between phases is quite smooth, making it challenging to define rigid boundaries. Relatively speaking, the Rising Phase has an average duration of approximately 10 hr, after which the Diminishing Phase persists. The exact duration of this transition varies across different cases. For the HSHDHF case, the Diminishing Phase starts earliest, at roughly 7 hr, whereas for the LSLDLF case, it takes the longest time, approximately 13 hr. For cases with relative tilt between CMEs, this transition period is slightly shorter compared to cases with no relative tilt.

The momentum gained by CME1, from the start of the interaction to the phase transition point and until the combined structure reaches 1 au, also varies significantly with the initial conditions of CME2. The eight cases shown in the subplots can be clustered into two groups: HD and LD cases. After 10 hr of interaction between CME1 and CME2, CME1 in HD cases gains approximately 5%–15% more momentum, whereas in LD cases, it gains 1%–5% more momentum. This momentum gain is slightly higher for cases with no relative tilt compared to cases with relative tilt. Moreover, the difference between the HD and LD cases is also more pronounced for the no-relative-tilt cases.

As the CME–CME structure reaches 1 au, starting from 25 hr, the difference in momentum gain between the HD and LD cases becomes more pronounced. Since shock strength is inversely correlated with shock propagation time, which is shorter for HD cases (see Table 2), these cases exhibit higher momentum gain. Additionally, CME2 in HD cases has more initial momentum to impart than in LD cases. After 30 hr of interaction, CME1 in HD cases demonstrates a 20%–40% increase in momentum compared to 7%–12% in LD scenarios. The gap between these two clusters of cases becomes 16% in the absence of relative tilt compared to 8% in the presence of tilt.

Since the total KE (KE = $0.5{\sum }_{j}{\rho }_{j},{v}_{{r}_{j}}^{2},{{dV}}_{j}$) of a radially evolving CME1 is related to its TRM, a similarity between the two is expected. The subplots (b1), (b2), (e1), and (e2) of Figure 5 depict the temporal evolution of the percentage change in the KE of CME1 for cases 1–16 relative to case 0. These subplots closely resemble the momentum plots but indicate a greater gain. On average, the percentage gain in KE is about two-thirds greater than the percentage gain in momentum of CME1 due to the interaction process. As with momentum, the eight cases can be clustered into two groups, based on their initial density. The difference between these two groups is more pronounced in the absence of relative tilt (up to 25%) compared to scenarios with nonzero tilt (up to 18%). After 30 hr of interaction, CME1 gains up to 60% more KE in the HSHDHF case and up to 20% more in the LSLDLF case, due to the interaction with CME2.

The evolution of the total ME of CME1 differs significantly from the changes observed in momentum and KE. The subplots (c1), (c2), (f1), and (f2) of Figure 5 showcase the temporal evolution of the percentage change in ME of CME1. These subplots reveal distinct behaviors based on the initial density of CME2. In LD cases, the interaction has little to no effect on the ME of CME1. In contrast, HD cases exhibit an initial increase in ME of approximately 1%–3%, followed by a subsequent decrease. This pattern is consistent in both the tilt and no-tilt scenarios, although the changes in ME occur more rapidly in the presence of relative tilt.

As the interaction begins, CME2 compresses CME1, potentially increasing the magnetic field strength and, consequently, the total ME if magnetic flux is conserved. Since CME1 and CME2 have the same chirality, this compression could also induce magnetic reconnection and other instabilities (e.g., tearing mode and plasmoid instabilities), dissipating magnetic fields and converting ME into kinetic and thermal energy. Despite significant compression even in the LSLDLF case, the gain in ME for LD scenarios is almost negligible, suggesting the dissipation of the magnetic field may happen even in weaker collisions.

In HD cases, there is a noticeable increase in ME, implying that the rate of magnetic field strength enhancement due to compression initially exceeds the rate of magnetic dissipation. However, after approximately 25 hr of interaction, the ME begins to decrease, indicating that the rate of compression diminishes while the rate of magnetic dissipation becomes dominant, leading to a reduction in the total ME of CME1.

Unlike momentum and KE profiles, ME profiles can be categorized into three groups: LD cases, HD high-flux (HDHF) cases, and HD low-flux (HDLF) cases. Given the insignificant changes in ME for LD cases, their slopes are excluded from subplots (c2) and (f2), to highlight the behavior of the remaining cases. Comparing the HSHDHF (red) and HSHDLF (orange) cases, the initial rise in ME is greater for the HSHDHF case, but the subsequent decrease is also more rapid. After 30 hr of interaction, the effective ME is lower for the HSHDHF case. This trend is consistent when comparing the LSHDHF to LSHDLF cases and across the tilt and no-tilt scenarios. This suggests that a higher initial magnetic flux results in a greater initial gain in ME due to compression, followed by a more substantial decrease in ME due to magnetic dissipation.

In the earlier sections, we demonstrated and discussed the nonuniform radial expansion of CME1 due to interactions with the SW from the front and CME2 from the rear (see Figures 1 and 2). We also noted in the last section that the only possible way of increasing the ME is the compression of CME1 due to CME2. Although multiple studies (M. Xiong et al. 2006; N. Lugaz et al. 2013, 2017) have examined the compression of the leading CME in CME–CME interactions, a quantitative analysis of the temporal evolution of this compression has not been performed.

Using the passive scalar tracing technique, we traced the CME1 structure and computed its radial compression along different longitudes. This method tracks plasma using the advection equation, allowing the study of transport and mixing without altering magnetic or fluid dynamics, as demonstrated in P. Mayank et al. (2023). By placing virtual spacecraft at the front and rear of CME1, we can measure the radial difference between them, thus determining the radial extent of CME1 along each longitude. Figure 6 shows the temporal evolution of the percentage change in CME1's radial width across all CME–CME interaction cases (1–16) versus the single-CME case, with maximum compression seen in the HSHDHF case and minimum in the LSLDLF case across −10°, 0°, and +10°.

Similar to the momentum, KE, and ME profiles, the compression features of the HD and LD cases are distinguishable in the subplots of Figure 6. Compression is consistently higher for HD cases from the onset, though no noticeable gap between the LSHDLF and HSLDHF cases appears until about 10 hr for the 45° tilt cases. At +10° (Figure 6(a1)), HD cases with 45° tilt (HD1) reach an 80% reduction over 30 hr, while LD1 shows 60% reduction. Along −10° (Figure 6(a2)), HD1 shows a 40% reduction and LD1 a 20% reduction. The absence of an SIR along this longitude means there is no significant obstruction to CME1's expansion from the front, resulting in less compression from the rear than at +10°, with a nearly constant decrease for all cases. At 0° (Figure 6(a3)), the HD1 cases show a 65% reduction compared to 45% for LD1, revealing that the presence of the SIR ahead of CME1 has an intermediate impact on this region, more than −10° but less than +10°.

Zoom In Zoom Out Reset image size

The abovementioned statistics pertain to cases with a 45° tilt between the CMEs. In cases with a 0° tilt (see Figure 6, panels (b1)–(b3)), the interaction begins slightly earlier, leading to faster and greater compression. On average, the compression in cases without tilt is approximately 5% greater than in cases with tilt. This suggests that the alignment of CME2 with CME1 can influence the efficiency of the compressive interaction. Due to this tilt-induced difference, an enhancement in the gap between the HD and LD cases is evident in Figure 6. When comparing the overall effect of CME2's initial properties, the HSHDHF case consistently shows the greatest compression, followed by the HSHDLF and LSHDHF cases.

To estimate the extent of the plasma mixing between the leading and trailing CMEs, we utilize the CME tracer described in earlier sections to define the boundary of the CME structure. This approach is analogous to methods used in studying the mixing of astrophysical jets (e.g., see S. Walg et al. 2013). To quantify this mixing, we define a mixing factor (${ \mathcal M }$), which represents the absolute mass fractions within a specific grid cell. We set ${ \mathcal M }=0$ for the case of no mixing, where only CME1 material is present, and ${ \mathcal M }=2$ for the case where only CME2 material is present. We set ${ \mathcal M }=1$ in the case of maximum absolute mixing, meaning equal amounts of CME1 and CME2 constituents are present within the grid cell. In this scenario, the mass fraction of CME1 is equal to the mass fraction of CME2.

At a given time (t) and distance ( r ), a tracer ${ \mathcal T }(t,{\boldsymbol{r}})$ is advected by the flow and obtains values within the range ${{ \mathcal T }}_{\min }\leqslant { \mathcal T }(t,{\boldsymbol{r}})\leqslant {{ \mathcal T }}_{\max }$. Here, ${{ \mathcal T }}_{\min }$ corresponds to the absence of that CME within the cell, while ${{ \mathcal T }}_{\max }$ corresponds to a cell purely containing the plasma of that CME. In this work, we have taken ${{ \mathcal T }}_{\min }$= 0.1 and ${{ \mathcal T }}_{\max }$= 1.0 for all CMEs. Given that a tracer value directly corresponds to the quantity of that CME, the mass fraction of CME1 (δ₁) within a grid cell can be expressed as follows:

$$\begin{eqnarray}&&{\delta }_{1}(t,r)=\left|\displaystyle \frac{{{ \mathcal T }}_{1}(t,r)-{{ \mathcal T }}_{\min }}{{{ \mathcal T }}_{\max }-{{ \mathcal T }}_{\min }}\right|.\end{eqnarray} \tag{ 4 }$$

Similarly, the mass fraction of CME2 (δ₂) can also be calculated. Furthermore, assuming that the mass fraction within a cell linearly scales with the amount of mixing, the mixing factor of the plasma of CME1 and CME2 in terms of their tracer values in a grid cell can be written as:

$$\begin{eqnarray}&&{ \mathcal M }=1-\left(\displaystyle \frac{{{ \mathcal T }}_{1}-{{ \mathcal T }}_{2}}{{{ \mathcal T }}_{1}+{{ \mathcal T }}_{2}}\right).\end{eqnarray} \tag{ 5 }$$

Figure 7 presents a 2D histogram of the calculated mixing factor values throughout the entire simulation domain for all ensemble cases at the time when the merged structure reaches 1 au. These values are computed for each grid cell within the computational domain and binned into eight uniform clusters, ranging from 0 to 2. A mixing factor of ${ \mathcal M }=1$ indicates the highest level of mixing (50%), meaning both CME1 and CME2 contribute equal amounts of plasma to that specific cell. The cases with and without relative tilt between the leading and trailing CMEs are plotted separately. The color of each block represents the number of grid cells associated with a specific amount of mixing (Y-axis) for each corresponding case (X-axis). For instance, in case 15, 261 grid cells have ${ \mathcal M }$ values between 0.75 and 1 (see Figure 7(b)).

Zoom In Zoom Out Reset image size

The most notable feature in these subplots corresponds to the HD cases (5–8 and 13–16), which exhibit a significantly larger volume of mixing. As discussed earlier, these scenarios also demonstrate stronger CME shocks, higher momentum exchange, and greater radial compression. This trend is consistent in terms of the extent of mixing, where ${ \mathcal M }$ values are several times higher in HD cases as compared to LD cases. The HS cases (9–16) also consistently exhibit greater mixing compared to LS cases. Although the difference is not as big as in HD versus LD cases, the trend is evident across all bins of mixing factor.

Apart from these global trends, we also observe nonuniform mixing across the CME2 front. Similar to the nonuniform interaction between the back of CME1 and the front of CME2 discussed in earlier sections, the bottom flank exhibited a higher percentage of mixing (100 × [1 − $| { \mathcal M }-1|$]) compared to the upper flank. Panels (c) and (d) of Figure 7 show the mean percentage of mixing along +10° (bottom flank), −10° (upper flank), and 0° longitudes for the 45° and 0° cases, respectively. In these plots, the +10° region consistently exhibits nonzero mixing, with the highest percentage in almost all the cases, suggesting that the bottom flank of CME2 experiences more intense interaction with CME1. In contrast, along −10°, only three cases with a 45° tilt and none with a 0° tilt show any mixing, while the 0° region exhibits mixing in a total of eight cases. This pattern supports the observed asymmetry in interaction strength, where the bottom flank experiences stronger compression, enhancing the mixing process. This again highlights the impact of the inhomogeneous ambient SW, which causes nonuniform interaction and mixing of the interacting CMEs.

The onset of the Dst main phase begins with the arrival of the first shock at 1 au. Figure 8 illustrates the complex variations among the 16 cases, highlighting the different starting and ending times of the main phase. The onset time variation primarily depends on the ambient SW conditions, which remain constant across all cases, resulting in similar onset times for different cases at the same longitude. However, in scenarios where shock–shock interactions occur, deviations from the initial onset times are observed. Since the shock first impacts the spacecraft at +10° longitude and no shock–shock mergers occur at +10° in any of the cases, the main-phase starting time is identical for the green profile across all cases. This timing aligns with the arrival of the first shock (∼20:00 hr on 2023 May 13).

Similar to +10° longitude, the two shocks arrive sequentially at 0° longitude at 1 au without merging in any of the cases. Both locations exhibit a two-dip profile, where the first dip is associated with the arrival of the first shock and the second dip begins with the arrival of the second shock. In the orange profile, the first Dst drop is very small (around 20:00 UT on May 13) and recovers almost fully to the pre-storm level within a few hours, just before the second shock arrives at about 12:00 UT on May 14. Depending on the strength of the second shock, the Dst value then drops again, reaching a minimum either more rapidly or gradually. A similar trend is observed for the blue profile, with the main difference being that it takes longer for Dst to recover to pre-storm levels before the arrival of the second shock. This variation is primarily due to the B_z profiles at these two locations, as shown in Figure 13. In the orange profile, the negative B_z quickly (∼4 hr) turns positive, whereas at 0° longitude, it remains negative for almost 7 hr.

Unlike the orange and blue profiles, the Dst main-phase profiles at +10° show only a single dip in all cases. In LD cases, this dip features two distinct slopes: an initial gradual decline followed by a sharper fall. For HD cases, there is a single sharp decline to the global minimum, representing the shock–shock merger and indicating that the CME–CME interaction at this longitude is in its fourth stage. Despite the sharpest decline occurring at this location, the minimum Dst values do not always correspond to this longitude. Interestingly, this is especially not true for HF cases without tilt, where the initial magnetic flux of CME2 is greater. In these cases, the minimum Dst value is seen at 0° (blue line) rather than at +10° (green line). Conversely, in LF cases without tilt, the minimum Dst value appears in the green (+10°) profile. This trend is intriguing, because the interaction is most prominent along the +10° longitude. Yet, an increase in the initial magnetic flux of CME2 leads to higher Dst values along this longitude.

There are also some peculiar trends in the main phase corresponding to the initial properties of CME2. The subplots of Figure 8(b) demonstrate clear differences between the HD and LD cases, showing a trend similar to that observed in Figure 5. The divergence between these cases begins around 15:00 hr on May 14, with the arrival of the trailing shock, imitating the differences shown in the strength of the shocks in Figure 3(b). Notably, by the time Dst reaches its global minimum, the temporal gap between the HD and LD cases further widens, with HD cases reaching their minimum Dst value before the single-CME case.

Additionally, cases with and without relative tilt exhibit significant differences. In most cases, both the cumulative Dst and minimum Dst values are higher when there is a relative tilt between the CMEs compared to when there is none. The primary reason for this is the change in orientation of the magnetic field of CME2, which should ideally enhance the effect along −10° and decrease it along +10°, while remaining the same at 0°. However, due to the presence of SW and prolonged interaction with CME1, this trend deviates slightly from the ideal scenario. Figure 13 shows the B_z values for all cases, where the minimum B_z values exhibit a mixed trend along −10°, 0°, and +10° longitudes with changes in tilt. However, the cumulative duration of negative B_z is consistently longer for +10°. Additionally, another noticeable difference is in the slope of the Dst fall, which is steeper for HD cases compared to LD cases.

In an ideal single-CME storm, the recovery phase is smooth and continuous, with the Dst index value generally increasing following an exponential trend and gradually returning to the pre-storm level (as in Equation (A3)). Any deviations or additional fluctuations from this ideal trend can indicate the influence of a second CME and subsequent SW. Figure 8(b) showcases that among the 16 ensemble cases, some exhibit significant deviations, others show minor variations, and a few have changes that are almost negligible.

The HD cases exhibit significant deviations from the idealized scenario of the Dst recovery phase. Notably, scenarios without relative tilt between CME1 and CME2 emphasize these deviations more clearly (see Figure 8). In cases 5, 7, 13, and 15, a discontinuity is observed in the Dst index increase across all three profiles (orange, blue, and green), resulting in the formation of a recovery phase plateau. These interruptions, while not affecting the minimum Dst value, significantly prolong the overall recovery time, delaying the return to pre-storm Dst levels. The mentioned HD cases demonstrate notably longer recovery phases compared to their LD counterparts. For instance, by May 16, the Dst index had almost recovered to 0 nT for cases 3 and 11 (LDHF0), whereas for cases 7 and 15 (HDHF0), the Dst index remained near −50 nT.

The recovery phase plateau can be attributed to fluctuations in the southward magnetic field, as southward interplanetary magnetic field (IMF) excursions drive the Dst down before recovery can resume. In the B_z profile (see Figure 13), such fluctuations are noticeable after approximately 18:00 hr on May 14 for HD cases without tilt. These fluctuations in the B_z profile are mirrored by similar, more pronounced fluctuations in the in situ speed profile. As discussed in Section 4.1.3, these fluctuations result from the reverse shock formed by the collision of CME2 with CME1. The propagation of this reverse shock through CME2 creates complex patterns of alternating compressed and rarefied regions, leading to the ripples observed in the in situ profiles. Consequently, these ripples extend the Dst recovery phase by forming a plateau, significantly delaying the return to pre-storm levels.

Apart from the major deviations from the idealized scenarios, many cases showcase minor or no deviations. For example, the blue profile in cases with tilt shows slight deviations; in cases 6 and 14, there is a break point, and in other cases, the profile does not increase smoothly but rather with some disruptions. Notably, there is no significant flipping in the B_z profile after reaching the global minimum in the Dst profile around May 15 for these cases. Interestingly, although there are a few break points around May 15 in the speed profile, the absence of flipping in the B_z profile prevents the formation of any plateau. In other cases, such as 3 and 11, the blue profile closely follows the idealized scenario. Similarly, the orange profile exhibits an ideal recovery phase in almost all tilt cases. In these scenarios, there is neither directional flipping in the B_z profile nor any break points in the speed profile, resulting in a smooth, continuous function.

Figure 9 illustrates the minimum Dst values for LD and HD cases across three different longitudes: −10°, 0°, and +10°. The percentage differences (nearest integer of $100\cdot \left[\tfrac{\mathrm{HD}-\mathrm{LD}}{\mathrm{LD}}\right]$) between paired cases highlight some significant trends. At 0° (panel (a2)) and +10° (panel (a3)), HD cases consistently show lower minimum Dst values compared to LD cases, with percentage differences ranging from 42% to 70% at 0° and 18% to 129% at +10°. This indicates that higher initial densities tend to exacerbate the severity of geomagnetic storms at these longitudes. This trend seems intuitive since a higher initial density in CME2 leads to a stronger trailing shock and higher in situ speed and density, which in turn would lead to a lower Dst index. However, at −10° (panel (a1)), the trend is reversed; HD cases exhibit higher minimum Dst values compared to LD cases, with percentage differences ranging from −8% to −23%. Although there is a significant increase in in situ speed at −10°, the in situ density does not show a major increment, possibly because the trailing shock has to cover a much larger radial distance in an overexpanded region of CME1. Additionally, when this second shock arrived, the B_z profile was not in the southward direction as it was along the 0° and +10° longitudes. Due to this combined in situ configuration, the effective Dst minimum became higher.

Zoom In Zoom Out Reset image size

At −10° (panel (b1)), the differences between the LS and HS cases are minimal, with percentage changes ranging from −11% to 2%. This suggests that the initial speed has a relatively minor impact on the severity of geomagnetic storms along this overexpanded CME1 region. In contrast, at 0° (panel (b2)), HS cases consistently exhibit lower minimum Dst values compared to LS cases, with percentage differences ranging from 7% to 15% and an average change of 11%. This potentially indicates that higher initial speeds of CME2 results in a more geo-effective configuration at this longitude. At +10° (panel (b3)), the trend is less consistent (in six out of eight cases), with percentage differences varying from −2% to 18% and an average change of 5%. Despite this variability, HS cases generally lead to lower minimum Dst values (in 18 out of 24 cases), suggesting a trend toward more severe storms with increased speeds. These findings demonstrate that while the influence of the initial speed is significant at 0° and somewhat at +10°, it is relatively minor at −10°.

Like the initial density and speed, the effect of the CME2 tilt on the minimum Dst values also shows distinct patterns across the three longitudes. With an average percentage difference of 3%, ranging from −2% to 6%, the Dst minimum was lower for six out of eight tilt cases along −10° longitude (panel (c1)). At 0° (panel (c2)), the tilt's effect is more varied, with percentage differences ranging from –5% to 28% and an average change of 7%, showing lower Dst minimum values for tilted cases in five out of eight instances. The most pronounced effect is observed at +10° (panel (c3)), where the differences are consistent and substantial, ranging from 32% to 128%, with an average decrease of 81%. This significant impact at +10° highlights how the tilt in CME2, which alters the magnetic field orientation, can greatly intensify the geomagnetic storm along one direction. Ideally, the Dst minimum should have consistently increased in the opposite direction (+10°), but the nonuniform deformation of CME1 due to inhomogeneous ambient SW causes the Dst values to deviate from this expected trend.

The increase in the initial magnetic flux of a CME results in an increase in its magnetic field strength. In an ideal isolated flux-rope scenario, this enhancement would typically lead to a lower Dst minimum upon the CME's interaction with Earth's magnetosphere. However, for CME–CME interactions in the presence of a realistic nonuniform ambient SW, the outcome can differ significantly. This is evident in Figure 9. At −10° (panel (d1)), the differences between the low- and high-magnetic-flux cases are relatively minor, with percentage changes ranging from 2% to 12% and an average difference of 5%. This consistent trend of lower Dst minimum for higher initial magnetic flux suggests that magnetic flux does impact the severity of geomagnetic storms, albeit slightly at this location. At 0° (panel (d2)), the impact of the magnetic flux is more pronounced, with percentage differences ranging from −2% to 36% and an average change of 13%. The trend is consistent, with one exception that has the lowest percentage change. Combining the results from −10° and 0°, the trend is clear: higher magnetic flux leads to lower Dst minimum values in 15 out of 16 cases. However, at +10° (panel (d3)), the trend is reversed in seven out of eight cases; high-magnetic-flux cases exhibit higher minimum Dst values compared to low-magnetic-flux cases, with percentage differences ranging from −46% to 2% and an average increase of −16%. It is important to emphasize that this trend reversal occurs along the longitude where the interaction has been strongest. This potentially indicates that higher magnetic flux in CME2 may lead to greater magnetic flux dissipation in strong CME–CME interaction regions, reducing the geo-effectiveness of CME–CME interactions at this longitude.

The histograms in Figure 10 illustrate the significant impact of the initial density on the cumulative Dst values. Panel (a1) shows the cumulative Dst index values for −10° longitude. The differences between the LD and HD cases are substantial, with percentage changes ranging from −79% to 107% and an average change of −29%. This trend of higher initial densities resulting in lower cumulative Dst values is observed in six out of eight cases, indicating a general reduction in geomagnetic impact at this longitude, consistent with the trend observed in the minimum Dst. At 0° (panel (a2)), HD cases exhibit consistently higher cumulative Dst values compared to LD cases, with percentage differences ranging from 33% to 246% and an average increase of 109%. At +10° (panel (a3)), the differences are even more pronounced, with percentage changes ranging from 28% to 2327% and an average increase of 731%. This increasing trend is observed in all eight cases at 0° and +10°, underscoring the substantial enhancement in storm duration due to the higher initial densities of CME2. These findings demonstrate that the initial density of CME2 plays a crucial role in determining the cumulative geomagnetic impact, with the most significant effects observed at 0° and +10°, and a complex behavior at −10°.

Zoom In Zoom Out Reset image size

At −10° (panel (b1)), the cumulative Dst values show a wide range of percentage changes, from −46% to 47%, with an average change of −16%. This variability highlights that higher initial speeds tend to reduce the cumulative geomagnetic impact at this location, though the trend is opposite to what was observed in the minimum Dst, and the effect is not uniform. In contrast, at 0° (panel (b2)), the cumulative Dst values for HS cases exceed those of LS cases, with percentage changes spanning from −7% to 51% and an average increase of 16%. Here, seven out of eight cases follow this pattern, consistent with the trend observed in the minimum Dst. At +10° (panel (b3)), the differences are more subtle, with percentage changes ranging from −9% to 23% and an average change of 9%. This pattern appears in six out of eight cases, suggesting a mild tendency for higher speeds to amplify the cumulative Dst values. These observations indicate that the impact of CME2's initial speed on the storm duration and intensity varies with longitude, showing enhancement at 0° and +10° but reduction at −10°.

The impact of tilt between CMEs shows a relatively more consistent trend. For −10°, the cumulative Dst values decrease, they show mixed behavior at 0°, and increase at +10°. At −10° (panel (c1)), the high-tilt cases have percentage changes ranging from −86% to 6%. This trend is observed in seven out of eight cases, with an average change of −36%, opposed to the general observation in minimum Dst. Moving to 0° (panel (c2)), the effect of tilt becomes inconsistent. Here, percentage changes vary from −16% to 170%, with an average change of 34%. The trend is followed in five out of eight cases, indicating a mixed influence of tilt on cumulative geomagnetic impact, and it is not as pronounced as at −10°. In contrast, at +10° (panel (c3)), the high-tilt cases dramatically increase the cumulative Dst values, with percentage changes ranging from 32% to 2402% and an average increase of 697%. This overwhelming trend is observed in all eight cases, highlighting a substantial enhancement in geomagnetic storm duration and intensity due to tilt at this longitude. This demonstrates that CME2's tilt has a variable influence on the cumulative geomagnetic impact, with a significant reduction at −10°, moderate variability at 0°, and immense enhancement at +10°.

Among the three longitudes, only the cumulative Dst values at 0° show a consistent trend, while the ±10° longitudes exhibit mixed behavior with the increasing magnetic flux of CME2. This inconsistency contrasts with the statistical trends observed in the minimum Dst. At −10° (panel (d1)), the percentage changes range from −46% to 362%. Decreases are observed in five out of eight cases, but the average change is +56%, reflecting minor decrements and significant increments. Similarly, at +10° (panel (d3)), the results are highly variable, with percentage changes ranging from −41% to 184% and an average increase of 22%. Decreases occur in five cases, while increases are seen in three out of eight cases. In contrast, at 0° (panel (d2)), the impact of higher magnetic flux is consistent and more pronounced, with percentage changes spanning from 15% to 156% and an average increase of 67%. This trend is observed in all eight cases, indicating a strong correlation between increased magnetic flux and enhanced geomagnetic storm duration and intensity.