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

, , , , , , and

Published 2019 February 12 © 2019. The American Astronomical Society. All rights reserved.
, , Citation Martin A. Reiss et al 2019 ApJS 240 35 DOI 10.3847/1538-4365/aaf8b3

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0067-0049/240/2/35

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 behavior 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 model (DCHB), 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 (RMSE ≈ 83 km s−1). 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.

Export citation and abstract BibTeX RIS

1. Introduction

The knowledge of the evolving ambient solar wind is an essential component of successful space weather forecasting (Owens & 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 (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 & 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 solar radii (R0) to best match observations (Hoeksema et al. 1982). To account for Ulysses measurements showing latitudinal invariance of the radial magnetic field component (Wang & 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 & 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 (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 R0 to 2.5 R0 (PFSS), or 20 R0 to 30 R0 (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 an inner boundary condition for the heliospheric model (20–30 R0 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 & 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 lines highlight models accessible at NASA's Community Coordinated Modeling Center (CCMC) online platform.

Figure 1.

Figure 1. Overview of coronal and heliospheric model combinations for forecasting the ambient solar wind conditions in the near-Earth environment. The orange lines highlight models hosted at NASA's Community Coordinated Modeling Center online platform (see, https://ccmc.gsfc.nasa.gov/models/).

Standard image High-resolution image

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 & 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 & Sheeley 1990; Riley et al. 2001; Arge et al. 2003) and the Heliospheric Upwind eXtrapolation (HUX) model (Riley & 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 R0 to 5 R0, 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.

Figure 2.

Figure 2. Breakdown structure of the framework implementation. The sections explaining the corresponding framework component are indicated.

Standard image High-resolution image

2.1. Input to the Models

We use full-disk photospheric magnetograms from the Global Oscillation Network Group (GONG) from the National Solar Observatory 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 × 360 grid points, where θ ∈ [0, π] and ϕ ∈ [0, 2π] are the latitude and longitude coordinates, respectively. They are available as near-real-time magnetic maps or full CR maps at the GONG online platform.4 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 θ, (ii) contain both poles at θ = 0 and θ = π, 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 Mi,j with Nθ = 180 and Nϕ = 360 grid points in the latitude and longitude, respectively. The values for the extra grid cells at the poles M0 and ${M}_{{N}_{\theta }+1}$ are given by

Equation (1)

and

Equation (2)

The latitude of the uniform θ grid is

Equation (3)

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

Equation (4)

where the index i is selected so that

Equation (5)

and

Equation (6)

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

Equation (7)

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.

Figure 3.

Figure 3. Comparison of the observed line-of-sight component of the photospheric field and the PFSS model solution as a function of the heliographic latitude and Carrington longitude for CR2077 (i.e., 2008 November 20–December 17). The magnetic field strength is saturated at ±5 G. (a) Remeshed GONG synoptic magnetogram with Nθ = 181; (b)–(d) Magnetic field solutions Br, Bθ, and Bϕ computed from the PFSS model using the spherical harmonics expansion up to the order of 80 spherical harmonics.

Standard image High-resolution image

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 (R1 = 2.5 R0) is radial only (Altschuler & Newkirk 1969; Schatten et al. 1969). As detailed in Appendix A, Equation (19) can be solved to derive the magnetic field components at any point in the coronal domain (R0 ≤ r ≤ R1). Figures 3(b)–(d) present the magnetic field components Br, Bθ, and Bϕ 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 & 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 R0 to a distance of R2 = 5 R0, 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 Br, Bθ, and Bϕ are reversed if Br < 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° 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 ≤ R1.

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 R0 but use the PFSS solution at 2.3 R0 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 fp 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 & 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 & 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 R0 ≤ r ≤ R2, we use a fourth-order Runge–Kutta method (RK4) and solve the following set of equations, expressed as

Equation (8)

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 = R0) 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 = R2), we label the footpoint as an "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 = R2) and trace the field lines down to the solar surface. By doing so, we compute the areal expansion factor fp as

Equation (9)

where Br 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 fp at the photosphere and the source surface, respectively.

Figure 4.

Figure 4. Illustration of the topology of the magnetic field between the photosphere and the source surface at 2.5 R0 computed from the described magnetic model of the solar corona for CR2077. Red and blue indicate positive and negative magnetic field lines, respectively. While the lower plane shows the great-circle angular distance from the nearest coronal hole boundary at the photosphere d, the upper plane shows the areal expansion factor fp, indicating how much the corresponding flux tube expands between the photosphere and the source surface at 2.5 R0.

Standard image High-resolution image

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 (∝ρv2). 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, fp)) at a reference sphere of 5 R0 (or 30 R0 for the MAS model). In the literature, these empirical relationships are known as the Wang–Sheeley (WS) model (Wang & 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 & Sheeley 1990), namely

Equation (10)

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 Equation (10), we use v0 = 250, v1 = 660, and α = 2/5 (see, Arge & 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

Equation (11)

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 v0. For a footpoint located deep inside a coronal hole, the solar wind speed is equal to v1. Hence, the farther away the footpoint is from a coronal hole boundary, the faster the expected solar wind speed. In this study, we use v0 = 350, v1 = 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

Equation (12)

where α, β, γ, and δ are additional model parameters. For the coefficients in Equation (12), we use v0 = 285, v1 = 910, α = 2/9, β = 1, γ = 0.8, w = 2, and δ = 3. Figure 5 shows a comparison between the different empirical relationships for specifying the solar wind speed v(fp, d) near the Sun.

Figure 5.

Figure 5. Comparison of empirical relationships for specifying the solar wind speed near the Sun on the example of CR2077. The magnetic field lines illustrate the magnetic connectivity between the computed coronal hole regions and the sub-Earth orbit at the outer boundary of the coronal model (5 R0). The gray colored pixels indicate closed magnetic field topology with negative (dark gray) and positive (light gray) polarity. The full black and blue lines indicate positive and negative magnetic field lines, respectively. (a) WS model; (b) DCHB model; (c) WSA model.

Standard image High-resolution image

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 & Lionello (2011) developed the HUX model. Riley & Lionello (2011) 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

Equation (13)

where the subscripts r and ϕ refer to radial and longitudinal grid cells with cell spacings Δr and Δϕ, and Ω denotes the equatorial angular rotation rate of the Sun (neglecting differential rotation; see, Riley & Lionello 2011). In order to match the spatial resolution of the coronal model, we use Δϕ = 2° and Δr = 1 R0, respectively. As the step size in Δ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 & 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 & Lionello 2011), the residual acceleration vacc expected beyond the coronal model is

Equation (14)

where ${v}_{{r}_{0}}$ is the initial speed at the outer boundary of the coronal model, r is the heliocentric distance, α is a factor by which ${v}_{{r}_{0}}$ is enhanced, and rh is the scale length for the acceleration (not crucial when mapping the solar wind to 1 au). For the coefficients in Equation (14), we use α = 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 Δr from 5 R0 to 215 R0 (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 & Lionello 2011). This is particularly useful to study a large number of initial conditions in the context of ensemble modeling.

Figure 6.

Figure 6. HUX model for mapping the WSA solar wind speed solution to 215 R0 (or 1 au). (a)–(b) Solar wind bulk speed as a function of the heliocentric distance. (c) The input solar wind speed at 5 R0 (blue), and the resulting solar wind speed at 215 R0 with (yellow) and without (red) the residual speed contribution vacc.

Standard image High-resolution image

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, 2013b). 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 θp around the sub-Earth orbit given by

Equation (15)

where ϕ is the corresponding Carrington longitude; θE is the sub-Earth latitude; θ1 is the amplitude of the deviation; n is the wavenumber, indicating how many oscillations the perturbed trajectory completes; and ϕ0 is the phase offset, indicating how two perturbed trajectories can be out of synchronicity with each other. For the coefficients in Equation (15), we select θ1 ∈ [0°, 15°] in 1° steps, n ∈ [0, 1] in 0.5 steps, and θ0 ∈ [0°, 330°] in 30° 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 ±2σ 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 R0 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 & Stephenson 2003). As discussed in Owens et al. (2017) and Henley & 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.

