Forecasting the Ambient Solar Wind with Numerical Models: I. On the Implementation of an Operational Framework

TL;DR

This paper presents an operational framework for forecasting the ambient solar wind using numerical models, combining magnetic field reconstruction, empirical relationships, and ensemble predictions to improve space weather forecasting accuracy.

Contribution

It introduces a comprehensive, validated framework integrating magnetic field models, empirical solar wind relations, and ensemble methods for operational solar wind forecasting.

Findings

WSA model most accurately predicts solar wind speed.

Ensemble forecasting slightly improves prediction accuracy.

Framework effectively studies coronal magnetic fields and solar wind properties.

Abstract

The ambient solar wind conditions in interplanetary space and in the near-Earth environment are determined by activity on the Sun. Steady solar wind streams modulate the propagation behaviour of interplanetary coronal mass ejections and are themselves an important driver of recurrent geomagnetic storm activity. The knowledge of the ambient solar wind flows and fields is thus an essential component of successful space weather forecasting. Here, we present an implementation of an operational framework for operating, validating and optimizing models of the ambient solar wind flow on the example of Carrington Rotation 2077. We reconstruct the global topology of the coronal magnetic field using the potential field source surface model (PFSS) and the Schatten current sheet model (SCS), and discuss three empirical relationships for specifying the solar wind conditions near the Sun, namely the Wang-Sheeley (WS) model, the distance from the coronal hole boundary (DCHB) model, and the Wang-Sheeley-Arge (WSA) model. By adding uncertainty in the latitude about the sub-Earth point, we select an ensemble of initial conditions and map the solutions to Earth by the Heliospheric Upwind eXtrapolation (HUX) model. We assess the forecasting performance from a continuous variable validation, and find that the WSA model most accurately predicts the solar wind speed time series. We note that the process of ensemble forecasting slightly improves the forecasting performance of all solar wind models investigated. We conclude that the implemented framework is well suited for studying the relationship between coronal magnetic fields and the properties of the ambient solar wind flow in the near-Earth environment.

Figures

Tables

Table 1. Table 1: The skill of model predictions of solar wind speed for CR2077 in terms of Arithmetic Mean (AM), Standard Deviation (SD), Mean Error (ME), Mean Absolute Error (MAE), and Root-Mean-Square Error (RMSE). A 4-day and 27-day persistence model of near-Earth solar wind conditions provides a baseline against which solar wind models can be compared.

Model	AM [km/s]	SD [km/s]	ME [km/s]	MAE [km/s]	RMSE [km/s]
WS	428.79	62.63	-27.90	74.09	85.27
DCHB	379.47	30.93	21.42	83.83	103.43
WSA	437.14	69.14	-36.25	68.54	82.62
Ensemble Median (WS)	435.55	51.84	-34.66	71.52	83.36
Ensemble Median (DCHB)	367.47	4.62	33.43	78.27	100.04
Ensemble Median (WSA)	437.48	63.64	-36.56	62.24	74.86
Persistence (4-days)	394.76	99.23	6.13	130.48	161.99
Persistence (27-days)	409.94	108.91	-9.05	66.54	78.86
Observation	400.89	96.12	-	-	-

Equations

$$M_0=\frac{1}{N_{\varphi }}\sum _{j=1}^{N_{\varphi }}{M}_{1,j},$$

$$M_{N_{\theta }+1}=\frac{1}{N_{\varphi }}\sum _{j=1}^{N_{\varphi }}{M}_{N_{\theta },j}.$$

$${\theta }_{i^{\prime }}^{\prime }=\pi \frac{i^{\prime }-1}{N_{\theta }^{\prime }-1},$$

$$M_{i^{\prime },j}^{\prime }=\alpha M_{i,j}+(1-\alpha )M_{i+1,j},$$

$$\theta \le {\theta }_{i^{\prime }}^{\prime }\le {\theta }_{i+1},$$

$$\alpha =\frac{{\theta }_{i+1}-{\theta }_{i^{\prime }}^{\prime }}{{\theta }_{i+1}-{\theta }_i}.$$

$$\min \left(\frac{2N_{\theta }^{\prime }}{3},\frac{N_{\varphi }}{3}\right).$$

$$\begin{array}{rl}\frac{dr}{ds}& =\frac{B_{\text{r}}}{B},\\ \frac{d\theta }{ds}& =\frac{1}{r}\frac{B_{\theta }}{B},\\ \frac{d\varphi }{ds}& =\frac{1}{r\sin \theta }\frac{B_{\varphi }}{B},\end{array}$$

$$f_{\text{p}}={\left(\frac{R_0}{R_1}\right)}^2\left|\frac{B_{\text{r}}(R_0,{\theta }_0,{\varphi }_0)}{B_{\text{r}}(R_1,{\theta }_1,{\varphi }_1)}\right|,$$

$$v_{\text{ws}}(f_{\text{p}})=v_0+\frac{(v_1-v_0)}{{f_{\text{p}}}^{\alpha }}.$$

$$v_{\text{dchb}}(d)=v_0+\frac{1}{2}\left(v_1-v_0\right)\left[1+\tanh \left(\frac{d-ϵ}{w}\right)\right],$$

$$v_{\text{wsa}}(f_{\text{p}},d)=v_0+\frac{\left(v_1-v_0\right)}{{\left(1+f_{\text{p}}\right)}^{\alpha }}{\left\{\beta -\gamma \exp \left[{-(d/w)}^{\delta }\right]\right\}}^3,$$

$$v_{r+1,\varphi }=v_{r,\varphi }+\frac{\Delta r\Omega }{v_{r,\varphi }}\left(\frac{v_{r,\varphi +1}-v_{r,\varphi }}{\Delta \varphi }\right),$$

$$v_{\text{acc}}(r)=\alpha v_{r_0}\left[1-\exp \left(\frac{-r}{r_h}\right)\right],$$

