1 Introduction

Variability in near-Earth solar-wind conditions can lead to a number of adverse effects on space- and ground-based technologies (Hapgood, 2011; Cannon et al., 2013). Forecasting ahead more than the approximately one hour afforded by the propagation time of solar wind from L1 to Earth requires prediction of the solar-wind conditions close to the Sun, which has yet to propagate to Earth. This is typically achieved using a coronal model (e.g. Linker et al., 1999; Arge et al., 2003; Toth et al., 2005) wherein the inner boundary conditions are determined by the observed photospheric magnetic field and the (open) outer boundary is set somewhere between 21 and 30 solar radii [R⊙], beyond the solar-wind Alfvén point. Solar-wind conditions are then propagated from 30 R⊙ to Earth using a heliospheric model (e.g. Riley, Linker, and Mikic, 2001; Odstrcil, 2003; Toth et al., 2005; Merkin et al., 2016; Pomoell and Poedts, 2018). Because the solar wind is both supersonic and super-Alfvénic, information only flows away from the Sun and no outer boundary conditions are required to generate any one model simulation run.

An example of a steady-state solar-wind solution is shown in Figure 1, using the Magnetohydrodynamics Algorithm outside a Sphere (MAS) global coronal model to compute the solar-wind speed at 30 R⊙ on the basis of the observed photospheric magnetic field, and the “HelioMAS” heliospheric model to propagate the solar wind from 30 R⊙ to Earth orbit (Linker et al., 1999; Riley et al., 2012). Carrington Rotation 1833 (spanning September 1990) is shown, as it generates a long-lived fast solar-wind stream at the Heliographic Equator. Figure 1a shows the radial solar-wind speed [\(V_{r}\)] at the corona/heliosphere model boundary of 30 R⊙. It is plotted as a latitude–time map for a fixed longitude in the frame of the Earth–Sun system (e.g. at Earth longitude). For steady-state solar wind conditions, this is an exact mirror image of a latitude–longitude map of \(V_{r}\). Figure 1c shows the same latitude–time map at 1 AU obtained from the HelioMAS model. Figures 1b and d show the time series of \(V_{r}\) at the Heliographic Equator at 30 R⊙ (e.g. at the sub-Earth point) and 1 AU (e.g. at Earth), respectively.

Figure 1
figure 1

An example of a three-dimensional numerical MHD solution of the solar wind, using the MAS coronal model and the HelioMAS heliospheric model. The photospheric magnetic field for Carrington Rotation 1833 (spanning September 1990) was used, as it results in both fast and slow wind in the equatorial plane. (a) Latitude–time plot at Earth longitude of \(V_{r}\) at the corona–heliosphere interface (30 R⊙) from MAS. (Note if the solar-wind structure is time stationary, this is an exact mirror image of a Carrington map.) The white line shows the Heliographic Equator. (b) The associated time series of \(V_{r}\) at the sub-Earth point. (c) Latitude–time plot at Earth longitude of \(V_{r}\) at Earth orbit (215 R⊙) from HelioMAS. (d) The associated time series of \(V_{r}\) at Earth from HelioMAS (black) and HUXt (red). See text for details of the HelioMAS and HUXt models.

Full size image

While no comparison with observed solar-wind speeds is shown here (the aim of this study is not to evaluate the existing solar-wind models), in general this approach provides a reasonable estimate of the steady-state component of the near-Earth solar wind (e.g. Owens et al., 2008; Jian et al., 2015, 2016; Reiss et al., 2016; MacNeice et al., 2018). However, the largest space-weather disturbances are driven not by the steady-state solar wind, but by dynamical structures associated with coronal mass ejections (CMEs) (e.g. Richardson, Cane, and Cliver, 2002). These disturbances are incorporated in the forecast models in an ad-hoc way. Coronagraph observations are used to characterize the trajectory, speed, and size of CMEs (Zhao, Plunkett, and Liu, 2002; Millward et al., 2013), which are then introduced at the inner boundary of heliospheric model (Odstrcil, Riley, and Zhao, 2004; Lee et al., 2013; Mays et al., 2015). While the internal magnetic-field structure of the CME is usually neglected, this approach enables a valuable estimate of CME arrival time at Earth (Riley et al., 2018).

