{'fraction_non_alphanumeric': 0.04293818061363919, 'fraction_numerical': 0.0165072529777925, 'mean_word_length': 4.525777913361806, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Predicting Solar Wind Streams from the Inner-Heliosphere to Earth via Shifted Operator Inference

Solar wind conditions are predominantly predicted via three-dimensional numerical magnetohydrodynamic (MHD) models. Despite their ability to produce highly accurate predictions, MHD models require computationally intensive high-dimensional simulations. This renders them inadequate for making time-sensitive predictions and for large-ensemble analysis required in uncertainty quantification. This paper presents a new data-driven reduced-order model (ROM) capability for forecasting heliospheric solar wind speeds. Traditional model reduction methods based on Galerkin projection have difficulties with advection-dominated systems -- such as solar winds -- since they require a large number of basis functions and can become unstable. A core contribution of this work addresses this challenge by extending the non-intrusive operator inference ROM framework to exploit the translational symmetries present in the solar wind caused by the Sun's rotation. The numerical results show that our method can adequately emulate the MHD simulations and outperforms a reduced-physics surrogate model, the Heliospheric Upwind Extrapolation model.


Introduction
Magnetohydrodynamic (MHD) modeling of coronal and interplanetary solar wind can significantly improve the prediction of catastrophic space weather events. Such space weather geomagnetic storms can have detrimental effects on spacecrafts, cause electric power outage, satellite collisions, telecommunication interruption, and expose astronauts to harmful radiation. Very recently, 40 out of 49 of SpaceX's starlink satellites failed to reach their low-Earth orbits presumably due to the effects of a geomagnetic storm around Feb. 2, 2022; the geomagnetic storm increased Earth's upper atmosphere density causing orbital drag [6]. This event further emphasises the need for real-time modeling of solar storms. The most substantial source of space weather events are coronal mass ejections, corotating interaction regions, and high-speed solar wind streams that reach the Earth's magnetosphere. Three-dimensional MHD solar wind models, such as the Magnetohydrodynamics Around a Sphere [55] model, Enlil [40], and the Space Weather Modeling Framework [63], can provide high-fidelity predictions. Apart from providing a global assessment of coronal and heliospheric properties, MHD modeling can connect in situ magnetic and plasma observations from one spacecraft to the other, providing crucial support to interplanetary missions [57]. Although MHD models are an important tool in understanding observed coronal and heliospheric dynamics, they require computationally intensive high-dimensional simulations. This renders them infeasible for time-sensitive predictions and large-ensemble methods such as quantifying forecast uncertainty and performing parameter sensitiv-ity analysis. Thus, there is a need for computationally efficient surrogate models that are capable of reproducing MHD results with sufficient fidelity.
The present work addresses this challenge by proposing a new method to derive a data-driven reduced-order model (ROM) for solar wind predictions that particularly focuses on issues that arise from the advection-dominated nature of the problem. The proposed method efficiently learns predictive ROMs from high-fidelity simulations (or other data) of solar wind models, while simultaneously producing an interpretable model and physical model form. We subsequently review the related literature, both focusing on the mathematical aspects of surrogate modeling for advection-dominated systems and the application domain of efficient heliospheric modeling.
Solar wind predictions are produced by a chain of coupled models in different parts of the Sun-Earth domain, i.e., the solar surface, corona, and heliosphere. In the heliospheric domain, there are mainly two classes of surrogate solar wind models: reduced-physics (white-box) and data-driven (black-box) models. The first approach is based on physical simplifications of the MHD equations. The simplest reduced-physics model is the ballistic approximation which assumes that each solar wind parcel maintains a constant radial speed as it propagates in the heliosphere. This approximation is mainly used to map solar wind streams for short radial distances [61]. An improved kinematic model that bridges the gap between the ballistic mapping and global three-dimensional MHD modeling is the Heliospheric Upwind Extrapolation model, where each parcel speed is dependent on its adjacent parcel speed [56,52,42,51]. A second approach is to build surrogate solar wind models via data-driven and statistical techniques [12]. Such methods mainly aim to forecast the solar wind at Earth's vicinity without computing the solar wind dynamics on the full heliospheric domain. Examples of data-driven models include an artificial neural network model [67], a gradient-boosting regression-based model [9], and a probability distribution function model based on past rotation solar wind observations [11]. As we will see later, our proposed reduced-order modeling methodology is data-driven, yet accounts for the physical properties of the solar wind rotation. It therefore serves as a hybrid gray-box approach that leverages available physical information while remaining computationally efficient.
From a methodological perspective, several ROM strategies have targeted advection-dominated scenarios. For background, traditional ROMs are derived via (Petrov-) Galerkin projection, where the full-order model (FOM) is projected onto a low-dimensional subspace, see [23,7,22]. This class of ROMs aims to identify a small set of basis functions that minimize a certain error metric. However, a well-reported issue with linear-subspace ROMs is that they fail to model advection-dominated problems due to poorly decaying Kolmogorov N-width [19,41,65] which results in a slow decay of the singular values, see, e.g., [59,24,39,50,34]. An accurate ROM would require a large number of basis functions, rendering it inefficient from a computational perspective. Additionally, a large number of global basis functions can lead to numerical instabilities. The efforts to address the challenge posed by advection-dominated systems can be roughly categorized into Lagrangian-based approaches and methods that leverage a transport-invariant coordinate frame. The first line of research leverages Lagrangian coordinate grids to build a ROM that propagates both the wave physics and the coordinate grid in time [39,34,33]. These methods work extremely well, yet require knowledge of the underlying equations to solve for the Lagrangian grid, limiting their range of applicability. The second line of research, which our work builds upon, is based on transforming the dynamics to a moving coordinate frame via a time-dependent shift that is added to the spatial coordinates. In the moving frame, the system dynamics are absent of advection. The shift function can be numerically learned in various ways. For instance, the shifted proper orthogonal decomposition method [50] proposes to detect the shift either through tracking peaks of the solution, or through an expansive SVD algorithm, where different candidate shifts are applied to the data matrix and then the SVD is computed. The shifts leading to the best singular value decay are then selected and a corresponding ROM is built via Galerkin projection. This strategy is computationally quite expensive due to the need for many SVDs of a large (yet rectangular) data matrix. The authors in [38] propose an implicit feature-tracking algorithm that is based on a minimal-residual ROM. This algorithm works well on complex geometries, however, finding the domain mapping can be quite expensive, too. Very recent work by [43] builds on the trend of machine learning to derive two separate neural networks, one for detecting nonlinear shifts in the transport velocity, and a second for interpolating a shifted solution back to the reference frame. The method is fully data-driven, yet does not propose a predictive ROM that integrates the shift detection with a projection framework. Most closely related to our work is [37], which proposes an unsupervised learning method to aid the identification and low-dimensional modeling of systems with translational symmetries. The data-based method uses sparse regression and ridge detection to identify models with non-constant wave speeds inherent to the data. The method performs well on several examples, including multiple waves travelling in opposite directions. While the examples all demonstrate interesting wave phenomena, the governing equations were always known, allowing to derive solid intuition about the wave-speed library. Other approaches that have been developed that do not fit precisely into these categories include the work by [44] that updates the ROM basis online to avoid slowly decaying Kolmogorov N-width, and the work by [24] that applies two separate mechanisms to deal with advection-diffusion systems, namely a representation of translational features via advection modes, and then the subsequent residual (features that are not purely advective) via global modes. Building on the philosophy of inducing time-dependent shifts into the ROM framework, we extend the non-intrusive projection-based operator inference ROM framework [46] towards advectiondominated systems by transforming the dynamics to a moving coordinate frame. Standard operator inference has successfully been applied to diverse applications such as combustion [62,36], chemical reactors [10], ocean flows [68], Hamiltonian systems [60], and general reaction systems in the presence of incomplete data [64]. Since operator inference learns the operators that would be obtained through intrusive Galerkin projection (which can be done exactly with additional data pre-processing, see [45]) it inherits problems that intrusive Galerkin ROMs face in the presence of strong advection.
We propose a new strategy for efficient data-driven heliospheric solar wind modeling. The method, shifted operator inference (sOpInf ), builds on standard operator inference and extends it towards the challenges faced in solar wind predictions. Our proposed method first determines a moving coordinate frame where the dynamics are absent of translation and rotation and subsequently transforms the system into the new time-dependent coordinates. Two methods for predicting the shift are proposed. It then performs model learning in the shifted coordinate systems and subsequently makes predictions with the sOpInf-ROM via interpretable ODE simulation. Our hypothesis aligns with the previously cited references, in that simple translational patterns in the data and model can (and should) be exploited in the ROM approach. Our proposed approach (1) speeds up the MHD simulation by several orders of magnitude, (2) preserves the solar wind spiral pattern created by the Sun's rotation, and (3) uncovers macroscopic coherent structures present in the evolution of solar wind streams by analyzing the velocity field modal decomposition. We present computationally-efficient data-driven ROMs for two heliospheric solar wind models: the MAS (Magnetohydrodynamics Around a Sphere) model and the HUX (Heliospheric Upwinding eXtrapolation) model. This paper is organized as follows. Section 2 describes the MAS and HUX heliospheric solar wind models and in Section 3 we present the proposed method, shifted operator inference. In Section 4 we demonstrate the performance of sOpInf on the MAS and HUX solar wind speed simulated data. Section 5 then offers conclusions and an outlook to future work.

Solar Wind Models
This section introduces the solar wind models considered in this study. Section 2.1 presents the MAS model. Section 2.2 discusses an approximation to that model with similar physical attributes, the HUX model.

Spherical Magnetohydrodynamics: The MAS Model
The MAS (Magnetohydrodynamics Around a Sphere) model is the primary MHD model in the CORHEL (CORona-HELiosphere) software and is publicly available at NASA's community-coordinated modeling center [3]. The MAS model solves the time-dependent resistive MHD equations and has been used to study coronal mass ejections [31], coronal dynamics [55], solar wind structure [54], and connect in-situ spacecraft observations [57]. Herein, we focus our effort on analyzing the MAS solar wind radial velocity results and therefore exclude the discussion of other plasma components such as the magnetic field, plasma temperature, density, pressure, etc. This is because many space weather operational forecast models for satellite control and Earth-based infrastructure are particularly interested in near-Earth solar wind speed. Predicting the solar wind speed is pivotal for assessing the risk of geomagnetic storms because (1) coronal mass ejections, which are the most fundamental source of space weather events, are modeled as perturbations to the ambient solar wind; (2) interaction regions between fast and slow solar wind, known as co-rotating interaction regions, mainly present during solar minimum, are a driver of moderate geomagnetic activity [54]; and (3) high-speed solar wind streams cause an additional acceleration of energetic electrons in the radiation belts [16].

Governing Equation
The MAS model solves a system of three-dimensional time-dependent resistive MHD equations in spherical coordinates (r, θ, φ), where r is the radial distance from the Sun, θ is Carrington latitude in heliographic (rotating) coordinate system (HG), φ is the Carrington longitude in the HG coordinate system. The governing equations are and the initial and boundary condition are described in [54,30,53]. Here, B(r, θ, φ, t) is the magnetic field, J(r, θ, φ, t) is the current density, E(r, θ, φ, t) is the electric field, v(r, θ, φ, t) is the plasma velocity, T (r, θ, φ, t) is the plasma temperature, ρ(r, θ, φ, t) is the plasma mass density, and p(r, θ, φ, t) is the plasma pressure, and p w (r, θ, φ, t) is the Alfvén wave pressure. The constant c denotes the speed of light in a vacuum and g(r) = − GMs r 2êr is the gravitational acceleration, whereê r is the unit vector in the radial direction, G is the universal gravitational constant, and M s = 1.99 × 10 30 kg is the solar mass. For the simulations used in this study, the constant resistivity is set to η = 4.6779 × 10 −5 s and the kinematic viscosity ν = 3.3503 × 10 16 m 2 s −1 . In the energy equation described in Eq. (6), the thermodynamic approximation sets the ratio of specific heats to γ = 5/3. Moreover, the energy source terms are denoted by S = S(T ); for more details about the thermodynamic energy source term see the MAS user guide [53]. The MAS boundary conditions exploit photospheric magnetic field observations (e.g. data from the Wilcox Solar Observatory, the Global Oscillation Network Group, and the Solar Dynamics Observatory spacecraft), see [1,53]. Here, we choose to analyze the scenario where the thermodynamic MAS results are driven by a synoptic map of the photospheric magnetic field as it reaches a dynamic steady state [55]. Most MHD models, e.g. MAS, solve for the ambient solar wind via time-dependent simulations and allow the solution to relax to steady state. While customary, this results in a large run time to compute the steady state. There are other MHD models, such as [47], that solve directly for the steady solar wind, yet they have their own set of numerical challenges as it is no longer straightforward to extract the physical variables from the flux variables.

Numerical Solver
The MAS model equations are numerically solved on a nonuniform logically-rectangular staggered grid using finite differences. The nonuniform mesh allows for adjustment of the grid point concentration based on transition and active regions. For more details about the numerical methods and their stability, see [32,14]. The MAS model divides its computational domain to two distinct regions: the corona and heliosphere. The corona is the region between 1R S to 30R S and the heliosphere is the region between 30R S to 1.1AU. The value R S denotes solar radii unit of distance which is 695, 700km, and 1/215th of an astronomical unit (AU), which is equal to the distance from the Sun to Earth. Numerical results and implementation details are discussed in Section 4.2.

Solar Wind Speed: The HUX Model
The HUX (Heliospheric Upwinding eXtrapolation) model developed by [56,52] is a two-dimensional time-stationary model that predicts the heliospheric solar wind speed. The HUX model has been incorporated into operational and ensemble-based space weather programs [29,5] as an MHD surrogate model to study retrospective time periods as well as real-time predictions. It has also been used to map streams directly from in-situ spacecraft observations (e.g. Helios A/B) to Earth [25]. The HUX model is based on simplified physical assumptions of the fluid momentum equation. In contrast to the MAS model, where the velocity field is solved via the MHD equations, the HUX model constructs a kinematic mapping where each plasma parcel speed is governed by its adjacent parcel's speed.
We introduce the HUX model as a reduced-physics surrogate solar wind model that is capable of capturing the solar wind speed as it propagates away from the Sun. The HUX model shares many similarities with the MAS model, such as advection-dominated solutions, and can therefore be used to test our proposed method, shifted operator inference. Moreover, since both sOpInf and HUX are surrogate models to the MAS model, we will compare their capability to approximate the MAS results. This section describes the HUX governing equations along with their spatial discretization and implementation.

Reduced-Physics Equation
The HUX model [56,61] is derived from the fluid momentum equation in the corotating frame of reference with the Sun by considering Eq. (5) in the absence of magnetic and viscous effects describing steady flow by replacing the time derivative ( ∂ ∂t ) with a spatial derivative (−Ω rot ∂ ∂φ ), i.e. the governing equations are where v = [v r (r, θ, φ), v θ (r, θ, φ), v φ (r, θ, φ)] is the solar wind proton velocity, ρ(r, θ, φ) is the plasma density, and g(r) is the gravitational acceleration specified in Section 2.1.1. The term Ω rot (θ) = 2π 25.38 − 2.77π 180 cos(θ − π 2 ) 2 is the angular frequency of the Sun's rotation, i.e., a function of latitude [52]. The analysis in [56,52,25] justifies neglecting the pressure gradient and gravity terms in Eq. (7) and only taking into account variations of the velocity in the radial direction. As a result, Eq. (7) reduces to the two-dimensional nonlinear scalar homogeneous time-stationary equation where the independent variables are r and φ and the dependent variable is the velocity in the radial direction v r (r, φ). The angular frequency of the Sun's rotation is evaluated at a constant Carrington latitudeθ; here we consider the Sun's equatorial plane (θ = 0) so that Ω rot (0) = 2π 25. 38 1/days at the solar equator. The initial-boundary value problem (IBVP) is subject to the initial condition v r (r 0 , φ) = v r0 (φ) and is defined on the periodic domain 0 ≤ φ ≤ 2π and r ≥ 30R S , where beyond 30R S , the solar wind travels along roughly radial trajectories, justifying the assumption of only considering the velocity in the radial direction. Additionally, to account for the residual acceleration present in the inner heliosphere, the authors in [56] suggested adding an acceleration boost to the initial velocity profile described by where v r0 (φ) is the initial radial velocity, α = 0.15 is the acceleration factor, and r h = 50R S is the radial location at which the acceleration ends. Hence, the acceleration boost, v acc (r 0 , v r0 (φ)), is added to the initial velocity profile v r0 (φ) prior to solving the HUX equation.

Discretization via the Upwind Scheme
This section describes the semi-discretization of Eq. (8) in longitude, which then results in a set of ODEs. To begin, we rewrite Eq. (8) in the hyperbolic conservation form where the physical flux function is We use the first-order conservative upwind method from [25], so to approximate the partial derivative of the flux function f with respect to φ by where n φ is the number of mesh points in longitude and j = 1, 2, . . . , n φ denotes the longitude grid index. We discretize the longitudinal direction uniformly with ∆φ mesh spacing and denote the discretized state vector as v(r) = [v r (r, φ (1) ), v r (r, φ (2) ), . . . , v r (r, φ (n φ ) )] ∈ R n φ . From here, we obtain the semi-discrete system of ordinary differential equations with the sparse matrix The initial condition v(r 0 ) ∈ R n φ is set as the MAS coronal solution, v MAS ∈ R n φ , evaluated at r 0 = 30R S along with adding the ad hoc acceleration boost described in Eq. (9), i.e.,

Shifted Operator Inference: A Non-intrusive Reduced-Order Model Approach for Advection-Dominated Systems
This section proposes shifted operator inference (sOpInf ), a non-intrusive data-driven modeling framework that expands standard operator inference [44] towards the challenge of modeling solar wind, and more generally, advection-dominated systems. The proposed method ensures that the learned ROM is able to (i) capture the dynamics of the translational systems with only a few modes, (ii) retain the translation and rotation properties of the physical system, and (iii) accurately predict the shift velocity in the testing regime. This is done by transforming the original coordinates to a moving coordinate frame where the dynamics are absent of translation and rotation. Section 3.1 describes the new method that leverages a coordinate shift to first transform the data and subsequently learn in the transformed coordinates. Section 3.2 proposes two alternative strategies for deriving this coordinate transformation. Lastly, Section 3.3 illustrates the sOpInf methodology on an introductory example: the one-dimensional inviscid Burgers' equation.

Shifted Operator Inference for Advection-Dominated Systems
To illustrate the proposed methodology, we consider a generic k−dimensional time-dependent partial differential equation (PDE) for the scalar function u(x 1 , x 2 , . . . , x k , t) of the form where x 1 , x 2 , . . . , x k ∈ R denote the spatial coordinates and t ∈ R + . One may think of t as time, or, as described in the previous section, we will also consider the independent variable to be the radial distance from the Sun, r. The function F defines the equations of motion of the system which include advective terms. For simplicity, we focused the illustration on a purely convective case, but our method can account for higher derivatives (e.g. diffusion terms) as well, see Section 4 where we consider a model with viscous forces. Our goal is to derive a data-driven ROM that can accurately predict the solutions to advectiondominated systems, i.e., that uses data of a semi-discretization of (15) and produces a predictive and efficient low-dimensional model. The proposed method proceeds in four steps as outlined next.
(I) Data collection and translation. The system (15) is in k dimensions and is typically solved via a spatial discretization scheme at fixed spatial locations which we refer to as the original coordinates. These are assembled into the following spatial grid where n x = Π k i=1 n i is the number of spatial grid points and n = k · n x . Depending on the context, using X in matrix form or x in vector form may be more preferred. We collect data (for instance, solar wind speed data) from the numerical solver in the original coordinates at instances t i , i.e., u i ≈ u(x, t i ) ∈ R n with 0 = t 0 < t 1 < · · · < t K = T . To account for the translational element in the data, we next shift each snapshot to a moving coordinate frame wherex(x, t i ) denotes the moving coordinate frame and c(t) ∈ R n represents the traveling wave speed. We evaluateũ i via piecewise linear interpolation, i.e., where P k i denotes the k > 1 dimensional piecewise linear interpolant of u i in the original grid x evaluated on the moving coordinate frame gridx(x, t i ). The piecewise linear interpolation is implemented in the Python package scipy under the function scipy.interpolate.LinearNDInterpolator(). Our method is not restricted to the type of interpolation, so higher order interpolation methods can be used as well (e.g. piecewise cubic interpolation). In the moving coordinate frame, the system no longer exhibits translational properties. Section 3.2 proposes two techniques to determine c(t). We then store the transformed data in the matrix where in the applications that we consider, n K, so the matrixŨ is tall and skinny.
(II) Data reduction via projection. Given high-dimensional data, we first identify the low-dimensional subspace in which to learn a ROM. In this work, we use the subspace spanned by the proper orthogonal decomposition (POD) modes [23], which is obtained by computing the economy-sized singular value decomposition of the snapshot matrix, i.e.,Ũ where , is given by the first columns of V. The basis dimensions can be chosen based on the cumulative energy criteria, e.g., where σ i = Σ ii are the singular values, and tol is commonly chosen to be tol = 0.95 or tol = 0.99 which encapsulate 95% and 99% of the energy in the data, respectively. Next, we project the state snapshot data onto the POD subspace spanned by the columns of V and obtain the reduced snapshot matrices where the columns of˙ U are computed from U using any time derivative approximation (see, e.g., [35,27,15]), or can be obtained-if available-by evaluating the right-hand-side of the governing equation (the residual) and projecting the resulting data.
(III) Model learning and prediction via operator inference. In this section we start with the assumption that the finite-dimensional data-generating model (a spatial discretization of Eq. (15) in the moving coordinate frame) is of the form of a polynomial nonlinear system of ODEs, written as with matrices A ∈ R n×n , H ∈ R n×n 2 and C ∈ R n×n 3 . The column-wise Kronecker product is denoted by ⊗. Boundary conditions can either be represented via the constant B ∈ R n or time-dependent BCs as Bη(t). The abbreviation "HOT" in Eq. (20) denotes higher-order terms, and represents terms that are quartic and higher order. For example, in the case of Burgers' equation in Section 3.3 the term u ∂u ∂x is quadratic in the PDE state u(x, t) and would therefore yield a discretized component H(ũ ⊗ũ).
Approximating the high-dimensional stateũ in a low-dimensional basis V ∈ R n× , with n, we writeũ ≈ V û. Using a Galerkin projection, this yields the ROM of Eq. (20) as where ⊗ is the compact Kronecker product, which removes redundant terms in the standard Kronecker product ⊗. For example, forû = [û 1 ,û 2 ] the standard Kronecker product yieldsû ⊗û = [û 2 1 ,û 1û2 ,û 2û1 ,û 2 2 ] and the compact Kronecker product yieldsû ⊗ û = [û 2 1 ,û 1û2 ,û 2 2 ] , which uses only unique terms. For here on, we only use the compact Kronecker product for learned sOpInf ROMs. Consequently, the ROM operators and their dimensions are We note again that projection preserves polynomial structure, that is, Eq. (21) has the same polynomial form as Eq. (20), but in the reduced subspace defined by V .
To simplify notation, we continue from now on with a quadratic system, but note that all results carry over directly to cubic, quartic and all higher-order polynomial terms. Nevertheless, we note that the number of elements in the ROM operators scales with 4 for the cubic operator, 5 for the quartic operator, etc., yet higher-order terms often exhibit significant block-sparsity that can be exploited in numerical implementations which limits the growth of computational cost to solve the ROM. For terms in the governing equations that are not in polynomial form, the introduction of variable transformations and auxiliary variables via the process of lifting [20,28,62,48] can convert these terms to polynomial form.
The goal at this stage is to learn a ROM that evolves the shifted and projected snapshots in time. Operator inference solves a least-squares problem to find the reduced operators that yield the ROM that best matches the projected snapshot data in a minimum residual sense. For a quadratic ROM (with C and HOT set to zero in Eq. (21)) operator inference solves the least-squares problem where 1 K ∈ R K is the length K row vector with all entries set to one. Note that this least-squares problem is linear in the coefficients of the unknown ROM operators A, H, and B. The appeal of the operator inference approach comes from the ability to compute the ROM operators A, H, and B directly from data without needing explicit access to the original high-dimensional operators A, H, and B. The unknown operators and known low-dimensional data are combined in the matrices respectively. The unknown operators are then obtained as a solution to the minimization problem For K > + 1 2 ( + 1) + 1 this overdetermined linear least-squares problem has a unique solution [18,Sec. 5.3]. It follows from linear algebra (and is noted in [46]) that Eq. (22) can be written as independent least-squares problems, each of the form and r i is a column of˙ U . This makes the operator inference approach efficient and scalable.
To avoid overfitting and prevent potential instability of the learned ROMs, regularization becomes necessary, see [36] for a detailed regularization study of operator inference. In this work, we use an Tikhonov regularization penalty so that the least-squares problem becomes is the diagonal matrix used for regularization. The parameter λ 1 is the regularization parameter of the operators B ∈ R and A ∈ R × and λ 2 regularizes the operator H ∈ R × 1 2 ( +1) . The regularization parameters λ 1 and λ 2 are problem specific and should be chosen accordingly. We provide details in Section 4 and refer to [62, Sec. IV.B] for more implementation details of operator inference.
Having learned the ROM in Eq. (21) from the shifted data, we can make efficient predictions in that low-dimensional subspace that go beyond the training data into the fully predictive regime. We thus simulate Eq. (21) to obtain a solutionû(t), which we then lift to n dimensions to get the approximate ROM solutionũ ROM (t) = V û(t) ≈ũ(t). We use the Operator Inference Python package version 1.2.1 [4] to implement the model learning and prediction step of sOpInf.
(IV) Re-shifting predicted ROM data. The predicted ROM solutionsũ ROM (t) will be in the moving coordinate frame and require reverse translation to the original coordinates. We shift the ROMpredicted solutions back to the original coordinate system via interpolation where P k i denotes the k > 1 dimensional (here: piecewise linear) interpolant ofũ ROM i in the moving coordinate frame gridx(x, t i ) evaluated on the original grid x; see part (I ) in Section 3.1 for more details about the interpolation implementation. After shifting back to the original coordinates, the predicted and reconstructed sOpInf snapshots are columns of the matrix U ROM ∈ R n×(K+m) , such that t K+m > t K , which is the final output of the algorithm.
The previous steps (I )-(IV ) are summarized in Algorithm 1, which is written for a quadratic system; yet extensions to cubic, quartic, and other polynomial systems are straightforward.

Determination of Spatial Shift Velocity
There are various ways in which the traveling wave speed, c(t), can be discovered. For instance, the authors in [37] use sparse regression and spectral clustering to uncover the function c(t). We present two methods: the method of characteristics discussed in Section 3.2.1, an analytic approach that requires knowledge of the underlying equations (but not the discretization or computer code), and the cross-correlation extrapolation method described in Section 3.2.2, a purely data-driven approach. In both cases, the shift function c(t) is learned from the batch of training data and extrapolated in the testing regime.

Method of Characteristics
The method of characteristics can be applied to quasi-linear partial differential equations, in which along the characteristic curves the PDE can be transformed to a set of coupled ODEs. We exploit the method of characteristics to find the scalar shift function c(t) for first-order PDEs describing the scalar quantity u( whose dynamics are described by the following quasi-linear PDE: subject to the initial condition u(x 1 , . . . , x k , t = 0) = u 0 (x 1 , . . . , x k ) with periodic boundary conditions at each spatial coordinate boundaries x i = a i and x i = b i for i = 1, 2, . . . , k. The transport speed , t] ∈ R is the source term. Then, by the method of characteristics, Eq. (24) can be written as a system of k + 1 ODEs, i.e.
We solve Eq. (25a) first, which results in where G [u 0 (x 1 (0), . . . , x k (0)), t] ∈ R describes the scalar quantity u along the characteristic curves. Given this ansatz, we proceed to obtain the characteristic curves by solving Eq. (25b) via separation of variables, such that For most equations that arise from conservation laws, the ODEs in Eq. (25a)-(25b) can be solved analytically via Eq. (26)- (27). In the case when Eq. (25a) can not be solved analytically, the coupled system of ODEs in Eq. (25a)-(25b) can be solved numerically on a discrete grid; or in the case when f i [G, t] is a nonelementary antiderivative, we can approximate the function f i using Taylor series, and integrate term-by-term. By following the characteristic curves described in Eq. (27) we are able to discover a moving coordinate frame absent of advection. In fluid dynamics, the characteristic paths are referred to as the Lagrangian specification of the flow field, whereby the fluid motion is observed following an individual fluid parcel. We construct a moving coordinate frame, x + c(t), that resembles the main direction of the Lagrangian frame of reference. The shift function c i (t) corresponding to the spatial coordinate x i can be obtained by computing the mean characteristic emerging from a certain spatial domain. Let the spatial domain of x i , for i = 1, 2, . . . , k be discretized on a grid with n i points such that x i = [x (1) i , x i ] with 1 ≤ p < q ≤ n i is given as The spatial interval [x i ] can be set to the entire spatial domain, yet it is usually set to be a specific region of interest.
For problems with shock formation, the shift function c(t) is computed via the mean characteristic curve only before shock formation, i.e. before the characteristic lines first intersect, and approximated via the shock curve after shock formation, i.e.
where a = 1 |x x 1 (0), . . . , x k (0)), t] , t] dt dx i . The time of shock formation is denoted by t s , and s i (t) is the shock trajectory in x i coordinate, see Section 3.3 for a one-dimensional inviscid Burgers' equation example. For one-dimensional problems where g = 0, the shock trajectory can be computed via the entropy (Rankine-Hugoniot) condition or the Whitham's geometric equal area rule in which multi-valued regions of the solution are replaced with a discontinuity that satisfies conservation [66]. The shock location is approximated by locating a vertical line that splits the multi-valued curve into two regions with equal area.
In problems where there are multiple shock curves, such as shown in Section 4.3, we suggest to approximate c(t) by tracking only one shock curve. The selection of the shock is problem dependent, e.g. in Section 4.3 we choose to follow the first shock emerging. Although, if the choice of shock curve or spatial interval [x (p) , x (q) ] is ambiguous, computing the shock curve is unfeasible or computationally expensive, or the initial condition is noisy, we recommend to use the cross-correlation extrapolation technique (Section 3.2.2) to find c(t). Additionally, it can be non-trivial to find such characteristics in the case of more complex hyperbolic PDEs that are not in the form of Eq. (24), e.g. coupled systems such as Eqs. (1)-(6) with multiple dependent variables resulting in more than one characteristic curve emanating from a single spatial point, in which we recommend employing the data-driven crosscorrelation extrapolation technique, which we present next.

