Forecasting the 2018 February 12th CME propagation with the P-DBM model: a fast warning procedure

The forecast of the time of arrival of a Coronal Mass Ejection (CME) to Earth is of critical importance for our high-technology society and for the Earth’s upper atmosphere status and LEO satellites. We realized a procedure based on the Drag-Based Model which uses probability distributions, rather than exact values, as input parameters, and allows the evaluation of the uncertainty on the forecast. We tested this approach using a set of CMEs whose transit times are known, obtaining extremely promising results. We realized a real-time implementation of this algorithm which ingests the outputs of automated CME tracking algorithms as inputs to provide early warning for those CME approaching Earth. We present the results of this real-time fast warning procedure for the case of the 2018 February 12 th CME.


I. INTRODUCTION
Coronal Mass Ejections (CMEs) are expulsion of plasma and magnetic field from the Sun with repercussions on the whole interplanetary space. Most of our knowledge about CMEs comes from coronagraphic imaging of the early stages of the CME propagation in the heliosphere. Their Plane of Sky (PoS) speed is among their characterizing attributes that we are able to measure while they are coronograph Field of View. During the launch phase, CME are accelerated (from a fraction to several km/s 2 ) near the Sun and then move with constant speed or slowly decelerate due to the drag force exerted by the ambient medium. Their speeds when they leave the coronagraph Field of View (at distance of ~30 Rs) range from 100 to >3000 km/s [Gopalswamy, 2016]. The CMEs that have higher speeds are the most energetic phenomena relevant to space weather [Gopalswamy, 2018] and primary space weather consequences of those CMEs are large energetic particle fluxes in the near-Earth environment and geomagnetic storms [e.g.: Gonzalez et al., 1994;Tsurutani et al., 1988]. More or less intense storms depend on the Interplanetary Coronal Mass Ejection (ICME) characteristics, such as its velocity, size and the associated magnetic field [Gopalswamy, 2006]. Strong geomagnetic storms are more likely when the incoming ICME is fast and its magnetic field orientation is opposite with respect to the Earth's magnetic field, as this configuration favors magnetic reconnection and the transfer of energy from the ICME to the terrestrial magnetosphere [Gonzalez et al., 1999;Koskinen and Huttunen, 2006;Tsurutani et al., 2003]. An empirical relation is known to represent how the storm strength, measured by the Disturbance Storm Time (Dst) index, is determined by the CME impact velocity V and the strength of its magnetic field component oriented in the direction opposite to the Earth's horizontal magnetic field, Bz [Gopalswamy, 2010]: (1) Our technological society, which strongly relies on inorbit infrastructures, has become more vulnerable to these events. For this reason, several space-weather centers have been created to continuously monitor the solar activity and possibly forecast its outcomes on Earth. The accurate prediction of the arrival time and velocity of ICMEs at Earth has always been a primary goal of the space-weather forecasting [e.g.: Daglis, 2001;Schrijver and Siscoe, 2010], and the space-weather centers employ different methods to model and predict the propagation of ICMEs. Among those, the Space Weather Prediction Center of the National Oceanic and Atmospheric Administration 1 , SWPC/NOAA, in USA, and the Met Office Space Weather Operations Centre 2 in UK, employ the WSA-ENLIL+Cone model [Odstrcil et al., 2004], while the Galileo Ionosphere Prediction Service [Albanese et al., 2018] and the SWERTO 3 regional service [Berrilli et al., 2018] in Italy employs the Probabilistic-Drag Based Model [P- DBM -Napoletano et al., 2018]. Those two ICME propagation models differ in approach, assumptions, and complexity. The WSA-ENLIL+Cone model is a 3-D numerical MHDbased model. It starts with a time-dependent simulation of the background solar-wind in the heliosphere where the CME is inserted as a hydrodynamic structure with an initial cone-shaped morphology. To return an accurate description of the ICME propagation, WSA-ENLIL+Cone needs a good knowledge of the background solar wind and of several parameters of the CME. In practice, multiple-view observations from the SOHO and STEREO spacecrafts, and from the GONG network are used to derive the inputs to the model, with the model becoming less and less accurate when there is no sufficient information to reproduce both a suitable background solar wind description and an appropriate CME replica. The P-DBM belongs to a family of models that apply a somewhat simplified description of the main interaction 1https://www.swpc.noaa.gov/products/wsa-enlil-solar-wind-prediction 2https://www.metoffice.gov.uk/public/weather/space-weather/ 3http://spaceweather.roma2.infn.it/ the ICME is subject to during its interplanetary journey. The drag-based models (DBMs) assume that beyond a certain heliospheric distance, a simple aerodynamic drag equation [Vrsnak et al., 2007] can describe the interaction between the ICME and the solar wind environment. Namely, the CME acceleration a can be evaluated by considering a drag effect according to: where v is the CME velocity, w is the solar wind velocity and γ is the drag parameter that conceals the details about the CME shape, mass, and in general all the physics related to the effectiveness of the drag effect. The implementation of the various DBMs may differ in the description of the CME morphology or other details [Vrsnak et al., 2013, 2014, Zic et al., 2015, but, due to their simplicity, all run in seconds, and can be made publicly available as on-line tools 4 . Notwithstanding the major difference in complexity, the DBMs performances are very close to those of the WSA-ENLIL+Cone [Vrsnak et al., 2014]. This indicates that, at present, the major limitation to an accurate ICME propagation forecast comes from the lack of reliable inputs, both in the form of large measure errors and of totally unknown parameters. Taking advantage of the DBM computational lightness, one can try to compensate for this lack of information applying an ensemble modeling approach. Ensemble modeling incorporate the intrinsic limitation of information due to measure errors or due to the lack of measure at all, in form of probability distributions. In practice, instead of a single run to forecast an ICME propagation, a set of runs, driven with input parameters extracted from suitable distributions are used to retrieve a distribution of output parameters. A similar approach was used a few years ago by  to evaluate the sensitivity of WSA-ENLIL+Cone to initial CME parameters. Their results indicated how the the predictions were heavily dependent on the CME initial parameters and their associated spread. In particular, they reported about the necessity of expanding the parameter space, since their "ensembles did not sample a wide enough spread in CME input parameters". Alas, expanding the parameter space was hardly a choice, since the WSA-ENLIL+Cone simulations are computationally very intensive and consequently need a relatively long time to complete. Recently, two ensemble approaches for DBMs, the P-DBM and the Drag Based Ensemble Model [DBEM 5 -Dumbovic et al., 2018], were independently developed, with the capacity to run thousands of single simulations in seconds and thus exploring thoroughly the parameter space. The model outputs are the most probable ICME travel time and velocity at arrival at Earth, as well as the 4http://oh.geof.unizg.hr/DBM/dbm.php 5http://oh.geof.unizg.hr/DBEM/dbem.php associated prediction uncertainties. In this paper, we resume the basic assumptions of the P-DBM and the outcomes of its validation using a dataset of past CMEs. Then, we focus on its real-time implementation and show the result of an actual application for forecasting the arrival time and speed of a CME at Earth.