In order for such forecasts to be useful, it is also necessary to have an estimate of their uncertainty. In its simplest form, this could be the average error in previous forecasts, although this treats all events the same and neglects the contextual information and knowledge that certain situations are inherently more predictable than others. The event-specific uncertainty can be estimated by using an ensemble of model runs with initial conditions perturbed to represent their uncertainty (Lee et al., 2013; Mays et al., 2015; Cash et al., 2015; Murray, 2018). The difficulty is the high dimensionality of the problem. Uncertainty in the steady-state solar wind likely represents a minimum of two degrees of freedom (to describe positional uncertainty, discussed later), and the CME properties typically consist of five parameters describing the CME speed, direction of propagation, width, and mass, bringing the total to at least seven dimensions. Modern three-dimensional magnetohydrodynamic (MHD) heliospheric models are extremely efficient, with the HelioMAS model (Riley et al., 2012) used in this study able to complete six days of simulation time in approximately 10 to 12 minutes on an NVIDIA TitanXP graphical processing unit (GPU). The ENLIL solar-wind model (Odstrcil, 2003) at low resolution requires similar resources, rising to approximately an hour on 36 processors at medium resolution. Fully exploring seven-dimensional parameter space is clearly prohibitive even at ten minutes per simulation run.

One solution is to adopt a simplified solar-wind model. The Wang–Sheeley–Arge (WSA) model (Arge et al., 2003) originally used a kinematic solar-wind propagation method that allows for stream interactions when propagating solar wind from the top of the corona to 1 AU. This produces \(V_{r}\) time series at 1 AU that show many of the observed properties (e.g. rapid rises in solar-wind speed, followed by more gradual declines). While this approach has been primarily used for steady-state solar-wind propagation, it could in principle be adapted for time-dependent boundary conditions and hence the study of CME propagation in large ensembles (Arge et al., 2004). However, this may not be ideal for many of the intended applications discussed below; it tracks plasma “particles” rather than solving the \(V_{r}\) on a discretized grid, making application of many standard data-assimilation (DA) methods problematic. Furthermore, the kinematic treatment is well tuned to match 1 AU observations, but it achieves this without explicit allowance for any “residual” solar-wind acceleration between the Sun and Earth (Schwenn, 1990). It will therefore likely overestimate solar-wind speeds inside 1 AU, which may subsequently affect CME dynamics.

A further level of abstraction is to assume \(V_{r}\) is prescribed and solve only for CME propagation, such as with a the one-dimensional drag-based model (DBM: Vrsnak and Gopalswamy, 2002; Cargill, 2004). While this typically assumes a uniform solar wind and provides no feedback between the solar wind and CME, it is very efficient and can be run in large ensembles (Dumbovic et al., 2018; Kay and Gopalswamy, 2018), with \(10^{5}\) ensemble members requiring only a few seconds on a modest desktop computer.

Here, we propose a compromise between the MHD and DBM approaches. The finite-difference model described below retains the structured solar wind and (some limited) feedback between the solar wind and CME, while making the necessary simplifications to give computation time comparable to the ensemble DBM approach. The goal is not to replace the MHD approach, but to build an adequate “surrogate” or “metamodel” (Blanning, 1975) of the MHD approach. This could be used to, e.g., further aid uncertainty quantification, explore parameter space in order to construct a suitable set of initial conditions for a (smaller) ensemble of full MHD runs, etc. While a surrogate could be constructed using a machine-learning approach (e.g. Lu and Ricciuto, 2019), we here adopt a more first-principles approach. The benefit is that this results in a simplified model that is also amenable to construction of an “adjoint” or inverse model, permitting data assimilation, a technique that shows great promise for improving the accuracy of solar-wind forecasting (Lang et al., 2017; Lang and Owens, 2019).

2 A Reduced Physics Approximation

In order to make the solar-wind model as efficient as possible, a large number of physical assumptions and approximations are required. While these individual approximations are partially justified below, the primary method by which we assess their collective validity is by making direct comparison with the full three-dimensional MHD (i.e. HelioMAS) solutions to the same inner boundary conditions. Note that the aim is not to construct a model that matches observations per se, but to construct a model that emulates the three-dimensional MHD with sufficient fidelity to act as a computationally efficient surrogate.