Cross-Correlation Extrapolation Method
Cross-correlation is a mathematical operation that is commonly used in signal processing and pattern recognition to measure similarity of two signals. For two finite discrete signals f , g ∈ C n , the univariate discrete circular cross-correlation is defined as where f denotes the complex conjugate of f , the bracket [j] denotes the jth element of the signal, and τ ∈ Z is the discrete displacement. The discrete circular cross-correlation can be extended to the multi-variate case, for snapshots with k ∈ Z variables and tensor-valued f , g ∈ C n1×n2×...×n k : where τ = [τ 1 , τ 2 , . . . , τ k ] ∈ Z k is the multi-variate discrete displacement. When the signals correlate, the value of f g is maximized.
We propose to find the optimal discrete displacement, τ * ∈ Z k , by maximizing the cross-correlation between the two discrete signals (or snapshots), which amounts to solving Once the shift is computed for all training snapshots by applying the circular discrete cross-correlation between each snapshot u i ∈ R n and the initial condition u 0 ∈ R n , we obtain the the multivariate discrete displacement for each time-step, i.e., Then, the shift function c j (t), ∀j ∈ {1, 2, . . . , k} is found via least squares polynomial curve fitting to the data points between the time increments and corresponding spatial location of the shift, such that the shift function is of the form where d is the degree of the polynomial approximation. To find the vector of real coefficients a = [a 0 , a 1 , . . . , a d ] ∈ R d+1 , we solve the minimization problem ] ∈ R K is the vector of corresponding spatial location of the shift, and the Vandermonde matrix T is defined as To make predictions outside the training interval, we approximate the shift by polynomial extrapolation of c(t) = [c 1 (t), c 2 (t), . . . , c k (t)] ∈ R k . As an illustration, Figure 1 shows the bi-variate discrete circular cross-correlation applied to the MAS CR2210 snapshot at 30R S (the initial condition) and at 1AU. Here, since the flow is steady, the independent variable is the radial distance from the Sun, r, (instead of time t). For this case, the convective shift is 45 • in longitude and 0 • in latitude. Figure 1 confirms that the translation in the solar wind is purely longitudinal due to the rotation of the Sun. The main idea behind the cross-correlation extrapolation technique is similar to the template-fitting technique studied in [59,58] where the data is periodic and the template is set to be the initial condition.

