Extended Drag-Based Model for better predicting the evolution of Coronal Mass Ejections
Mattia Rossi Sabrina Guastavino Michele Piana Anna Maria Massone




empty 



Extended Drag-Based Model for better predicting the evolution of Coronal Mass Ejections



Mattia Rossi^1, Sabrina Guastavino^1,2, Michele Piana^1,2 and Anna Maria Massone^1 


^1MIDA, Dipartimento di Matematica Università di Genova, via Dodecaneso 35 16146 Genova, Italy

^2Osservatorio Astrofisico di Torino, Istituto Nazionale di Astrofisica, Strada Osservatorio 20 10025, Pino Torinese, Italy


empty 

§ ABSTRACT 
The solar wind drag-based model is a widely used framework for predicting the propagation of Coronal Mass Ejections (CMEs) through interplanetary space. This model primarily considers the aerodynamic drag exerted by the solar wind on CMEs. However, factors like magnetic forces, pressure gradients, and the internal dynamics within CMEs justify the need of introducing an additional small-scale acceleration term in the game. Indeed, by accounting for this extra acceleration, the extended drag-based model is shown to offer improved accuracy in describing the evolution of CMEs through the heliosphere and, in turn, in forecasting CME trajectories and arrival times at Earth. This enhancement is crucial for better predicting Space Weather events and mitigating their potential impacts on space-based and terrestrial technologies.


Keywords: coronal mass ejections – interplanetary propagation – drag-based model – accelerated dynamics – spacecraft alignment



§ INTRODUCTION 
 
Coronal Mass Ejections <cit.> are massive outbursts of magnetized plasma from the solar corona into the interplanetary space. When directed toward Earth, they cause severe geomagnetic disturbances <cit.> and can pose a persistent hazard as harmful radiation to space and ground-based facilities, and human health. Therefore, predicting the CMEs' arrival time and impact speed to Earth is essential in the context of the Space Weather forecasting science <cit.>.

One of the most popular and commonly used approaches to predict the transit time of a CME and its speed to Earth is known as the Drag-Based Model (DBM) <cit.>. This model assumes that the kinematics of the CME is governed by its dynamic interaction with the Parker spiral-shaped interplanetary structures (i.e., high- and low-speed streams) where it propagates, via the magnetohydrodynamic (MHD) equivalent of the aerodynamic drag force. The model, which mathematically reduces to a rather simple equation of motion, thus essentially predicts that the speed of the CME will balance that of the ambient solar wind in which it is expanding. Recent efforts have also been devoted to incorporating the physics of aerodynamic drag into methodologies based on Artificial Intelligence (AI) techniques, paving the way for innovative (hybrid) approaches known as physics-driven AI models <cit.>.

Although the DBM has been subject to continuous refinements <cit.> and is now a well-established approach, the several drawbacks associated with its intrinsic approximations are evident. Indeed, it is clear that the complex dynamical interaction of the CME with its surroundings cannot be properly described solely by the drag force. Other important physical processes are certainly at play in the evolution of CMEs: these include CME rotation, reconfiguration, deformation, deflection, erosion along with any other magnetic reconnection-driven processes <cit.>, resulting in additional accelerations beyond that predicted by the trivial DBM acting on the CME as it travels through the heliosphere.

<cit.> recently pointed out that the DBM is often quite ineffective in describing the proper propagation of CMEs in interplanetary space. In this study, a CME was observed by two radially aligned probes separated by a distance of just 0.13 AU. Although the model predicted that the CME would decelerate, the velocity profiles measured by the two spacecraft instead revealed a residual acceleration, pointing to an additional force to the drag that overpowered its braking effect, and thus resulting in an increase in velocity. This work presents a more refined and realistic drag-based model, with the aim to overcome the limitations of current versions by introducing into the equation of motion describing the dynamic interaction of the CME with the solar wind an extra acceleration, representing any other forces involved. 

After obtaining and discussing the mathematical solutions of the resulting new equations of motion (<ref>.), the updated version is applied to the observation of the same CME already studied in <cit.>, showing that it satisfactorily succeeds in describing its dynamic evolution, and thus becoming a significant breakthrough in the prediction of CME travel time in Space Weather studies (<ref>.). Interpretation of which physical process(es) the additional acceleration is due to, is tentatively given in <ref>., where our conclusions are also offered. Computational details to derive formulae in <ref>. are summarized in the Appendix.