We begin by assuming that the solar wind can be approximated as a hydrodynamic flow and hence neglect magnetic forces, i.e. the plasma beta is large, as is expected above the mid-corona (Gary, 2001) and observed in the inner heliosphere by the Helios spacecraft as well as at 1 AU (Tu, Marsch, and Qin, 2004). Thus the fluid momentum equation becomes

$$ \frac{\partial \boldsymbol{V}}{\partial t} + (\boldsymbol{V}\cdot \nabla ) \boldsymbol{V} = - \frac{1}{\rho } \nabla P - \frac{\mathrm{G} \mathrm{M_{\odot }}}{r^{2}}, $$
(1)

where \(\boldsymbol{V}\) is the solar-wind velocity, \(\rho \) is the plasma mass density, \(P\) is the plasma pressure, \(\mathrm{G}\) is the universal gravitation constant, \(\mathrm{M_{\odot }}\) is the solar mass, and \(r\) is the radial distance from the Sun. In the heliosphere, pressure-gradient and gravitation terms are assumed to be small compared with the flow momentum and hence neglected. We further only consider variations in the radial direction. While non-radial flows are known to be significant in the vicinity of fast CMEs (Owens and Cargill, 2004) and CMEs observed to undergo non-radial deflections (Kay, Opher, and Evans, 2015), these effects are largest within 30 R⊙. Additionally, it has been argued that, to first order, CMEs in the heliosphere behave as individual radial elements due to the high expansion speeds involved (Owens, Lockwood, and Barnard, 2017). These approximations reduce the momentum equation to

$$ \frac{\partial {V_{r}}}{\partial t} + V_{r} \frac{\partial {V_{r}}}{\partial r} = 0, $$
(2)

where \(V_{r}\) is the radial velocity. This is the inviscid Burgers’ equation. Pizzo (1978) and Riley and Lionello (2011) further assumed time-stationary flows in the frame corotating with the Sun, in order to equate heliographic longitude with time and hence translate the time dependency [\(\partial / \partial t\)] into spatial dependency [\(\partial / \partial \phi \), where \(\phi \) is heliographic longitude]. This time-stationary method, referred to as Heliospheric Upwind eXtrapolation (HUX), was recently used to quantify uncertainty in steady-state solar-wind conditions (Owens and Riley, 2017; Reiss et al., 2019). Here, however, we retain the explicit time dependence and hence refer to the model as HUXt.

By 30 R⊙, the solar wind is highly supersonic, but observations and theoretical considerations indicate that the solar wind is still accelerating (Schwenn, 1990). In MHD simulations, this “residual acceleration” within the heliospheric domain is incorporated through the energy equation. Here, we explicitly add an additional velocity as a function of \(r\): \(accV (r)\). Thus in a uniform solar wind the velocity at a given \(r\) is

$$ V(r) = V_{0} + accV (r), $$
(3)

where \(V_{0}\) is the solar-wind speed that the plasma parcel had at the reference height \(r_{0}\), in this case the 30 R⊙ inner boundary. The \(r\)-subscript has been dropped for convenience.

On the basis of previous MHD simulation results, Riley and Lionello (2011) proposed \(accV\) of the form

$$ accV (r) = \alpha V_{0} \left [ 1 - \exp \left ( \frac{-(r-r_{0})}{r_{\text{H}}} \right ) \right ], $$
(4)

where \(\alpha \) is the acceleration factor and \(r_{\text{H}}\) is the scale height over which it applies. Riley and Lionello (2011) used \(\alpha = 0.15\) and \(r_{\text{H}} = 50~\text{R}_{\odot }\), which was found to match the HelioMAS \(V_{r}\) profile and has been verified within the numerical scheme presented here.

Substituting this expression into Equation 3 and rearranging gives:

$$ V_{0}= \frac{V(r)}{1 + \alpha \left [ 1 - \exp (\frac{r_{0}-r}{r_{\text{H}}}) \right ] }. $$
(5)

(Note that this estimate of \(V_{0}\) on the basis of a given \(V(r)\) is unique within uniform, unstructured solar wind. However, when time-dependent inner boundary conditions and solar-wind stream interactions are considered, a given value of \(V(r)\) will not have a unique value of \(V_{0}\). But just as the value of \(V(r)\) is effectively a weighted sum of the solar-wind speeds that have interacted within a distance \(r\), so too will be the value of \(V_{0}\).)