Illustrative Example: Shifted Operator Inference for the Inviscid Burgers' Equation
The one-dimensional inviscid Burgers' equation is of the form of Eq. (24) describing the scalar quantity u(x, t) : subject to the initial condition u(x, t = 0) = u 0 (x) with appropriate boundary conditions at x = a and x = b. In the moving coordinate frame defined bỹ x(x, t) = x + c(t) and u(x, t) =ũ(x(x, t), t) Burgers' equation (35) becomes This can be written in conservative form as Then, by the conservative first-order upwind scheme, we approximate the spatial derivative by where j = 1, 2, . . . , n denotes the grid index inx. Based on the discretization scheme in Eq. (39), we can write the dynamics of Eq. (38) in vector form as where ⊗ denotes the Kronecker product andũ(t) = [ũ(x (1) , t),ũ(x (2) , t), . . . ,ũ(x (n) , t)] ∈ R n denotes the state vector discretized over n spatial points at time t. Here, H ∈ R n×n 2 is the quadratic operator that corresponds to the discrete term 1 2∆x ũ(x (j) , t) 2 −ũ(x (j−1) , t) 2 from Eq. (39). Moreover, A(t) ∈ R n×n is the linear time-dependent operator corresponding to the discrete term 1 ∆x dc dt ũ(x (j) , t) −ũ(x (j−1) , t) in the moving coordinate frame. In the case where c(t) ∈ R is a linear function (so dc dt =const.), meaning the wave is traveling at constant speed, the linear operator is time independent, i.e., A(t) ≡ A.
The traveling wave speed c(t) can be estimated by the method of characteristics (see Section 3.2.1), where we can rewrite Eq. (35) as the two coupled ODEs Hence, along the characteristic lines the quantity u(x(t), t) remains constant, which can be verified by Then, by integrating Eq. (41) the characteristic curves are linear before shock formation. Let the spatial domain of x be discretized on a uniform grid with n points such that x = [x (1) , x (2) , . . . , x (n) ] ∈ R n . From here, we can approximate the shift function via Eq. (30) as a piecewise continuous function: where t s ∈ R is the time of shock formation, a = 1 q−p q j=p u(x (j) , t = 0)t s , and s(t) ∈ R is the shock trajectory.
To demonstrate the sOpInf method, we now consider a specific spatial domain x ∈ [0, 3] and time domain t ∈ [0, 2] along with a Gaussian initial condition u 0 (x) = 0.8 + 0.5 exp (−(x − 1) 2 /0.1) and periodic boundary conditions. Equation (35) in the regular coordinates is solved numerically via the forward Euler method in time and the conservative upwind scheme in space on an equidistant computational grid with 500 discretization points in space, and 1000 points in time. With this choice, the CFL condition u(x, t) ∆t ∆x ≤ 1 is satisfied, where ∆x and ∆t denote the grid spacing in x and t, respectively. The training dataset consists of 80% of the snapshots and the testing dataset consists of 20% of the snapshots. The shock is formed at t s = min  resulting in c(t) = 1.05t. Whereas, via the method of characteristics (Section 3.2.1) the shift function, c(t) described in Eq. (43), is linear only before shock formation, where we compute the mean of characteristics emanating from the spatial interval [1.15, 1.4] before shock formation, see the green characteristic curves in Figure 2e. By the Rankine-Hugoniot jump condition [66], the shock curve s(t), and hence the shift function after shock formation, is non-linear, see the red curve in Figure 2e. However, the absolute difference between the shift function c(t) obtained by the method of characteristics vs. the cross-correlation extrapolation is less than 10 −2 for all t ∈ [0, 2], indicating there is not a significant difference between the two methods, see Figure 2c for a visual comparison. With that, we found that using the method of characteristics to find c(t) resulted in up to 2% more accurate ROM results in comparison to the results obtained by the cross-correlation extrapolation linear shift function. By aligning the wave discontinuity at one spatial point, the ROM modes can better approximate the shock. We continue the Burgers' ROM analysis using the shift function computed by the method of characteristics. Figure 2a shows the inviscid Burgers' equation results on the original coordinate system (x, t) while Figure 2b shows the inviscid Burgers' equation results on the shifted coordinates (x, t). It is apparent from Figure 2b that the dynamics are largely absent of translation in the shifted coordinate frame. The evolution of the initial wave depends solely on shock formation as the discontinuity in the wave sharpens. Figure 2d compares the singular value decay of the solution data matrix U on the original coordinates and ofŨ on the shifted coordinates. For the former, the singular value decay is very slow due to the translation properties in the system; a ROM in those coordinates would require a large number of basis functions. In contrast, once the dynamics are absent of translation on the shifted coordinates the system can be approximated using only a few modes. Figure 3 illustrates the predicted solutions from sOpInf compared to the FOM results for the inviscid Burgers' equation. The sOpInf model is of the formu =Âû +Ĥ(û ⊗ û) with = 9 modes. The regularization coefficients of the operatorsÂ ∈ R × andĤ ∈ R × 1 2 ( +1) are λ 1 = 1 and λ 2 = 1, respectively. The sOpInf results in 0.9999958 Pearson correlation coefficient (PCC) and 1.204 × 10 −4 mean relative error (MRE) in comparison to the inviscid Burgers' results, indicating that the sOpInf model is able to capture well the evolution of the high-fidelity Burgers' simulation with only a few modes. On the shifted coordinates we are able to construct a ROM with only 9 modes, meanwhile, in the original coordinates, we would need 89 modes to achieve the same projection error. This demonstrates that the shifting procedure produces significant computational speedups.
We next assess how the amount of training data affects the accuracy of the resulting ROM by training sOpInf on 50%, 60%, 70%, and 80% of the total snapshots, see Figure 4. The numerical results show that-as anticipated-increasing the amount of training data improves the ROMs accuracy. Additionally, the relative error measured in the L 2 -norm is bounded by 10% for all four simulations, resulting in an adequate ROM using as low as 50% of the total snapshots for training. To analyze the framework's robustness to noise, we tested the sOpInf methodology on the inviscid Burgers' simulated data with added noise. The results are presented in Appendix A and demonstrate that the method seems to be robust to the addition of moderate levels of Gaussian noise.