$${\theta }_p(\varphi )={\theta }_E+{\theta }_1\sin \left(n\varphi +{\varphi }_0\right),$$

$$\mathbf{B}=-\nabla \Psi .$$

$${\nabla }^2\Psi =0.$$

$$\begin{array}{rl}\Psi (r,\theta ,\varphi )& =\left[R_0{\left(\frac{R_0}{r}\right)}^{n+1}-R_1{\left(\frac{R_0}{R_1}\right)}^{n+2}{\left(\frac{r}{R_1}\right)}^n\right]\\ & \times \sum _{n=1}^{\infty }\sum _{m=0}^n\left(g_n^m\cos m\varphi +h_n^m\sin m\varphi \right)P_n^m(\cos \theta ),\end{array}$$

$$\mathbf{B}=(B_{\text{r}},B_{\theta },B_{\varphi })=\left(-\frac{\partial \psi }{\partial r},-\frac{1}{r}\frac{\partial \psi }{\partial \varphi },-\frac{1}{r\sin \theta }\frac{\partial \psi }{\partial \varphi }\right).$$

$$\begin{array}{rl}{B}_{\text{r}}& =\left[(n+1){\left(\frac{R_0}{r}\right)}^{n+2}+{\left(\frac{R_0}{R_1}\right)}^{n+2}{\left(\frac{r}{R_1}\right)}^{n-1}\right]\\ & \times \sum _{n=1}^{\infty }\sum _{m=0}^n\left(g_n^m\cos m\varphi +h_n^m\sin m\varphi \right)P_n^m(\cos \theta ),\end{array}$$

$$\begin{array}{rl}{B}_{\theta }& =-\left[{\left(\frac{R_0}{r}\right)}^{n+2}-{\left(\frac{R_0}{R_1}\right)}^{n+2}{\left(\frac{r}{R_1}\right)}^{n-1}\right]\\ & \times \sum _{n=1}^{\infty }\sum _{m=0}^n\left(g_n^m\cos m\varphi +h_n^m\sin m\varphi \right)P_n^m(\cos \theta ),\end{array}$$

$$\begin{array}{rl}{B}_{\varphi }& =\left[{\left(\frac{R_0}{r}\right)}^{n+2}-{\left(\frac{R_0}{R_1}\right)}^{n+2}{\left(\frac{r}{R_1}\right)}^{n-1}\right]\\ & \times \sum _{n=1}^{\infty }\sum _{m=0}^n\frac{m}{\sin \theta }\left(g_n^m\sin m\varphi -h_n^m\cos m\varphi \right)P_n^m(\cos \theta ).\end{array}$$

$$\frac{1}{4\pi }{\int }_0^{\pi }{\int }_0^{2\pi }{P}_n^m(\theta )\left\{\begin{array}{c}\cos m\varphi \\ \sin m\varphi \end{array}\right\}{P}_{n^{\prime }}^{m^{\prime }}(\theta )\left\{\begin{array}{c}\cos m^{\prime }\varphi \\ \sin m^{\prime }\varphi \end{array}\right\}\sin \theta d\theta d\varphi =\frac{1}{2n+1}{\delta }_n^{n^{\prime }}{\delta }_m^{m^{\prime }}.$$

$$\left\{\begin{array}{c}{g}_n^m\\ h_n^m\end{array}\right\}=\frac{2n+1}{4\pi \left(n+1+n{\left(\frac{R_0}{R_1}\right)}^{2n+1}\right)}{\int }_0^{\pi }d\theta \sin \theta P_n^m(\theta ){\int }_0^{2\pi }d\varphi B_{\text{r}}(R_0,\theta ,\varphi )\left\{\begin{array}{c}\cos m\varphi \\ \sin m\varphi \end{array}\right\}.$$

$${\int }_0^{\pi }d\theta \sin \theta F(\theta )\approx \sum _{i=1}^{N_{\theta }}{ϵ}_i{w}_{i}F({\theta }_i),$$

$$w_i=-\frac{2}{H}\sum _{k=0}^H\frac{{ϵ}_k^{\prime }}{4k^2-1}\cos \left(\frac{\pi k(i-1)}{H}\right)$$

$$\left\{\begin{array}{c}{g}_n^m\\ h_n^m\end{array}\right\}=\frac{2n+1}{4\pi \left(n+1+n{\left(\frac{R_0}{R_1}\right)}^{2n+1}\right)}\frac{2\pi }{N_{\varphi }}\sum _{i=1}^{N_{\theta }}\sum _{j=1}^{N_{\varphi }}{ϵ}_i{w}_i{P}_n^m({\theta }_i)B_{\text{r}}(R_0,{\theta }_i,{\varphi }_j)\left\{\begin{array}{c}\cos m{\varphi }_j\\ \sin m{\varphi }_j\end{array}\right\}.$$

$$\Psi =R_0\sum _{n=0}^{\infty }\sum _{m=0}^n\left[{\left(\frac{R_0}{r}\right)}^{n+1}\left(g_n^m\cos m\varphi +h_n^m\sin m\varphi \right)P_n^m(\varphi )\right].$$

$$\begin{array}{rl}{B}_{\text{r}}& =-\frac{d\Psi }{dr}=\sum _{n=1}^{\infty }(n+1){\left(\frac{R_1}{r}\right)}^{n+2}\sum _{m=0}^n\left(g_n^m\cos m\varphi +h_n^m\sin m\varphi \right)P_n^m(\cos \theta ),\\ B_{\theta }& =-\frac{1}{r}\frac{d\Psi }{d\theta }=-\sum _{n=1}^{\infty }{\left(\frac{R_1}{r}\right)}^{n+2}\sum _{m=0}^n\left(g_n^m\cos m\varphi +h_n^m\sin m\varphi \right)\frac{d{P}_n^m(\cos \theta )}{d\theta },\\ B_{\varphi }& =-\frac{1}{r\sin \theta }\frac{d\Psi }{d\varphi }=\sum _{n=1}^{\infty }{\left(\frac{R_1}{r}\right)}^{n+2}\sum _{m=0}^n\frac{m}{\sin \theta }\left(g_n^m\sin m\varphi -h_n^m\cos m\varphi \right)P_n^m(\cos \theta ).\end{array}$$