We take a finite-difference approach to these equations in order that they can be solved numerically. Adopting a first-order upwind scheme (Press et al., 1989), the time-dependent Burgers’ equation becomes

$$\begin{aligned} \begin{aligned}[b] V^{n+1}_{j} &= V^{n}_{j} - V^{n}_{j} \frac{\Delta t}{\Delta r} \left [ V^{n}_{j} - V^{n}_{j-1} \right ] \\ & \quad {} + V^{n}_{j-1} \frac{\Delta t}{\Delta r} \left [ accV^{n}_{j} - accV^{n}_{j-1} \right ], \end{aligned} \end{aligned}$$
(6)

where \(V^{n}_{j}\) is the radial speed at radial coordinate \(j\) and temporal coordinate \(n\). \(\Delta r\) and \(\Delta t\) are the radial and temporal grid steps, respectively. That is, \(\Delta r = r_{j} - r_{j-1}\) and \(\Delta t = t^{n+1}-t^{n}\), where \(r_{j}\) is the radial distance at radial coordinate \(j\) and \(t^{n}\) is the time at temporal coordinate \(n\).

Substituting for \(accV\) gives the full HUXt equation:

$$\begin{aligned} V^{n+1}_{j} =& V^{n}_{j} - V^{n}_{j} \frac{\Delta t}{\Delta r} \left [ V^{n}_{j} - V^{n}_{j-1} \right ] \\ &{} + V^{n}_{j-1} \frac{\Delta t}{\Delta r} \alpha V^{n}_{j} \left [ \frac{A_{j}}{1+\alpha A_{j}} \right ] \\ &{} - V^{n}_{j-1} \frac{\Delta t}{\Delta r} \alpha V^{n}_{j-1} \left [ \frac{A_{j-1}}{1+\alpha A_{j-1}} \right ], \end{aligned}$$
(7)

where

$$ A_{j} = 1-\exp \left [ \frac{r_{0} - r_{j}}{r_{\text{H}}} \right ]. $$
(8)

In order to make direct comparison with HelioMAS, we use the same radial grid, namely 140 uniformly spaced grid cells between an inner boundary of 30 R⊙ and an outer boundary of 236 R⊙. As HUXt is a one-dimensional model, there is no intrinsic grid in the azimuthal and meridional directions, but where equatorial and meridional plane cuts are shown, we use radial columns at each longitude and latitude point in the HelioMAS grid, namely 128 cells equally spaced in longitude and 111 cells equally spaced in sine latitude.

The time step [\(\Delta t\)] is chosen to meet the Courant–Friedrichs–Lewy (CFL) condition for numerical stability, i.e. that \(\Delta t < \Delta r / V_{\mathrm{MAX}}\), where \(V_{\mathrm{MAX}}\) is the maximum speed in the model domain. As HUXt contains a residual acceleration term, we conservatively set \(V_{\mathrm{MAX}}\) to be 50% larger than fastest speed at the inner boundary. Decreasing \(\Delta t\) by a factor of ten does not significantly affect the results, suggesting convergence has been reached.

HUXt is initialized with a solar-wind speed of 400 km s−1 at all radial distances. The time-dependent boundary conditions are used for five days before the model “spin up” is considered completed. Model results from the spin-up period are discarded.

3 Performance Evaluation: Ambient Solar Wind

A typical example of the level of agreement between the HUXt and the full three-dimensional MHD model (HelioMAS) solutions is shown in the bottom panel of Figure 1. Here, the photospheric magnetic field for CR1833 is used to generate a steady-state solar-wind solution with the MAS coronal model. The computed \(V_{r}\) at 30 R⊙ is used to construct a 27.27-day time series at the Heliographic Equator, shown in Figure 1b, which serves as input to HUXt. Note that 27 days of simulation time for HelioMAS would require somewhere between 30 and 60 minutes on a GPU. For HUXt, the computation time is 0.05 seconds on a modest desktop CPU.

It can be seen that at 1 AU, HUXt matches the general structure produced by HelioMAS, with a similar timing of \(V_{r}\)-structures and speed contrast between fast and slow solar wind. The equatorial plane cuts shown in Figure 2 show this agreement is present over all radial distances. The mean absolute error (MAE) between HUXt and HelioMAS over the whole equatorial plane is 15.7 km s−1, which is approximately 4%. The most obvious difference is at the slow-to-fast solar-wind transition. The upwind nature of HUXt means the transition is much steeper than in the MHD solution, where there is a more gradual rise. In this study, we are assessing how well HUXt can act as a surrogate for the full MHD solution and thus this difference is regarded as an error in the HUXt approach. But it is nevertheless important to remember that the MHD model should not be considered “perfect”, and one of the issues of numerical MHD is the diffusive nature, which means that it does not readily capture velocity sharp gradients present in observations (e.g. Jian et al., 2015; Shen et al., 2018).

