STEREOSCOPIC STUDY OF THE KINEMATIC EVOLUTION OF A CORONAL MASS EJECTION AND ITS DRIVEN SHOCK FROM THE SUN TO THE EARTH AND THE PREDICTION OF THEIR ARRIVAL TIMES

On board the STEREO spacecraft, SECCHI contains five separate instruments: the Extreme Ultraviolet Imager measuring EUV radiation at four wavelengths; two coronagraphs, COR1 and COR2, which observe white light ranging from about 1.4–4 R_☉ and about 2.5–15 R_☉, respectively; and two Heliospheric Imagers, HI-1 and HI-2, which have ranges of about 15–84 R_☉ and 66–318 R_☉, respectively. The Heliospheric Imagers, like the coronagraphs, observe Thomson scattered white light. By combining observations from each instrument, a CME can be tracked continuously from its eruption in the low corona into the heliosphere (Howard et al. 2008).

The primary data sets used for this work are the COR2, HI-1, and HI-2 instruments on board the two STEREO spacecraft, as there is an SEP related data gap in SOHO LASCO data for much of the propagation. In situ data was obtained from the ACE, a spacecraft at the L1 point that has instruments for measuring both the magnetic field (Magnetometer) and the properties of the plasma (Solar Wind Electron-Proton and Alpha Monitor) of the solar wind.

Different methods and models for extracting CME height data from STEREO have been developed (Tappin & Howard 2009; Lugaz et al. 2009; Thernisien et al. 2009; Liu et al. 2010; Mishra & Srivastava 2013). For the purpose of studying the evolution of the shock and CME ejecta fronts in this paper, two different methods were used.

The first is a forward modeling technique that uses the graduated cylindrical shell (GCS), or croissant geometric model, utilizing the raytrace code developed by Thernisien et al. (2006, 2009). A demonstration of the model is shown in Figure 1. This model involves taking a cylindrical croissant defined by six free parameters and superimposing it upon images from multiple spacecraft at once; the parameters are chosen such that the model image best matches the observed images from multiple view points. These parameters are the direction of propagation (longitude and latitude), the leading edge distance of the tracked feature from the Sun, the half angle width of the shell, the tilt angle of the flux rope central axis orientation, and the ratio between the major and minor radius of the flux rope (aspect ratio). Making the assumption that the other five parameters are roughly constant as the CME propagates, the leading edge height can be used to determine a height profile of the CME with time. Performing a numerical derivative of this data can then yield a velocity.

The geometric structures utilized in the raytrace code are idealized structures, and for some events it can be difficult to match the model to the data. For the July 12 event, there is a bright feature that is not being fit as part of the ejecta, as seen in Figure 1. As can be seen more clearly in Figure 2, this feature is not part of the cavity structure of the CME ejecta, but instead is likely a streamer disturbed by the CME shock sheath region.

Zoom In Zoom Out Reset image size

In addition to the croissant-shaped flux rope, there is also a spheroid shock model in the code, with six free parameters as well. Some of these, also in the flux rope model, are propagation direction (longitude and latitude), tilt, and height. The other two parameters, ε and κ, control the major and minor axis of the spheroid.

As this work is focusing mostly on the shock and ejecta leading edges, less attention was paid to ε and κ, which are more important to the detailed three-dimensional shock geometry. However, both features can be tracked in such a way that the heights along the leading edge can be tracked, leading to an understanding of the independent evolution of each front, as well as the size of the sheath region between them.

In order to track the different features into the heliosphere, it is important to process image data in such a way that both features are distinguishable. In white light observations, the CME flux rope is best seen in direct images, which are made through a base-ratio technique where an average background is used as the denominator, while a particular image is numerator in the ratio; a direct image allows the bright features of the flux rope and the internal structure to be seen. On the other hand, to make the shock front and sheath more clear, a running-difference image subtracting the previous image from the current one is more useful, as the outermost features corresponding to the shock are often faint and diffuse. An example of both types of images and the different features is seen in Figure 2.