$$F=\sum _{i=1}^{N_i}\sum _{j=1}^{N_j}\sum _{k=1}^3{\left[B_k({\theta }_i,{\varphi }_j)-\sum _{n=0}^{N_s}\sum _{m=0}^n\left(g_n^m{\alpha }_{nmk}({\theta }_i,{\varphi }_j)+h_n^m{\beta }_{nmk}({\theta }_i,{\varphi }_j)\right)\right]}^2,$$

Full text

K****eywords Solar wind  $\cdot$ Solar-terrestrial relations  $\cdot$ Sun: heliosphere  $\cdot$ Sun: magnetic fields

1 Introduction

The knowledge of the evolving ambient solar wind is an essential component of successful space weather forecasting (Owens and Forsyth, 2013). The ambient solar wind flows and fields play an important role in understanding the propagation of coronal mass ejections (CMEs) and are themselves an important driver of recurrent geomagnetic activity. More specifically, high-speed solar wind streams contribute about $70\,\%$ of geomagnetic activity outside of solar maximum and about $30\,\%$ at solar maximum Richardson et al. (2000). Even when the occurrence rate of CMEs increases from $0.3$ per day during solar minimum to about $4$–$5$ per day during solar maximum St. Cyr et al. (1999), key properties in interplanetary space, such as the solar wind bulk speed, magnetic field strength, and orientation, are determined by the ambient solar wind flow Luhmann et al. (2002).

The topology of open magnetic field lines along which ambient solar wind flows accelerate to supersonic speeds plays a fundamental role in understanding phenomena that drive our evolving space weather. To date, there are no routine measurements of the coronal magnetic field. Magnetic models for the solar corona, therefore, have to rely on extrapolations calculated from the observed line-of-sight component of the photospheric field. The most widely applied extrapolation technique to reconstruct global solutions for the coronal magnetic field is the Potential Field Source Surface model (PFSS; Altschuler and Newkirk, 1969; Schatten et al., 1969). Using the current-free (or potential field) assumption for regions above the photosphere, solutions for the magnetic field can be expressed as the gradient of a magnetic scalar potential. Since potential fields give closed magnetic fields, a spherical source surface, where the magnetic field is assumed to be only radial is added as an outer boundary condition. The radius of the spherical source surface is an adjustable parameter typically set to a reference height of $2.5\,R_{0}$ to best match observations Hoeksema et al. (1982). To account for Ulysses measurements showing latitudinal invariance of the radial magnetic field component Wang and Sheeley (1995), the so-called Schatten current sheet (SCS) model Schatten (1971) is typically added beyond the PFSS model to create a more uniform radial field strength solution.

The effect of the solar wind is to drag and distort magnetic field lines, and thus to distort the coronal field from the assumed current-free configuration. Ideally, a model should account for the complex dynamics of the solar wind flow by solving a set of nonlinear partial differential equations of magnetohydrodynamics (MHD) (e.g., Altschuler and Newkirk, 1969). As such, the Magnetohydrodynamics Algorithm outside a Sphere model (MAS; Linker et al., 1999; Mikić et al., 1999) and the Space Weather Modeling Framework (SWMF; Tóth et al., 2005) are three-dimensional MHD models. The PFSS solutions are usually used as an initial condition, and the MHD equations are integrated in time until the plasma and magnetic fields relax into equilibrium. The final steady-state solution is characterized by closed magnetic fields that confine the solar wind plasma and open magnetic field lines along which solar wind flows accelerate to supersonic speeds.

The state-of-the-art framework for forecasting the ambient solar wind couples models of the corona with those of the inner heliosphere Riley et al. (2001); Lee et al. (2008). The coronal part spans the range from 1 solar radii ($R_{0}$) to $2.5\,R_{0}$ (PFSS), or $20\,R_{0}$ to $30\,R_{0}$ (MHD). The inner boundary condition of the heliospheric part either results directly from the coronal model or is computed from the topology of the coronal magnetic field. The outer boundary condition computed from the coronal model is used as inner boundary condition for the heliospheric model ($20$–$30\,R_{0}$ to 1 au). Typically, the heliospheric model is based on the MAS model Linker et al. (1999); Mikić et al. (1999), Enlil Odstrcil (2003), or the recently developed European Heliospheric Forecasting Information Asset (EUHFORIA; Pomoell and Poedts, 2018), each of which are three-dimensional numerical MHD models that derive stationary solutions for the ambient solar wind in interplanetary space. Figure 1 shows some different coronal and heliospheric model combinations, together with other approaches for forecasting the ambient solar wind based on empirical relationships (e.g., Robbins et al., 2006; Vršnak et al., 2007; Shugay et al., 2011; Reiss et al., 2016), and statistics or machine learning techniques (e.g., Owens et al., 2013; Riley et al., 2017). The orange colored lines highlight models accessible online at NASA’s Community Coordinated Modeling Center (CCMC) online platform.

Many recent studies (e.g., Owens et al., 2008; Jian et al., 2011; Devos et al., 2014; Reiss et al., 2016) have assessed the performance of operational frameworks for forecasting the ambient solar wind and reported typical uncertainties of about 1 day in the arrival time of high-speed streams. Furthermore, it is now well established that the performance of models of the ambient solar wind is, if at all, only slightly better than a 27 day persistence model (e.g., Owens et al., 2013; Kohutova et al., 2016), assuming that the near-Earth solar wind conditions will replicate after each Carrington Rotation (CR). Overall, these results highlight the fact that forecasting the conditions in the ambient solar wind in interplanetary space and in the near-Earth environment is challenging, even during times of low solar activity when the large-scale interplanetary magnetic field configuration evolves less rapidly and disturbances due to CMEs are less frequent Owens and Forsyth (2013). As outlined in the space weather roadmap for the years 2015–2025 (see, Schrijver et al., 2015), continued efforts are needed to improve our capabilities for forecasting the conditions in the ambient solar wind.