Figure 2
figure 2

Radial solar-wind speed in km s−1 in the equatorial plane between 30 and 236 R⊙ for Carrington Rotation 1833 using (a) HelioMAS and (b) HUXt. (c) Absolute error in \(V_{r}\) between HelioMAS and HUXt.

Full size image

We now expand this same analysis to the whole 40+ year period for which steady-state HelioMAS solutions are available (i.e. the data set used by Owens, Lockwood, and Riley, 2017). For each of the 578 Carrington Rotations, HelioMAS and HUXt are computed for the same \(V_{r}\)-maps at 30 R⊙, produced by the MAS model. Figure 3 shows the MAE between HUXt and HelioMAS, averaged over the whole equatorial plane, as a function of time. The average MAE over the interval is 25.6 km s−1. As a percentage error, this is 6.4%.

Figure 3
figure 3

Comparison of HUXt and HelioMAS \(V_{r}\) in the equatorial plane averaged over the 30 to 236 R⊙ domain. Top: Carrington Rotation averages of \(V_{r}\) from HelioMAS (black) and HUXt (red). Middle: Carrington Rotation averages of \(V_{r}\) MAE between HUXt and HelioMAS. Bottom: MAE as a percentage error.

Full size image

To give a better indication of what this level of error means in terms of 1-AU solar-wind structure, the left-hand panels of Figure 4 show the three Carrington Rotations (CRs) with the highest percentage MAE in the whole 578-CR set. For all three CRs, there are numerous strong fast streams in the equatorial plane. The general solar-wind structures are in close agreement, with the largest differences being at the slow-to-fast wind interfaces, as illustrated for CR1833. The general trend for HUXt to overestimate the maximum speed at the stream interaction front may be a result of the one-dimensional assumption, as solar wind is unable to be deflected non-radially, as in the full three-dimensional HelioMAS (Pizzo, 1978). For the three best CRs, shown in the right-hand panels, there is a general absence of fast wind (although CR1833 shows that good agreement is also present when fast wind is present). Whether this level of agreement is adequate will, of course, depend on the application.

Figure 4
figure 4

Examples of the solar-wind speed at 215 R⊙ from HelioMAS (black) and HUXt (red) using the same input conditions. The left-hand side shows the three worst CRs out of the 578 considered. The right-hand side shows the three best.

Full size image

4 “Cone-Model” CME Example

We now compare HelioMAS and HUXt results for truly time-dependent boundary conditions. The aim is not to simulate a specific CME, but to illustrate the degree to which HUXt approximates HelioMAS. Both models are first “spun up” for the ambient solar-wind structure of CR2214, which spans mid-February to mid-March 2019. As can be seen in Figure 5, this is a typical solar-minimum solar-wind-speed structure with fast wind over the poles and mid-latitudes, and a narrow band of slow wind near the equatorial plane. As demonstrated above, the general ambient solar-wind-speed structure is in good agreement between HUXt and HelioMAS.

Figure 5
figure 5

An example of a fast CME-like disturbance into a solar-minimum solar wind. The left-hand panels show the equatorial plane at 0.45, 1.00, and 1.81 days after the leading edge of the CME-like disturbance passes the 30 R⊙ inner boundary. The position of Earth is shown by the black cross. The right-hand panels show the meridional plane at Earth longitude. The color scale shows solar-wind speed in km s−1. Top: The HelioMAS solution. Bottom: The HUXt solution.

Full size image