§ THE EXTENDED DRAG-BASED MODEL
 

Let us consider a generalization of the DBM where the total net acceleration acting on the CME in the interplanetary phase is made of two contributions:

    r̈=a_drag(r,t)+a_extra(r) ,

where a_drag= -γ r^-α|ṙ - w(r,t)|^k(ṙ-w(r,t)), k∈ 2ℕ+1, and a_extra=ar^-β, for appropriate exponents α,β>0 and coefficients γ>0, a≠0; r=r(t) and v(t)=ṙ(t) are the CME’s instantaneous radial position and speed (typically the CME front distance and front speed); w(r,t) is the background solar wind speed given as a known function of position and time. Physically, the model describes the same
dynamics of the DBM perturbed by an extra (e.g., magneto-gravitational) force acting on the CME along the motion, altogether exponentially damped over distance.

In general, equation (<ref>) does not admit an analytical solution, which hampers the computation of a_extra from time and space measurements, i.e., by solving a boundary value problem. A closed-form time solution of (<ref>) is possible by assuming α=β=0 and a constant w(r,t)≡ w, for any fixed odd integer k. Therefore, we introduce the Extended Drag-Based Model (EDBM hereafter) as the equation

    r̈=-γ|ṙ-w|(ṙ-w)+a ,

in which k=1 and which corresponds to a straightforward perturbation of the simplest form of the DBM studied in <cit.>. The sign of a≠0 in (<ref>) establishes the form of the solution and the properties of the associated dynamical system. 

We start from the equilibria, i.e.:
    
    
  * if a<0, v(t)≡ w-√(-a/γ) is an asymptotically stable (in the future) constant solution of (<ref>); thus, for v_0>w-√(-a/γ) (v_0<w-√(-a/γ)) the CME monotonically decelerates (accelerates) for positive times;
   
  * if a>0, v(t)≡ w+√(a/γ) is an asymptotically stable (in the future) constant solution of (<ref>); thus, for v_0>w+√(a/γ) (v_0<w+√(a/γ)) the CME monotonically decelerates (accelerates) for positive times.
   
These assertions clarify the role of the acceleration term a≠0: it shifts the asymptotic solution from v=w (standard DBM) to v=w±√(± a/γ) (EDBM). In contrast to the case a=0, this means that for a>0 (a<0) initial speeds below (above) the wind speed can increase (decrease) up (down) to w and beyond. A schematic of the dynamics around the equilibrium points is provided in Figure <ref>.
    