Numerical Results for MAS and HUX Solar Wind Models
We apply the sOpInf methodology to learn low-dimensional models for the ambient solar wind radial velocity predicted by the HUX and MAS models described in Section 2. We consider a specific event of interest and give some background on that event in Section 4.1. In Section 4.2, we give details on the data and implementation. Section 4.3 and Section 4.4 present the sOpInf twodimensional steady results trained on HUX and MAS equatorial plane data, respectively. Section 4.5 presents sOpInf three-dimensional steady results trained on the full-Sun MAS results. Section 4.6 compares sOpInf and HUX reconstruction of MAS equatorial plane streamlines. Lastly, Section 4.7 discusses the choice of ROM model form. The public repository https://github.com/opaliss/ Space-Weather-ROM contains a collection of Jupyter notebooks in Python 3.9 containing the code and data used in this study.  The sensitivity of sOpInf ROM to the amount of training data with 50%, 60%, 70%, and 80% of the total snapshots used for training. The relative error in the L 2 -norm is shown as a function of time. The numerical results show that the sOpInf methodology is able to learn adequate ROMs with as low as 50% of the total snapshots utilized for training. The solar wind radial velocity at the equator has three main peaks, originating from coronal holes, which make this period a great candidate to study large-scale structure in the solar wind.

Physical Relevance of Data
We focus on Carrington Rotation (CR) 2210, which occurred from 26 October to 23 November 2018, during a solar minimum. Figure 5a shows the heliospheric solar wind radial velocity MAS results for CR2210 in the latitude region of θ ∈ [−30 • , 30 • ] and Figure 5b presents the equatorial solar wind radial velocity profile. As observed in Figure 5, during CR2210, the solar wind radial velocity at the equator has three main peaks, essentially making this period a great candidate to study large-scale structure in the solar wind. Although not shown here, the origin of the fast wind (ranging from 500 to 700 km/s) is from equatorial coronal holes located at approximately 180 • − 210 • , 280 • − 300 • , and 330 • − 360 • in longitude at the source surface (see [8, Fig. 8(a)] for a synoptic view at 2R S of the coronal hole regions). Although fast streams commonly originate in large coronal holes, slow streams come from various coronal sources (e.g., coronal hole boundaries, coronal loops, etc.). Another reason to consider this data is that this Carrington period is well-studied in literature due to Parker Solar Probe (PSP) reaching its first perihelion pass of 35.7R S after its first Venus gravity assist on 6 November 2018, as it broke records by becoming the closest spacecraft to the Sun. Moreover, the authors in [57] showed that the thermospheric MAS solar-wind speed results highly match the observations made by PSP during CR2210. Thus, while we do not present a comparison with in-situ solar-wind observations we can treat the MAS results as a physically meaningful representation of the solar wind.