A circular-cross-section velocity perturbation, intended to approximate a CME, is added to the steady-state solar-wind inner boundary conditions. This approach is often referred to as a “cone-model” CME (Odstrcil, Riley, and Zhao, 2004), with the CME dimensions and velocity vector specified by coronagraph observations (Zhao, Plunkett, and Liu, 2002; Millward et al., 2013). Of course, this CME approximation neglects the CME’s internal magnetic field and thus is quasi-hydrodynamic even when used within an MHD solar-wind model such as HelioMAS, which is discussed further in Section 6. In this demonstration, the CME proprieties are chosen arbitrarily. Approximately three days into the simulation run, a CME-like structure is introduced at 30 R⊙ by raising the solar-wind speed to 1000 km s−1 in a circular-cross-section front spanning 30∘ about the Earth–Sun line. This elevated solar-wind speed is maintained for 12 hours, equating to a CME “thickness” of approximately ×3.5 the CME radius. Increasing this CME thickness reduces the rate at which it decelerates to the ambient solar-wind speed, effectively increasing its momentum. This is also often achieved by increasing the plasma density within the CME-like perturbation (Odstrcil, Riley, and Zhao, 2004).

In the equatorial plane, this CME-like perturbation propagates through HUXt and HelioMAS solar-wind solutions in a similar manner, with the front distorting asymmetrically due to the ambient solar-wind-speed gradient in longitude. The most obvious difference is again that HUXt produces a sharp front to the perturbation, whereas HelioMAS shows a more gradual speed gradient. In the meridional plane, the differences are greater, with the HUXt CME-like disturbance maintaining constant angular width, but the HelioMAS CME-like disturbance expanding meridionally into the fast wind.

Figure 6 shows the resulting \(V_{r}\)-time series in near-Earth space. The CME-like disturbance appears to arrive earlier in HelioMAS than HUXt. Taking the time of maximum \(V_{r}\)-gradient to define the CME arrival, the HelioMAS disturbance arrives 0.2 days (4.8 hours) ahead of the HUXt equivalent. We note, however, that the time of peak \(V_{r}\) agrees to within an hour. A more systematic study of many CME-like disturbances in MHD and HUXt, as well as with comparison with observations, will form the basis of future work. Here, we proceed to demonstrating the value of large ensembles of HUXt CME-like runs.

Figure 6
figure 6

The \(V_{r}\)-time series in near-Earth space for HelioMAS (black) and HUXt (red) simulations of the CME-like disturbance displayed in Figure 5. The HUXt ambient solar-wind solution (i.e. with no CME-like disturbance) is shown as the blue-dashed line.

Full size image

5 CME Arrival Time Sensitivity

To investigate the sensitivity of CME arrival time to initial condition-uncertainties, we again use the ambient solar-wind structure from Carrington Rotation 1833, which produced long-lived, fast solar-wind streams in the Ecliptic plane. Two CME cases are considered, as shown in Figure 7. The CMEs are identical except for the time into the Carrington Rotation at which they pass 30 R⊙, being five days past the start of the CR for Case A, and ten days for Case B. Both CMEs are directed along the Earth–Sun line, and they have a velocity of 1000 km s−1, an angular half-width with respect to the Sun of 30∘, and a CME thickness of ×0.5 the CME radius.

Figure 7
figure 7

CME-like disturbances inserted into the ambient solar wind for CR1833 (top) five days and (bottom) ten days after the start of the CR. The CMEs are otherwise identical. The left-hand panels show \(V_{r}\) in the equatorial plane for one, two, and three days after the CME insertion. The position of Earth is shown by the black cross. The black dots indicate the boundary of the CME-like disturbance. The right-hand panels show transit time from 30 R⊙ to 1 AU. The red line shows the value for the case shown, while the black curve shows the probability distribution function of 100,000 ensemble runs with perturbed input conditions, as described in the text.

Full size image

For Case A, the CME-like disturbance spans the fast–slow solar-wind boundary, distorting the CME shape significantly from its originally circular shape. Along the Earth–Sun line, the ambient solar wind encountered by the leading edge of the CME is primarily slow wind, and thus the transit time of the CME from 30 R⊙ to 1 AU is approximately 2.82 days. For Case B, the CME-like disturbance is primarily into fast wind, leading to a much reduced transit time of 1.72 days. Such strong control of CME transit time by the ambient solar-wind conditions is also found in three-dimensional MHD models (e.g. Case et al., 2008).

As HUXt is a computationally efficient code, it can be used in large ensembles to efficiently investigate the effect on the CME transit time of uncertainty in the initial CME parameters and ambient solar wind. Here, the uncertainty in each initial condition parameter is assumed to take the form of a Gaussian distribution centred on the best-guess value. Correctly specifying the width of the Gaussian (and thus the uncertainty in the parameter) requires significant further study (Kay and Gopalswamy, 2018). Here, we use arbitrary values to demonstrate the principle. Table 1 shows the Gaussian parameters from which initial conditions are randomly selected.

