Ensemble modeling of CMEs using the WSA-ENLIL+Cone model
Ensemble modeling of coronal mass ejections (CMEs) provides a probabilistic forecast of CME arrival time which includes an estimation of arrival time uncertainty from the spread and distribution of predictions and forecast confidence in the likelihood of CME arrival. The real-time ensemble modeling of CME propagation uses the Wang-Sheeley-Arge (WSA)-ENLIL+Cone model installed at the Community Coordinated Modeling Center (CCMC) and executed in real-time at the CCMC/Space Weather Research Center. The current implementation of this ensemble modeling method evaluates the sensitivity of WSA-ENLIL+Cone model simulations of CME propagation to initial CME parameters. We discuss the results of real-time ensemble simulations for a total of 35 CME events which occurred between January 2013 - July 2014. For the 17 events where the CME was predicted to arrive at Earth, the mean absolute arrival time prediction error was 12.3 hours, which is comparable to the errors reported in other studies. For predictions of CME arrival at Earth the correct rejection rate is 62%, the false-alarm rate is 38%, the correct alarm ratio is 77%, and false alarm ratio is 23%. The arrival time was within the range of the ensemble arrival predictions for 8 out of 17 events. The Brier Score for CME arrival predictions is 0.15 (where a score of 0 on a range of 0 to 1 is a perfect forecast), which indicates that on average, the predicted probability, or likelihood, of CME arrival is fairly accurate. The reliability of ensemble CME arrival predictions is heavily dependent on the initial distribution of CME input parameters (e.g. speed, direction, and width), particularly the median and spread. Preliminary analysis of the probabilistic forecasts suggests undervariability, indicating that these ensembles do not sample a wide enough spread in CME input parameters. Prediction errors can also arise from ambient model parameters, the accuracy of the solar wind background derived from coronal maps, or other model limitations. Finally, predictions of the geomagnetic index differ from observed values by less than one for 11 out of 17 of the ensembles and prediction errors computed from the mean predicted show a mean absolute error of 1.3.
Ensemble modeling has been employed in weather forecasting in order to quantify prediction uncertainties and determine forecast confidence Sivillo, Ahlquist, and Toth (1997). Individual forecasts which constitute an ensemble forecast represent possible scenarios which approximate a probability distribution that reflects forecasting uncertainties. Such uncertainties which when considered as a group include those associated with initial conditions (such as observational uncertainties), techniques and models. Different forecasts in the ensemble can start from different initial conditions and/or be based on different forecasting models/procedures. In the simplest application, the ensemble mean or a weighted mean can be taken as a single forecast. The ensemble mean should perform better than individual ensemble members by emphasizing systematic features found in all members. However, an ensemble also contains additional information about possible scenarios and their probabilities and thus provides a probabilistic forecast. For example, ensemble modeling provides a quantitative description of the forecast probability that an event will occur by giving event occurrence predictions as a percentage of ensemble size. This conveys the level of uncertainty in a given forecast in contrast to a categorical yes/no forecast. Additionally, all ensemble forecast members can be plotted together to allow visualization of the uncertainty among ensemble members, and their clustering distribution. An example of such a visualization is hurricane track “plume” maps in weather forecasting. Regions where members tend to coincide/cluster can be taken to have a higher forecast confidence.
To understand the uncertainties in space weather forecasting, ensemble coronal mass ejection (CME) forecasting efforts have now begun in space weather models of the heliosphere. \inlinecitefry2003, \inlinecitemckenna2006, and \inlinecitesmith2009 compared the performance of real-time shock arrival time forecasts following solar events (since 1997) from the three “Fearless Forecast” models: Shock Time of Arrival (STOA)Dryer (1974), Interplanetary Shock Propagation Model (ISPM)Smith and Dryer (1990), and Hakamada-Akasofu-Fry (HAFv.2)Dryer et al. (2001). While there are many models predicting the evolution of CMEs (see Zhao and Dryer (2014) and references therein), only the Wang-Sheeley-Arge (WSA) coronal model Arge and Pizzo (2000); Arge et al. (2004) coupled with the global heliospheric ENLIL solar wind model Odstrčil (2003) has been used extensively in space weather operations world-wide. The first effort in utilizing this model for ensemble forecasting of CME propagation was reported by \inlinecitepulkkinen2011. \inlineciteemmons2013 performed WSA-ENLIL ensemble CME modeling using 100 ensemble members for 15 historical events with automatically determined cone model CME parameters Pulkkinen, Oates, and Taktakishvili (2010). They found that the observed CME arrival was within the ensemble prediction spread for 8 out of the 15 events. \inlinecitelee2013 discuss ensemble modeling of CME propagation with WSA-ENLIL for an event study using eight ensemble members and various synoptic background maps. Differences found in the predicted arrival time of each individual simulation were mostly due to CME initial speed and the time at which the CME was inserted at the WSA-ENLIL inner boundary, resulting in propagation through a different background solar wind. They used National Solar Observatory Global Oscillation Network Group (GONG) Harvey et al. (1996) synoptic magnetograms and Air Force Data Assimilative Photospheric flux Transport (ADAPT) maps Arge et al. (2010); Henney et al. (2012). For their CME event they show that when using ADAPT maps, the WSA-ENLIL model values were in better agreement with in-situ observations, and the arrival time predictions were improved due to the more accurate background solar wind representation. However, the overall spread in CME arrival times did not change significantly.
This paper describes the WSA-ENLIL+Cone ensemble modeling system installed at the Community Coordinated Modeling Center (CCMC) and results from the past 1.5 years of real-time execution at the CCMC/Space Weather Research Center. This is the first ensemble space weather prediction system for CME propagation of its kind employed in a real-time environment. The current version of the system evaluates the sensitivity of CME arrival time predictions from the WSA-ENLIL+Cone model to initial CME parameters. The CCMC, located at NASA Goddard Space Flight Center, is an interagency partnership to facilitate community research and accelerate implementation of progress in research into space weather operations. The SWRC is a CCMC sub-team which provides space weather services to NASA robotic mission operators and science campaigns, and prototypes new models, forecasting techniques and procedures. The CCMC also serves the CME Scoreboard website111http://kauai.ccmc.gsfc.nasa.gov/CMEscoreboard to the research community who may submit CME arrival time predictions in real-time for a variety of forecasting methods. The website facilitates model validation under real-time conditions and enables collaboration. For every CME event table on the site, the average of all submitted forecasts is automatically computed, thus itself providing a world-wide ensemble mean CME arrival time forecast from a variety of models/methods.
In Section 2 a brief description of the WSA-ENLIL+Cone model is given. The triangulation algorithm for determining CME parameters for the ENLIL model is described in Section 3. The real-time ensemble modeling methodology is explained in Section 4 followed by an example of an ensemble simulation given in Section 5. Results and the evaluation of the first 1.5 years of simulations are described in Section 6. In Section 7 we discuss a parametric event case study of the sensitivity of the CME arrival time prediction to model free parameters for the CME and ambient solar wind. Finally, a summary and discussion are presented in Section 8.
2 WSA-ENLIL+Cone Model Description
The global 3D MHD WSA-ENLIL model provides a time-dependent description of the background solar wind plasma and magnetic field into which a CME can be inserted Odstrčil, Smith, and Dryer (1996); Odstrčil and Pizzo (1999a, b); Odstrčil (2003); Odstrčil, Riley, and Zhao (2004). This modeling system does not simulate CME initiation but uses kinematic properties of CMEs inferred from coronagraphs to launch a CME-like hydrodynamic structure into the solar wind and interplanetary magnetic field computed from the WSA coronal model Arge and Pizzo (2000); Arge et al. (2004). A common method to estimate the 3D CME kinematic and geometric parameters is to assume that the geometrical CME properties are approximated by the Cone model Zhao, Plunkett, and Liu (2002); Xie, Ofman, and Lawrence (2004) which assumes isotropic expansion, radial propagation, and constant CME cone angular width. Generally, a CME disturbance is inserted in the WSA-ENLIL model as slices of a homogeneous spherical plasma cloud with uniform velocity, density, and temperature as a time-dependent inner boundary condition at 21.5 solar radii () with an unchanged background magnetic field. While the simplest geometrical case is employed in this work, the WSA-ENLIL+Cone model can also support an elliptical geometry including tilt, an elongated spheroid or ellipsoid, and leading and trailing edge velocities. Measurements derived from coronagraphs (described in Section 3.1) determine the cloud velocity, location, and width. The CME cloud density (dcld) is a free parameter which by default is 4 times larger than typical mean values in the ambient fast wind providing a pressure of four times larger than that in the ambient fast wind. The cloud temperature is taken to be equal to the ambient fast wind temperature. Another ENLIL free CME parameter is the cavity ratio which allows the CME to be represented by a spherical shell of plasma and is based on coronagraph observations of CME cavities. The cavity ratio radcav is defined as the ratio of the radial CME cavity width to the CME width, with the default being no cavity (radcav=0).
WSA-ENLIL+Cone runs performed for research and operations have shown that accurate descriptions of the heliosphere and transients are achieved only when the background solar wind is well-reproduced and if coronagraph observations from multiple views, for example from the SOlar and Heliospheric Observatory (SOHO) spacecraft near the Earth Domingo, Fleck, and Poland (1995) and the Solar TErrestrial RElations Observatory (STEREO) spacecraft Kaiser et al. (2008), are used to derive CME parameters Lee et al. (2013); Millward et al. (2013). WSA coronal maps provide the magnetic field and solar wind speed at the boundary between the coronal and heliospheric models, usually at 21.5 , and they are generated from synoptic magnetograms. Small latitudinal shifts in the magnetogram-derived coronal maps caused by inaccuracies in solar magnetic field observations, particularly in the polar regions, can cause large longitudinal shifts in the solar wind structure, for example in characterizing high speed stream arrival times (e.g. \inlinecitemacneice2009; \inlinecitejian2011). Other coronal models, such as MAS (MHD around a Sphere) Riley, Linker, and Mikić (2001) or heliospheric tomography from interplanetary scintillation (IPS) Jackson et al. (2011) can also provide the background solar wind and have been coupled with ENLIL heliospheric simulations.
CCMC/SWRC has been carrying out routine WSA-ENLIL+Cone simulations for several years using solar magnetic synoptic maps and CME geometric and kinematic properties inferred from coronagraph observations Zheng et al. (2013). Each ENLIL run uses a WSA model synoptic map computed from the single GONG daily-updated synoptic magnetogram (see e.g. \inlinecitearge2000) closest to the time the simulation is executed. These low 4 resolution real-time simulations complete in 20 minutes running on 2 nodes with 16 processors/node on a spherical grid size of 2563090 () with 5-10 minute output cadence at locations of interest. The simulation range is 0.1 to 2 AU in radius , -60 to +60 in latitude , and 0 to 360 in longitude . CME parameters are derived using real-time coronagraph observations from spacecraft and a geometric triangulation algorithm. The measurements are an approximation of the true 3D speed and width of the CME at 21.5 (ENLIL inner boundary). However, often the coronagraph derived measurements are inferred from just a few data points, and some CMEs may be missed due to real-time data gaps. CME parameters derived in real-time and simulation graphical outputs are publicly available from the CCMC Space Weather Database Of Notifications, Knowledge, Information (DONKI) database222http://kauai.ccmc.gsfc.nasa.gov/DONKI.
3 Ensemble CME Parameters
3.1 Stereocat Triangulation Algorithm for Determining Cme Parameters
CME parameters are determined using the Stereoscopic CME Analysis Tool (StereoCAT), developed by the CCMC for real-time CME analysis carried out by the CCMC/SWRC forecasting team. The goal was to develop a tool that can be used quickly, yet reliably in a real-time environment with any possible combination of spacecraft available for analysis. It was also required that the tool was intuitive and simple enough to be employed by a wide variety of users such as space weather forecasters, scientists, students, and citizen scientists. The basic methodology of the tool, i.e., tracking of CME kinematic properties from two different fields-of-view, is similar to that of the NOAA Space Weather Prediction Center CME Analysis Tool (CAT) by \inlinecitemillward2013 and the geometric localization developed by \inlinecitepizzo2004. However, StereoCAT does not attempt to capture the volumetric structure of CMEs but is based on tracking specific CME features. The algorithm is most similar to the CME geometric triangulation method of \inlineciteliu2010a. For a more detailed discussion of different CME analysis techniques in the context of cone model-based CME simulations, see Pulkkinen, Oates, and Taktakishvili (2010); Millward et al. (2013). Other stereoscopic methods for determining the kinematic properties of CMEs include those by \inlinecitethernisien2006, \inlinecitelugaz2010, and \inlinecitedavies2013.
StereoCAT is based on triangulation of transient CME features from two different coronagraph fields-of-view or planes-of-sky. We will call these planes-of-sky and , which may designate, for example, fields-of-view of the SECCHI/COR2 instruments onboard the STEREO A and STEREO B spacecraft Howard et al. (2008). The tool is used to manually identify the same CME features in two consecutive images which are then used to calculate the plane-of-sky velocities for and , and , respectively. Note that these velocities are in local plane-of-sky coordinates indicated by and . These data need to be brought into the same coordinate system (Heliospheric Earth Equatorial (HEEQ) coordinates in this case), which can be accomplished by rotations:
where operators and carry out transformations from planes-of-sky and coordinates into a common base such as HEEQ, respectively.
We then define two projection matrices as
where is a 3 3 identity matrix. The unit vectors normal to the planes-of-sky of coronagraphs and are defined as and , where is the transpose of matrix . The matrices and project any vector to plane-of-skies and , respectively. Consequently, plane-of-sky speeds can be expressed as
from which we can solve
Importantly, exists as long as planes-of-sky and are different, i.e. when and are not co-linear (parallel to each other). Therefore large triangulation errors occur when the spacecraft separation angle is very small or around 180.
A similar approach can be used to track the three-dimensional location of a feature from plane-of-sky measurements and as
Often the time stamps of coronagraph imagery from spacecraft and do not match exactly. This is handled in StereoCAT by propagating the tracked feature in with speed to a new that matches the time stamp. Consequently, matching time stamps are used for and in Eq. (9).
The angular size of a CME is estimated in StereoCAT simply by manually selecting the two outer edges of the CME. These two lines that connect through the center of the Sun are then used to compute the opening angle of the CME. It is noted that this process does not take into account projection of the outer CME edges to the spacecraft plane-of-sky, and is therefore a measurement of the projected CME width. While this is not an issue if the CME propagation direction is not too far away from the plane-of-sky of the spacecraft which is used to measure the opening angle, one needs to be very careful with events with propagation directions substantially away from the plane-of-sky, as in such cases the opening angle can be overestimated. This issue will be addressed in the future versions of StereoCAT.
Other limitations of StereoCAT arise from the user’s ability to reliably identify the same structures in images from both spacecraft due to ambiguities from the different viewing angles. It may at times be difficult or impossible to track the same structure since different sections of the CME contribute most strongly to images in different planes-of-sky Howard and DeForest (2012). Consequently, StereoCAT is not suitable to use with coronagraph data in which the CME appears as a halo, since the CME leading edge is not visible.
3.2 Performing Cme Measurements With Stereocat
StereoCAT has three modes: “two-timepoint”, “ensemble”, and “frame series”, and is available online via a web interface333http://ccmc.gsfc.nasa.gov/analysis/stereo LaSota (2013). Available coronagraphs include the LASCO C2 and C3 instruments on board the SOHO spacecraft Brueckner et al. (1995), and the SECCHI/COR2 instruments on the STEREO A and B spacecraft. All three modes are based on the same triangulation algorithm, described in section 3.1. In the basic “two-timepoint” mode the user manually measures the CME leading edge height for two different times in each coronagraph image for two different coronagraph viewpoints. The plane-of-sky sky speed for each viewpoint is calculated, from which the triangulated speed and direction is computed using the algorithm described in section 3.1. The user also manually measures the CME opening angle in each coronagraph view. Because this is a projected width measurement, both widths and their average are displayed for the user.
In “ensemble” mode the user manually repeats the same procedure as for the “two-timepoint” mode, by measuring the same feature for the same pair of coronagraphs at two different times. Between each “two-timepoint” measurement, the display is fully reset such that the user is forced to carefully remeasure the CME leading edge height and opening angle. This series of repeated measurements leads to a range of CME parameters which can be used to initialize an ensemble simulation. For every “two-timepoint” measurements made, ensemble CME parameter members are automatically generated by combining different spacecraft measurement pairs. For example, for “two-timepoint” measurements, there are ways to combine the first and second time step height measurements in viewpoints and to triangulate the CME. Since the two projected width measurements made for each measurement are not triangulated, they are randomly assigned to each ensemble member. An example screenshot of =6 “two-timepoint” measurements performed in “ensemble” mode using StereoCAT is shown in Figure 1. Two image pairs are shown from the SECCHI/COR2 instruments from STEREO B (top row) and STEREO A (bottom row), for two different times separated by 30 minutes in the left and right columns. The white circles indicate the 6 individual “two-timepoint” plane-of-sky leading edge height measurements (near the center of the CME front) and the width measurements are marked by the green circles (near the CME edges). The green lines in panel c of Figure 1 illustrate the CME opening angle measurements for one of the coronograph images. In this example the 6 individual “two-timepoint” measurements were combined by the algorithm to create ensemble members.
After completing the measurements the user may inspect histograms of their CME parameters. The web interface allows the user to remove any ensemble members, and add any “custom” members. Generally, members are removed when they have nearly identical parameters, or members for which the triangulation appears unreliable. Custom members can be measurements from different image time pairs, from plane of sky estimates which incorporate the source location, or from any other CME measurement technique. The same procedure can be applied to create individual ensemble measurements for CMEs for a series of events which are then combined one-to-one to be simulated together such that there are ensemble members containing CMEs each.
In “frame series” mode the user can measure a series of different frames (times) for each spacecraft, which are then triangulated to create a CME height-time profile. The user selects a range of time, and steps through the images available from each instrument, measuring the CME in as many images as they choose. The software chooses time pairs of measurements for triangulation based on a user-specified maximum allowed time difference. From these measurements, plane-of-sky and triangulated height-time, velocity, latitude, and longitude profiles of the CME are generated. Triangulations made with different spacecraft pairs are shown as separate height-time profiles. Several methods are used to calculate the CME speed, acceleration, the time the CME passes 21.5 (ENLIL inner boundary), and the time it erupts from the Sun. These include least-squares linear and quadratic fits, averages over selected data points, and averages from only the first and last data points. Results for each method are reported separately, allowing the user to choose the most appropriate fitting technique depending on the acceleration profile of the CME. Plane-of-sky values are also reported, which can be used when coronagraph projection effects make this triangulation method unreliable. This can occur if the CME is very wide, appears as a halo, or is heavily projected in the coronagraph data. In these cases the user will not be able to identify the same CME leading edge feature in the data from two coronagraphs. The user can inspect the triangulated height values directly on the height-time plot to evaluate triangulation accuracy in these cases.
4 Ensemble Modeling with WSA-ENLIL+Cone
The current implementation of this ensemble modeling method evaluates the sensitivity of WSA-ENLIL+Cone model simulations of CME propagation to initial CME parameters. As described in Section 3.1, StereoCAT is used to create an ensemble of CME parameters which are used as input to WSA-ENLIL+Cone simulations. We have observed that 36 to 48 provides an adequate spread of input parameters, but this number can be increased if necessary. For =48 a typical run takes 130 minutes to complete on 24 nodes with 4 processors/node on the initial development system. We estimate that the same run will take 80 minutes on the CCMC production system that has 16 processors/node.
The simulations provide profiles of MHD quantities (density, velocity, temperature, and magnetic field components) and a distribution of predicted arrival times at locations of interest within the computational domain. Currently, ensemble modeling is performed for spacecraft at the following locations: Mercury (MESSENGER), Venus (VEX), Earth (ACE, Wind, SOHO, and orbiting spacecraft), Mars (MSL, MAVEN, MEX), Spitzer Space Telescope, STEREO-A and B. The CME-associated disturbance/shock arrival time is then automatically computed in post-processing from any sharp increases in the modeled solar wind dynamic pressure at a given location. In this work, we focus on the ensemble results of the Earth-directed events.
For Earth-directed CMEs, the CCMC/SWRC also computes estimates of the geomagnetic index using the WSA-ENLIL+Cone model plasma parameters at Earth. The geomagnetic three-hour planetary index, , is a measure of general planetary wide geomagnetic disturbances at mid-latitudes based on ground-based magnetic observations Bartels, Heck, and Johnston (1939); Rostoker (1972); Menvielle and Berthelier (1991). The index is created from standardized indices from individual stations, which measure the magnitude of horizontal geomagnetic field disturbances (not including daily variations). is a quasi-logarithmic index ranging from 0 to 9. Real-time estimated planetary indices are available from NOAA using real-time data from a limited number of geomagnetic observatories, and the final definitive is from the Helmholtz Center Potsdam GFZ German Research Centre for Geosciences.
The predicted estimate is made by using the \inlinecitenewell2007 coupling function arising from their correlation of 20 candidate coupling functions with geomagnetic indices. The function which represents the rate of magnetic flux opening at the magnetopause and correlated best with 9 out of 10 indices is given as
where is the bulk solar wind speed, the interplanetary magnetic field (IMF) clock angle is given by , and the perpendicular component of the magnetic field is given by (in GSM coordinates). An exponential fit to the correlation of this coupling function with the index yields the following relation used for the estimate
emmons2013 showed for their sample of 15 events, that predictions using Eq. 11 computed directly from in-situ solar wind observations had a mean absolute error of 0.5. Because ENLIL modeled CMEs do not contain an internal magnetic field and the magnetic field amplification is caused mostly by plasma compression, only the magnetic field magnitude is used and three magnetic field clock angles scenarios of 90 (westward), 135(south-westward), 180 (southward) are assumed. This provides a simplistic estimate of three possible maximum values which the index might reach following arrival of the predicted CME shock/sheath. For the forecast, the estimates are rounded to the nearest whole number.
Another commonly used activity index is the (disturbance storm time) index, which is a measure of magnetosphere storm activity primarily from the strength of the ring current. The index is obtained from the measurement of the perturbations in the horizontal component of the Earth’s magnetic field from ground-based observatories that are sufficiently distant from the auroral and equatorial electrojets, are located at approximately geomagnetic latitude, and are evenly distributed in longitude Sugiura (1964). Although the ring current makes the largest contribution to the , all magnetospheric current systems contribute, such as the Chapman-Ferraro magnetopause current which is strengthened during sudden storm commencement (SSC) and increases the Earth’s surface field and gives a sudden positive jump in . Currently, ENLIL model results are not used to predict the , however in principle this can be computed in a similar manner to the index by using the \inlinecitenewell2007 relation.
5 Example Ensemble: 18 April 2014 CME
In this section we describe the real-time ensemble modeling of an Earth-directed partial halo CME that was first observed at 13:09 UT on 18 April 2014 by by SECCHI/COR2A. Figure 2 shows this CME as viewed from SOHO LASCO/C2 and C3, STEREO SECCHI/COR2 A, and B near 14:50 UT. This CME was associated with an M7.3 class solar flare from Active Region (AR) 12036 located at S18W29 with peak at 13:03 UT. The eruption and a coronal wave were visible south of the active region in SDO/AIA 193Å and a nearby filament eruption was visible in AIA 304Å. Subsequently starting at 13:35 UT, an increase in solar energetic particle proton flux above 0.1 pfu/MeV (1 pfu = 1 particle cm sr s) was observed by the GOES-13 EPEAD instrument in Earth orbit.
Figure 1 shows StereoCAT measurements for the 18 April 2014 CME. As discussed above, the central white circles indicate the individual leading edge measurements and the green outer circles near the CME edges are the projected width measurements. The six leading edge measurements are combined together using the triangulation algorithm discussed in Sections 3.1-3.2 to generate ensemble members. The distribution of the resulting CME parameters which are used as initial conditions for =36 WSA-ENLIL+Cone simulations is shown in Figure 3 in the (a) equatorial plane (latitude=0) and (b) meridional plane (longitude=0). The plots show the CME velocity vectors in spherical HEEQ coordinates with the grids showing the degrees longitude (a) and latitude (b), and the radial coordinate showing the speed in km/s. The Sun-Earth line is along 0 longitude and latitude. The arrow directions on the grid indicate the CME central longitude and latitude respectively, with CME half width indicated by the color of the vector. The arrow lengths correspond to the CME speed. CME propagation directions are clustered between -30 to -40 latitude, and around 10 west of the Sun-Earth line in longitude, while CME speeds range from 1300 to 1600 km/s. Median CME parameters are: speed of 1394 km/s, direction of 9 longitude, -35 latitude, and a half-width of 46.
Model results for the 36-member ensemble WSA-ENLIL+Cone run for this CME are shown in Figures 4-5. For the ensemble member with median CME input parameters, Figure 4 shows a scaled velocity contour plot for the (a) constant Earth latitude plane, (b) meridional plane of Earth, and (c) 1 AU sphere in cylindrical projection on 20 April at 06:00 UT. Panel (d) shows the measured (red) and simulated (blue) radial velocity profiles at Earth, with the simulated CME duration shown in yellow. This simulation figure shows the northeastern portion of the CME impacting Earth. Figure 5 shows the modeled magnetic field, velocity, density, and temperature profiles at Earth plotted as color traces for all 36 ensemble members, along with the observed in-situ L1 observations from ACE, plotted in black. The model traces are color coded by CME input speed such that slow to faster input speeds are colored from light green to dark blue. The arrival of the CME-associated shock was observed by Wind and ACE on 20 April 2014 at around 10:20 UT, and energetic storm particles were observed by ACE. The provisional SYM-H index (1 minute ) shows a sudden storm commencement of +25 nT at 11:01 UT. The observations in Figure 5 show clear signatures of the arrival of an Interplanetary Coronal Mass Ejection (ICME), including a leading shock (abrupt increase in all the solar wind parameters at around 10:20 UT) with enhanced post-shock temperatures, enhanced magnetic field with rotations in direction, and declining solar wind speed. This CME was predicted to arrive at Earth and also at Mars for all of the 36 runs. The mean predicted arrival at Earth was on 20 April 2014 at 05:07 UT with arrival times from individual runs ranging from 20 April 2014 at 01:08 to 11:16 UT. A histogram showing the distribution of arrival times at Earth is shown in Figure 6 with individual arrivals marked by the blue arrows. This figure shows a normal distribution with 50% of the predicted arrivals within one hour of the mean. The prediction error for the mean predicted CME arrival time was -5.2 hours and the observed arrival time was just within the ensemble predicted spread. The spread in ensemble member predictions can also be seen in Figure 5 compared to the observations, showing that most of the predictions are earlier than the observed arrival with a few after. From the CME input parameters plotted in Figure 3 the ensemble members with arrival times closest to the observed time had CME input speeds in the range of 1200-1400 km/s, latitudes near -40 and half widths around 35-40. This suggests that the early arrival time predictions for this event could be due to overestimations of the CME input speed and half width.
The NOAA real-time observed index (and the Potsdam final Kp) reached 5 during the synoptic period 12:00-15:00 UT on 20 April associated with the CME shock arrival. The reached a minimum of -24 nT at 15:00 UT on 21 April and thus based on , this CME only resulted in very weak geomagnetic activity. As discussed in Section 4, Eq. 11 can be used to forecast the maximum index from maximum ENLIL predicted quantities at CME shock/sheath arrival at Earth (colored traces shown in Figure 5). Figure 7 shows the predicted probability distribution of for three clock angle scenarios (green), 135 (purple), 180 (orange). The figure also shows the overall forecast probability distribution calculated for all three angles combined 90-180, assuming each scenario is equally likely, in black. The standard deviation of the overall forecast probability distribution is 1.1, with 84% of the forecasts falling between = 5 to 7. The most likely forecast is for =7 at 41%, followed by =5 at 27% and =6 at 16% likelihood of occurrence. Using the most likely forecast of =7, the prediction error for this event is (overprediction). The overprediction of may be related to the overestimation of the CME input speed. In Sections 6.1-6.2 and 7 we discuss various factors which can contribute to early arrival time predictions and Kp overpredictions.
6 Real-time Ensemble Modeling: First Results
For 35 Earth-directed CME events from January 2013 through June 2014, real-time ensemble modeling was carried out by the CCMC/SWRC team following the methods described in Sections 4-5. In Table 1 we list a summary of the ensemble simulation results for these 35 CME events. The first and second columns give the CME onset date and time based on the first appearance in C2 or COR2. Generally, if two CMEs occur within a day of each other they will both be included in the same simulation as separate CMEs which may or may not merge during their propagation. A few of the ensemble simulations listed in the table contain two CMEs as part of a single run. In these cases, CMEs that were simulated together with the CME listed on the previous row are indicated by . The third column lists (for 2013) the 2nd order plane-of-the-sky (POS) speed at 20 reported in the SOHO LASCO CDAW CME catalog444http://cdaw.gsfc.nasa.gov/CME_list Yashiro et al. (2004); Gopalswamy et al. (2009). If measurements were not made to 20 , the 2nd order POS speed at the time of last observation is used. The next four columns provide the median ensemble CME input parameters of , latitude, longitude (HEEQ), and half-width () measured using StereoCAT. In columns 8, 9 and 10 we list the mean predicted arrival time of all ensemble members, followed by the spread in arrival times in hours relative to the mean. The next column (11) shows , the number of ensemble members out of , the total number of ensemble members, that predict that the CME will arrive at Earth. This ratio gives a forecast probability and conveys the forecast uncertainty about the likelihood that the CME will arrive. Columns 12, 13, and 14, list the actual arrival time of the CME-associated shock or disturbance observed in-situ at the Wind spacecraft, followed by the total in-situ observed CME transit time relative to the CME start time. In the last column the prediction error is calculated for predictions indicating hits. The prediction error is defined as , which is negative when ENLIL predictions are earlier than the observed CME arrival time, and late predictions are positive. When possible, ICME and magnetic cloud catalogs were used to help assess whether the CME did arrive at Earth. These included the \inlineciterichardson2010 ICME catalog555http://www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm, and the Wind ICME catalog666http://wind.nasa.gov/index˙WI˙ICME˙list.htm with circular flux rope model fitting (based on \inlinecitehidalgo2000). Shocks identified by the SOHO CELIAS/MTOF/PM “shockspotter” program were also used in arrival time assessment. Determining the measured in-situ arrival time of the CME-associated shock or disturbance can be subjective and therefore can be a source of error in the prediction error calculation. Taking this into consideration, in-situ signatures which could not be unambiguously identified as the arrival of the CME-related disturbance are indicated by and these five ensembles are not included in the following forecast verification. This reduces the sample size from 35 to 30 ensembles.
In the following subsections we discuss ensemble CME arrival and forecast verification inspired by methods used in ensemble weather forecasting and applied here for the first time.
6.1 Cme Arrival Forecast Verification
To begin with a simple forecast evaluation of CME arrival time, the ensemble mean can be taken as a single forecast. Using the prediction error (last column of Table 1, the mean absolute error (MAE) is 12.3 hours, the Root Mean Square Error (RMSE) is 13.9 hours, and the mean error (ME) is -5.8 hours (early) for all 17 ensembles containing hits. Considering the sample size in this study, these errors are comparable to CME arrival time prediction errors (a RMSE of 10 hours) reported by others Millward et al. (2013); Romano et al. (2013); Vršnak et al. (2014); Mays et al. (2014). Similarly, \inlinecitecolaninno2013 used a variety of methods to evaluate CME arrival time predictions (non real-time) based on imaging data analysis only, and found an error 6 hours for 78% of their sample, and 13 hours for their full sample of 9 CMEs. The CME arrival time prediction error is inevitably related to the CME propagation speed, thus it is useful to consider the input speed and in-situ observed transit time relative to the prediction error. For this sample, the average in-situ observed transit time was 66 hours. In Figure 8a the CME arrival time prediction error is plotted against the CME input speed, and in Figure 8b the prediction error as a percentage of the CME transit time is plotted against the CME input speed. The error bars are computed using the predicted ensemble range as listed in column 10 of Table 1. The dashed horizontal line indicates the mean arrival time prediction error (a) and mean of the prediction error/transit time percentage (b). These figures show a nearly consistent negative prediction error for fast CMEs above 1000 km/s such that these fast CMEs are generally predicted to arrive earlier than they are observed. This could be a sign of the modeled CME having too much momentum as defined by a combination of the input speed and half width (which is related to the modeled CME mass). The overestimation of the modeled CME velocity compared to in-situ observed values is also due to the modeled CME having a lower magnetic pressure than is observed in typical magnetic clouds.
|CME Onset||SOHO||Median CME parameters||Mean Predicted Arrival||In-situ Arrival Transit|
|Date Time||Lat||Lon||Date Time Spread||Date Time Time|
|(yyyy-mm-dd) (UT)||(km/s)||(km/s)||()||()||()||(UT)||(h)||(UT) (h)||(h)|
|CMEs was simulated together with the CME listed on the previous row as part of a single ensemble.|
|In-situ signature could not be unambiguously identified as arrival of CME-related disturbance and is not included in forecast verification.|
|CME arrival forecast|
|Observation||Will occur||Will not occur|
|Occurs||Hit (17)||Miss (0)|
|Does not occur||False alarm (5)||Correct rejection (8)|
Ensemble modeling produces a probabilistic forecast of the likelihood of CME arrival for each ensemble , but we begin with a more simple forecast evaluation by binning the probability into a categorical yes/no forecast. Categorical forecasts only have two probabilities, zero and one. Therefore we start by binning the probability forecast into two categories: “yes” the CME will arrive, and “no” the CME will not arrive. In the signal detection theory model of weather forecasting, event forecasting performance can be evaluated in terms of a 22 contingency table, as shown in Table 2 Harvey et al. (1992); Weigel et al. (2006); Jolliffe and Stephenson (2011). For CME arrival prediction, the “event” is taken as the “CME arrival”. Hits are then defined as CME arrivals which were both predicted and observed to occur. Misses are defined as CME arrivals which were not predicted, but were observed to occur. False alarms (FA) are defined as CME arrivals that were predicted to occur, but were observed not to occur. And correct rejections (CR) are CME arrivals that were not predicted, and were observed not to occur. To bin each ensemble’s probabilistic forecast, correct rejections were identified when the criterion of the forecast probability 15% was met; i.e., that less than 15% of the total predictions in the ensemble indicated CME arrival. Similarly, the inverse criterion is used to identify hits. Table 2 shows the contingency table definitions and values for this 30 event sample: 17 hits, 8 correct rejections, 5 false alarms, and 0 misses (see Table 1 for specific CR and FA events). For this sample zero misses indicates that there were no ensemble simulations which did not predict CME arrivals which were observed to occur. There were 8 out of 30 correct rejections and 5 false alarms for events that were not observed in-situ, giving a correct rejection and false-alarm rate of 62% (8/ 13) and 38% (5/ 13) respectively. The correct alarm ratio, defined as the number of hits over the number of hits and false alarms, is 77% and the false alarm ratio is 23%.
Let us now consider a more nuanced technique to evaluate the probabilistic forecast without partitioning it into a categorical forecast with only two probabilities as described above. A method defining the magnitude of probability forecast errors is the Brier Score () Brier (1950); Murphy (1973); Wilks (1995), defined as
where is the number of events, is the forecast probability of occurrence for event , and is 1 if the event was observed to occur and 0 if it did not occur. For CME arrival prediction, the “event” here is taken as the “CME arrival” and is listed in column 11 of Table 1 for each ensemble. This score is a probability mean square error which weights larger errors more than small ones and ranges from 0 to 1, with 0 being a perfect forecast. The computed from all = 30 ensemble CME arrival probabilities (Table 1, column 11) is 0.15, which indicates that in this sample, on average, the probability of the CME arriving is fairly accurate. However, such verification scores reduce the problem to a single measure which can only consider one dimension, whereas there are many dimensions to the system. For example, consider the aspect of forecast reliability. Reliable forecasts are those where the observed frequencies of events are in agreement with the forecast probabilities. To evaluate the reliability of probabilistic ensemble forecasts, a set of probabilistic forecasts must be evaluated using observations that demonstrate that those events either occurred or did not occur. Multiple forecasts must be evaluated because a single probabilistic forecast cannot be simply assessed as “right” or “wrong” e.g. if a forecast suggests a 30% chance of CME arrival, and the CME does arrive, the forecast is not clearly either “right” or “wrong”. Therefore, to provide forecast verification for a =30% chance of CME arrival one would need to compile the statistics of observed CME arrivals for a set of forecasts that predicted a 30% chance of arrival. In this way a reliability diagram can be constructed to determine how well the predicted probabilities of an event correspond to their observed frequencies Wilks (1995); Jolliffe and Stephenson (2011). Figure 9a shows the reliability diagram of the likelihood of CME arrival forecast for the 30 event sample, with the reliability for this sample shown as the black line with points and the perfect reliability diagonal as a dotted line. The line of perfect reliability is diagonal because, for example, when a 60% probability forecast is made, it is considered perfectly reliable if the event is observed to occur 60% of the time over multiple ensemble forecasts. The number of events used in each calculation is shown next to each point, and the sample size is smaller than needed for a robust diagram. Nevertheless, the diagram shows that overall ensemble modeling is underforecasting in the forecast bins between 20-80%, and slightly overforecasting in the 1-20% and 80-100% forecast bins. Overforecasting is when the forecast chance of CME arrival (forecast probability) is higher than is actually observed; i.e., the CME is observed to arrive less often than is predicted. Similarly, underforecasting is when the chance of CME arrival is lower than is actually observed; i.e., the CME is observed to arrive more often than is predicted.
Another aspect of forecast reliability is to assess how well the ensemble spread of the forecast represents the true variability of the observations. For 8 out of 17 of the ensemble runs containing hits the observed CME arrival was within the spread of ensemble arrival time predictions. This indicates that roughly half of the observations fall outside of the extremes of the predicted ensemble spread. However, one aspect of a reliable forecast is that the set of ensemble member forecast values for a given event and observations should be considered as random samples from the same probability distribution. This reliability then implies that if an member ensemble and the observation are sorted from earliest to latest arrival times, the observation is equally likely to occur in each of the possible “ranks”. Therefore a histogram of the rank of the observation, “rank histogram”, tallied over many events should show be uniform (flat) Anderson (1996); Hamill and Colucci (1997); Talagrand, Vautard, and Strauss (1997). While more samples would be desirable, it is still instructive to examine the rank histogram for the CME arrival time predictions from the 17 ensembles containing hits in this sample, shown in Figure 9b. Since each ensemble run in our sample does not have the same number of members, the rank has been normalized to 10 (9 member ensemble). To construct this rank histogram the CME arrival time predictions of each ensemble are sorted from earliest to latest and the rank of where the observed arrival falls among the predicted times is noted. For example, an ensemble with a rank of 8 has the meaning that 7 arrival time predictions fall before the observed arrival, a rank of 10 would mean that all 9 predictions occur before the observation, and a rank of 1 means that the observation occurs before all of the predictions. The non-uniform U-shape of this histogram partly illustrates that roughly half of the observed arrivals are outside the spread of predictions (ranks 1 and 10), with a tendency for an overall early spread of predictions (rank=10) compared to observations (also quantified by mean arrival time error of -7.0 hours). U-shaped rank histograms can indicate lack of variability in the ensemble, but can also be a sign of a combination of conditional biases in the model Hamill (2001). However, when evaluating the WSA-ENLIL+Cone model in this sample of ensembles, and 70 regular runs containing hits performed by SWRC Romano et al. (2013), an overall negative bias (early predictions) was found, with less bias for CME input speeds below 1000 km/s. Therefore, it is unlikely that a combination of positive and negative model biases within the ensembles contributed to the U-shaped rank histogram for our sample. Most likely, the U-shape suggests undervariability, indicating that these ensembles to not sample a wide enough spread in CME input parameters.
6.2 Forecast Verification
For each event for which a hit is predicted in Table 1, ensemble modeling provides a probabilistic forecast (see Section 4) for three magnetic field clock angles scenarios of 90 (westward), 135 (south-westward), 180 (southward). An overall probabilistic forecast can then be obtained by making the simple assumption that each clock angle is equally likely to occur. Table 3 lists the overall probabilistic forecast for each bin (e.g. the distribution shown in Figure 7 in black) for these 17 events. The observed , sudden storm commencement (SSC) and minimum indices are also shown. The mean predicted is listed in column 12, along with the overall predicted spread (using plus or minus notation). Underlined probabilities indicate that the NOAA real-time observation falls within this bin, and the final definitive values are listed in column 13. The values are from the real-time (quicklook) index provided by the World Data Center for Geomagnetism in Kyoto, Japan. In order to estimate the reliability of the probabilistic forecast, the Brier Score is calculated for each bin and listed on the last line of the table.
To evaluate forecast performance, a single categorical predicted forecast can be derived from the probabilistic forecast distribution. For example, the single categorical forecast can be taken as the mean predicted , or the most probable value. This allows a prediction error to be computed as for each ensemble, where positive values of indicate an over prediction of the index and negative values indicate that has been under predicted. If the categorical is taken as the bin which has the highest likelihood in the probabilistic forecast for each ensemble, the prediction errors are calculated to give a mean absolute error (MAE) of 1.9, Root Mean Square Error (RMSE) of 2.5, and mean error (ME) of +1.4. However, if the categorical is taken as the mean predicted in each ensemble (last column of Table 3) these errors are reduced to MAE=1.5, RMSE=2.0, and ME=+0.6. Consequently, utilizing the ensemble mean yields a more accurate forecast in this sample, however both forecast choices show an overall tendency for the overprediction of . Given that the modeled CMEs do not have an internal magnetic field structure, the \inlinecitenewell2007 coupling function using ENLIL results as input performs surprisingly well. For comparison, using ACE solar wind data as input to the coupling function for this sample gives prediction errors of MAE=0.67, RMSE=0.77, and ME=+0.22.
In Figure 11 the prediction error (from the ensemble mean ) is compared to the CME input speed; the error bars show the ensemble prediction spread. This figure shows that is usually overpredicted when CME input speeds are above 1000 km/s. This bias is also apparent in the predictions made from a sample of 70 regular WSA-ENLIL+Cone runs reported by \inlineciteromano2013. The overprediction is most likely due to an overestimation of the CME dynamic pressure at Earth by the WSA-ENLIL+Cone model, due to the CME having a lower magnetic pressure than is observed in typical magnetic clouds. Also since the CME dynamic pressure is linearly related to the density and the square of the velocity, this quantity will be in particular more sensitive to higher CME input speeds, and produce higher in-situ speeds than those measured. Another factor in the higher CME dynamic pressure can arise from the approximation of the CME as a cloud with homogeneous density in the model.
|CME Onset||Binned Probabilistic Forecast (%)||Mean||Obs.||(nT)|
|Date Time (UT)||1||2||3||4||5||6||7||8||9||& spread||SSC||min.|
Other factors contributing to overprediction may include the magnetic field direction - two out of the three field configurations assumed produce persistent southward fields (135 and 180), so there is a bias towards geoeffective field configurations. Examining the distribution of north-south magnetic fields associated with the ICMEs of \inlineciterichardson2010 and the associated sheaths, in only 2% of cases are southward fields completely absent, so the bias towards geoeffective field configurations is consistent with observations. However, both small and large maximum southward fields are observed relatively infrequently (e.g., maximum southward fields are 4 nT in 17% of events, and 15 nT in 16%), suggesting that the weighting of 90 and 180 clock angles should be reduced. In particular, reducing the 180 clock angle weight would be expected to reduce the Kp overprediction.
The last line of Table 3 lists Brier Score calculated for each bin. Here, the is a measure of the magnitude of error in the probability forecast (how likely a given bin will occur) in each bin. The values indicate that in this sample, the probability forecast is reliable for the =5 and 6 bins (=0.17 for both), and less so for the =3 and 4 bins (=0.27 and 0.19). Although the scores also indicate that the forecast is most reliable for the smallest and largest bins, most of the observations in this sample did not fall in these extreme bins, hence a larger sample is needed to verify forecast reliability for these bins. Figure 10a shows the overall observed distribution and the forecast probability distribution for the events in Table 3 used to calculate the .
To further evaluate probability forecast reliability we compare the observed to the spread in ensemble predictions. For most (12 out of 17) of the ensembles, the observed was within the overall predicted spread (column 11). The observed was also within the predicted mean 1 for 11 out of 17 of the ensembles. A rank histogram was also constructed for the predictions for all ensembles and is shown in Figure 10b, again normalized to an ensemble size of 9. To construct this rank histogram the predictions are sorted from smallest to largest and the rank of where the observed value falls among the predicted values is noted. For example, an ensemble with a rank of 6 has the meaning that 5 predictions are less than the observed value, a rank of 10 would mean that all 9 of the predictions are less than the observed (underprediction), and a rank of 1 means that the observed value is less than all of the predictions (overprediction). The histogram has an overall flat shape, with more occurrence at rank 1 (the observed was less than the predicted range) and less occurrence in the higher ranks which shows the bias for overprediction (mean error=+0.6). Note, that the rank histogram does not indicate how “good” forecasts are but only measures whether the observed probability distribution is well represented by the ensemble. Therefore, a uniform, flat rank histogram is a necessary but not sufficient condition for determining the reliability of ensembles Hamill (2001).
7 ENLIL Parameter Sensitivity: 11 April 2013 Event Case Study
In the current configuration, other than the measured CME speed, direction, and size, the real-time WSA-ENLIL+Cone ensemble simulations use the default values for the model CME free parameters. In this section, we present a case study which examines the effect of changing these model free parameters on the ensemble modeling. The CME starting on 11 April 2013 at 07:24 UT was chosen for this study due to the large early arrival time prediction error obtained for all members of the model ensemble. \inlinecitetaktak2010 studied the dependence of arrival time predictions on the uncertainty in CME input parameters (speed, width, density ratio) for three Earth directed CME events of varying speeds. A similar procedure was adopted for this case study, and by employing the ensemble modeling technique, the parameter space can be sampled more systematically.
The original set of simulations performed in real-time were chosen as the “base ensemble” (ensemble I). Subsequently, ten ensemble runs (ensembles II-XI), each containing 36 members for 360 total simulations, were performed to assess the sensitivity of the CME arrival time prediction to changes in the model free parameters and ambient solar wind model, while keeping the CME speed and direction input parameters fixed. The ENLIL model free parameters considered in this study include the CME half width, CME density ratio, CME cavity ratio, and ambient solar wind reduction factor. The CME density ratio (dcld) is a free parameter which is by default is a set factor of 4 times larger than typical mean values in the ambient fast wind providing a pressure of four times larger than that of the ambient fast wind. The cavity ratio radcav is defined as the ratio of the radial CME cavity width to the CME width, with the default being no cavity radcav=0. The ambient speed reduction factor vred reduces the solar wind speed provided by the WSA coronal map in order to account for expansion of the solar wind from the WSA boundary to 1 AU since WSA is calibrated against 1 AU in-situ observations.
Figure 12 shows the CME starting on 11 April 2013 at 07:24 UT as viewed from SOHO LASCO C3, STEREO A and B SECCHI/COR2 near the time of 09:55 UT. On this date the STEREO B spacecraft was located at -142 and STEREO A was at 133 in HEEQ coordinates. This CME was associated with an M6.5 class flare from AR 11719 located at N07E13 with peak intensity at 07:16 UT. The eruption, coronal dimming and wave were visible mostly southeast of the active region in SDO/AIA 193Å. Additionally, an increase in solar energetic particle proton flux was observed starting at around 07:40 UT by the SOHO COSTEP (reaching 1 pfu/MeV, in the 16-40 MeV energy range), ACE EPAM (100 pfu/MeV, 1.22-4.94 MeV), and GOES-13 EPEAD (5 pfu/MeV, 15-40 MeV energy range) instruments starting at 08:00 UT, and by the IMPACT HET instruments on STEREO B (5 pfu/MeV, 24-41 MeV energy range) and A (0.001 pfu/MeV, 24-41 MeV energy range). This solar energetic particle event and its longitudinal extent is studied in detail by \inlinecitecohen2014 and \inlinecitelario2014.
Due to the lack of availability of real-time concurrent coronagraph images, triangulation of CME parameters with the StereoCAT “ensemble” mode method was not possible for this CME. Therefore, the ensemble was composed of “custom” members. The CME parameters for each member were derived from plane-of-sky CME speed measurements combined with the source location at the Sun. The distribution CME input parameters for 36 ensemble members are shown in Figure 13. Median CME parameters are: speed of 1000 km/s, direction of -15 longitude, 0 latitude, and a half-width of 55. This figure shows that custom ensemble members were chosen with speeds of 850, 900, 1000, 1100, and 1200 km/s, between 10 latitude, -10 to -25 longitude with a half width of 55. Subsequent re-analysis of the CME height-time evolution gives average plane-of-sky speeds of 800 km/s and 700 km/s for SECCHI COR2B and LASCO C3 respectively, yielding a triangulated speed of 850200 km/s, -55 latitude, -1510 longitude, 505 half width, which is represented within the ensemble members derived in real-time.
The WSA-ENLIL+Cone model scaled velocity contour plot is shown in Figure 14 on 13 April at 06:00 UT for the ensemble member with median CME input parameters. This simulation figure shows a nearly direct CME impact at Earth, slightly eastward. Figure 15 shows the base ensemble WSA-ENLIL+Cone modeled quantities for all 36 ensemble members (color traces) at Earth along with in-situ ACE (black) and Wind (grey) observations (when there are ACE datagaps). The model traces are color coded by CME input speed such that slow to faster input speeds are colored from light green to dark blue. All 36 of the ensemble members predicted that the CME would arrive (100%) and the mean predicted arrival at Earth was 13 April 06:14 UT (range from 13 April 00:47 to 12:20 UT). The histogram of the distribution of arrival times is shown in Figure 16. The clustering of predicted arrival times in this histogram (and also in Figure 15) reflects the limited number of discrete CME input speeds represented in the ensemble (see Figure 13), with faster CMEs arriving first. The CME-associated shock was observed to arrive at ACE and Wind on 13 April at 22:13 UT, giving an average prediction error of -16 hours. Clear ICME signatures including enhanced low variability magnetic field, declining solar wind speed, and low proton temperatures, start at around 16:45 UT on 14 April through about 18:30 UT on 15 April. The overall spread in arrival time predictions of all of the members in the base ensemble (including the clustering by CME input speed) can also be seen in Figure 15 as the color traces increase ahead of the observed arrival. The traces also show that the velocity, density, and temperature are overpredicted, while the maximum magnetic field strength is similar to that actually observed. The passage of this CME did not produce a geomagnetic storm due to an almost persistently northward magnetic field, shown in red in the top panel of Figure 15. The NOAA real-time observed index reached 3 during the synoptic period 21-24:00 UT on 13 April, while the Potsdam final was 3+. The index shows a sudden storm commencement of +21 nT at 23:00 UT on 13 April, and reached a minimum of only -7 nT at 11:00 UT on 15 April.
In Figure 17 the arrival time prediction error () for the members in the base ensemble is plotted against the CME input speed for different CME input propagation directions (gray scale coded) and a fixed half width of 55 (full angular width of 110). On 11 April 2013 Earth was located at -5.9 latitude and 0 longitude in HEEQ coordinates, thus the input propagation direction of -10 latitude and -10 longitude (black, and dark blue in subsequent figures) represent the members with the most direct impact. This Figure shows that the arrival time prediction error ranges from -9.9 hours to -21.4 hours and increases with initial CME speed.
Considering that a source of the prediction error could be due to uncertainty in the CME width, an identical ensemble (II) simulation was performed with the same input conditions but decreasing the half width by 10 (full angular width decreased from 110 to 90). Figure 18 shows the difference from the original predicted arrival times from the base ensemble against the CME input speed for different propagation directions (as shown in Figure 17) when the full angular width decreased from 110 to 90. In this figure (and those subsequent), the new CME arrival time prediction error is shown in hours relative to the original “base ensemble” prediction error. Compared to the original arrival time estimates, the overall prediction error decreases by 0.2 to 1.8 hours with increasing initial CME speed. Since all of the predictions in the base ensemble were too early (negative prediction error), a decrease in prediction error means that the new predictions are shifted to later times, closer to the observed arrival time. Nevertheless, the improvement is small compared to the prediction error.
Next, the dependence of prediction on the input CME density ratio dcld was considered. Two ensembles were performed (III and IV) for which all parameters of the base ensemble were held fixed but the CME density ratio was adjusted from four (default), to two, and three. The results are shown in Figure 19 which shows the difference from the predicted arrival time for the base ensemble as a function of CME speed for the two different density ratios. The prediction error decreases by 3.3 to 4.3 hrs for a CME density ratio of two and by 1.3 to 1.7 hrs for a density ratio of three, as a function of increasing initial CME speed. Hence, reducing the density ratio from four to two or three improves the arrival time prediction by around 3.5 or 1.5 hours, respectively.
Another ENLIL CME parameter is the cavity ratio which allows the CME to be represented by spherical shell of plasma, based on coronagraph observations of CME cavities. The cavity ratio radcav is defined as the ratio of the radial CME cavity width to the CME width, with the default being no cavity radcav=0. Figure 20 shows the results of five ensembles (V-IX) with the CME cavity ratio adjusted to 0.1, 0.3, 0.5. 0.6, and 0.7, i.e., the CME is modeled as a progressively thinner shell as the ratio increases, using the base ensemble for all other parameters fixed. Specifically, the differences from the arrival times obtained for the base ensemble are plotted as a function of CME speed and direction (indicated by the symbol/line type) for each of these ensembles (indicated by the line color). For a cavity ratio of 0.1 the prediction remains largely unchanged compared to the base ensemble (with 0.15 hours). For the other cavity ratios, increasing differences in Figure 20 as the cavity ratio increases correspond to the predicted arrival moving to later times, reducing the prediction error. Furthermore, for each cavity ratio there is a spread (1 to 3 hours) in prediction time difference (compared to the base ensemble) for the different CME input directions with the more Earth directed inputs showing the largest difference from the base ensemble. Overall the prediction error decreases by between 0-1.6 hrs, 0.9-3.9 hrs, 1.8-4.7 hrs, 2.4-5.6 hrs for cavity ratios of 0.3, 0.5, 0.5, 0.6 respectively.
Considering now the influence of the ENLIL ambient solar wind solution, the in-situ data-model comparison (Figure 15) for the base ensemble indicates that the modeled background solar wind speed is 125 km/s higher than the observed in-situ values, whereas the default value of ambient solar wind reduction factor vred is 25 km/s . To examine the role of the speed reduction factor in the prediction error, vred was increased to 50 km/s and 75 km/s for two ensembles (X and XI). This factor reduces the speed provided by the WSA coronal map in order to account for expansion of the solar wind from the WSA boundary to 1AU since WSA is calibrated against 1 AU in-situ observations. Figure 21 shows the prediction time difference from the base ensemble for these two ensembles which show differences of 1-1.3 hrs and 2.2-2.8 hrs for vred of 50 km/s and 75 km/s respectively. Since the differences are positive, this indicates that the predicted arrival times are moved later, reducing the error relative to the observed arrival time. As might be expected, the modeled CME propagates more slowly when the ambient solar wind is slower. Figure 22 illustrates how the modeled background solar wind speed better matches the observed speed prior to CME arrival when the ambient speed reduction factor vred is increased to 75 km/s in ensemble XI.
Overall, this parametric case study shows that after the CME input speed, the cavity ratio and density ratio assumed in ENLIL have the greatest effects on the predicted CME arrival time, each changing this time by about 3 hours on average. Their effect is also more noticeable with higher CME input speeds. The CME input speed, cavity and density ratios define the CME momentum which defines the CME deceleration. In the addition to the using the default values, new ensemble runs could be performed with changes to the CME cavity ratio and density ratio as informed by coronagraph measurements of the CME. Here, we have only examined the effect of changing the ad hoc ambient speed reduction factor in ENLIL, but we could also produce an ensemble of ambient solar wind WSA-ENLIL simulations using different ambient speed reduction factors which can be compared to in-situ measurements to determine “best” factor to use in subsequent CME simulations. An ensemble forecast reflecting uncertainties in the background solar wind could also be produced by using a variety of magnetograms (from different observatories or processed using different techniques) as input to the WSA or WSA-ADAPT models.
8 Summary and Discussion
This study evaluates the first ensemble CME prediction system of its kind employed in a real-time environment, providing unique space weather information for NASA users. The ensemble prediction approach provides a probabilistic forecast which includes an estimation of arrival time uncertainty from the spread in predictions and a forecast confidence in the likelihood of CME arrival. The current implementation explores the sensitivity of CME arrival time predictions from the WSA-ENLIL+Cone model to initial CME parameters. First results give a mean absolute arrival time error of 12.3 hours, RMSE of 13.9 hours, and mean error of -5.8 hours (early bias), based on a sample of 30 CME events for which ensemble simulations were performed. The arrival time is generally based on the arrival of the CME-generated shock at the Earth. The ensemble mean absolute error and RMSE are both comparable with other CME arrival time prediction errors reported in the literature.
When considering the overall performance of CME arrival prediction, it was found that the correct rejection rate is 62%, the false-alarm rate is 38%, correct alarm ratio is 77%, and false alarm ratio is 23%. Each ensemble CME arrival time forecast includes a forecast probability , which conveys a forecast uncertainty about the likelihood that the CME will arrive, which can be compared with observations to determine forecast reliability. The Brier Score () of 0.15 for all 30 ensemble CME arrival probabilities indicates that in this sample, on average, the predicted probability of the CME arriving is fairly accurate. (A of 0 on a range of 0 to 1 is a perfect forecast.) However, the reliability diagram (Figure 9a) shows that the ensemble simulations are underforecasting the likelihood that the CME will arrive in the forecast bins between 20-80%, and slightly overforecasting in the 1-20% and 80-100% forecast bins. Overforecasting is when the forecast chance of CME arrival is higher than is actually observed; i.e., the CME is observed to arrive less often than is predicted. More ensemble simulations are needed for a more robust forecast verification of these probabilistic CME arrival time forecasts.
For 8 out of 17 of the ensemble runs containing hits, the observed CME arrival was within the spread of ensemble arrival time predictions. The initial distribution of CME input parameters was shown to be an important influence on the accuracy of CME arrival time predictions. Particularly, the median and spread of the input distribution should accurately represent the range of CME parameters derived from observations. This is evidenced by the rank histogram (Figure 9b) which illustrates that roughly half of the observed arrivals are outside the spread of predictions, and also suggests undervariability in initial conditions; i.e., these ensembles do not sample a wide enough spread in CME input parameters.
Each set of ensemble simulations also provides a probabilistic forecast for each bin which can be compared with observations to determine forecast reliability. The Brier Score () for the probabilistic forecast bins show reliability for the =5 and 6 bins (=0.17 for both), and less so for the =3 and 4 bins (=0.27 and 0.19). If choosing a single categorical forecast value, the mean predicted was found to have smaller prediction errors compared to using the bin with the highest likelihood from the probabilistic forecast. The observed was within 1 of the predicted mean for 11 out of 17 of the ensembles. The prediction errors computed from the mean predicted show a mean absolute error of 1.4, RMSE of 1.8, and mean error +0.4. There is a known overall tendency for the overprediction of , generally found for CME input speeds above 800-1000 km/s. Again, more ensemble simulations are needed to provide better forecast verification and to calibrate the forecast.
This paper focuses on the forecast verification of the ensemble modeling aspect of CME arrival and predictions. More events, as well as comparison of results using different CME propagation models, are needed for more comprehensive forecast verification. These aspects are being investigated in a separate verification study which evaluates 400 single WSA-ENLIL+Cone simulations (of which there are 70 simulations containing CME arrivals) performed at the CCMC/SWRC.
The parameter sensitivity studied discussed in Section 7 suggests future directions for this ensemble system. In the addition to the using the default model values, new ensemble runs could be performed with changes to the CME cavity ratio and density ratio as informed by coronagraph measurements of the CME. As discussed in Section 2 an accurate representation of the background solar wind is necessary for simulating transients, and prediction errors arising from background characterization and other model limitations should be considered. An ensemble forecast reflecting uncertainties in the background solar wind could be produced by using a variety of magnetograms (from different observatories or processed using different techniques) as input to the WSA or WSA-ADAPT models. From these results one can produce an ensemble of ambient solar wind WSA-ENLIL model outputs which can be compared to in-situ measurements to determine “best” coronal maps/model instance. These sub-selected WSA or WSA-ADAPT maps could then be used for a series of ensemble WSA-ENLIL+Cone CME simulations. Such an improved ensemble forecast would produce predictions which also reflect uncertainties in the WSA-ENLIL modeled background solar wind in addition to the uncertainties in CME input parameters (as considered in this work).
Another improvement could be the use of real-time interplanetary scintillation (IPS) observations by the Ooty Radio Telescope Manoharan (2006). These data can provide crucial information about the CME propagation and interaction with the surrounding solar wind which could be used to provide updated information on CME parameters as the CME moves out from the Sun. This information could then be used to refine model predictions of the propagation of the CME. The STEREO Heliospheric Imagers also provide CME propagation information out to 1AU. However, it is not always possible to extract this information from real-time data, and the imagers do not always have an optimal viewing angle for Earth-directed CMEs. Comparisons of CME propagation from WSA-ENLIL with near real-time observations of the CME location inferred from IPS, the STEREO heliospheric imagers, or some other source, can be used to select ensemble members with the best agreement using quantitative and visual inspection employing advanced visualization techniques such as “3D volumetric rendering” Bock et al. (2014).
Finally, the forecasting of CME arrival would benefit from the use of other propagation models, in addition to WSA-ENLIL, each with its own set of independently assessed input parameters, leading to a community-wide ensemble prediction capability. A first step to such a capability is provided by the CME Scoreboard, described in Section 1, where anyone is invited to post their estimate of the arrival time of a recently observed CME in real-time.
The work was carried out as a part of NASA’s Game Changing Development Program Advanced Radiation Protection Integrated Solar Energetic Proton (ISEP) project. LKJ acknowledges the support of NSF grants AGS 1242798 and 1321493. MLM thanks T. Nieves-Chinchilla and B.J. Thompson for useful discussions. We gratefully acknowledge the participants of the CME Arrival Time Scoreboard http://kauai.ccmc.gsfc.nasa.gov/CMEscoreboard. The ACE and Wind solar wind plasma and magnetic field data were obtained at NASA’s CDAWeb (http://cdaweb.gsfc.nasa.gov). OMNI data was obtained from NASA’s COHOWeb (http://omniweb.gsfc.nasa.gov/coho). The geomagnetic index was obtained from the World Data Center for Geomagnetism in Kyoto, Japan. Estimated real-time planetary indices are from NOAA and the NGDC, and final definitive indices are from the Helmholtz Center Potsdam GFZ German Research Centre for Geosciences. The SOHO LASCO CME catalog is generated and maintained at the CDAW Data Center by NASA and the Catholic University of America in cooperation with the Naval Research Laboratory. SOHO is a mission of international cooperation between the European Space Agency and NASA. The STEREO/SECCHI data are produced by an international consortium of the NRL, LMSAL and NASA GSFC (USA), RAL and University of Birmingham (UK), MPS (Germany), CSL (Belgium), IOTA and IAS (France). Some figure colors based on ColorBrewer.org.
- Anderson (1996) Anderson, J.L.: 1996, A Method for Producing and Evaluating Probabilistic Forecasts from Ensemble Model Integrations. Journal of Climate 9, 1518 – 1530.
- Arge and Pizzo (2000) Arge, C.N., Pizzo, V.J.: 2000, Improvement in the prediction of solar wind conditions using near-real time solar magnetic field updates. J. Geophys. Res. 105, 10465 – 10480. doi:10.1029/1999JA000262.
- Arge et al. (2004) Arge, C.N., Luhmann, J.G., Odstrčil, D., Schrijver, C.J., Li, Y.: 2004, Stream structure and coronal sources of the solar wind during the May 12th, 1997 CME. Journal of Atmospheric and Solar-Terrestrial Physics 66, 1295 – 1309. doi:10.1016/j.jastp.2004.03.018.
- Arge et al. (2010) Arge, C.N., Henney, C.J., Koller, J., Compeau, C.R., Young, S., MacKenzie, D., Fay, A., Harvey, J.W.: 2010, Air Force Data Assimilative Photospheric Flux Transport (ADAPT) Model. Twelfth International Solar Wind Conference 1216, 343 – 346. doi:10.1063/1.3395870.
- Bartels, Heck, and Johnston (1939) Bartels, J., Heck, N.H., Johnston, H.F.: 1939, The three-hour-range index measuring geomagnetic activity. Terrestrial Magnetism and Atmospheric Electricity (Journal of Geophysical Research) 44, 411. doi:10.1029/TE044i004p00411.
- Bock et al. (2014) Bock, A., Mays, M.L., Rastaetter, L., Ynnerman, A., Ropinski, T.: 2014, VCMass: A Framework for Verification of Coronal Mass Ejection Ensemble Simulations. IEEE Scientific Visualization Conference Abstracts.
- Brier (1950) Brier, G.W.: 1950, Verification of Forecasts Expressed in Terms of Probability. Monthly Weather Review 78, 1.
- Brueckner et al. (1995) Brueckner, G.E., Howard, R.A., Koomen, M.J., Korendyke, C.M., Michels, D.J., Moses, J.D., et al.: 1995,. Solar Phys. 162, 357 – 402.
- Cohen et al. (2014) Cohen, C.M.S., Mason, G.M., Mewaldt, R.A., Wiedenbeck, M.E.: 2014, The Longitudinal Dependence of Heavy-ion Composition in the 2013 April 11 Solar Energetic Particle Event. Astrophys. J. 793, 35. doi:10.1088/0004-637X/793/1/35.
- Colaninno, Vourlidas, and Wu (2013) Colaninno, R.C., Vourlidas, A., Wu, C.C.: 2013, Quantitative comparison of methods for predicting the arrival of coronal mass ejections at Earth based on multiview imaging. Journal of Geophysical Research (Space Physics) 118, 6866 – 6879. doi:10.1002/2013JA019205.
- Davies et al. (2013) Davies, J.A., Perry, C.H., Trines, R.M.G.M., Harrison, R.A., Lugaz, N., Möstl, C., Liu, Y.D., Steed, K.: 2013, Establishing a Stereoscopic Technique for Determining the Kinematic Properties of Solar Wind Transients based on a Generalized Self-similarly Expanding Circular Geometry. Astrophys. J. 777, 167. doi:10.1088/0004-637X/777/2/167.
- Domingo, Fleck, and Poland (1995) Domingo, V., Fleck, B., Poland, A.I.: 1995, The SOHO Mission: an Overview. Solar Phys. 162, 1 – 37.
- Dryer (1974) Dryer, M.: 1974, Interplanetary Shock Waves Generated by Solar Flares. Space Sci. Rev. 51, 403 – 468.
- Dryer et al. (2001) Dryer, M., Fry, C.D., Sun, W., Deehr, C., Smith, Z., Akasofu, S.-I., Andrews, M.D.: 2001, Prediction in Real Time of the 2000 July 14 Heliospheric Shock Wave and its Companions During the ‘Bastille’ Epoch. Solar Phys. 204, 265 – 284. doi:10.1023/A:1014200719867.
- Emmons et al. (2013) Emmons, D., Acebal, A., Pulkkinen, A., Taktakishvili, A., MacNeice, P., Odstrčil, D.: 2013, Ensemble forecasting of coronal mass ejections using the WSA-ENLIL with CONED Model. Space Weather 11, 95 – 106. doi:10.1002/swe.20019.
- Fry et al. (2003) Fry, C.D., Dryer, M., Smith, Z., Sun, W., Deehr, C.S., Akasofu, S.-I.: 2003, Forecasting solar wind structures and shock arrival times using an ensemble of models. Journal of Geophysical Research (Space Physics) 108, 1070. doi:10.1029/2002JA009474.
- Gopalswamy et al. (2009) Gopalswamy, N., Yashiro, S., Michalek, G., Stenborg, G., Vourlidas, A., Freeland, S., Howard, R.: 2009, The SOHO/LASCO CME Catalog. Earth Moon and Planets 104, 295 – 313. doi:10.1007/s11038-008-9282-7.
- Hamill (2001) Hamill, T.M.: 2001, Interpretation of Rank Histograms for Verifying Ensemble Forecasts. Monthly Weather Review 129, 550.
- Hamill and Colucci (1997) Hamill, T.M., Colucci, S.J.: 1997, Verification of Eta RSM Short-Range Ensemble Forecasts. Monthly Weather Review 125, 1312.
- Harvey et al. (1996) Harvey, J.W., Hill, F., Hubbard, R.P., Kennedy, J.R., Leibacher, J.W., Pintar, J.A., Gilman, P.A., Noyes, R.W., Title, A.M., Toomre, J., Ulrich, R.K., Bhatnagar, A., Kennewell, J.A., Marquette, W., Patron, J., Saa, O., Yasukawa, E.: 1996, The Global Oscillation Network Group (GONG) Project. Science 272, 1284 – 1286. doi:10.1126/science.272.5266.1284.
- Harvey et al. (1992) Harvey, L.O., Hammond, K.R., Lusk, C.M., Mross, E.F.: 1992, The Application of Signal Detection Theory to Weather Forecasting Behavior. Monthly Weather Review 120, 863.
- Henney et al. (2012) Henney, C.J., Toussaint, W.A., White, S.M., Arge, C.N.: 2012, Forecasting F with solar magnetic flux transport modeling. Space Weather 10, 2011. doi:10.1029/2011SW000748.
- Hidalgo et al. (2000) Hidalgo, M.A., Cid, C., Medina, J., Viñas, A.F.: 2000, A new model for the topology of magnetic clouds in the solar wind. Solar Phys. 194, 165 – 174. doi:10.1023/A:1005206107017.
- Howard et al. (2008) Howard, R.A., Moses, J.D., Vourlidas, A., Newmark, J.S., Socker, D.G., Plunkett, S.P., et al.: 2008, Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI). Space Science Reviews 136, 67 – 115. doi:10.1007/s11214-008-9341-4.
- Howard and DeForest (2012) Howard, T.A., DeForest, C.E.: 2012, The Thomson Surface. I. Reality and Myth. Astrophys. J. 752, 130. doi:10.1088/0004-637X/752/2/130.
- Jackson et al. (2011) Jackson, B.V., Hick, P.P., Buffington, A., Bisi, M.M., Clover, J.M., Tokumaru, M., Kojima, M., Fujiki, K.: 2011, Three-dimensional reconstruction of heliospheric structure using iterative tomography: A review. Journal of Atmospheric and Solar-Terrestrial Physics 73, 1214 – 1227. doi:10.1016/j.jastp.2010.10.007.
- Jian et al. (2011) Jian, L.K., Russell, C.T., Luhmann, J.G., MacNeice, P.J., Odstrčil, D., Riley, P., Linker, J.A., Skoug, R.M., Steinberg, J.T.: 2011, Comparison of Observations at ACE and Ulysses with Enlil Model Results: Stream Interaction Regions During Carrington Rotations 2016 - 2018. Solar Phys. 273, 179 – 203. doi:10.1007/s11207-011-9858-7.
- Jolliffe and Stephenson (2011) Jolliffe, I.T., Stephenson, D.B. (eds.): 2011, Forecast Verification: A Practioner’s Guide in Atmospheric Science, 2nd edn. Wiley, New Jersey, USA.
- Kaiser et al. (2008) Kaiser, M.L., Kucera, T.A., Davila, J.M., St. Cyr, O.C., Guhathakurta, M., Christian, E.: 2008, The STEREO Mission: An Introduction. Space Science Reviews 136, 5 – 16. doi:10.1007/s11214-007-9277-0.
- Lario et al. (2014) Lario, D., Raouafi, N.E., Kwon, R.-Y., Zhang, J., Gómez-Herrero, R., Dresing, N., Riley, P.: 2014, The Solar Energetic Particle Event on 2013 April 11: An Investigation of its Solar Origin and Longitudinal Spread. Astrophys. J. 797, 8. doi:10.1088/0004-637X/797/1/8.
- LaSota (2013) LaSota, J.A.: 2013, STEREO Analysis. Undergraduate Honors Thesis, University of Alaska Fairbanks.
- Lee et al. (2013) Lee, C.O., Arge, C.N., Odstrčil, D., Millward, G., Pizzo, V., Quinn, J.M., Henney, C.J.: 2013, Ensemble Modeling of CME Propagation. Solar Phys. 285, 349 – 368. doi:10.1007/s11207-012-9980-1.
- Liu et al. (2010) Liu, Y., Davies, J.A., Luhmann, J.G., Vourlidas, A., Bale, S.D., Lin, R.P.: 2010, Geometric Triangulation of Imaging Observations to Track Coronal Mass Ejections Continuously Out to 1 AU. Astrophys. J. Lett. 710, L82 – L87. doi:10.1088/2041-8205/710/1/L82.
- Lugaz et al. (2010) Lugaz, N., Hernandez-Charpak, J.N., Roussev, I.I., Davis, C.J., Vourlidas, A., Davies, J.A.: 2010, Determining the Azimuthal Properties of Coronal Mass Ejections from Multi-Spacecraft Remote-Sensing Observations with STEREO SECCHI. Astrophys. J. 715, 493 – 499. doi:10.1088/0004-637X/715/1/493.
- MacNeice (2009) MacNeice, P.: 2009, Validation of community models: 2. Development of a baseline using the Wang-Sheeley-Arge model. Space Weather 7, 12002. doi:10.1029/2009SW000489.
- Manoharan (2006) Manoharan, P.K.: 2006, Evolution of Coronal Mass Ejections in the Inner Heliosphere: A Study Using White-Light and Scintillation Images. Solar Phys. 235, 345 – 368. doi:10.1007/s11207-006-0100-y.
- Mays et al. (2014) Mays, M.L., Taktakishvili, A., Romano, M., MacNeice, P.J., Zheng, Y., Pulkkinen, A.A., Kuznetsova, M.M., Odstrčil, D.: 2014, Validation of Real-time Modeling of Coronal Mass Ejections Using the WSA-ENLIL+Cone Heliospheric Model. Space Weather. In preparation.
- McKenna-Lawlor et al. (2006) McKenna-Lawlor, S.M.P., Dryer, M., Kartalev, M.D., Smith, Z., Fry, C.D., Sun, W., Deehr, C.S., Kecskemety, K., Kudela, K.: 2006, Near real-time predictions of the arrival at Earth of flare-related shocks during Solar Cycle 23. Journal of Geophysical Research (Space Physics) 111, 11103. doi:10.1029/2005JA011162.
- Menvielle and Berthelier (1991) Menvielle, M., Berthelier, A.: 1991, The K-derived planetary indices - Description and availability. Reviews of Geophysics 29, 415 – 432. doi:10.1029/91RG00994.
- Millward et al. (2013) Millward, G., Biesecker, D., Pizzo, V., Koning, C.A.: 2013, An operational software tool for the analysis of coronagraph images: Determining CME parameters for input into the WSA-Enlil heliospheric model. Space Weather 11, 57 – 68. doi:10.1002/swe.20024.
- Müller et al. (2009) Müller, D., Dimitoglou, G., Caplins, B., Ireland, J., Wamsler, B., Hughitt, K., Agheksanterian, D. A. Amadigwe: 2009, Jhelioviewer - Visualizing large sets of solar data using JPEG 2000. Comput. Sci. Eng. 11 38.
- Murphy (1973) Murphy, A.H.: 1973, A New Vector Partition of the Probability Score. Journal of Applied Meteorology 12, 595 – 600.
- Newell et al. (2007) Newell, P.T., Sotirelis, T., Liou, K., Meng, C.-I., Rich, F.J.: 2007, A nearly universal solar wind-magnetosphere coupling function inferred from 10 magnetospheric state variables. Journal of Geophysical Research (Space Physics) 112, 1206. doi:10.1029/2006JA012015.
- Odstrčil (2003) Odstrčil, D.: 2003, Modeling 3-D solar wind structure. Advances in Space Research 32, 497 – 506. doi:10.1016/S0273-1177(03)00332-6.
- Odstrčil and Pizzo (1999a) Odstrčil, D., Pizzo, V.J.: 1999a, Three-dimensional propagation of CMEs in a structured solar wind flow: 1. CME launched within the streamer belt. J. Geophys. Res. 104, 483 – 492. doi:10.1029/1998JA900019.
- Odstrčil and Pizzo (1999b) Odstrčil, D., Pizzo, V.J.: 1999b, Three-dimensional propagation of coronal mass ejections in a structured solar wind flow 2. CME launched adjacent to the streamer belt. J. Geophys. Res. 104, 493 – 504. doi:10.1029/1998JA900038.
- Odstrčil, Riley, and Zhao (2004) Odstrčil, D., Riley, P., Zhao, X.P.: 2004, Numerical simulation of the 12 May 1997 interplanetary CME event. J. Geophys. Res. 109, 2116. doi:10.1029/2003JA010135.
- Odstrčil, Smith, and Dryer (1996) Odstrčil, D., Smith, Z., Dryer, M.: 1996, Distortion of the heliospheric plasma sheet by interplanetary shocks. Geophys. Res. Lett. 23, 2521 – 2524. doi:10.1029/96GL00159.
- Pizzo and Biesecker (2004) Pizzo, V.J., Biesecker, D.A.: 2004, Geometric localization of STEREO CMEs. Geophys. Res. Lett. 31, 21802. doi:10.1029/2004GL021141.
- Pulkkinen, Oates, and Taktakishvili (2010) Pulkkinen, A., Oates, T., Taktakishvili, A.: 2010, Automatic Determination of the Conic Coronal Mass Ejection Model Parameters. Solar Phys. 261, 115 – 126. doi:10.1007/s11207-009-9473-z.
- Pulkkinen et al. (2011) Pulkkinen, A.A., Taktakishvili, A., Odstrčil, D., MacNeice, P.J.: 2011, Ensemble forecasting of coronal mass ejection propagation in the interplanetary medium. NOAA Space Weather Workshop Abstracts.
- Richardson and Cane (2010) Richardson, I.G., Cane, H.V.: 2010, Near-Earth Interplanetary Coronal Mass Ejections During Solar Cycle 23 (1996 - 2009): Catalog and Summary of Properties. Solar Phys. 264, 189 – 237. doi:10.1007/s11207-010-9568-6.
- Riley, Linker, and Mikić (2001) Riley, P., Linker, J.A., Mikić, Z.: 2001, An empirically-driven global MHD model of the solar corona and inner heliosphere. J. Geophys. Res. 106, 15889 – 15902. doi:10.1029/2000JA000121.
- Romano et al. (2013) Romano, M., Mays, M.L., Taktakishvili, A., MacNeice, P.J., Zheng, Y., Pulkkinen, A.A., Kuznetsova, M.M., Odstrčil, D.: 2013, Validation of Real-time Modeling of Coronal Mass Ejections Using the WSA-ENLIL+Cone Heliospheric Model. AGU Fall Meeting Abstracts, A2156.
- Rostoker (1972) Rostoker, G.: 1972, Geomagnetic indices. Reviews of Geophysics and Space Physics 10, 935 – 950. doi:10.1029/RG010i004p00935.
- Sivillo, Ahlquist, and Toth (1997) Sivillo, J.K., Ahlquist, J.E., Toth, Z.: 1997, An Ensemble Forecasting Primer. Weather and Forecasting 12, 809 – 818.
- Smith and Dryer (1990) Smith, Z., Dryer, M.: 1990, Mhd study of temporal and spatial evolution of simulated interplanetary shocks in the ecliptic plane within 1 au. Solar Physics 129(2), 387 – 405. doi:10.1007/BF00159049. http://dx.doi.org/10.1007/BF00159049.
- Smith et al. (2009) Smith, Z.K., Dryer, M., McKenna-Lawlor, S.M.P., Fry, C.D., Deehr, C.S., Sun, W.: 2009, Operational validation of HAFv2’s predictions of interplanetary shock arrivals at Earth: Declining phase of Solar Cycle 23. Journal of Geophysical Research (Space Physics) 114, 5106. doi:10.1029/2008JA013836.
- Sugiura (1964) Sugiura, M.: 1964, Hourly values of equatorial Dst for the IGY. Ann. Int. Geophys. Year 35(9), 945.
- Taktakishvili, MacNeice, and Odstrčil (2010) Taktakishvili, A., MacNeice, P., Odstrčil, D.: 2010, Model uncertainties in predictions of arrival of coronal mass ejections at Earth orbit. Space Weather 8, 6007. doi:10.1029/2009SW000543.
- Talagrand, Vautard, and Strauss (1997) Talagrand, O., Vautard, R., Strauss, B.: 1997, Evaluation of probabilistic prediction systems. In: Proceedings, ECMWF Workshop on Predictability, 1 – 25.
- Thernisien, Howard, and Vourlidas (2006) Thernisien, A.F.R., Howard, R.A., Vourlidas, A.: 2006, Modeling of Flux Rope Coronal Mass Ejections. Astrophys. J. 652, 763 – 773. doi:10.1086/508254.
- Vršnak et al. (2014) Vršnak, B., Temmer, M., Žic, T., Taktakishvili, A., Dumbović, M., Möstl, C., Veronig, A.M., Mays, M.L., Odstrčil, D.: 2014, Heliospheric Propagation of Coronal Mass Ejections: Comparison of Numerical WSA-ENLIL+Cone Model and Analytical Drag-based Model. Astrophys. J. Supp. 213, 21. doi:10.1088/0067-0049/213/2/21.
- Weigel et al. (2006) Weigel, R.S., Detman, T., Rigler, E.J., Baker, D.N.: 2006, Decision theory and the analysis of rare event space weather forecasts. Space Weather 4, 5002. doi:10.1029/2005SW000157.
- Wilks (1995) Wilks, D.S.: 1995, Statistical Methods in Atmospheric Sciences: An Introduction, Academic Press, Massachusetts, USA.
- Xie, Ofman, and Lawrence (2004) Xie, H., Ofman, L., Lawrence, G.: 2004, Cone model for halo CMEs: Application to space weather forecasting. J. Geophys. Res. 109, 3109. doi:10.1029/2003JA010226.
- Yashiro et al. (2004) Yashiro, S., Gopalswamy, N., Michalek, G., St. Cyr, O.C., Plunkett, S.P., Rich, N.B., Howard, R.A.: 2004, A catalog of white light coronal mass ejections observed by the SOHO spacecraft. Journal of Geophysical Research (Space Physics) 109, 7105. doi:10.1029/2003JA010282.
- Zhao and Dryer (2014) Zhao, X., Dryer, M.: 2014, Current status of CME/shock arrival time prediction. Space Weather 12, 448 – 469. doi:10.1002/2014SW001060.
- Zhao, Plunkett, and Liu (2002) Zhao, X.P., Plunkett, S.P., Liu, W.: 2002, Determination of geometrical and kinematical properties of halo coronal mass ejections using the cone model. J. Geophys. Res. 107, 1223. doi:10.1029/2001JA009143.
- Zheng et al. (2013) Zheng, Y., Macneice, P., Odstrčil, D., Mays, M.L., Rastaetter, L., Pulkkinen, A., Taktakishvili, A., Hesse, M., Masha Kuznetsova, M., Lee, H., Chulaki, A.: 2013, Forecasting propagation and evolution of CMEs in an operational setting: What has been learned. Space Weather 11, 557 – 574. doi:10.1002/swe.20096.