Data and Implementation Details 4.2.1. MAS Data
The MAS model is described in detail in Section 2.1. The MAS model solar wind velocity results are time-stationary in spherical coordinates. The data covers the entire domain in longitude 0 • ≤ φ ≤ 360 • , latitude −90 • ≤ θ ≤ 90 • , and radial axis (a.k.a. heliocentric distance) 0.14AU ≤ r ≤ 1.1AU. The MAS simulation results are on a rectangular grid with n φ = 128 uniformly spaced points in Carrington longitude, n θ = 111 uniformly spaced points in heliographic latitude, and n r = 140 points on a non-uniformly spaced grid in the radial axis. Hence, the snapshot data dimension is n x = n φ × n θ = 14, 208 where we treat r as the independent variable. A convergence check with a higher-resolution grid found that this resolution is sufficient for physical accuracy. The MAS coronal and heliospheric models (implemented in FORTRAN) were run with medium and high resolution by the authors of [1], which took approximately four and 28 hours of wall-clock time using four NVIDIA RTX 2080Ti GPUs, respectively [13]. We derive data-driven ROMs from the medium-resolution MAS heliospheric simulation which is publicly available at PSI's web page [2]. The MAS training datasets contain 70% of the snapshots, with 98 and 42 snapshots for training and testing, respectively. The training domain is from 0.14 to 0.82AU and the testing domain is from 0.82 to 1.1AU.