In this paper, we present an implementation of a numerical framework for operating, validating, and optimizing models of the evolving ambient solar wind. To study a large number of initial conditions, we rely on the coupling between well-established coronal models Wang and Sheeley (1990); Riley et al. (2001); Arge et al. (2003) and the Heliospheric Upwind eXtrapolation (HUX) model Riley and Lionello (2011) that bridges the gap of ballistic mapping and MHD modeling. While this paper outlines the breakdown structure and the mathematical foundation of the operational framework, a subsequent paper will be devoted to the validation and optimization of the numerical models by examining the relationship between the coronal magnetic fields and properties of the ambient solar wind. This paper is organized into three sections. In Section 2, we discuss the numerical framework for forecasting the ambient solar wind, including the remeshing of magnetograms (Section 2.1), the magnetic model of the solar corona (Section 2.2), the features of the coronal field solution (Section 2.3), the empirical relationships for specifying the solar wind conditions near the Sun (Section 2.4), the mapping of the solar wind solutions to Earth (Section 2.5), the application of a sensitivity analysis (Section 2.6), and the quantification of the forecasting performance (Section 2.7); and in Section 3 we conclude with a summary of the results and outline future applications of the framework.

2 Description of the Numerical Framework

We present an implementation of a numerical framework for operating, validating, and optimizing models of the ambient solar wind. Figure 2 illustrates the breakdown structure of the current framework implementation. The coronal domain spans the range from $1\,R_{0}$ to $5\,R_{0}$, where the outer boundary condition computed from the coronal part is used as an inner boundary condition for the heliospheric domain. The model in the heliospheric domain where the solar wind flow is supersonic then uses the boundary condition as an input for propagating the solar wind solutions near the Sun to 1 au. The subsequent sections are concerned with explaining the components of the framework in detail.

2.1 Input to the Models

We use full-disk photospheric magnetograms from the Global Oscillation Network Group (GONG) from the National Solar Observatory (NSO) as an inner boundary condition for the coronal model. The global maps of the solar magnetic field measured in Gauss (G) are given on the $\left(\sin\theta,\phi\right)$ grid with $180\times 360$ grid points, where $\theta\in[0,\pi]$ and $\phi\in[0,2\pi]$ are the latitude and longitude coordinates, respectively. They are available as near-real-time magnetic maps or full CR maps at the GONG online platform111http://gong.nso.edu. Throughout this paper, we illustrate our model solutions on the example of CR2077 (2008 November 20–December 17), that is, during solar minimum when the global magnetic field is dominated by its dipolar component and pronounced polar coronal holes are observable.

We use the full-disk magnetic maps as an inner boundary condition for the coronal model. The coronal model is based on a spherical harmonic decomposition of the input magnetogram. When the raw magnetogram on the $\left(\sin\theta,\phi\right)$ grid is used as an input and the order of spherical harmonics expansion is large compared to the resolution of the magnetic map, inaccurate results, especially at the crucial polar regions, are expected. Tóth et al. (2011) concludes that the information from magnetograms is used more efficiently when the input magnetograms (i) are remeshed to a grid that is uniform in the latitude $\theta$, (ii) contain both poles at $\theta=0$, and $\theta=\pi$, and (iii) have an odd number of grid points.

We remesh the magnetograms according to the linear interpolation procedure discussed in Tóth et al. (2011). To do so, we begin with adding two additional grid cells at the north and south poles of the input magnetic map $M_{i,j}$ with $N_{\theta}=180$ and $N_{\phi}=360$ grid points in the latitude and longitude, respectively. The values for the extra grid cells at the poles $M_{0}$ and $M_{N_{\theta}+1}$ are given by

[TABLE]

and

[TABLE]

The latitude of the uniform $\theta$ grid is

[TABLE]

where $N^{\prime}_{\theta}$ is the number of grid points at the uniform $\theta$ grid, and the index $i^{\prime}=1\dots N^{\prime}_{\theta}$. Finally, we interpolate the grid points from the raw magnetogram mesh to the uniform $\theta$ mesh using a linear interpolation relation (no magnetic flux conservation) of the form

[TABLE]

where the index $i$ is selected so that

[TABLE]

and

[TABLE]

Notice that the maximum degree for the expansion of spherical harmonics is now limited only by the anti-alias limit (see, Tóth et al., 2011) given by

[TABLE]

By utilizing remeshed magnetograms as an inner boundary condition for the coronal model, we can, therefore, expect accurate solutions up to the order of 120 spherical harmonics Tóth et al. (2011). Figure 3(a) depicts the computed remeshed magnetogram as a function of the heliographic latitude and Carrington longitude for CR2077.

2.2 Magnetic Model of the Solar Corona

The magnetic model of the solar corona couples the PFSS model and the SCS model to reconstruct the coronal magnetic field. The PFSS model is based on the assumptions that the region above the photosphere is current free, and that the magnetic field at an imaginary reference sphere, called the source surface with radius $R_{1}$ is radial only Altschuler and Newkirk (1969); Schatten et al. (1969). As detailed in Appendix A, Eq.(19) can be solved to derive the magnetic field components at any point in the coronal domain ($R_{0}\leq r\leq R_{1}$). Figure 3 (b)–(d) present the magnetic field components $B_{\text{r}},B_{\theta},andB_{\phi}$ derived from the PFSS model at the solar surface as a function of the latitude and Carrington longitude for CR2077. We note that we used a color scheme for all illustrations that is color-blind friendly (see, Light and Bartlein, 2004). From the comparison in Figure 3, it is apparent that the PFSS solution using the spherical harmonics expansion up to $n=80$, is in reasonable agreement with the observed line-of-sight component of the photospheric magnetic field.