While success has been shown using the raytrace model to track CMEs (Thernisien et al. 2009; Poomvises et al. 2010; Nieves-Chinchilla et al. 2012), another method was sought as a complement to the raytrace measurements. The height measurements can also be made using so-called jmaps (Sheeley et al. 1999), which are stack plots of image slices along a given position angle from one spacecraft that provides plots in terms of elongation angle versus time. However, due to various observational effects (Howard & Tappin 2009; Lugaz et al. 2009), getting the true de-projected heights in distance from the Sun from jmaps is not as easy as just converting the elongation angles to a distance. There are different methods for deriving true distances. The simplest of these, so-called fixed-φ is given by

$$\begin{equation} R=D_{{\rm obs}} \frac{{\rm sin}(\varepsilon (t))}{{\rm sin}(\varepsilon (t)+\varphi)}, \end{equation} \tag{ 1 }$$

where ε is the elongation angle of the object being tracked, as determined from the jmap, and φ is the direction of propagation of the CME relative to the observer, and D_obs is the distance from the observer to the Sun.

Fixed φ has been shown to be useful near the Sun, but as the CME propagates out it begins to greatly overestimate the leading edge position (Rollett et al. 2012). For this reason, the harmonic mean method was also used. Of all methods utilizing single spacecraft observations, Mishra & Srivastava (2013) found the harmonic mean method to be the most successful. The equation for this method is

$$\begin{equation} R=D_{{\rm obs}} \frac{2 {\rm sin}(\varepsilon (t))}{1+{\rm sin}(\varepsilon (t)+\varphi)}. \end{equation} \tag{ 2 }$$

This equation represents the harmonic mean between fixed ϕ and the Point-P method, given by

$$\begin{equation} R=D_{{\rm obs}} {\rm sin}(\varepsilon (t)), \end{equation} \tag{ 3 }$$

which is the simplest, and most inaccurate, method of deriving true heights from single-spacecraft observations. While it is inaccurate, as Point-P tends to underestimate the true heights, it helps provide a correction for the geometric effects of fixed ϕ (Lugaz et al. 2009).

While most work in the past has attempted to fit the data to provide the best estimation of the direction of propagation, this study is different in that the propagation direction from raytrace was used, allowing for a straight-forward calculation of leading edge distances. Using the same direction of propagation as raytrace could introduce bias in the jmap data, however this should not be a significant issue because the direction of propagation has a relatively small error in the stereoscopic-based raytrace method. The jmaps used along the leading edge position angle are presented in Figure 3, along with the measurements used for the shock.

Zoom In Zoom Out Reset image size

The jmaps from STEREO A and B shown in Figure 3 both show a bright feature beginning in the COR2 field of fiew (FOV), which remains into the HI-1 FOV. The leading front of these two was taken to be the shock front, while the brighter front trails behind. Judging from the appearance of the front in the HI-1 FOV, these fronts are part of the same overall structure, which continually becomes more diffuse. For this reason, the outermost front is the one that was tracked.

These methods solved the problem of how to convert measurements from jmap elongation angles into heights, but there was still the issue of how to interpret the jmap data for the purpose of studying the separate fronts of the CME. Rouillard (2011) relates bright features in jmap data to density enhancements in situ, which makes practical physical sense because the higher the density of a region, the more light will be Thomson scattered and the brighter these areas will be in white light images. Therefore, following this logic, the bright features in the white light images will be taken to be the sheath region, which will have the highest density. The shock front will be the leading edge of the bright feature in the jmap.

In the coronagraph FOV near the Sun, the interpretation of the bright features is different. Due to the observational fact that the CME ejecta is brighter than the shock near the Sun, the assumption in the previous paragraph may be invalid. The HI-1 data is where the two fronts begin to separate themselves more distinctly, so this is more important, and in the majority of the HI-1 data the sheath gains more optical significance.

While this general assumption of the fronts should be valid, there is one more significant issue that should be considered with this interpretation of jmap data. Jmaps are constructed using running difference images. In the case of the shock front, this should not matter, as the shock front is generally traveling into an empty region of space and can be well seen in a running difference image. However, if the sheath region is larger than the distance traveled by the shock front between two images, the shock region that gets subtracted out from the previous image will remove the back of the sheath and leave a dark void. In a jmap it is possible that this void would be interpreted as the ejecta, when it is actually an observational artifact. Therefore with the jmap data, the focus will remain on the shock front.