II. THE P-DBM MODEL
The P-DBM makes use of probability distributions, rather than exact values, as input parameters, to model the CME propagation in the interplanetary space, allowing us to forecast ICME arrival times and velocities and to evaluate also the uncertainty on the forecast. Details on the implementation and an extended description are in Napoletano et al. (2018). In the following, we briefly introduce some concept and present an example of the validation outputs.

The model
Considering constant the solar wind speed w and the drag parameter γ, eq. (2) can be solved explicitly to return both the ICME distance and speed as functions of time. Consequently, the DBM needs four quantities [r0, v0, w, γ] (r0 and v0 are the initial position and speed, respectively) to compute the ICME distance and speed at any time.
With the further simplification of a spherical CME front, we can compute the transit time and the velocity at 1AU, for any CME that has been launched Earthwards and whose [r0, v0, w, γ] are measured or known. However, these quantities are either affected by measure errors or simply not known. Whenever those errors are large or that knowledge is lacking, as is the case for real-time ICME forecast implementation, the incertitude on the results may be large. The P-DBM addresses this problem by returning a statistical evaluation of the time of arrival and velocity of the ICME at a chosen distance from the Sun, transforming the probability or error distribution functions associated to the input parameters into probability distribution functions (PDFs) for the output results. This allows computing the best estimates and the errors for both the time of arrival and the velocity of a given ICME. For each ICME whose r0 and v0 are measured with an associated uncertainty, we will generate N different [r0, v0, w, γ] initial condition sets, randomly chosen from the relative PDFs, to compute the transit time and the velocity at 1AU, for example. This process generates the PDFs associated to t1AU and v1AU, which can be used to estimate the ICME most probable time of arrival and velocity and their associated uncertainties at 1AU.