Equation (<ref>) can be integrated from 0 to t>0 to obtain explicit formulae for v(t) and r(t), given the initial conditions v(0)=v_0, r(0)=r_0. Depending on the choice of v_0, the solutions may be (differentiably) piecewise-defined for positive or negative values of a due to the presence of the absolute value term |ṙ-w|, and present obvious symmetries in the form. Specifically,

    
  Case a>0.


    
    
  * if v_0≤ w, then
    
    v(t)=
             w+√(a/γ)tan(√(aγ)t-σ_+) ,    for  0≤ t≤1/√(aγ)σ_+
     w+√(a/γ)e^2(√(aγ)t-σ_+)-1/e^2(√(aγ)t-σ_+)+1 ,    for   t>1/√(aγ)σ_+
             ,

     
    r(t)=
             wt+r_0-1/γln(S_+cos(√(aγ)t-σ_+)) ,    for  0≤ t≤1/√(aγ)σ_+
    (w-√(a/γ))t+r_0+1/γ(ln(e^2(√(aγ)t-σ_+)+1/2S_+)+σ_+) ,    for   t>1/√(aγ)σ_+
             ,

    where σ_+arctan(√(γ/a)(w-v_0)) and S_+√((a+γ(v_0-w)^2)/a);

 
  * if v_0>w, then
     
    v(t)=
            w+√(a/γ)A_+e^2√(aγ)t+B_+/A_+e^2√(aγ)t-B_+ ,  for    t≥ 0 ,

    
    r(t)=(w-√(a/γ))t+r_0
                +1/γln(A_+e^2√(aγ)t-B_+/2√(a)) ,  for    t≥ 0 ,

    where A_+√(γ)(v_0-w)+√(a) and B_+√(γ)(v_0-w)-√(a).
    


  Case a<0.


    
    
  * if v_0≤ w, then
    
    v(t)=w+√(-a/γ)A_-e^-2√(-aγ)t+B_-/A_-e^-2√(-aγ)t-B_- ,  for    t≥0 ,

    
    r(t)=(w-√(-a/γ))t+r_0-1/γln(A_-e^-2√(-aγ)t-B_-/2√(-a)) ,  for    t≥0 ,

    where A_-√(γ)(v_0-w)+√(-a) and B_-√(γ)(v_0-w)-√(-a);    
    
  * if v_0> w, then
    
    v(t)= w-√(-a/γ)tan(√(-aγ)t-σ_-) ,    for   0≤ t≤1/√(-aγ)σ_-
     w+√(-a/γ)e^-2(√(-a γ)t-σ_-)-1/e^-2(√(-a γ)t-σ_-)+1 ,    for   t>1/√(-aγ)σ_-
             ,

    
    r(t)= wt+r_0+1/γln(S_-cos(√(-aγ)t-σ_-)) ,    for    0≤ t≤1/√(-aγ)σ_-
    (w-√(-a/γ))t+r_0-1/γ(ln(e^-2(√(-aγ)t-σ_-)+1/2S_-)-σ_-) ,    for    t>1/√(-aγ)σ_-
             ,

    where σ_-arctan(√(-γ/a)(v_0-w)) and S_-√((a-γ(v_0-w)^2)/a).   
    
    
    
The expressions for the CME's speed reflect the dynamical behavior of Figure <ref>: in (<ref>), v(t)≤ w in the former expression while v(t)>w in the latter, and v→ w+√(a/γ) as t→+∞; in (<ref>), v(t)>w with v→ w+√(a/γ) as t→+∞, and v(t) is never smaller than or equal to w for positive times; in (<ref>), v(t)≤ w with v→ w-√(-a/γ) as t→+∞, and v(t) is never larger than or equal to w for positive times; finally, in (<ref>), v(t)≥ w in the former expression while v(t)<w in the latter, and v→ w-√(-a/γ) as t→+∞.

The derivation of (<ref>)–(<ref>) requires standard calculus techniques, whose details are given in Appendix <ref>.



§ VALIDATION OF THE EDBM: THE NOVEMBER 3RD – 5TH 2021 EVENT


As discussed in <ref>., the closely spaced SolO-Wind detections of a CME of early November 2021 provided a reliable test-bed to assess the effectiveness of the EDBM. Indeed, in agreement with the analysis of <cit.>, the wind speed profiles measured by SolO and Wind (Figure <ref>) suggested an acceleration of the CME Magnetic Cloud (MC) front from SolO to Wind (bottom panel) rather than the expected deceleration due to the drag force induced by the background solar wind. Although the physical reasons for this local behavior remain unclear (plausible interpretations are discussed in <ref>.), in the following we applied the extended model described in <ref>. against data collected at SolO and Wind locations r_SolO and r_Wind to estimate the additive acceleration a>0 between the two instruments. Specifically, given the measured MC front speeds v_SolO,v_Wind at times t_SolO,t_Wind, respectively, we focused on the difference between the mean acceleration 

    a_mean=Δ v/Δ t=v_Wind-v_SolO/t_Wind-t_SolO  ,

which is an indicator of the approximate total measured acceleration exerted on the CME between the two spacecraft, and the model-dependent acceleration contributions a_drag(SolO)+a and a_drag(Wind)+a, where

    a_drag(SolO)=-γ|v_SolO-w|(v_SolO-w) ,
    
        a_drag(Wind)=-γ|v_Wind-w|(v_Wind-w) .



The main task, therefore, was to determine, from the solutions v(t),r(t) in <ref>., the values of the extra-acceleration term a that are compatible with the set of boundary values

    (v_SolO,v_Wind,r_SolO,r_Wind)=(690.86 km/s,705.97 km/s,0.85 AU,0.98 AU)  ,