In addition, errors are also introduced in jmaps because generating a clear and coherent image of the different images stacked together requires pre-processing and smoothing of the images. While this creates an image that is easier to interpret, it is possible that data will be smoothed over. Between these issues and the fact that a jmap only takes data from one vantage point, as opposed to the three views in raytrace, the raytrace is considered the more accurate measurement (but limited to within HI-1 FOV), although all data sets are useful.

With both the GCS model and jmap interpretation to get height time data, a method is needed to obtain velocity profiles for each front. As shown in Byrne et al. (2013), simply taking a numerical derivative of the measurements is not a reliable way to get an accurate result for the CME velocity. Typically, to get data that makes sense using a numerical derivative such as this, smoothing and other manipulations must be performed, introducing more errors. Error is also caused by the non-uniform time cadences of the different instruments and fields of view. Given all this information, it is typically more reliable to choose a certain functional form to fit the data in the height time space and analytically derive this function to obtain a velocity and acceleration.

For this work, the model used to fit the height time data is the drag model (Vršnak et al. 2013), which is a semi-empirical model that assumes that after the initial CME acceleration, magnetic forces are no longer a significant factor in governing CME evolution and aerodynamic drag is the primary force acting on the CME. The equations of motion in the drag model begin with an assumption of a quadratic acceleration and can be analytically integrated to produce velocity and height as a function of time. The equations of motion are

$$\begin{equation} \ a(t)=-\gamma (v(t)-v_{{\rm sw}})|v(t)-v_{{\rm sw}}|,\qquad\quad \end{equation} \tag{ 4 }$$

$$\begin{equation} v(t)=\frac{v_0-v_{{\rm sw}}}{1+\gamma (v_0-v_{{\rm sw}})t}+v_{{\rm sw}},\qquad\qquad \end{equation} \tag{ 5 }$$

$$\begin{equation} r(t)=\frac{1}{\gamma }{\rm ln}[1+\gamma (v_0-v_{{\rm sw}})t]+v_{{\rm sw}}t+r_0, \end{equation} \tag{ 6 }$$

where v₀ is the initial CME speed, v_sw is the upstream solar wind speed, r₀ is the initial CME height, and γ is the drag parameter that controls the rate at which the CME reaches the solar wind speed. For our work with these equations, we are assuming that both γ and v_sw are fixed throughout the CME propagation. Both of these parameters are actually functions of radial distance, especially lower in the corona as the solar wind gets accelerated and the CME begins to expand (Vršnak et al. 2004). This also assumes the entire CME is encountering one static solar wind regime, which may be very incorrect for some events. Approximating γ and v_sw as constants allows us to use the equations above, and fit them as a static parameter.

For theses fittings and with our assumptions, most of the parameters can be determined with just a little information. The initial CME speed and height can be well constrained from observations. v_sw is obtained by taking the average upstream speed of the solar wind eight hours ahead of the arrival of the CME at ACE. For events that experience drastically different solar wind regimes, the assumption that this solar wind speed is representative of the entire upstream regime may introduce significant error, but for most events it should work. This leaves the only totally unknown parameter in the model as γ, the drag parameter. Assuming this is a fixed value throughout propagation, it can be easily fit using the height measurements to get the height time profile.

Applying these methods to the 2012 July 12 CME, both the shock front and ejecta front were tracked. Extreme ultraviolet images and GOES X-ray data show an eruptive X1.4 flare originating from AR 11520 (S17° W08°) with an X-ray peak around 16:45 UT on the 12th. This same active region produced a number of strong flares and CMEs, including the July 23 backside CME, which was considered the fastest CME observed in the STEREO era (Liu et al. 2014). The CME was well observed by COR2 and HI-1 in both STEREO A and B and can also be seen in HI-2, as shown in Figure 4. However the associated SEPs limit the SOHO LASCO data to just a few frames between the eruption and the CME passing beyond the C3 FOV. Both heliospheric imagers observe the event, allowing it to be tracked out to about 80 R_☉ before the different structures become too diffuse and faint to be clearly identified.