HUX Data
The HUX model is described in detail in Section 2.2. To simulate the HUX model at the heliographic solar equatorial plane we numerically solve the ODE in Eq. (12) via the forward explicit Euler's method in time and setting the initial condition to the MAS velocity profile results at 30R S ≈ 0.14AU. The HUX dataset is two-dimensional with n φ = 128 and n r = 140 (satisfying the CFL condition) on the same domain as listed above for the MAS results. On the equatorial plane (θ e.p = 0), the angular frequency of the Sun is Ω rot (θ e.p ) = 2π 25.38 1/days. The HUX dataset took 0.03s to simulate on a MacBook Pro 2.3 GHz Quad-Core Intel Core i7 processor with 16 GB RAM. We note that the HUX model is already computationally efficient. We train a ROM for HUX merely as a introductory example since it produces the same physics (advection-dominated solutions) that we expect in much more expensive solvers, and not to show any computational improvements. The HUX training and testing domains are identical to the MAS equatorial training domain, see details in Section 4.2.1.

ROM Implementation
The MAS and HUX numerical results are on a non-uniformly spaced grid in the radial axis, consequently, we approximate the derivatives of the training data with respect to r, i.e. d drû in Eq. (19), via a second-order accurate central difference in the interior points and second-order accurate onesided (forward and backward) difference at the boundaries. For the case of uniform grid meshing, e.g. Burgers' equation ROM presented in Section 3.3, we used a sixth-order finite difference scheme. We simulate the ROM via an implicit multi-step variable method based on a backward differentiation formula using the scipy.integrate.solve_ivp() Python function. The ROM regularization coefficients, λ 1 and λ 2 , are chosen from the logarithmically spaced set, i.e. {10 0 , 10 1 , 10 2 , . . . , 10 10 }, such that the best coefficients minimizes the relative error measured via the L ∞ -norm over the training regime.

HUX Equatorial Plane Numerical Results
We apply the sOpInf framework to learn a ROM of the HUX CR2210 equatorial plane radial velocity. As derived in Section 2.2.1, the HUX underlying equation is where r, φ are the independent variables. To begin, we shift the HUX dynamics to a moving coordinate frame defined byφ (r, φ) = φ + c(r) and v r (r, φ) =ṽ r (φ(r, φ), r).
The shift function c(r) can be learned via either the method of characteristics (Section 3.2.1) or the cross-correlation extrapolation method (Section 3.2.2), in particular, the circular uni-variate crosscorrelation method described in Eq. (31). The linear-fit cross-correlation shift function resulted in c(r) = −51.711 • r + 7.032 • . Our numerical studies found that there is no substantial differences between the numerical results of the two methods. Following the steps described in Section 3.2.1, the HUX characteristic curves are derived by the following two coupled ODEs: Then, by integration of Eq. (44), the characteristics before shock formation are straight lines described by where r 0 = 0.14AU. The HUX characteristic curves are also called the ballistic approximation, which assumes that each spiral field line of plasma continues at a constant speed throughout the heliosphere [52]. After obtaining the characteristic curves, we are able to approximate the shift function by if r > r s so that r s is the radial position where the characteristics first intersect, s(r) is the shock trajectory, and a = 1 . The indices p, q can include the whole spatial domain or instead a bounded spatial interval to track specific regions of the initial wave. For CR2210, we limit p, q to include the characteristics emanating from a main equatorial high solar wind peak, originating from an equatorial coronal hole, on the longitudinal interval [180 • , 260 • ]. There are three main shock curves, and we choose to follow the first shock curve that emerged at r s = 0.344AU and φ s = 211.737 • . The learned sOpInf ROM is of the forṁ and it is able to sufficiently model the dynamics of the HUX model with only = 4 modes. A practical implementation of operator inference requires regularization, and for the least-squares fitting we found λ = 10 3 to give good results. The comparison between sOpInf and HUX velocity profiles is provided in Figure 6. The figure shows that the advective solutions are well approximated both in the training regime until r = 0.82AU and in the testing regime, where the ROM is fully predictive. This conclusion is also supported by Table 1 where we provide the mean/median/maximum relative error and the Pearson correlation coefficient comparing the HUX solutions with the ROM solutions, both in the testing and training regime. The error measures show that a sOpInf ROM can sufficiently predict the HUX dynamics while reducing the dimensionality of the problem from n φ = 128 to = 4, i.e., a factor of 32 reduction of state-space dimension.