obtained from the data time series at initial time t_SolO=0 and final time t_Wind=17820 s, and the parameters (γ,w). More specifically, we considered several experiments by choosing w∈[400,800] km/s with incremental step Δ w=50 km/s, and we used the same value γ=0.24×10^-7 km^-1 as in <cit.>, compatible with the CME erupted on November 2nd 2021 at 02:48 UT detected by the Solar and Heliospheric Observatory (SOHO) LASCO C2 coronograph (see, e.g., <cit.>). Through a standard root-finding Newton-Raphson method <cit.>, for every w one has to search for the solution(s) (if any) of

    f_v(a) v(a;t_SolO,t_Wind,v_SolO,w,γ)-v_Wind=0

or

    f_r(a) r(a;t_SolO,t_Wind,v_SolO,r_SolO,w,γ)-r_Wind=0

using formulae (<ref>)–(<ref>) (case a>0). We initialized the root-finding algorithm with an initial guess of approximately the same order of magnitude of |a_drag(SolO)|,|a_drag(Wind)|, and iterate until convergence to a local positive value (did not the scheme converge, we would set a=0).

For the sake of simplicity, we applied this scheme to 
(<ref>) and, since formula (<ref>) is case-defined and the time intervals depend on the unknown a, we eventually checked that for each experiment the corresponding time condition was fulfilled once a is found (did not the time condition apply, we would reject the solution).

Figure <ref> contains the results of this analysis. Specifically, in the top left panel a positive extra-acceleration a was obtained for each choice of w (cyan curve), as opposed to a drag deceleration at SolO and Wind until w=700 km/s (orange and green curves). The sum of the corresponding contributions provided two profiles that symmetrically fit the constant value a_mean in (<ref>) with a notable degree of accuracy, independently of the ambient solar wind (from the magnification in the top right panel, the maximum error committed is about 0.1 m/s^2 attained at w=400 km/s). Furthemore, it is worth mentioning the optimality reached at w=700 km/s, with almost an exact match between the three accelerations: indeed, this is the value for the solar wind closest to v_SolO and v_Wind.


The bottom panels of the same figure describe the outcomes of two further tests. First, we generated two sets containing ten values of the extra-acceleration a computed for two sets of ten random realizations of the initial speed in the range [v_SolO - 50,v_SolO + 10 ] km/s and of the final speed in the range [v_Wind-10,v_Wind+50] km/s, respectively (the reason for this choice of ranges is two-fold: it guarantees that v_SolO<v_Wind, and a maximum error of 50 km/s is plausible while accounting for the uncertainty on the temporal location of the MC boundary). The left panel contains average values ⟨ a⟩ and the corresponding standard deviations σ_v(SolO),σ_v(Wind) computed over the two sets (these standard deviations stabilize after ten random realizations of the inital/final speeds). Note that σ_v(SolO)≈σ_v(Wind)≈ 1 m/s^2 independently of w, and we coherently re-obtained the best agreement between the two profiles for w=700 km/s. Second, in the bottom right panel of Figure <ref>, we computed the absolute error ε=|f_r(a)| from (<ref>) at Wind location for the solutions a obtained as in the top left panel using (<ref>). The overall error as a function of w did not exceed ε=0.046805 AU (relative error ≈ 5%), attained at w=700 km/s. Note that for this value of the wind speed, we found, at the same time, the best outcome as far as a is concerned, though the largest error on r. This suggests to rely on a trade-off strategy when fitting the real data either with the v(t) model (equation (<ref>)) or the r(t) model (equation (<ref>)).

In this respect, we infer that the EDBM can accurately describe the dynamics of a vast sample of interplanetary CMEs, especially of those excluded by the simplest form of the DBM, like the ones propelled beyond the solar wind speed (case with w=700 km/s in Figure <ref>). Indeed, when w=700 km/s is assumed, the intermediate condition v_SolO<w<v_Wind holds, which cannot be modelled using the classical DBM (see <ref>). 



§ DISCUSSION AND CONCLUSIONS