The ICME was observed in situ by ACE as seen in Figure 5. The shock is seen as a spike in the magnetic field, velocity, and temperature at approximately 17:15 UT on July 14. Curiously, the density does not show a strong jump indicative of a shock. The sheath lasts for almost 14 hr until around 07:00 UT on the 15th, when the magnetic cloud structure of the flux rope arrives. As is often the case, the boundaries of the magnetic cloud are more ambiguous than those of the more clearly defined shock. In this case, the main signature used to define the flux rope as seen in Figure 5 was the lowered proton temperature, compared to expected temperature (Lopez & Freeman 1986; Richardson & Cane 1995) and the drop in plasma-β. There is certainly room for the interpretation of the arrival of the flux rope, but changing the exact arrival slightly should not have a significant impact on the results.

Zoom In Zoom Out Reset image size

The flux rope initially shows a strong southward magnetic field with a maximum strength of 16.7 nT. The total field has a peak of 28.4 nT, which gradually declines as the flux rope passes the spacecraft. The x, y, and z directions of the magnetic field all show rotational signatures representative of the helically twisted field lines of a flux rope. The velocity shows gradual decline, a sign of an expanding flux rope and as expected, low density and temperature are also shown in the magnetic cloud. The ICME causes a strong geomagnetic storm with a Dst peak of −127 nT on the 15th that lasts for many days before Dst recovers.

To keep the raytrace data consistent, only one of the six parameters that was altered from time step to time step was the height, the other five were held constant for each measurement. While this certainly introduces errors, as the CME will change in time due to internal forces and interaction with the ambient solar wind, it is assumed that these changes are not significant over the short timeframe (less than one day) for which the raytrace measurements were carried out. For this event, the five constant parameters are a latitude of −89, a heliographic longitude of 03, tilt angle of 581, half angle width of 254, and aspect ratio of 0.37. The croissant shape used is comparable to the same model presented in Möstl et al. (2014), which also fit the initial propagation longitude to be within 1 of the Sun–Earth line.

The propagation direction was chosen by first locating the source region of the event in EUV coronal images. All the parameters were then varied until they provided the best match to the images in the coronagraph/heliospheric imager data for as long as the CME is clearly visible. Since the entire CME structure is not visible in the HI-1 FOV, most of the parameters must be determined while the CME is in the COR2 FOV, and by the time the CME passes beyond this FOV, only the height can be determined.

The propagation direction for this event indicates that is a nearly perfect event for comparing remote-sensing observations with in situ data, as both longitude and latitude are within 10° of the Sun–Earth line. The jmaps used for measurement were taken along the same position angle as the leading edge of the CME, so a direct comparison with the height profile from the raytrace method could be made. For predictive purposes, the CME evolution along the equatorial plane is more relevant than that along the CME leading edge, but for this particular event we found that the propagation direction of the CME was close enough to the Sun–Earth line that this distinction was not significant.