To extend the magnetic field solution from the outer boundary of the PFSS model at $2.5\,R_{0}$ to a distance of $R_{2}=5\,R_{0}$, we employ the SCS model Schatten (1971). The SCS model uses the PFSS solutions at the source surface as an inner boundary condition. The SCS model is similar to the PFSS model but involves solving an additional Laplace equation with different boundary conditions. In the first step, the magnetic field on the source surface is oriented to point away from the Sun. In other words, the signs of the magnetic field components $B_{\text{r}}$, $B_{\theta}$, and $B_{\phi}$ are reversed if $B_{\text{r}}<0$ at the source surface. While the magnetic field solution is defined to vanish at infinity, the non-negative magnetic field values as model input ensure that a thin current sheet is retained in the field solution. The underlying idea of the SCS model is to extend the magnetic field solution in a nearly radial way. To retain the resolution, we match the grid of the SCS model with the PFSS model with equal steps of $2\degree$ in latitude and longitude. As discussed in Appendix B, we use estimates of the coefficients $g_{n}^{m}$ and $h_{n}^{m}$ from the least mean square fit to compute the magnetic field above the source surface. Finally, the polarity needs to be corrected to match the polarities in the region $r\leq R_{1}$.

The imposed boundary conditions of the PFSS model and the SCS model are not compatible, causing discontinuities in the tangential component of the magnetic field across the source surface. When coupling the PFSS model and SCS model solutions, discontinuities in the form of kinks in the magnetic field topology are observable. In order to account for this effect, McGregor et al. (2008) proposed a more flexible coupling between the two models. The authors concluded that their approach leads to more accurate forecast results. Following the procedure of McGregor et al. (2008), we set the radius of the source surface to $2.5\,R_{0}$ but use the PFSS solution at $2.3\,R_{0}$ as an inner boundary condition for the SCS model. By doing so, we couple the PFSS and SCS model solutions and reconstruct the global topology of the coronal magnetic field.

2.3 Computing Features of the Coronal Field Solution

Models for specifying the solar wind conditions near the Sun rely on the areal expansion factor $f_{\text{p}}$ and the great-circle angular distance from the nearest open-closed boundary $d$, both of which are derived from the coronal magnetic field solution. It is noteworthy that both features are linked to different fundamental theories on the origin of slow solar wind (see, Riley and Luhmann, 2012; Riley et al., 2015). While the expansion factor assumes that the solar wind flow along open field lines that diverge the most leads to the slow solar wind Wang and Sheeley (1990), the distance to the coronal hole boundary is more related to the boundary layer idea of interchange reconnection for the origin of the slow solar wind Antiochos et al. (2011).

To trace magnetic field lines and reconstruct the magnetic field topology in the coronal model domain $R_{0}\leq r\leq R_{2}$, we use a fourth-order Runge-Kutta method (RK4) and solve the following set of equations, expressed as

[TABLE]

where $ds$ is a segment along the magnetic field line.

More specifically, our approach for constructing the large-scale topology of the coronal field is two-fold. First, we start from the base of the model at the solar surface ($r=R_{0}$) and trace magnetic field lines upwards. When the field line returns to the solar surface, we label the corresponding footpoint at the solar surface as a “closed field”. In contrast, when the field line reaches the upper boundary of the coronal model ($r=R_{2}$) we label the footpoint as a “open field”. In this way, we compute coronal hole regions, i.e., magnetically open regions expected to guide high-speed solar wind streams out into the heliosphere. Using the location of coronal holes at the solar surface, we run a perimeter tracing algorithm to compute coronal hole boundaries. Then, we compute the great-circle angular distance $d$ between footpoints of open field lines and their nearest coronal hole boundary. Second, we start from the outer boundary of the coronal model ($r=R_{2}$) and trace the field lines down to the solar surface. By doing so, we compute the areal expansion factor $f_{\text{p}}$ as

[TABLE]

where $B_{\text{r}}$ is the radial magnetic field component at a given reference height. The areal expansion factor is the amount by which a flux tube expands from the solar surface to another reference height in the corona. Figure 4 shows the topology of the coronal magnetic field and illustrates the features $d$ and $f_{\text{p}}$ at the photosphere and the source surface, respectively.

2.4 Computing Solar Wind Speed Maps

The structure in the solar wind flow is governed by the dynamic pressure term in the momentum equation ($\propto\rho\,v^{2}$). This indicates that the near-Sun solar wind solution used as an inner boundary condition for heliospheric models is highly sensitive to errors in the computed solar wind bulk speed Riley et al. (2015). Recently, Riley et al. (2015) highlighted three empirical relationships for specifying the solar wind speed $(v=v(d,f_{\text{p}}))$ at a reference sphere of $5\,R_{0}$ (or $30\,R_{0}$ for the MAS model). In the literature, these empirical relationships are known as the Wang-Sheeley (WS) model Wang and Sheeley (1990), the Distance from the Coronal Hole Boundary (DCHB) model Riley et al. (2001), and the Wang-Sheeley-Arge (WSA) model Arge et al. (2003).

The WS model is based on the inverse relationship between the solar wind speed and the magnetic field expansion factor Wang and Sheeley (1990), namely

[TABLE]

Previous research has established that low magnetic field expansion between the photosphere and some reference height in the corona is correlated with a fast solar wind speed (e.g., Levine et al., 1977). For the coefficients in Eq.(10), we use $v_{0}=250$, $v_{1}=660$, and $\alpha=2/5$ (see, Arge and Pizzo, 2000).