MAS Equatorial Plane Numerical Results
We apply the sOpInf framework to learn a ROM for the MAS CR2210 equatorial plane velocity field. Since the MAS Eqs. (1)- (6) are not in the form of Eq. (24), the method of characteristics can not be applied to approximate the shift function c(r); instead, we use the cross-correlation extrapolation method, which resulted in c(r) = −54.98 • r + 7.39 • , where r is measured in AU. Figure 7 shows the MAS equatorial plane heliospheric results on the original and shifted polar coordinates along with the cumulative singular value energy criteria described in Eq. (18) in each coordinate system. Evidently, we see that shifting the snapshots to a moving coordinate frame creates a faster singular value decay, leading to an accurate representation of the shifted data with far fewer modes. This observation is in line with the results we showed for Burgers' equation in Figure 2d above.
The results shown in this section are for the ROM model forṁ v =Ĥ(v ⊗ v),Ĥ ∈ R × 1 2 ( +1) Figure 6: The HUX solar wind radial velocity results at the heliographic equatorial plane for CR2210 along with sOpInf quadratic ROM results. The sOpInf ROM aligns very well with the data past the training interval; the snapshots at r = 0.91AU and r = 1.1AU are testing data where the ROM is fully predictive. with = 9 modes. This model form provided the best overall results in training and extrapolation, compared to other model combinations including linear and constant terms, see Section 4.7 for further discussion. The regularization coefficient for computing the operatorĤ is λ = 10 5 . Figure 8 visually demonstrates that sOpInf is capable of accurately approximating the MAS equatorial results, where the snapshots at r = 0.91AU and r = 1.1AU are in the fully predictive regime of the ROM. We highlight that the MAS data contains more complex features than the HUX data, specifically in the region 0 • ≤ φ 120 • where small localized wave structures exist. Nevertheless, sOpInf covers those equally well as the HUX data in the previous section. Figure 9 presents are more qualitative assessment of the relative error between the sOpInf ROM and MAS velocity fields, which shows that the relative error is less than 8.3% in the entire domain. The mean/median/maximum relative error and PCC in the testing and training regime are again provided in Table 1 above.

3-D Full-Sun MAS Numerical Results
We showcase sOpInf trained on MAS three-dimensional steady-state velocity results. As mentioned previously, the MAS Eqs. (1)- (6) are not in the form of Eq. (24), hence, the method of characteristics is not a valid choice in approximating the shift function c(r) ∈ R 2 . Therefore, we turn to the crosscorrelation extrapolation method described in Section 3.2.2. More specifically, the full-Sun MAS   snapshots in matrix form, v(r) ∈ R n φ ×n θ , are bi-variate in Carrington longitude (φ) and latitude (θ), with n φ , n θ mesh points in longitude and latitude, respectively. The shift function c(r) ∈ R 2 is found via the bi-variate circular cross-correlation defined by Eqs. (32)- (33) with k = 2. The bivariate circular cross-correlation is applied between each snapshot and the velocity profile at the initial condition (30R S ). Since the translation in the MAS solar wind results is due to the Sun's rotation, the shift function can be expressed as c(r) = c(r)ê φ , whereê φ is the unit vector in longitude direction, and c(r) ∈ R is the shift in longitude, where we computed c(r) = −52.44 • r + 6.94 • with r measured in AU. The results in this section are for a ROM of the forṁ v =Âv +Ĥ(v ⊗ v) +B whereB ∈ R ,Â ∈ R × ,Ĥ ∈ R × 1 2 ( +1) with only = 8 modes. The regularization coefficient for computing the operatorB andÂ is λ 1 = 10 4 and the regularization coefficient forĤ is λ 2 = 10 8 . Figure 10a shows the relative error between the two velocity fields and the mean/median/max relative error in the testing and training regime is shown in Table 1. For a visual comparison, Figure 10b shows a comparison between the sOpInf reconstructed and predicted two-dimensional full-Sun snapshots; the visual comparison and error estimates indicate that sOpInf can successfully reproduce the high fidelity full-Sun MAS dataset, where n x = n φ × n θ = 14, 208 with only = 8 modes, leading to a substantial reduction in the model's dimensionality.
In practice, the sOpInf framework can be employed to speed up the MAS computational time by setting the MAS heliospheric outer boundary condition to 0.82AU (which is 70% of the current computational domain) and run the time-dependent MAS simulation until it relaxes to steady-state. Then, the steady-state MAS snapshots are used as training data for sOpInf, which takes 0.355 seconds to simulate from 30R s up to 1.1AU on a MacBook Pro 2.3 GHz Quad-Core Intel Core i7 processor with 16 GB RAM (for 3-D full-Sun simulation). If we assume that running MAS on 70% of the computational domain would take 70% of the MAS current computational time, then we would be able to speed up the MAS computational time by approximately 1.2 hours and 8.4 hours for the medium and high resolution runs, respectively (i.e. save 30% of the MAS computational time). Another important sOpInf speed-up contribution can be in the case when one is interested in studying the solar wind dynamics much further in the heliosphere (e.g. conditions in the vicinity of Mars or Jupiter). In this case, one can use the MAS simulation up to 1.1AU for training sOpInf and use the reduced model to predict up to 5AU, resulting in a more significant speed up. It is important to mention that the MAS model solves for several plasma flow quantities, i.e. Eqs. (1)- (6), and sOpInf is currently only solving for the radial velocity component. Thus, a direct comparison of their run-time can be equivocal.