Table 1 Table of initial condition perturbations. Values are randomly drawn from a Gaussian distribution of mean \(\mu \) and width \(\sigma \). A minimum CME thickness of 0 is imposed.
Full size table

The ambient solar-wind structure is perturbed in the same manner as by Owens and Riley (2017) and Reiss et al. (2019). In summary, poorly observed polar photospheric magnetic fields are assumed to result in spatial errors in \(V_{r}\) at the 30 R⊙ inner boundary. Thus the \(V_{r}(30~\text{R}_{\odot })\) is rotated in latitude by a small angle [\(\mathrm{\theta _{SW}}\)] with the pole descending through a longitude \(\phi _{\text{SW}}\). We consider a Gaussian distribution of \(\theta _{\text{SW}}\), centered on 0∘, while \(\phi _{\text{SW}}\) is equally likely to take any value between 0 and 360∘ (i.e. the latitudinal rotation can proceed through any longitude). Once \(V_{r}(30~\text{R}_{\odot })\) has been rotated, the input to HUXt is again extracted along the Heliographic Equator.

Taking 100,000 random samples from the Gaussian initial condition distributions results in the probability distribution function (PDF) of transit times shown on the right-hand side of Figure 7. The best-guess transit time, in red, is reasonably close to the mode of the distribution in both cases (however, this need not always be the case). For Case A, the PDF is approximately symmetric; the mode is 2.85 days, with the half-maximum of the PDF spanning transit times of 2.71 to 3.03 days. The tails of the PDF are similarly symmetric about the mode, with minimum and maximum values of 2.37 and 3.48 days, respectively. For Case B, the PDF is more asymmetric. The mode is 1.74 days, with the half-maximum of the PDF spanning 1.68 to 1.86 days. The minimum and maximum values are 1.56 and 3.00 days, respectively.

Figures 8 and 9 show the sensitivity of CME transit time to different initial conditions for Case A and B, respectively. For Case A, which was primarily propagating into slow solar wind, perturbation of the ambient solar-wind conditions does not exhibit much influence on the CME transit time, although very large latitudinal rotations of the solar wind structure (\(>10^{\circ }\)) are associated with slightly reduced transit times. In general, this is to be expected, as faster wind is present at higher latitudes, which is brought down to the Equator, reducing the effective drag on fast CMEs in the equatorial plane. Surprisingly, the CME speed does not have a large effect on the transit time. The most influential factors are the CME width and thickness, both of which affect the effective momentum of the CME, reducing how quickly the CME is decelerated by the ambient solar wind. For Case B, which propagates into fast wind, the sensitivities to initial conditions are very different. The latitudinal perturbation of the ambient solar wind has a strong effect on the CME transit time, enabling both small reductions in transit time and large increases. While not displaying a strong correlation, the CME speed does place a hard lower limit on the transit time.

Figure 8
figure 8

Two-dimensional histograms of CME transit time as a function of initial conditions for the 100,000 perturbed initial condition runs of Case A.

Full size image
Figure 9
figure 9

Two-dimensional histograms of CME transit time as a function of initial conditions for the 100,000 perturbed initial condition runs of Case B.

Full size image

For simulated CMEs in a forecast scenario, quantitative measures of sensitivity to input conditions could be computed from such HUXt ensembles (e.g. Saltelli and Tarantola, 2002; Paruolo, Saisana, and Saltelli, 2013). This could enable more accurate and efficient seeding of the (smaller) initial condition ensembles for exploration with a full three-dimensional MHD model.

6 Discussion

In this study we have described and demonstrated HUXt; a computationally efficient solar-wind model that can be used for time-dependent solar-wind conditions, such as coronal mass ejections. For 40+ years of ambient solar-wind solutions, HUXt was shown to adequately emulate the full three-dimensional magnetohydrodynamic solution to the same 30 R⊙ boundary conditions, with solar-wind speeds throughout the simulation domain agreeing to within 6%. The biggest difference between the two solutions is that HUXt, by adopting an upwind numerical approach, produces much sharper boundaries between slow and fast wind than the MHD solutions.