The DCHB model relates the speed at the photosphere with the distance of an open field footpoint from the coronal hole boundary and maps the calculated solar wind speed solution out along the magnetic field lines to a given reference sphere Riley et al. (2001). The DCHB relation is defined as

[TABLE]

where $\epsilon$ is a measure for the thickness of the slow flow band, and $w$ denotes the width over which the solar wind reaches coronal hole values. For an open field footpoint at the coronal hole boundary, the solar wind speed is equal to $v_{0}$. For a footpoint located deep inside a coronal hole, the solar wind speed is equal to $v_{1}$. Hence, the farther away the footpoint is from a coronal hole boundary, the faster the expected solar wind speed. In this study, we use $v_{0}=350$, $v_{1}=750$, $\epsilon=0.1$, and $w=0.05$ (see, Riley et al., 2001).

Finally, the WSA model is a hybrid of the WS model and the DCHB model, combining aspects of the expansion factor computed from the topology of the magnetic field and the distance to the coronal hole boundary Arge et al. (2003). The WSA model for specifying solar wind speed is given by

[TABLE]

where $\alpha$, $\beta$, $\gamma$, and $\delta$ are additional model parameters. For the coefficients in Eq.(12) we use $v_{0}=285$, $v_{1}=910$, $\alpha=2/9$, $\beta=1$, $\gamma=0.8$, $w=2$, and $\delta=3$. Figure 5 shows a comparison between the different empirical relationships for specifying a solar wind speed of $v(f_{\text{p}},d)$ near the Sun.

2.5 Mapping of the Solar Wind to Earth

Numerous models dealing with the mapping of solar wind solutions near the Sun to Earth have been developed. The spectrum extends from a simple ballistic approximation where each parcel of plasma is assumed to travel with a constant speed through the heliosphere to global heliospheric MHD models that aim to cover all relevant dynamical processes (e.g., Riley et al., 2011; Odstrcil, 2003). In an attempt to bridge the gap, Riley and Lionello (2011) developed the Heliospheric Upwind eXtrapolation (HUX) model. Riley and Lionello (2011) have simulated the kinematic of the solar wind flow numerically by simplifying the MHD equations as much as possible. By neglecting the pressure gradient and the gravity terms in the fluid momentum equation (e.g., Pizzo, 1978) the solar wind solution on a discretized grid may be written as

[TABLE]

where the subscripts $r$ and $\phi$ refer to radial and longitudinal grid cells with cell spacings $\Delta r$ and $\Delta\phi$, and $\Omega$ denotes the equatorial angular rotation rate of the Sun (neglecting differential rotation) (see, Riley and Lionello, 2011). In order to match the spatial resolution of the coronal model, we use $\Delta\phi=2\degree$ and $\Delta r=1\,R_{0}$, respectively. As the step size in $\Delta r$ is increased by two orders of magnitude, the solutions are not significantly different, indicating that the Courant-Friedrichs-Lewy condition for numerical convergence is well met.

Furthermore, Riley and Lionello (2011) concluded that it is convenient to account for the residual wind acceleration beyond the coronal model. According to previous model results (e.g., Riley et al., 2001; Riley and Lionello, 2011), the residual acceleration $v_{\text{acc}}$ expected beyond the coronal model is

[TABLE]

where $v_{r_{0}}$ is the initial speed at the outer boundary of the coronal model, $r$ is the heliocentric distance, $\alpha$ is a factor by which $v_{r_{0}}$ is enhanced, and $r_{h}$ is the scale length for the acceleration (not crucial when mapping the solar wind to 1 au). For the coefficients in Eq.(14) we use $\alpha=0.15$ and $v_{r_{0}}=50\,R_{0}$.

Figure 6 illustrates the HUX model for mapping the solar wind solution near the Sun to Earth. While the top panel shows the propagation of the solar wind solutions in steps of $\Delta r$ from $5\,R_{0}$ to $215\,R_{0}$ (or 1 au), the bottom panel shows the propagation in spherical coordinates (left) and the calculated speed time series at Earth with and without the residual speed contribution (right). It is also noteworthy that computing the solar wind solution for other heliospheric distances is straightforward. The benefit of the HUX model is that it can match the dynamical evolution explored by global heliospheric MHD models that demand only low computational requirements Riley and Lionello (2011). This is particularly useful to study a large number of initial conditions in the context of ensemble modeling.

2.6 Creating Ensembles of Initial Conditions for Solar Wind Forecasting

The knowledge on the sensitivity of the model performance to initial conditions and model parameter settings is crucial for models of the ambient solar wind. In recent years, the concept of ensemble forecasting has been successfully applied to a number of applications in the space weather forecasting regime. As an example, the sensitivity of PFSS and MHD coronal models on magnetic maps from different observatories has been studied by Riley et al. (2013a, b). More details on recent activities can be found in the reviews by Knipp (2016) and Murray (2018). Ensemble forecasting is a technique based on the use of a sample of possible future states to forecast a future state. Intuitively, one would expect that the forecast of the ambient solar wind conditions is very uncertain when the ensemble members (e.g., different solar wind models) are very different from one another. The strength of ensemble forecasting is its capability to deduce confidence bounds by quantifying the uncertainty in the ensemble of possible future states.

Owens et al. (2017), for instance, has pointed out that the forecast uncertainty in solar wind models can be studied by adding uncertainty in the sub-Earth orbit at which the coronal solutions near the Sun are sampled. Using the approach of Owens et al. (2017), we perform a sensitivity analysis that considers ensemble members at perturbed latitudes $\theta_{p}$ around the sub-Earth orbit given by

[TABLE]