The validation
To perform a validation of the P-DBM, we employed the sample of CMEs from [Shi et al., 2015]. That work provides the travel times from the Sun to Earth, allow the estimation of the initial position and velocity and the associated errors for 21 CMEs, extracted from the larger sample of the GMU CME/ICME catalogue compiled by Phillip Hess and Jie Zhang (http://solar.gmu.edu/heliophysics/index.php/GMU_CME /ICME_List), spanning from 2007 to 2017 from SOHO/LASCO [Domingo et al., 1995;Brueckner et al., 1995] and STEREO/SECCHI [Kaiser et al., 2008;Howard et al., 2008] data. [Shi et al., 2015] selected only those CMEs whose front could be unambiguously identified. For each event, they fitted a Graduated Cylindrical Shell [Thernisien, 2011] model and employed the J map technique [Sheeley et al., 1999] to extract the CME launch parameters and ensuring relatively small errors in their estimation. Among these CMEs, the authors highlighted six that were likely to have undergone interactions with non-radial forces (e.g.: a possible CME-CME interactions) and were therefore not suitable to be described by the DBM. The CME arrival times were obtained from the in situ data by WIND [Acuña et al. 1995], however, in some cases, the travel time had to be corrected, since the authors gave the arrival time of the shock preceding the arrival of the ICME itself by a few hours. The initial positions and the velocities and their relative errors have been used to generate Gaussian PDFs for r0 and v0. The γ PDF has been modeled with a Log-Normal function (μ = 0.70 and σ = 1.01). The w PDF has been chosen considering the slow or fast solar wind condition measured at 1AU in the hours before the ICME arrival at Earth, using a Gaussian function centered at 400 km/s with σ = 33 km/s for the slow solar wind, and a Gaussian function centered at 600 km/s with σ = 66 km/s for the fast solar wind. From this data, we realized N=50000 different propagations for each ICME, we were able to compute the forecast travel times, the velocity at 1AU and their associated statistical ( ±1σ ) errors, of the ICMEs in the sample. In Fig. 1 we plot the comparison between the forecast travel times and the measured travel times as green dots, with the error bars to represent the error associated to our forecast travel times. In the plot, we shaded in light blue the perfect match zone, where the forecast and measured travel times are the same within the measure errors. For 10 out of 15 CMEs of this sample, the absolute difference between measure and forecast travel time is smaller than the forecast error. The average value of this quantity is 8.5±5 hours. Since, in principle, we may not know at launch time whether a CME would interact with another CME or not, in this plot we also show (red symbols) those six CMEs that were not suitable to be described by the P-DBM, and that were excluded from the whole sample by Shi et al. [2015]. Considering also those events, the average of the absolute residual times grows to 12±13 hours, largely dominated by those two CMEs that arrived with large advance with respect to the forecast.

III. THE 2018 FEBRUARY 12 TH CME
To highlight the capabilities of the P-DBM for real-time application, we present one case study, where we compared the forecast provided by the P-DBM fed by SOHO/LASCO data and a simple de-projection program with the measured ICME characteristics at its arrival at Earth. We stress that in the real-time implementation, only information from SDO and SOHO is used to estimate the CME characteristics. This leads to larger errors in the characteristics estimation than when using also information from STEREO, but it is a nice test of the operational capabilities of the P-DBM as forecast tool when STEREO data would not be available any more.
The CME detection A CME launch was detected on the SOHO/LASCO coronagraphic images on February 12th, 2018 (Fig. 2) The CME direction assessment To compute the whole 3D propagation vector of the CME, we have to complement the SOHO/LASCO information with a possible source on the solar surface. With the source on the Sun, the hypotheses of a propagation vector normal to the solar surface and of a conical shape for the CME, we can de-project the velocity and estimate the true width. By intersecting the CME cone with the ecliptic, we can estimate the probability of the CME hitting Earth (or any other target). To find the most probable CME source on the Sun, we developed a source finding algorithm that make use of the Heliophysics Events Knowledgebase [Hurlburt et al., 2012] to check which solar features, that may have caused the CME (e.g.: flares, filaments, etc.), are within the area compatible with the CME launch parameters. The algorithm returns the coordinate of the CME source and the Line of Sight (LOS) velocity of the CME. This algorithm also returns the position of the coronal holes at the CME onset, allowing us to use a fast (CME source close or within a coronal hole) or slow (CME source far from a coronal hole) solar wind PDF for the successive ICME propagation computation. Fig. 3 shows the graphical output of the CME source finding algorithm, visualizing the characteristics of the CME over an SDO/AIA [Lemen et al., 2012] image of the Sun: the green arrow shows the CME POS angle, the yellow sector shows the angle subtended by the CME POS width, the white boxes and the red dots shows the ARs and recent flare locations, respectively, taken into account to estimate the CME onset location (shown as a black circle). The gray polygons show the position of coronal holes at the time of the CME onset.
Heliographic longitude: 13°±1° Heliographic latitude: -1°±1° LOS velocity: 847±367 km/s Solar wind: slow The ICME propagation with P-DBM Using the estimate of r0 and its associated error from the CACTUS data, the estimate of v0 and its error from the de-projection procedure, and the solar wind state assessment, we can generate the input PDFs and compute the output PDFs at 1AU for the ICME under consideration. Also for this real-time implementation, we used N=50000 CME realizations. The PDFs of the computed t1AU and v1AU are shown in Figg. 4 and 5.
In particular, using the PDF of the travel times, we can estimate the most probable t1AU and compute the inner solar system planet positions at that moment. Fig. 6 shows the planet positions at the most probable arrival time at 1AU of the CME. Considering Gaussian errors on longitude of the CME origin site on the Sun and on the CME width, and the error in the travel time forecast, we define four zones of probability about the CME position and extension, which are represented by different shades of colors, from red (impact probability ~90%), to dark orange (impact probability ~40%), to light orange (impact probability ~10%), to yellow (impact probability ~5%). For this particular CME, the probability of hitting Earth was greater than 90%. This forecast was available within one hour from the CME detection by CACTUS.

The measured ICME characteristics
The ICME under investigation actually hit Earth. In Fig.  7, we plot several relevant plasma and magnetic field quantities measured by the Magnetometer and Faraday Cup instruments on-board the DSCOVR satellite [Szabo, 2012] at L1 from 06 to 12 UTC of the 2018 February 15th. From these data, there is no conclusive signature of the difference between the preceding shock and the actual ICME magnetic structure arrivals. Therefore, it has been agreed [Mays, 2018] that the CME Arrival Time-stamp at Earth is: 2018-02-15T07:38 UTC, by summing 1 hour of time for the CME to travel from the DSCOVR satellite position at L1 to Earth. It is worth to note that this particular CME had a large POS velocity uncertainty from the initial CACTUS estimate. This large error propagated in the algorithm, leading to a relatively large uncertainty on both t1AU and v1AU. This notwithstanding, the forecast arrival time differs by a mere 1.5h from the measured arrival time and the ICME velocity is partially compatible with the measured value. ICME arrival date: 2018 Feb 15 07:38 UTC ICME arrival velocity: ~350 km/s

IV. CONCLUSIONS
We developed a statistical approach for the computation of kinematic propagation of ICME, making use of PDF to cope with our lack of knowledge on CME or IMF parameters. We tested this approach using a set of CMEs whose transit times are known, and obtained extremely promising results. We presented a case study where this approach has been used to forecast in real-time the propagation and arrival at Earth of a single CME launched on February 12 th , 2018. The forecast was available within one hour from the CME detection and both the forecast arrival time and the velocity at Earth are compatible with the measured ICME values.   The 2018 February 12th CME as imaged by the LASCO C2 coronagraph on-board SOHO. The image shown is a contrast enhanced difference of two images acquired few minutes apart. This allows the visualization of the CME front few hours after the onset time. Fig. 3: Graphical output of the CME source finding algorithm. We visualize several characteristics of the CME over a SDO/AIA 0304 image of the Sun: the green arrow shows the CME POS angle, the yellow sector shows the angle subtended by the CME POS width, the white boxes and the red dots shows the ARs and recent flare locations, respectively, taken into account to estimate the CME onset location (shown as a black circle). The gray polygons show the position of coronal holes at the time of the CME onset.   : Graphical output of the CME propagation algorithm. The inner part of the solar system is drawn at the most probable arrival time at 1AU of the CME. The estimated CME propagation path is represented by dashed lines, and the CME itself is represented by different shades of colors, representing different confidence levels of the CME extension: yellow=10% probability, light orange=50% probability, dark orange=90% probability. Inner heliosphere planets are represented by colored symbols: Mercury= purple star, Venus= green triangle, Earth= blue circle; Mars= red pentagon. Fig. 7: Relevant plasma and magnetic field quantities measured by the DSCOVR satellite at L1 from the 6 to 12 UTC of the 2018 February 15 th . From top to bottom: the total strength (in blue) and the component perpendicular to the ecliptic (in orange) of the interplanetary magnetic field; the solar wind proton density; the solar wind proton velocity. The evident discontinuity (highlighted with the red vertical dashed line) at 06:38 UTC marks the arrival of the ICME structure at L1.