The data in Figure 6 show the raytrace measurements as well as the in situ arrival times with the drag model fittings plotted with them. The associated velocity and acceleration profiles are shown in Figure 7. The parameters of the fitting are shown in full in Columns 1–4 of Table 1. The top two rows show the flux rope and shock parameters, respectively, for the full fitting with the in situ data point included. Rows 3 and 4 show the parameters for a fitting without the inclusion of any in situ data that would be more predictively useful. Also presented in Table 1 is the average in situ velocity of each feature (Column 5), the velocity of the drag model at the point it reaches 1 AU (Column 6), and the difference in time between the actual and modeled in situ arrivals (Column 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The drag model, using a final asymptotic speed for both fronts of 353.7 km s⁻¹, is still able to reproduce the speed of the shock and the flux rope as observed in situ quite well, even though at 1 AU both fronts were still in the process of decelerating. This is a sign that the drag model is able to accurately capture both the flux rope and shock propagation using height-time measurements near the Sun and one in situ data at L1. More events will now have to be studied to determine the physical significance of the drag parameter, and how it may be predicted for the creation of an empirical space weather forecasting tool, with a more physical basis than other purely empirical models (Gopalswamy et al. 2005; Kim et al. 2010; Wood et al. 2012).

Our drag parameter, γ, for each front, 1.47 × 10⁻⁸ and 3.09 × 10⁻⁸ km⁻¹, is similar to those derived theoretically in Vršnak et al. (2013), although more work must be done to understand the parameters in the data that might be used to successfully predict this important value. However, the great difference in the two separate γ values is an illustration that the two fronts, while following a similar evolution, are still separate and independent.

To perform an analysis of the error, we assume that the GCS method for the flux rope in COR2 has an error of ±.5 R_☉, with the total 1 R_☉ error representing about 3% of the total FOV of the instrument. Since the shock is often diffuse and difficult to see in COR2, we double this to ±1 R_☉. For both fronts in HI-1, we use an error of ±2.5 R_☉, which is a total error of about 7% of the HI-1 FOV. The increase in error in HI-1 is due to the increase in the ambiguity of the CME fronts further from the Sun, as well as the errors introduced by the image processing techniques required to bring the CME to the forefront as it propagates outward and expands and becomes less intense in the observations. For the drag-model fitting of the ejecta on the 2012 July 12 CME, this yielded a χ²-value of 12.5. With 23 degrees of freedom in the fit, this gives a p-value for the fit of.96. For the shock, the χ²-value was 28.4 with 22 degrees of freedom, for a p-value of just.16. This still indicates a decent fit. Examining the shock fit in Figure 6, it shows a good agreement further from the Sun, but under 30 R_☉ the shock height is ahead of the model fit.

There are a number of possible explanations for this difference, with the possibility that the drag model is more physically descriptive for the flux rope traveling in the solar wind than it is for shock propagation. Still, the fit is good enough to provide an approximation of the shock propagation, and with the study of more events the model may be more adaptable to the shock front.

To confirm the raytrace measurements, they were compared with the jmap data, specifically along the shock, as seen in Figure 8. The ejecta was not compared due to the difficulty in reliably locating it in the jmaps. Both the fixed-φ and harmonic mean are included, rather than attempting to determine a point at which the latter begins to outperform the former. The data sets show qualitative agreement, at least before the fixed φ data begins to blow up, but there are still some significant discrepancies. The jmap data frequently under measures the shock compared to raytrace in the low HI-1 FOV (50 R_☉), but later overestimates it. Instead of representing the drag behavior of the raytrace measurements where the CME starts fast and decelerates to the ambient solar wind speed, the jmap data shows a CME with a lower initial velocity that has a more linear profile.

Zoom In Zoom Out Reset image size

The data are still close enough comparing the jmap with the raytrace measurements, thus it is fair to say they are generally consistent. However, we have much more trust in the raytrace measurement due to the number of different errors introduced in the jmap methods. Not only is the jmap just one two-dimensional projection rather than a stereoscopic reconstruction like raytrace, but there are also a number of image processing errors introduced, as well as the errors from determining where the shock boundary is on the jmap and the assumptions needed to convert elongations into height time measurements. The raytrace model also seems to align better with the in situ arrival point, which means we can be confident with our measurements.

Figure 9 shows the standoff distance, or the distance between the leading edge of the shock front and the leading edge of the ejecta front from the raytrace measurements. The standoff distance (black diamonds) is calculated by subtracting the measured ejecta front from the measured shock front. The line is the standoff distance from the models, which is taken by subtracting the ejecta drag fitting from the shock drag fitting. The solid vertical line shows the arrival of the shock at ACE and the dashed vertical line shows the arrival of the flux rope. The star symbol in the middle of the sheath timings at ACE shows the average in situ sheath size, as determined by integrating the ACE velocity data with time during the observed duration of the sheath. Error bars were obtained by doubling the error of the raytrace measurement in each FOV, since the standoff distance is a combination of two separate measurements. The model fitting matches reasonably well to the raytrace data in the HI-1 range, but there is significant deviation in the COR2 FOV. This is partially due to boundary effects of the parametric fitting. It might also be related to the fact that, nearer the Sun, the size of the standoff distance is small, compared to the error of the distance measurements. The model shows a standoff distance that has an exponential form in the early evolution, but for much of the later propagation, is roughly linear. This is supported by the observed sheath—integrating the velocity data at ACE over the sheath region, an average sheath size during the propagation of the sheath through the satellite can be obtained. The average in situ sheath was 43.5 R_☉, or about 4–12 R_☉ higher than the modeled sheath. This error could be an effect of either the model failing to fully capture the shock evolution, or an effect of the sheath being larger away from the CME nose; the location of the sheath from the ACE measurement is about 10 degrees from the nose of the CME. The model is representing the standoff distance right at the CME nose, whereas the in situ data is not.

Zoom In Zoom Out Reset image size

Our model standoff distance results are in good agreement with Maloney & Gallagher (2011), who found a standoff distance of about 20 R_☉ at 0.5 AU. Likewise, a model of standoff distance based on extrapolating fits to low corona measurements by Temmer et al. (2013) shows a standoff distance profile in the heliosphere that is very similar to our results, which also agree with typical in situ sheath sizes (Russell & Mulligan 2002).

While the previous fittings with the drag model are useful for demonstrating the ability of the drag model to fit the CME ejecta and shock, those fittings are not particularly useful for a predictive purpose, because they include the in situ arrival of the fronts in the fitting. In order to use the model for prediction, the fittings must only contain information that can be obtained early in CME propagation. This means that the in situ arrival time and solar wind speed must be removed, which allows for testing the method as a predictive tool. Furthermore, even without prior knowledge of the solar wind speed and drag parameter, γ, predictions can still be improved by including limits on these values based on studies of different events.

Vršnak et al. (2013) show that a good range for γ is 1 × 10⁻⁸ to 1 × 10⁻⁷ km⁻¹ and a solar wind speed of 300 to 500 km s⁻¹ for a typical slow solar wind upstream of the CME, are useful for prediction. Using these limits on the parameters and fitting the raytrace measurements from eruption until 06:49 UT on the 13th, predicted arrival times could be calculated. The parameters are included in Table 1 for comparison to the parameters with the in situ arrival included.

Both the shock and ejecta front skew toward the maximum allowed solar wind speed of 500 km s⁻¹, which is not surprising. The data being fit is all from within 80 R_☉, at which point the CME is still in the process of slowing down and is traveling well above 500 km s⁻¹. Without the in situ data point being included, no information on the propagation beyond this point is going into the model, so there is no way for the fitting to approach a more realistic final speed. The derived initial speeds are also higher for each front. Even with these higher initial and final speeds, the model still predicts the shock to arrive later than the actual in situ signatures. This is because of a significant increase in γ in both the shock and ejecta in the predictive fitting.

While the initial and final velocities of both structures are higher in the predictive fitting, with the higher amount of drag, the speed declines more quickly and the average speed of the propagation will be closer to the final speed. In the fittings including the in situ data, the drag is more gradual, leading to a higher overall average speed.

That the predictive fitting, which includes only data up to 80 R_☉, has more drag is unsurprising. This is consistent with what was found by Poomvises et al. (2010) that most drag occurs within 50 R_☉. This indicates that rather than using one static γ value, using a decreasing value with a stronger drag early in the propagation could be more physically accurate.

In terms of the success of the predictions, the shock prediction is about four hours late, and the ejecta prediction is about four hours early. Accordingly, the predicted shock velocity is lower than what is observed in situ, and the flux rope velocity is higher. This is an encouraging result for space weather prediction. It had an error of about 18 hr predicting the shock arrival for this same event, when a purely empirical model and far less data are used (Gopalswamy et al. 2013). The prediction also improves upon most of the methods used for this same event in Möstl et al. (2014) for the arrival of the shock, and compares favorably to the most successful prediction therein, the SSEF corrected model.

With the advantage of using data much further out into the heliosphere, we can greatly improve the accuracy of the prediction, with the caveat that even if the data could be accessed, processed, measured, and fit instantly, it would still provide a lead time of just about 30 hr on the prediction for this particular event. As more events are studied, it is our goal to have greater insight on the drag parameter, thus allowing us to predict the drag parameter with less data and have a model that combines the accuracy of the drag model predictions with the lead time of a more empirical model.