where $\phi$ is the corresponding Carrington longitude; $\theta_{E}$ is the sub-Earth latitude; $\theta_{1}$ is the amplitude of the deviation; $n$ is the wave number, indicating how many oscillations the perturbed trajectory completes; and $\phi_{0}$ is the phase offset, indicating how two perturbed trajectories can be out of synchronicity with each other. For the coefficients in Eq.(14), we select $\theta_{1}\in[0\degree,15\degree]$ in $1\degree$ steps, $n\in[0,1]$ in $0.5$ steps, and $\theta_{0}\in[0\degree,330\degree]$ in $30\degree$ steps, which ensures that all latitudes are well represented at each longitude (see, Owens et al., 2017).

Figure 7(a) shows three individual trajectories together with the median value and the $\pm 2\sigma$ quantiles of the ensemble members. In this way, we obtain information from a sampling of 576 initial speed solutions, each of which are mapped to Earth using the HUX model as outlined in Section 2.5. Figure 7(b) compares the resulting ensemble median at $215\,R_{0}$ and the observed solar wind speed in the near-Earth environment. In this study, the ensemble median is the preferred average measure as the ensembles of solar wind solutions around the sub-Earth latitude are often highly skewed, such that the ensemble mean can yield very biased measures of the ensemble average. Further, we note that ensemble averaging of possible future states is not necessarily expected to provide a systematic improvement of deterministic forecasts (see, Jolliffe and Stephenson, 2003). As discussed in Owens et al. (2017) and Henley and Pope (2017), it is the assessment of forecast uncertainty that is the key advantage of the ensemble approach rather than the ensemble average providing an improved deterministic forecast.

2.7 Assessing the Forecast Performance

The information from the implemented models is used most efficiently when the uncertainty of their results is constantly validated. To enable that, we measure the performance of the framework by the Operational Solar Wind Evaluation Algorithm (OSEA; Reiss et al., 2016). OSEA222https://bitbucket.org/reissmar/solar-wind-forecast-verification runs various validation procedures to compare the forecasts and the observations to which they pertain. Traditionally, the relationship between forecast and observation can be studied in terms of continuous variables and binary variables. While the former can take on any real values, the latter is restricted to two possible values such as event/non-event. In the context of solar wind forecasting, the solar wind speed time series can be interpreted in terms of both aspects. The forecasting performance can either be evaluated in terms of an average error or in terms of the capability of forecasting events of an enhanced solar wind speed Owens et al. (2005); MacNeice (2009a, b). OSEA is capable of quantifying both aspects, i.e., a continuous variable evaluation that uses simple point-to-point comparison metrics and an event-based validation analysis that assesses uncertainty of the arrival time of high-speed solar wind streams at Earth.

As an example, Table 1 lists the results obtained from the continuous variable validation of different solar wind models for CR2077 in terms of Arithmetic Mean (AM), Standard Deviation (SD), Mean Error (ME), Mean Absolute Error (MAE), and Root-Mean-Square Error (RMSE). A 4 day and 27 day persistence model of near-Earth solar wind conditions provides a baseline against which the performance can be compared. We find that the RMSE for the WS model, DCHB model, and WSA model is 85.27 km/s, 103.43 km/s, and 82.62 km/s, respectively. The 27 day persistence model has the same statistics as the measurements and greatly benefits from the quasi-steady and recurrent nature of the evolving ambient solar wind in the solar minimum phase and thus is expected to score very high in all measures (e.g., $\text{RMSE}=78.86$ km/s). We conclude that the WSA model gives the best forecast results in terms of a continuous variable validation and that the process of ensemble forecasting slightly improves the performance of all solar wind models in this study. It is important to note that the present analysis is indented to illustrate the application of the implemented framework components, which means that our results apply only to CR2077 and that our conclusions are not reliable as a general guide.

3 Discussion

We present a numerical framework for forecasting the evolving ambient solar wind that uses magnetic maps as an input for the PFSS and SCS model to reconstruct the global topology of the coronal magnetic field, specifies the solar wind speed using different established empirical relationships (WS, DCHB, and WSA models) based on the areal expansion factor and the distance to the nearest coronal hole boundary, maps the near-Sun solar wind solution outward to Earth by the HUX model, creates an ensemble of initial conditions by adding uncertainty in the latitude about the sub-Earth point, and uses an automated forecast validation module to quantitatively assess the forecasting skill. The framework relies on established models of the ambient solar wind and is conceptually very similar to already existing numerical frameworks. We carefully compared our model solutions to existing frameworks (e.g., Nikolic et al., 2014; Arge et al., 2003), and found that our modular framework implementation in the C++ and Matlab programming language using tools from the Armadillo library Sanderson and Curtin (2016) is both robust and fast.

The coronal part of the framework relies on empirical relationships between the magnetic field topology and the near-Sun solar wind conditions Wang and Sheeley (1990); Riley et al. (2001); Arge et al. (2003). Ideally, one would prefer the application of a physics-based MHD model for the corona to capture the complex dynamics at solar wind stream interaction regions, which are not included in the described model approach. Nevertheless, recent studies have shown that the forecast skill of empirical and full physics-based coupled corona-heliosphere models in terms of established metrics is very similar. As an example, Owens et al. (2008) studied the performance of different forecast models (WSA, WSA-Enlil, and MAS-Enlil) over an 8 yr period and concluded that the coupled empirical approach currently gives the best forecast results in terms of the mean square error. Considering the trade-off between accuracy and computational requirements of a full MHD code, it thus seems reasonable to follow the described methodology in the context of solar wind forecasting. We conclude that the efficient implementation of the framework using the heliospheric HUX model is well suited for studying the long-term relationship between coronal magnetic fields and the properties of the ambient solar wind.

