Speeding Up Solar Wind Forecasts with Reduced-order Modeling

By Opal Issan and Boris Kramer

The high-speed solar wind streams that cause the mesmerizing northern and southern lights are also a major source of geomagnetic storms. Such space weather events can have significant ramifications for satellites, electric power grids, aviation, and signal accuracy of the Global Positioning System; they also expose astronauts to harmful corpuscular radiation [4]. A moderate geomagnetic storm on February 2, 2022, increased the density of Earth’s upper atmosphere, which in turn increased the drag of low-Earth orbit satellites. As a result, 40 of the 49 SpaceX Starlink internet satellites that launched the following day failed to reach their orbits and burned up as they reentered Earth’s atmosphere [1].

The solar wind is a continuous flow of charged particles—mostly protons, helium nuclei, and electrons—that escape the Sun’s gravitational pull and travel supersonically into interplanetary space. The solar wind moves at speeds that range from 400 to 800 kilometers per second and typically takes between two to four days to travel from the surface of the Sun to Earth. Its topology resembles the spiral pattern of a spinning sprinkler, in that the solar wind streams (analogous to water ejected from a sprinkler) travel radially outwards as the Sun (or sprinkler’s head) rotates. The solar wind is also tied to the Sun’s magnetic field lines via magnetohydrodynamic (MHD) effects, which bend the magnetic field lines and result in an advection-dominated system. 

Solar wind streams are primarily predicted via three-dimensional numerical MHD models. Although MHD models can reproduce solar wind’s structure and properties with high fidelity, they are computationally expensive. For example, a recent study [3] took about 7,000 central processing unit hours to simulate a single solar rotation (which is approximately 25 days). Due to their high computational costs, MHD models cannot feasibly make time-sensitive predictions and generate the large ensembles that are required for uncertainty quantification. There is hence a need for computationally efficient surrogate models that accurately approximate MHD simulations.

To address this challenge, we developed shifted operator inference (sOpInf): a data-driven, projection-based, reduced-order model that predicts the solar wind near Earth [2]. The proposed method is computationally efficient and preserves the solar wind’s Archimedean spiral structure. Here, we briefly describe the sOpInf framework and demonstrate its ability to approximate the three-dimensional simulations of the Magnetohydrodynamics Algorithm outside a Sphere (MAS) model.

sOpInf: Non-intrusive Reduced-order Modeling for Advection-dominated Systems

Operator inference is a data-driven, non-intrusive, reduced-order modeling framework that aims to infer operators from spatiotemporal data, such as complex high-fidelity simulations [5]. sOpInf extends operator inference towards advection-dominated systems that are described on a periodic domain, such as the solar wind. In order to preserve the translational and rotational properties of the physical system while remaining low dimensional, this method transforms the original coordinates to a moving coordinate frame where the dynamics are devoid of translation and rotation.

To illustrate the sOpInf framework, we consider a generic spatially \(k\)-dimensional time-dependent partial differential equation for the scalar function \(v(x_1, x_2, \ldots, x_k, t)\) of the form

\[\mathcal{F}\left (v,x_{1},\ldots, x_{k}, \frac{\partial v}{\partial t}, \frac{\partial v}{\partial x_{1}}, \ldots , \frac{\partial v}{\partial x_{k}}, t \right ) = 0,\tag1\]

where \(x_1, x_2, \ldots, x_k \in \mathbb{R}^{k}\) denote the spatial coordinates and \(t \in \mathbb{R}^+ \) denotes time. The main steps of sOpInf are as follows.

(I) Data Collection and Shifting to Advection-Invariant Coordinates. The first step is to numerically simulate \((1)\). The numerical solution is described on a discretized mesh, which we call the original coordinates and denote with \(\mathbf{x} \in \mathbb{R}^n\). Next, we interpolate the snapshots \(\mathbf{v}_i\) at instances \(t_i\) with \(0 = t_0 < t_1 < \dots < t_K = T\) onto the shifted coordinates \(\tilde{\mathbf{x}} (\mathbf{x}, t_i)\):

\[\mathbf{v}_i \approx v(\mathbf{x}, t_i) \mapsto \tilde{v} (\tilde{\mathbf{x}}(\mathbf{x}, t_i), t_i) \approx \tilde{\mathbf{v}}_i \qquad \textrm{with} \qquad \tilde{\mathbf{x}}(\mathbf{x}, t) = \mathbf{x} + \mathbf{c}(t).\tag2\]

Here, we can find \(\mathbf{c}(t) \in \mathbb{R}^n\)—which represents the advection speed—using analytic and fully data-driven approaches [2]. As Figure 1a illustrates, the solar wind speed in the heliographic coordinate system (original coordinates) exhibits advection in the opposite direction of the Sun’s rotation; however, in Figure 1b, the solar wind dynamics in the shifted coordinates are nearly free of advection.