A Comparison of Surrogate Model Accuracy via Equatorial Plane Streamlines
Given the two surrogate models of MAS, namely HUX (a reduced-physics approximation) and sOpInf (a data-driven ROM trained on MAS), we are interested in comparing their accuracy as surrogates of MAS. We do so by drawing the equatorial streamlines for each model. The streamlines of a flow field are curves that are tangential to the local velocity vector. In Figure 11(a-c), the streamlines are mapped from the inner-heliosphere at 0.14AU to 1.1AU. The shape of the streamlines depends on the velocity field, such that the fast solar wind results in less tightly wound lines than the slow solar wind. At regions where the streamlines interact, a compression wave is formed, whereas, at regions where the streamlines are distant there is a rarefaction wave. At these regions of compression and rarefaction, the solar wind streams go through substantial changes in density and flow speed [26]. For a better understanding of the mapped streamline accuracy, Figure 11(e-f) presents a histogram of the mapped streamlines longitude difference at 1.1AU along with plotting the cumulative distribution function (CDF) in Figure 11d of the streamlines longitude difference at 1.1AU for each surrogate model: HUX and sOpInf. The MAS in comparison to HUX and sOpInf trained on MAS mean/median/maximum and standard deviation (SD) of the streamline longitude error are presented in Table 2. The streamline mean longitude absolute error at 1.1AU of the HUX model is a factor of 12 larger than the ones associated with sOpInf. The numerical results show that the sOpInf model is   The choice of polynomial ROM form for the HUX and MAS dataset is an approximation of the governing equations (unlike the inviscid Burgers' example in Section 3.3) since both models have nonpolynomial terms. The above Sections 4.3-4.5 showcase the ROM model form that performed the best in the testing regime, i.e. purely-linear, purely-quadratic, or a combination thereof. Here, we present the results of a detailed investigation of how other polynomial ROM forms performed on each dataset. Figure 12 compares the state error in the fully predictive (testing) regime for each model: (12a) shows the HUX-2D equatorial plane results, (12b) shows the MAS-2D equatorial plane results, and (12c) shows the MAS-3D full-Sun results. In each case, we consider strictly-linear, strictlyquadratic, and linear-quadratic plus a constant term model forms in our analysis. Figure 12a shows a comparison of three different model forms for HUX dataset, in which strictly-quadratic and linearquadratic plus constant term ROMs perform better than the strictly-linear ROM. There is not a significant difference between the two quadratic forms, yet since the strictly-quadratic ROM has fewer model parameters and resulted in a slightly better relative error in the testing regime, we choose to employ a strictly-quadratic ROM form. For the MAS-2D dataset, the results in Figure 12b show that the strictly-quadratic model outperformed the other two model forms in the testing regime. Lastly, for the MAS-3D dataset, Figure 12c shows that the quadratic model with linear and constant terms outperforms the strictly-linear and strictly-quadratic ROMs with the same amount of modes.

Comparing Models with Similar Computational Cost
The main question we seek to answer: Is it better to have a linear ROM with larger or a quadratic ROM with smaller ? The cost of simulating the sOpInf ROM depends on the number of ROM parameters in the matrices on the right-hand side of Eq. (21). That number of parameters is determined by both the model form and the reduced basis dimension and determines the models computational cost. As mentioned in Section 3.1(III), we use the compact Kronecker product in the sOpInf ROM. Consequently, the number of model parameters d( ) (i.e., parameters in the system matrices) is for the different ROM forms we investigated. We thus compare the quadratic model forms chosen in Section 4.3-4.5 to a linear ROM with a comparable number of model parameters d( ) in Figure 13 The above numerical evidence highlights the importance of including quadratic nonlinearity in the sOpInf ROM. From another perspective, adding a quadratic term to the model allows us to have a lower ROM dimension than if only linear terms were present.

Conclusion
We proposed a reduced-order modeling strategy that uses simulated data to learn low-dimensional models for efficient solar wind predictions. The method leverages physical knowledge in that it first seeks to detect a spatial shift in the data/model (arising from advection) either through the method of characteristics or the fully data-based cross-correlation method. Given that shift, the system is then transformed into a moving coordinate frame, where a ROM can efficiently be learned via operator inference. The numerical results showed that for the full-Sun MAS simulations, a ROM with = 8 modes was sufficient to accurately predict the solar winds, yet produced significant computational speedup compared to the full-order model. From a surrogate modeling perspective, we investigated and compared the accuracy of two surrogates for the MAS model: a reduced-physics approximation (HUX) and the proposed sOpInf ROM approximation. We found that the latter is a much more accurate model than HUX; therefore, it is worth investigating ROM approaches for solar physics applications. Although we developed the sOpInf methodology to learn ROMs from simulated solar wind data, we found that the sOpInf framework is robust to moderate levels of noise, which is promising as one hopes to use sOpInf for noisy observational data. Additionally, when considering which sOpInf ROM model form should be chosen, we found that smaller quadratic ROMs perform better than larger, purely-linear ROMs in the testing regime (with comparable number of model parameters), highlighting the importance of quadratic terms in the ROM. While applied to use cases where solar wind velocities are most relevant, our methodology is applicable to forecasting additional solar wind quantities such as the density, pressure, etc. These models are our focus of future work. Moreover, while most quadratic forms of the ROM were sufficient to represent the physics with good accuracy, we expect more nonlinearly behaving systems to benefit from additional variable transformations (and lifting approaches similar to [28,49]). A long-term goal of our project is to solve uncertainty quantification problems, by efficiently assessing the impact of model input uncertainties, such as uncertain model coefficients, boundary conditions, and initial conditions, for which we anticipate our ROMs to be very useful. The ROMs can substantially accelerate ensemble methods, e.g. the most direct approach of Monte Carlo simulations and Bayesian inference, which are highly valuable in space weather operational forecasting.

Acknowledgement
This research was partially supported by the National Science Foundation under Award 2028125 for "SWQU: Composable Next Generation Software Framework for Space Weather Data Assimilation and Uncertainty Quantification".

Appendix A. Extending Shifted Operator Inference to Applications with Noisy Data
We developed the sOpInf methodology to predict the ambient solar wind from simulated data, yet the sOpInf framework extends to a wide class of advection-dominated systems described on a periodic domain and can be trained directly from observational data which is always noise-corrupted. We thus tested the sOpInf sensitivity to noise on the inviscid Burgers' equation example described in Section 3.3 with Gaussian noise added to each snapshot entry. Following the work by [21], the noise is drawn from a N (0, ν 2 ) with ν = ζ · (max x u(x, 0) − min x u(x, 0)) = ζ · (1.3 − 0.8) = ζ/2, where the coefficient ζ is the noise level. The same Gaussian initial condition, periodic boundary conditions, and first-order finite difference numerical solver described in Section 3.3 are used to generate the training and testing snapshots, where 80% of the total snapshots are used for training.
We generate noisy snapshots with various levels of noise ζ = 1, 2, . . . , 20%, and for each noise level compute the shift function c(t) via the cross-correlation extrapolation technique fitting a linear polynomial (see Section 3.2.2). Figure A.14a shows the cross-correlation extrapolation linear shift function c(t) slope as we vary the level of noise ζ = 1, 2, . . . , 20%. We found that the cross-correlation technique is robust to noise as the shift function remained within c(t) = (1.05 ± 0.012)t for ζ = 1, 2, . . . , 20%. Note that c(t) = 1.05t is the linear shift function for the noiseless data as mentioned in Section 3.3. Figure A.14b shows the singular values of the noisy and noiseless Burgers' training snapshots on the shifted and original coordinate frame for ζ = 2%. In agreement with Figure 2d, we see in Figure A.14b that the singular values decay faster in the shifted coordinates (in comparison to the original coordinates). Next, we compute the POD basis V = [v 1 , . . . , v ] ∈ R n× , shown in Figure A.14c. Inspecting the POD modes, we choose to keep = 5 POD basis functions, since after the fifth mode the POD modes are polluted with noise. This choice of ROM dimension = 5 is consistent with the singular values of the noisy data in the shifted coordinates, which plateau after = 5. The projected snapshotsÛ = V Ũ , where the rows ofÛ ∈ R ×K are denoted byÛ i,: ∈ R K for i = 1, 2, . . . , are the temporal ROM coefficients and are shown in Figure A.14d. Since the temporal coefficients are polluted with noise, using a uniform sixth-order finite difference scheme as we did in Section 3.3 will be highly inaccurate. Instead, the time derivativeU of the reduced states is computed via a simple factor method based on bimodal kernels of [17]. Figure A.15 shows the noisy snapshots for noise level ζ = 2%. Correspondingly, the sOpInf ROM results with = 5, λ 1 = 1 and λ 2 = 10 4 trained on noisy (noise level ζ = 2%) snapshots are shown in Figure A. 16. The numerical results illustrate that the sOpInf framework successfully reconstructs and predicts the dynamics of the inviscid Burgers' equation from noisy snapshots. The results show small oscillations near the shock in the testing regime. We suspect this behavior is due to the low number of modes ( = 5). The mean relative error is 1.1416 × 10 −3 and the Pearson correlation coefficient is 0.99987 in comparison to the noiseless snapshots. The numerical results show that the sOpInf methodology is robust to noise and can be potentially extended to applications with noisy observational data. The results show that the cross-correlation extrapolation technique is robust to noise as the shift function slope remains relatively close to 1.05 as we increase the noise level. Graphic (b) shows the singular value decay of the noisy (with ζ = 2%) and noiseless Burgers' equation snapshots on the original and shifted coordinates. Graphic (c) shows the first eight POD modes normalized between [−1, 1] of the noisy (noise level ζ = 2%) and noiseless inviscid Burgers' training snapshots. Of those, the first five POD modes of the noisy dataset are able to filter most of the noise, suggesting to learn a ROM with = 5. Graphic (d) shows the first five temporal ROM coefficients (normalized to [−1, 1]) in the training regime for noisy and noiseless snapshots.