In the future, we shall work on several topics to try to improve the forecasting performances. While this study presents the implementation of the numerical framework, a subsequent paper will be devoted to the validation and optimization of the present solar wind models. By studying the coupling between magnetic models of the corona and those of the inner heliosphere, we aim toward optimizing the empirical relationships for specifying solar wind properties to advance further their predictive capabilities. In context, a reliable forecast of the ambient solar wind might be beneficial for forecasting the arrival of CMEs. Here, we plan to combine the ambient solar wind framework with the ELlipse Evolution model based on Heliospheric Imager observations (ELEvoHI; Rollett et al., 2016; Amerstorfer et al., 2018). The kinematics of the elliptical-shaped CME front in the ELEvoHI model is governed by the ambient solar wind flow. In the current version of ELEvoHI, the solar wind speed is assumed to be constant during the propagation of the CME in the inner heliosphere. Therefore, we would also like to include information on the complex dynamics of the evolving ambient solar wind to simulate the dynamic deformation of the CME front during the propagation phase. To make the model runs accessible to the space weather community, we also plan to install a later version of the solar wind framework at NASA’s CCMC online platform. The release of this online resource is scheduled for 2019 July.

Acknowledgments

The work utilizes data obtained by the Global Oscillation Network Group (GONG) Program, managed by the National Solar Observatory, which is operated by AURA, Inc. under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory. L.N. performed this work as part of Natural Resources Canada’s Public Safety Geoscience program. M.A.R., C.M. and T.A. acknowledge the Austrian Science Fund (FWF): J4160-N27, P26174-N27 and P31265-N27.

Appendix A PFSS Model

The PFSS model is based on the assumption that the magnetic field $\mathbf{B}$ above the photosphere is current free ($\nabla\times\mathbf{B}=0$). With this approximation, the magnetic field can be expressed as the gradient of a scalar potential $\Psi$,

[TABLE]

Using $\nabla\cdot\mathbf{B}=0$, which expresses the fact that there are no magnetic monopoles, we can write the Laplace equation for the potential $\Psi$ as

[TABLE]

The solution of the Laplace equation in spherical coordinates with $\theta\in[0,\pi]$ and $\phi\in[0,2\pi]$ in the region $R_{0}\leq r\leq R_{1}$ can be expressed as a infinite series of spherical harmonics expressed as

[TABLE]

where $R_{1}=2.5\ R_{0}$ is the source surface radius, $P_{n}^{m}(\cos\theta)$ are the associated Legendre polynomials, and $g_{n}^{m}$ and $h_{n}^{m}$ are coefficients that are computed from the input magnetograms. Using Eq.(A1), the solution for the magnetic field components is given by

[TABLE]

The magnetic field components are computed as follows:

[TABLE]

[TABLE]

[TABLE]

To compute the coefficients $g_{n}^{m}$ and $h_{n}^{m}$, we multiply Eq.(A3) for $r=R_{0}$ with $P_{n^{\prime}}^{m^{\prime}}\cos m^{\prime}\phi$ and $P_{n^{\prime}}^{m^{\prime}}\sin m^{\prime}\phi$, respectively. Using the orthogonality of the Schmidt normalized Legendre polynomials and integrating over the spherical surface, we find

[TABLE]

This yields the solution of the coefficients $g_{n}^{m}$ and $h_{n}^{m}$ in the form

[TABLE]

For the use of remeshed magnetograms as described in Section 2.1, we modify this relation to a discrete representation. To do so, we use the Clenshaw-Curtis quadrature rule given by

[TABLE]

where $\epsilon_{i}$ is $1/2$ for $i=0$ or $i=H$, and $1$ elsewhere. The weights $w_{i}$ are given by

[TABLE]

where $H=(H_{\theta}-1)/2$, $\epsilon^{\prime}_{k}$ is $1/2$ for $i=0$ or $i=H$, and 1 elsewhere. Using this equation, we write Eq.(A9) as

[TABLE]

Using Eq.(A4) and Eq.(A10), we compute the magnetic field components at any point in the region between the solar surface and the source surface (e.g., Altschuler and Newkirk, 1969; Nikolic, 2017).

Appendix B SCS Model

The solution for the Laplace equation not bounded by the spherical source surface is given by

[TABLE]

In the SCS model, solutions for the PFSS model for the magnetic field on the source surface are first oriented to point away from the Sun. This means that if $B_{\text{r}}<0$ on the source surface, the sign of all components $B_{\text{r}}$, $B_{\theta}$, and $B_{\phi}$ are reversed. In this study, we match the grid of the PFSS model to the SCS model with equal steps of $2\degree$ in latitude and longitude. The components of the magnetic field beyond the source surface are

[TABLE]

In this way, we ensure that no closed magnetic field exists beyond the source surface. We use the least-squares approach to best fit the vector field on the source surface. In order to minimize the sum of squared residuals, we write

[TABLE]

where $B_{k}(\theta_{i},\phi_{j})$ is the reoriented field of the source surface, and the index $k=1,2$, and 3 refers to the radial, latitudinal, and azimuthal fields component of the grid point $(\theta_{i},\phi_{j})$ of the reoriented magnetic field. Further, $\alpha_{nmk}$ and $\beta_{nmk}$ are

[TABLE]

We compute the derivations to minimize $F$ given by $\partial F/\partial g_{n}^{m}$, and $\partial F/\partial h_{n}^{m}$. For each $(n,m)$ we can write

[TABLE]

The same expression in matrix form is

[TABLE]

with

[TABLE]

Here, the dimension of $\widehat{B}$ is $3IJ\times 1$, the dimension of $\widehat{GH}$ is $(N_{s}+1)^{2}\times 1$, and the dimension of $\widehat{\alpha\beta}$ is $(N_{s}+1)^{2}\times 3IJ$. By choosing $\widehat{AB}=\widehat{\alpha\beta}\cdot\widehat{\alpha\beta}^{\top}$ we find that the solution for $GH$ is of the form

[TABLE]

This means that the solution of $GH$ requires the inversion of the square matrix $\widehat{AB}$. Finally, we use the estimates of $g_{n}^{m}$ and $h_{n}^{m}$ from the least mean square fit to compute the magnetic field solution above the source surface. We note that the polarity needs to be refined to match the polarities for $r\leq R_{1}$.