Figure 1. The importance of coordinate transformation in model reduction for advection-dominated systems. 1a. The Magnetohydrodynamics Algorithm outside a Sphere (MAS) model solar wind radial velocity results at the equatorial plane for Carrington Rotation 2210 from the inner heliosphere to Earth. 1b. The solar wind data in shifted coordinates eliminate the translational properties that are caused by the Sun’s rotation. 1c. The singular value cumulative energy of the data in the original and shifted coordinates as obtained by (3). The singular values decay more rapidly in the shifted coordinates, indicating that the ROM will require fewer modes in those coordinates. Figure courtesy of [2].

(II) Data Reduction via Projection. This step projects the dynamics onto a low-dimensional subspace that captures the dominant structure and patterns of the snapshot data. To do so, we use the subspace that is spanned by the proper orthogonal decomposition (POD) modes, which construct the basis that optimally represents the solution data in the \(L_2\) sense. We obtain the POD modes by computing the economy-sized singular value decomposition of the snapshot matrix

\[\tilde{\mathbf{V}} = [\tilde{\mathbf{v}}_1 \quad \dots \quad  \tilde{\mathbf{v}}_K] = \boldsymbol{\Psi}\boldsymbol{\Sigma}\mathbf{W}^{\top},\]

where \(\boldsymbol{\Psi}\in \mathbb{R}^{n \times K}\), \(\boldsymbol{\Sigma} \in \mathbb{R}^{K \times K}\), and \(\mathbf{W} \in \mathbb{R}^{K \times K}\). The \(\ell \ll n\) dimensional POD basis \(\boldsymbol{\Psi}_\ell = [\mathbf{v}_1,\ldots,\mathbf{v}_\ell]\) is given by the first \(\ell\) columns of \(\boldsymbol{\Psi}\). We typically choose the truncation based on the cumulative energy criteria

\[\ell = \underset{\hat{\ell}}{\operatorname{arg min}} \frac{\sum_{i=1}^{\hat{\ell}}\sigma_{i}}{\sum_{i=1}^{K}\sigma_{i}} > \epsilon_{\text{tol}},\tag3\]

where \(\sigma_{i}=\boldsymbol{\Sigma}_{ii}\) are the singular values and \(\epsilon_{\text{tol}}\) is often \(\epsilon_{\text{tol}} =0.99\), which encapsulates 99 percent of the energy in the data. Figure 1c demonstrates that the MHD simulation results on the shifted coordinates can capture 99 percent of the system’s cumulative energy with only three modes; in contrast, those same results require many more modes on the original coordinates. This outcome emphasizes the need to first shift the data to advection-invariant coordinates in order to find the optimal low-dimensional basis. 

Lastly, we project the state snapshot data onto the POD subspace that is spanned by the columns of \(\boldsymbol{\Psi}_\ell\) and obtain

\[\widehat{\mathbf{V}} = \boldsymbol{\Psi}^{\top}_\ell \tilde{\mathbf{V}} = [\widehat{\mathbf{v}}_1 \quad \dots \quad \widehat{\mathbf{v}}_K] \in \mathbb{R}^{\ell \times K}
\qquad \text{and} \qquad 
\dot{\widehat{\mathbf{V}}} = [ \dot{\widehat{\mathbf{v}}}_1\quad \dots \quad  \dot{\widehat{\mathbf{v}}}_K]\in \mathbb{R}^{\ell \times K}.\]

Here, we use a finite-differencing time derivative approximation to compute the columns of \(\dot{\widehat{\mathbf{V}}}\) from \(\widehat{\mathbf{V}}\).

(III) Model Learning via Operator Inference. Next, we learn the low-dimensional operators of a polynomial reduced-order model. For example, consider the quadratic reduced-order model

\[\frac{\mathrm{d} \widehat{\mathbf{v}}} {\mathrm{d}t} = \widehat{\mathbf{A}} \widehat{\mathbf{v}}+ \widehat{\mathbf{H}} (\widehat{\mathbf{v}}\otimes \widehat{\mathbf{v}}) + \widehat{\mathbf{B}}, \qquad \widehat{\mathbf{v}} \in \mathbb{R}^\ell,\]

where \(\otimes\) is the Kronecker product. How do we infer the low-dimensional operators \(\widehat{\mathbf{A}}\in \mathbb{R}^{\ell \times \ell}\), \(\widehat{\mathbf{H}} \in \mathbb{R}^{\ell \times \ell^2}\), and \(\widehat{\mathbf{B}}\in \mathbb{R}^{\ell}\)? At the heart of operator inference, we solve a regression problem to find the low-dimensional operators:

\[\underset{ \widehat{\mathbf{A}} \in \mathbb{R}^{\ell \times \ell} , \widehat{\mathbf{H}} \in \mathbb{R}^{\ell \times \ell^2}, \widehat{\mathbf{B}} \in \mathbb{R}^\ell}{\min} \sum_{j=0}^{K-1}
\left \Vert \widehat{\mathbf{A}} \widehat{\mathbf{v}}_{j} + \widehat{\mathbf{H}}(\widehat{\mathbf{v}}_{j} \otimes \widehat{\mathbf{v}}_{j})  + \widehat{\mathbf{B}} - \dot{\widehat{\mathbf{v}}}_{j} \right \Vert^2_{2}.\]

Once we learn the operators, we can simulate the reduced-order model on the shifted coordinates. We then solve the minimization problem via Tikhonov regularization, which handles ill-posed problems and avoids overfitting.

(IV) Re-shifting to Original Coordinates. Finally, we shift the reduced-order model’s predicted dynamics back to the original coordinates via interpolation, essentially reversing the coordinate transformation in \((2)\).

Shifted Operator Inference Learns Accurate Reduced-order Models

We tested the methodology with the MAS model for Carrington Rotation 2210, which occurred during a solar minimum from October 26 to November 23, 2018. Figure 2 compares the sOpInf reconstructed and predicted solar wind radial velocities. The visual comparison and error estimates indicate that sOpInf can successfully reproduce the high-fidelity MAS simulations—where \(n=14{,}208\) grid points—with only \(\ell = 8\) modes, leading to a substantial reduction in the model’s dimensionality.

Figure 2. The shifted operator inference (sOpInf) model is able to successfully approximate the complex Magnetohydrodynamics Algorithm outside a Sphere (MAS) simulations. 2a. sOpInf model predictions (left column) with \(\ell=8\) modes compared to the MAS solar wind velocities (right column). The training ends at \(0.82\) astronomical units (AU), so the velocity profile at \(r=1.1\) AU (bottom row) is in the purely predictive regime. 2b. The relative error of the learned full-Sun sOpInf reduced-order model versus MAS. The sOpInf model results at \(r=0.968\) AU and \(r=1.1\) AU (bottom row) are in the purely predictive regime. The reduced-order model demonstrates good agreement with the MAS solution but can be evaluated at a much lower cost. Figure courtesy of [2].

The Need for Reduced-order Modeling in Space Weather Operational Forecasting

This work highlights the potential of reduced-order modeling approaches in space weather applications. But the story does not end here; in the future, we plan to use reduced-order modeling to solve space weather uncertainty quantification problems with techniques such as parametric global sensitivity analysis and Bayesian inference. These types of methods directly address the field’s ongoing challenge to estimate appropriate values for input parameters, initial conditions, and boundary conditions in the presence of sparse, heterogeneous, and noisy measurements. Reduced-order modeling can play a key role in accelerating uncertainty quantification tasks and enabling more informed decisions in operational settings.

SIAM Student Travel Awards are made possible in part by the generous support of our community. To make a gift to the Student Travel Fund, visit the SIAM website. 

References
[1] Andrews, R.G. (2022, February 9). Solar storm destroys 40 new SpaceX satellites in orbit. The New York Times. Retrieved from https://www.nytimes.com/2022/02/09/science/spacex-satellites-storm.html.
[2] Issan, O., & Kramer, B. (2023). Predicting solar wind streams from the inner-heliosphere to Earth via shifted operator inference. J. Comput. Phys., 473, 111689.
[3] Jivani, A., Sachdeva, N., Huang, Z., Chen, Y., van der Holst, B., Manchester, W., … Toth, G. (2023). Global sensitivity analysis and uncertainty quantification for background solar wind using the Alfvén Wave Solar atmosphere Model. Space Weather, 21(1), e2022SW003262.
[4] Moldwin, M. (2008). An introduction to space weather. Cambridge, U.K.: Cambridge University Press.
[5] Peherstorfer, B., & Willcox, K. (2016). Data-driven operator inference for nonintrusive projection-based model reduction. Comput. Methods Appl. Mech. Eng., 306(1), 196-215.

	Opal Issan is a Ph.D. student in the Department of Mechanical and Aerospace Engineering at the University of California, San Diego. She develops reduced-order models and Bayesian inference methods for uncertainty quantification in space weather applications.
	Boris Kramer is an assistant professor of mechanical and aerospace engineering at the University of California, San Diego. He develops reduced-order models to enable the control, design, optimization, and uncertainty quantification of complex and large-scale systems.