As CMEs travel through interplanetary space, they can experience residual acceleration during their expansion due to several factors beyond solar wind drag. These include:


	
  * Magnetic forces: CMEs are highly magnetized plasma structures, carrying their own magnetic field. As they expand into interplanetary space, their magnetic field interacts with the Sun’s Interplanetary Magnetic Field (IMF). This interaction can generate magnetic forces that can lead to residual acceleration, depending on the alignment and strength of the magnetic fields. The magnetic pressure from the Sun’s field, which decreases with distance, may provide a residual push on the CME as it expands.
	
  * Pressure gradients: as CMEs move away from the Sun, they encounter regions of lower density and pressure. The difference between the internal pressure of the CME and the external pressure of the surrounding solar wind can cause the CME to continue expanding and accelerate. If the internal pressure of the CME remains higher than the external pressure for an extended period, this imbalance can drive residual acceleration during the CME’s expansion.
	
  * Internal magnetic reconfiguration: the internal dynamics of a CME, including its magnetic tension forces and plasma flows, can also contribute to residual acceleration. Indeed, CMEs contain complex magnetic structures that can undergo reconfiguration or magnetic reconnection as they expand. These internal processes can release energy, contributing to the acceleration of the CME. For example, if magnetic loops within the CME reconnect, the release of magnetic energy could provide a push that accelerates the CME further into space.
	
  * Gravitational forces: although relatively weak at large distances from the Sun, gravitational forces from the Sun can still play a role. Near the Sun, gravity decelerates the CME, but as it moves farther away, the influence of gravity decreases. If the CME has not yet reached a terminal velocity, the reduction in gravitational influence can result in a relative acceleration as the opposing force weakens.
	
  * Plasma and magnetic pressure balance: the expansion of the CME involves the balance between plasma pressure and magnetic pressure within the CME and in the surrounding solar wind. As the CME expands and its internal pressure decreases, the balance between these pressures can change, leading to further acceleration. If the magnetic pressure within the CME remains relatively high, it could continue to push the CME outward.


In general, these factors certainly contribute to the complex dynamics of CMEs as they travel through space, influencing their speed and trajectory beyond the initial influence of the solar wind. More specifically, the present study showed that the combination of these effects can be modelled by an extra-acceleration term that, when added to the drag force, contributes to explain the observations of a CME performed by SolO and Wind much more reliably than the standard DBM.

Understanding the processes itemized above is essential for predicting the behavior of CMEs and their potential impact on space weather and Earth’s environment. Disentangling which of these processes is currently at work in the evolution of CME under studio needs further analysis, which, as it is beyond the scope of the present work, is devoted to a future paper. Finally, it is also worth noting the interesting possibility of combining this extended drag-based model in the neural network developed by <cit.>, so as to also refine their AI-based model and potentially make it even more predictive.



§ ACKNOWLEDGMENTS

SG was supported by the Programma Operativo Nazionale (PON) “Ricerca e Innovazione” 2014–2020. All authors acknowledge the support of the Fondazione Compagnia di San Paolo within the framework of the Artificial Intelligence Call for Proposals, AIxtreme project (ID Rol: 71708). AMM is also grateful to the HORIZON Europe ARCAFF Project, Grant No. 101082164. SG, MP and AMM are also grateful to the Gruppo Nazionale per il Calcolo Scientifico
- Istituto Nazionale di Alta Matematica (GNCS -
INdAM). MR is also grateful to the Gruppo Nazionale per la Fisica Matematica
- Istituto Nazionale di Alta Matematica (GNFM -
INdAM).



§ APPENDIX




 §.§ Computation of the solutions of the EDBM





  §.§.§ v_0,v≤ w


Assume to integrate (<ref>) in [0,t] such that v_0,v(t)≤ w:

    ∫_v_0^vdv'/a+γ(v'-w)^2=∫_0^tdt' .

Distinguishing between a>0 and a<0, we obtain two different primitives for the left-hand side and get:

    t=
    1/√(aγ)(arctan(√(γ/a)(v-w))-arctan(√(γ/a)(v_0-w))) ,   a>0
     -1/2√(-aγ)ln((√(γ)(v-w)+√(-a))(√(γ)(v_0-w)-√(-a))/(√(γ)(v_0-w)+√(-a))(√(γ)(v-w)-√(-a))) ,   a<0
     ;

