Coronal Mass Ejections (CMEs) are the primary driver of major geomagnetic storms, and forecasting their arrival time at Earth hinges on knowing how fast they travel and how strongly the ambient solar wind decelerates (or accelerates) them. In operational space weather forecasting, the workhorse for this is the Drag-Based Model (DBM), which reduces the messy magnetohydrodynamics of a CME’s propagation to a single aerodynamic drag equation. Given a handful of noisy height-time measurements from coronagraphs or heliospheric imagers, we can invert this model to recover the CME’s initial speed, the ambient solar wind speed, and the drag parameter — a classic nonlinear parameter estimation problem.
This post walks through a complete, self-contained example: simulating synthetic CME tracking data, fitting the DBM to it with nonlinear least squares, and visualizing both the fit quality and the shape of the underlying cost landscape in 3D.
The Drag-Based Model
The DBM assumes the CME leading edge experiences an aerodynamic drag force proportional to the square of its relative speed with respect to the solar wind:
$$
\frac{dv}{dt} = -\gamma (v - w),|v - w|
$$
where $v$ is the CME speed, $w$ is the ambient solar wind speed, and $\gamma$ is a drag parameter (units of $\text{km}^{-1}$) that lumps together the CME’s cross-sectional area, mass, and the ambient density.
For the deceleration branch ($v_0 > w$, with constant $\gamma$), this equation integrates analytically:
$$
v(t) = \frac{v_0 - w}{1 + \gamma (v_0 - w) t} + w
$$
$$
r(t) = r_0 + w t + \frac{1}{\gamma}\ln!\big(1 + \gamma (v_0 - w) t\big)
$$
Given a time series of observed heights $r(t)$, the inverse problem is to recover $(v_0, w, \gamma)$ by nonlinear least squares. Because the forward model is fully analytic, this fit is extremely cheap computationally — no numerical ODE integration is required, which is what keeps the code below fast even when we sweep a dense grid for the cost surface.
Full Colab Source
1 | # ========================================================== |
Code Walkthrough
Sections 1–2 — the physical model. dbm_height and dbm_speed implement the closed-form solution of the drag equation directly, with np.log1p used instead of np.log(1 + x) for better numerical stability when the argument is small. Because $v_0 > w$ throughout this example, the argument of the logarithm stays strictly positive and no domain errors can occur.
Section 3 — synthetic observations. Real CME height-time data comes from manually or automatically tracking the leading edge in a sequence of coronagraph (LASCO) or heliospheric imager (STEREO/HI) images, which is inherently noisy. Here we generate 24 “tracked” points over 60 hours from known ground-truth parameters and add 1.5% multiplicative Gaussian noise, mimicking that measurement uncertainty.
Section 4 — the inversion. scipy.optimize.curve_fit performs bounded nonlinear least squares (the Trust Region Reflective algorithm is selected automatically once bounds are supplied). The bounds keep $\gamma$ strictly positive and $v_0, w$ within physically reasonable ranges, which both prevents the optimizer from wandering into the singular region near $\gamma = 0$ and speeds up convergence. pcov gives the parameter covariance matrix, from which we extract 1-sigma uncertainties via its diagonal.
Section 5 — the cost surface, vectorized. A naive implementation would loop over every $(w, \gamma)$ grid point and every time sample. Instead, t_sec, W, and G are reshaped so NumPy’s broadcasting rules evaluate the entire $24 \times 120 \times 120$ tensor of model heights in one call, then reduce over the time axis. This is the “fast” version of what would otherwise be a triple-nested loop, and it finishes in well under a second even on Colab’s default CPU runtime.
Section 6 — the plots. The left panel overlays the noisy synthetic observations with the fitted curve; the middle panel shows residuals to check that no systematic trend remains (a good fit should look like scattered noise around zero); the right panel is a 3D log-cost surface over $w$ and $\gamma$ with the best-fit point marked, which visually confirms the fit landed in the basin of the minimum rather than a spurious local optimum.
===== CME Drag-Based Model: Parameter Estimation Results ===== v0 : true= 1450.00 km/s fit= 1393.57 +/- 96.73 km/s w : true= 400.00 km/s fit= 384.41 +/- 36.20 km/s gamma : true=3.500e-08 1/km fit=3.099e-08 +/- 9.541e-09 1/km RMSE (height fit) : 1.7015 solar radii Speed at final observation - true: 517.5 km/s, fit: 514.5 km/s
Interpreting the Results
If the fit is behaving well, the recovered $v_0$, $w$, and $\gamma$ should sit close to their true values, with the reported 1-sigma uncertainties reflecting how well-constrained each parameter is by the observation cadence and noise level — in practice, $\gamma$ is typically the hardest parameter to pin down because its effect on the height-time curve is subtle compared to $v_0$ and $w$. The 3D cost surface makes this concrete: a shallow, elongated valley along the $\gamma$ axis indicates that many $(w, \gamma)$ combinations produce nearly indistinguishable trajectories, which is exactly the kind of parameter degeneracy that makes real-world CME arrival-time forecasting genuinely difficult, even when the underlying physical model is this simple.