Whereas an MHD solution to a (27-day) Carrington rotation requires between 30–60 minutes computing time on an NVIDIA TitanXP GPU, HUXt requires approximately 0.05 seconds on a modest desktop CPU. We stress that HUXt is not intended to replace the more complete physics-based approaches of three-dimensional magnetohydrodynamic solar-wind models such as ENLIL, SWMF, HelioMAS, HelioLFM, and EUHFORIA (Odstrcil, 2003; Toth et al., 2005; Riley et al., 2012; Merkin et al., 2016; Pomoell and Poedts, 2018). Instead, we suggest HUXt can act as a surrogate for these models when very large ensembles are required, which would otherwise be computationally prohibitive. This complementary “reduced-physics” approach, wherein a simpler model is used in a large ensemble alongside a single realization of the full-physics model, is well used in meteorology for both ensemble data assimilation and more generally for sensitivity analysis (Paruolo, Saisana, and Saltelli, 2013). Indeed, in future work, we intend to use HUXt for similar data-assimilation experiments, acting as an effective “adjoint” model for the three-dimensional MHD models, enabling variational data-assimilation methods (e.g. Lang and Owens, 2019) without the need to modify the MHD codes themselves.

We have also demonstrated that the main features of a three-dimensional MHD simulation of a CME-like disturbance are replicated by HUXt. While beyond the scope of this initial model description study, further validation of HUXt representation of CME-like disturbances is obviously required, with respect to both MHD models and observations. The notable difference in the case considered was the non-radial expansion of the CME-like disturbance in the latitudinal direction in the MHD simulation, which was absent in HUXt. Despite this, we note that the estimated CME transit times agreed to within zero to four hours (depending on how precisely the CME arrival is defined). It should be noted, however, that MHD simulations are capable of incorporating the internal CME magnetic field (e.g. Torok et al., 2018), although it is difficult to implement in a forecasting situation and is not currently used operationally. The effect of non-hydrodynamic effects on CME arrival times has yet to be quantified, although it is expected to be relatively small. Regardless, it is likely to be a systematic effect that could be parameterized within HUXt.

For operational CME arrival-time forecasting, the UK Met Office currently uses an ensemble of 24 ENLIL runs with slightly different CME parameters, such as speed, propagation direction, timing, width, and density. This ensemble size is typical for MHD modeling (Lee et al., 2013; Emmons et al., 2013; Cash et al., 2015). Obviously, 24 (or even 100) ensemble members is not adequate to sample seven-dimensional parameter space. Dimensionality will increase further when accounting for uncertainty in the ambient solar wind (which is known to have a large effect: Case et al., 2008), or the initial conditions of multiple interacting CMEs (Lee et al., 2015). Thus seeding the ensemble with the correct initial conditions to either provide an adequate measure of forecast uncertainty, or the full range of possible values from “perfect storm” conditions, is difficult and requires careful attention (Pizzo et al., 2015). We suggest that HUXt could be used for such purposes in a more brute-force manner.

Sensitivity of CME travel time to different initial conditions has been demonstrated here using HUXt to model two identical CME-like disturbances, one into slow wind, the other into fast wind. The resultant transit times from 30 R⊙ to 1 AU differ by more than a day. Using a large (100,000 members) ensemble of perturbed initial conditions shows that the transit times are also sensitive to different controlling parameters in these two cases. In these particular instances, for the slow-wind CME its size is the primary controlling factor. For the fast-wind CME, perturbations to the ambient solar wind and CME initial speed have the most effect on the transit time. These examples serve to highlight that the optimal initial conditions are context dependent and need to be evaluated on a case-by-case basis (see also Pizzo et al., 2015).

Finally, we note a number of other possible future uses for HUXt. Uncertainty in (and hence range of perturbation of) the ambient solar wind and CME parameter input conditions can and should be specified by direct analysis of those observations and models. But they could also be estimated by HUXt ensemble runs of a large number of observed CMEs (e.g. Riley et al., 2018) and the uncertainties calibrated so that the model probabilities match the observed occurrence frequencies (Johnson and Bowler, 2009). HUXt could also be used to quantify the uncertainty in CME arrival times using ambient solar-wind conditions estimated from different photospheric observations, coronal models, and realizations of photospheric evolution (Arge et al., 2013; Hickmann et al., 2015).