now solving for v in both the expressions and setting v(t)≤ w yield formulae (<ref>) in the case 0≤ t≤arctan(√(γ/a)(w-v_0))/√(aγ) and (<ref>). 

As regards r(t), for a>0 a second integration provides

    r(t)   =∫_0^t v(t')dt'=wt+r_0+√(a/γ)∫_0^ttan(√(aγ)t'+arctan(√(γ/a)(v_0-w)))dt'
       =wt+r_0-1/γln(|cos(√(aγ)t+arctan(√(γ/a)(v_0-w)))|√(a+γ(v_0-w)^2/a)) ,

using |cos x|=1/√(1+tan^2x). In the interval [0,arctan(√(γ/a)(w-v_0))/√(aγ)] the cosine is positive, so we can remove the absolute value and obtain equation (<ref>) (first case).

For a<0, we conveniently set A√(γ)(v_0-w)+√(-a), B√(γ)(v_0-w)-√(-a) and C-2√(-aγ). We have

    r(t)=∫_0^tv(t')dt'=wt+∫_0^tAe^Ct'+B/Ae^Ct'-B√(-a/γ)dt'+r_0=(w-√(-a/γ))t+r_0+2/C√(-a/γ)ln|Ae^Ct-B/A-B| ,

where the integral is first computed by splitting the fraction as

    Ae^Ct'+B/Ae^Ct'-B=1+2B/Ae^Ct'-B ,

and then performing the two subsequent (monotonic) changes of variable u=Aexp(Ct')-B and U=1+B/u. Again, we can disregard the absolute value for t≥0, and replacing back the values of A,B,C we get (<ref>).



  §.§.§ v_0,v≥ w


The procedure is the same as in Appendix <ref>. Since now the integration of (<ref>) reads

    ∫_v_0^vdv'/a-γ(v'-w)^2=∫_0^tdt' ,

formulas derived for v(t) and r(t) are simply swapped for a≶0. Upon substituting a↦-a, √(aγ)↦-√(-aγ) or the other way around depending on the sign of a, an analogous argument for time intervals, absolute values, and a corresponding re-definition of constants A,B,C hold. This leads to equations (<ref>) (a>0), (<ref>) (first case, a<0) for v(t), and (<ref>) (a>0), (<ref>) (first case, a<0) for r(t).



  §.§.§ v_0<w<v


This time the integration of the EDBM gives rise to two contributions:

    ∫_v_0^wdv'/a+γ(v'-w)^2+∫_w^vdv'/a-γ(v'-w)^2=t ;

in addition, from Figure <ref>, the forward dynamics rules out the case a<0 (it is possible only backward in time) and requires t>t_*arctan(√(γ/a)(w-v_0))/√(aγ) to arrive at v(t)>w (cf. Appendix <ref>). Then, we have

    ∫_w^vdv'/a-γ(v'-w)^2=t-t_* ,

which is the case of Appendix <ref> with lower limit of integration equal to w. So, upon solving for v, we obtain expression (<ref>) in the case t>t_*.

Concerning r(t), we need to integrate equation (<ref>) (second case) from t_* to t:

    r(t)=w(t-t_*)+∫_t_*^tAe^Ct'+B/Ae^Ct'-B√(a/γ)dt'+r_* ,

with r_*=r(t_*), Aexp(2arctan(√(γ/a)(v_0-w))), B -1, C2√(aγ). This relationship is formally identical to the one of r(t) in Appendix <ref>, case a<0. We find

    r(t)=(w-√(a/γ))(t-t_*)+r_*+1/γln(Ae^Ct-B/Ae^Ct_*-B) .

Lastly, we determine r_* by enforcing continuity at t=t_* with the former expression in (<ref>):

    r_*=wt_*+r_0-1/γln√(a+γ(v_0-w)^2/a) ;

hence, replacing back in the expression of r(t), we retrieve the latter of (<ref>).



  §.§.§ v<w<v_0



Following Appendix <ref>, an analogous reasoning is applied in the case a<0. The corresponding sign adaptation of quantities involving a (see Appendix <ref>) and the continuity requirement with the first relationship of equation (<ref>) produce formulae (<ref>), (<ref>) for t>arctan(√(-γ/a)(v_0-w))/√(-aγ). 


aasjournal