Figure 7.

Figure 7. Illustration of the ensemble of initial conditions and the process of mapping the solar wind solutions to Earth by the HUX model. (a) Solar wind solution for CR2077 from the WSA model at 5 R0 overlaid with 3 out of 576 individual trajectories together with the median value (dashed red line) and the ±2σ quantiles (red line). (b) A comparison of the observed and WSA model solar wind speed time series at 215 R0. The ensemble median of solar wind solutions at Earth is indicated by the dashed blue line and the ±2σ quantiles are shown in blue.

Standard image High-resolution image

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). OSEA5 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, 2009b). 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 the arithmetic mean (AM), the standard deviation (SD), the mean error (ME), the mean absolute error (MAE), and the 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 RMSEs for the WS model, DCHB model, and WSA model are 85.27 km s−1, 103.43 km s−1, and 82.62 km s−1, 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., RMSE = 78.86 km s−1). 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.

Table 1.  The Skill of Model Predictions of the Solar Wind Speed for CR2077 in Terms of the AM, the SD, the ME, the MAE, and the RMSE

Model AM (km s−1) SD (km s−1) ME (km s−1) MAE (km s−1) RMSE (km s−1)
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

Note. 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.

Download table as:  ASCIITypeset image

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 model) 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., Arge et al. 2003; Nikolic et al. 2014) and found that our modular implementation in the C++ and Matlab programming language using tools from the Armadillo library (Sanderson & 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 & 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.

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 B above the photosphere is current free (∇ × B = 0). With this approximation, the magnetic field can be expressed as the gradient of a scalar potential Ψ:

Equation (16)

Using ∇ · B = 0, which expresses the fact that there are no magnetic monopoles, we can write the Laplace equation for the potential Ψ as

Equation (17)

The solution of the Laplace equation in spherical coordinates with θ ∈ [0, π] and ϕ ∈ [0, 2π] in the region R0 ≤ r ≤ R1 can be expressed as a infinite series of spherical harmonics expressed as

Equation (18)

where R1 = 2.5 R0 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 Equation (16), the solution for the magnetic field components is given by

Equation (19)

The magnetic field components are computed as follows:

Equation (20)

Equation (21)

Equation (22)

To compute the coefficients ${g}_{n}^{m}$ and ${h}_{n}^{m}$, we multiply Equation (18) for r = R0 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

Equation (23)

This yields the solution of the coefficients ${g}_{n}^{m}$ and ${h}_{n}^{m}$ in the form

Equation (24)

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

Equation (25)

where epsiloni is 1/2 for i = 0 or i = H, and 1 elsewhere. The weights wi are given by

Equation (26)

where H = (Hθ − 1)/2, ${\epsilon }_{k}^{{\prime} }$ is 1/2 for i = 0 or i = H, and 1 elsewhere. Using this equation, we write Equation (24) as

Equation (27)

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

Appendix B: SCS Model

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

Equation (28)

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 Br < 0 on the source surface, the sign of all components Br, Bθ, and Bϕ are reversed. In this study, we match the grid of the PFSS model to the SCS model with equal steps of 2° in latitude and longitude. The components of the magnetic field beyond the source surface are

Equation (29)

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

Equation (30)

where Bk (θi, ϕ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 (θi, ϕj) of the reoriented magnetic field. Further, αnmk and βnmk are

Equation (31)

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

Equation (32)

The same expression in matrix form is

Equation (33)

with

Equation (34)

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

Equation (35)

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 ≤ R1.

Footnotes

Please wait… references are loading.
10.3847/1538-4365/aaf8b3