© K.A. Bunting & H. Morgan, Published by EDP Sciences 2022

1 Introduction

Rapid changes in solar wind conditions have a direct impact on Earth’s magnetosphere (e.g., Meziane et al., 2014), ionosphere (e.g., Milan et al., 2007), and both ground-based and space-based technology (e.g., Baker et al., 2004; Doherty et al., 2004; Imken et al., 2018). Estimates of the economic impact of space weather on European power grids alone range from €10s to 100s billion (e.g., Eastwood et al., 2018). Risk can be mitigated considerably by early warnings of impending space weather, currently undertaken by organisations such as the Meteorological Office in the UK. We believe that improvements in space weather forecasting depend primarily on three related categories: (1) Better observations of the Sun, corona, and solar wind; (2) Greater understanding of the physical processes operating in the corona and solar wind; (3) Improvements in data analysis and forecasting methods. This work presents a novel boundary condition to solar wind models based on tomography maps created from coronagraph observations of the solar atmosphere, thus an advancement that belongs to the third category and can contribute to the second.

The use of tomography maps as an inner boundary condition to a solar wind model, combined with an ensemble approach, plays a key role in the Coronal Tomography (CorTom) module to the Space Weather Empirical Ensemble Package (SWEEP) project. SWEEP is funded by the UK government’s Space Weather Instrumentation, Measurement, Modelling, and Risk (SWIMMR) scheme and will provide an operational space weather prediction package to the UK Meteorological Office by 2023. The SWEEP package operates a robust, complementary framework of multiple models using different boundary conditions, including the tomography described in this paper and both simple and more sophisticated magnetic models (Weinzierl et al., 2016; Yeates et al., 2018; Gonzi et al., 2021).

Approaches to space weather forecasting can be broadly placed in two groups: simulations that use observations of the Sun’s photosphere to drive a model of the solar wind (e.g., Wang–Sheeley–Arge (WSA)/ENLIL: Arge & Pizzo, 2000; Wang & Sheeley, 1990) and a persistence-based approach which extrapolates the future behaviour of the solar wind based on its past behaviour over various timescales (Owens, 2018). A persistence-based approach assumes that the ambient solar wind does not evolve drastically over a solar rotation. This is shown in observational data by a weak recurrence in geomagnetic activity and solar wind conditions (Diego et al., 2010) over a ~27-day period. Hence a persistence approach can provide a baseline for comparison of solar wind forecasting models (Owens et al., 2013). The heliospheric simulations are primarily based on remote solar observations and depend on photospheric magnetic field observations to build a model of the corona (e.g., MAS: Linker et al., 1999; Riley et al., 2012, AWSom: der Holst et al., 2014, etc.). It is the modelled conditions of the outer corona that forms the inner boundary for heliospheric solar wind models (e.g., ENLIL: Odstrcil, 2003, EUHFORIA: Pomoell & Poedts, 2018, HUXt: Owens et al., 2020).

The Magnetohydrodynamic Algorithm outside a Sphere (MAS) coronal model uses photospheric magnetic field observational data in order to gain the magnetic field configuration at the solar wind base (Linker et al., 1999; Riley et al., 2012). MAS derives an inner boundary condition at 30 R_⊙ that can be used in solar wind heliospheric models to predict solar wind conditions at 1AU (Owens & Riley, 2017). In this work, the MAS inner boundary condition is used to benchmark the results of the tomography-derived inner boundary condition. MAS adopts a simplified coronal “polytropic” model, which assumes the thermal pressure is greater or equal to the magnetic pressure at distances closer to the Sun (i.e., β ≥ 1). This is used to empirically convert the pressure, and density gained via simplified physical laws into a solar wind solution (Parker et al., 1964). Magnetohydrodynamic (MHD) equations are used and integrated forward in time until the solar wind parameters reach a steady state and give a full three-dimensional state of the solar wind at heights greater than the solar wind Alfvén point (Linker et al., 1999; Riley et al., 2012). MHD heliospheric models use this information at 30 R_⊙ as inner boundary conditions in order to propagate the solar wind conditions to 1AU (Odstrcil, 2003; Riley & Lionello, 2011; Owens et al., 2020). Three-dimensional MHD models offer a comprehensive physical model of the solar wind at large scales but are computationally expensive.

Riley & Lionello (2011) proposed that the magnetic, gravitational, and pressure gradient forces of the solar wind plasma can be neglected at distances greater than 30 R_⊙ and that a purely radial flow of the ambient solar wind plasma can be assumed, thus vastly reducing the complexity of the MHD equations. Heliospheric models that use this reduced physical approach have a greatly increased computational efficiency at the expense of reduced physics while still yielding results comparable to full 3D MHD models (e.g., Owens et al., 2020). Owens et al. (2020) have adapted this reduced physical model for the time domain, namely the “Time-dependent Heliospheric Upwind Extrapolation” (HUXt) model. The model complexity is reduced further when limited to only radial components of the equatorial plane, allowing the model to run on a standard desktop computer in seconds, even with a moderate angular (~2.8°), radial (1.5 R_⊙), and time (~4 h) resolution.

While space weather forecasting has developed enormously over the past few decades, there is still room for further improvement. One of the main constraints on the accuracy of space weather forecasting is the lack of adequate observational data that constrains direct empirical models between the photosphere and Earth. According to the concluding sentences of MacNeice et al. (2018): “the pace of [physical model] development has outstripped the pace of improvements in the quality of the input data which they consume, and until this is remedied, these models will not achieve their full forecasting potential”. This statement is a strong argument for new instruments and missions focused on operational space weather to provide higher quality data. It is also an argument for full exploitation of current resources: new or improved constraints on coronal structure and density at the coronal-heliosphere boundary are therefore important. Recent efforts focus on improving the diagnostics from photospheric observations through advanced magnetic modelling or using alternative observations such as radio scintillation (Gonzi et al., 2021). Our efforts are focused on using coronagraph observations through coronal rotation tomography to provide a direct constraint on solar wind models. A less direct approach developed by Poirier et al. (2021) is to use coronagraph observations to constrain magnetic models without resolving the line of sight (LOS).

There are ample, routine observations of the coronal-heliospheric boundary region by space-based coronagraphs (e.g., LASCO: Brueckner et al., 1995, COR2: Howard et al., 2008). These are largely neglected in the context of ambient solar wind modelling due mainly to the LOS problem: given the complex spatial distribution of high- and low-density streams along the LOS, it is impossible, from a single observation, to estimate this distribution. Coronal rotational tomography techniques aim to find a distribution of electron density in a 3D corona which best satisfies a set of coronagraphic polarised brightness (pB) observations made over half a solar rotation (half a rotation since both east and west limbs are observed), subject to some reasonable assumptions such as the smoothness of the reconstruction, thus helping to rectify the LOS problem. A comprehensive review is given by Aschwanden (2011). An iterative regularised least-square fitting method was presented by Frazin (2000) and developed and applied to other cases (e.g., Butala et al., 2005). A similar method has been applied to very low heights in the corona (Kramar et al., 2014) and expanded to include a time-dependency (e.g., Vibert et al., 2016). A novel method for creating qualitative maps of the distribution of coronal structure was introduced by Morgan et al. (2009), resulting in a comprehensive study of coronal structure over a solar cycle (Morgan & Habbal, 2010) and measurements of coronal rotation rates (Morgan, 2011a). Machine learning approaches are currently being developed (Jang et al., 2021). A new quantitative method based on spherical harmonics is presented by Morgan (2019), utilising the calibration and processing methods of Morgan (2015). Further refinements to the method, and initial results, are presented in Morgan & Cook (2020), and a study of coronal rotation rates based on tomography is made by Edwards et al. (2022). The method is based on a spherical harmonic model of the coronal electron density and gives reconstructions at distances of 4–10 R_⊙. At these heights and above, the corona can be assumed to flow radially outward. The tomography results when, compared at different heights, confirm this radial structure (Morgan & Cook, 2020). A desirable goal over the next decade would be an unified approach, where solar wind models are driven by coronal models based on as many empirical constraints as possible, including magnetic field extrapolations, coronal density estimations, and any other routine empirical constraints such as outflow velocity estimations. A crucial advancement would be the inclusion of a time-dependent inner boundary condition based on time-dependent magnetic models (e.g. Yeates et al., 2008; Weinzierl et al., 2016) and time-dependent tomography (e.g. Vibert et al., 2016; Morgan, 2021). Morgan (2021) states that the time-dependent method requires further development for routine use; thus, our current work uses the static tomography method of Morgan & Cook (2020).

This work aims to improve the accuracy of predictions of the ambient solar wind velocity at Earth and to investigate the relationship between coronal electron density and solar wind velocity at a distance of 8 R_⊙. Although iterative tomography methods have been coupled with MHD models to forecast space weather before (e.g., HELTOMO: Jackson et al., 2020), these models and studies are primarily focused on coronal mass ejections (CME) (Jackson et al., 2010) and use an iterative integrated kinematic model in order to predict ambient solar wind conditions (Jackson et al., 2013). Such heliospheric models incorporate CME models such as the cone model (Odstrcil et al., 2004), which use observational constraints on CME characteristics in order to model CME propagation throughout the heliosphere. CMEs are not considered further in this work, although our method and results are relevant to CME propagation and arrival time predictions.

This work proposes to use the tomography results as a new inner boundary condition for heliospheric solar wind models (HUXt), describes an initial implementation, and presents initial results compared to boundary conditions based on MAS as proof of concept. The methodology used in this study, as well as the simple empirical relationship used to derive solar wind velocity from density at a distance of 8 R_⊙, are described in Section 2. The iterative method used to improve the model’s match to in situ data through adjustment of input parameters is given in Section 3.1. MAS inner boundary conditions, when coupled with 3D MHD heliospheric models, have produced results that compare well with in situ (OMNI) bulk solar wind velocity data (e.g., Owens & Riley, 2017). The optimised model output is compared to results gained using MAS-derived inner boundary conditions and in situ data obtained via Operating Missions as Nodes on the Internet (OMNI) satellite network in Section 3.2. The tomography-based model is then applied to dates at different stages of the solar cycle 24 in Section 3.3. Section 3.4 explores the operational feasibility of the tomography inner boundary condition by adopting a persistence-based approach. Section 4 presents an implementation of an ensemble model that demonstrates how uncertainties can be quantified by the system. Conclusions are given in Section 5.

2 Method

This section gives a brief overview of the tomography method (Sect. 2.1), outlines the method to generate the inner boundary velocity condition from the tomography density maps (Sect. 2.2), gives an overview of the heliospheric solar wind model (Sect. 2.3), and an overview of the dynamic time warping (DTW) algorithm and how this will be exploited in the context of this study (Sect. 2.4).

2.1 Tomography maps of the solar corona

The COR2 coronagraphs are part of the Sun–Earth Connection Coronal and Heliospheric Investigation (SECCHI: Howard et al., 2008) suite of instruments aboard the twin Solar Terrestrial Relations Observatory (STEREO A & B: Kaiser, 2005). Density maps are calculated from COR2A data in 3 main steps:

• Calibration is applied to COR2 polarised brightness (pB) observations over a period of half a Carrington rotation (±1 week from the required date) using the procedures of Morgan (2015). The processing includes a method to reduce the signal from coronal mass ejections (CMEs) (Morgan et al., 2012). Following calibration, the data is remapped in solar polar coordinates, and an annular slice at the required distance from the sun centre is extracted over several hundred observations taken over time. Figure 1a shows an example of data used as input for the tomography. For this example, the central date is 2018/11/11 12:00, the data spans from 2018/11/04 00:00 to 2018/11/17 23:08, and the distance chosen is 8 R_⊙. This date range corresponds approximately to the mid-date of Carrington rotation (CR) 2210.

  Fig. 1 (a) The observed coronal polarised brightness at a distance of 8 R⊙, plotted as a function of position angle (measured counter-clockwise from solar north) and time. This is used as input to the tomography method. (b) The reconstructed polarised brightness, or synthetic observations gained from the tomographical density distribution. (c) The electron density at a distance of 8 R⊙ mapped in Carrington longitude and latitude.

• Tomography is applied to the data using the regularised spherical-harmonic optimisation approach of Morgan (2019). The method is based on a spherical harmonic distribution of density at the height of reconstruction, with an r^−2.2 decrease in density above this inner height; thus, a purely radial density structure is prescribed. Line-of-sight integrations are made of the spherical-harmonic based densities; thus, a set of brightnesses are gained, one for each order, with the lines-of-sight corresponding to the COR2A observations. The problem is then reduced to finding the coefficients of each order based on a regularised least-squares fitting between the observed brightness and the spherical harmonic brightnesses. In this work, a 22nd order spherical harmonic basis results in a density reconstruction with 540 longitude and 270 latitude bins for a chosen distance (restricted to between 4 and ≈10 R_⊙, with this range limited by the instrument’s useful field of view). In this work, we use only the reconstructions for a distance of 8 R_⊙. Note that the tomography reconstruction is static and has no time dependence.

• High-density streamers are narrowed, and a correction for “excess” density, possibly F-corona contamination (Morgan & Habbal, 2007c), is applied according to the method of Morgan & Cook (2020). The narrowing method is applied based on the gradients within the initial tomography density map, with the degree of narrowing controlled by a single parameter. This parameter is adjusted to find an optimal fitting to the observations. The “excess density” is estimated based on an analysis of densities within large low-density regions (coronal holes) between tomography maps made at a range of distances (4–8 R_⊙) and consideration of mass flux.

Figure 1a shows the observed coronal polarised brightness. This distribution is similar to the reconstructed polarised brightness shown in Figure 1b. The reconstructed brightness results from synthetic observations that use the electron density distribution. The tomography density map resulting from the above steps is shown in Figure 1c. This visualises the distribution of electron density at a distance of 8 R_⊙ resulting from the steps described above. Note that this is a static reconstruction since it gives a non-time-dependent density distribution that is spatially smooth and best satisfies the data.

The distribution of the reconstruction shows a good agreement with the observed, and all the large-scale streamer features are well reconstructed. However, the reconstruction lacks the fine-scale detail of the observed, including a general “fuzziness” of the brightness and temporal changes over timescales of less than a day. To reconstruct fine-scale detail, it is necessary to apply the time-dependent tomography described by Morgan (2021). There is considerable small-scale variation in the nascent slow solar wind (e.g., Alzate et al., 2021, and references within), and the evolution of this variation in the heliospheric wind is an active field of research. For operational forecasting of the ambient solar wind, our immediate problem is to model the large spatial scale and longer-timescale variations, and the static tomography reconstruction is adequate for this purpose.

2.2 Generating the inner boundary condition

We wish to use the tomography density map as an inner boundary condition to the HUXt model and to compare the resulting solar wind velocities near Earth with the results of using a MAS-based inner boundary. HUXt allows the user to set the heliocentric distance of the inner boundary. In this work, a distance of 8 R_⊙ is used for the tomography-based inner boundary and 30 R_⊙ for the MAS-based inner boundary. For the purpose of solar wind modelling in the Sun-Earth equatorial plane, a slice of density is extracted from the tomography map at Earth’s Carrington latitude at the initiation of the Carrington rotation (this latitude is 4.8° for the 2018 November example). The heliospheric model requires only 128 solar wind velocity values at equally spaced longitudes with increments of ≈2.8°. The extracted tomography longitudinal profile is rebinned to this size. Figure 2 shows the density profile at this latitude as a blue curve.

Fig. 2 Coronal electron density (blue, right axis) taken from the tomographical map of Figure 1c at Earth’s latitude, and solar wind velocity at 8 R⊙ (black, left axis) gained from the density using the simple inverse relationship of equation (1). For this example, Vmax = 480 km s−1 and Vmin = 220 km s−1.

The solar wind model requires a set of radial outward velocities as an inner boundary condition. We adopt a simple linear inverse relationship between densities and velocities, approximately consistent with the general properties of the solar wind as revealed by both in situ and remote observations (e.g., Habbal et al., 1997; Schwenn, 2006; Allen et al., 2020), given by:

$$V=V_{\max }-\left[\left(\frac{n-n_{\min }}{n_{\max }-n_{\min }}\right)\left(V_{\max }-V_{\min }\right)\right],$$(1)

where V and n are the solar wind velocity and electron particle density, respectively, n_max and n_min are the maximum and minimum densities in the equatorial plane, respectively, and V_max and V_min are model parameters specifying the maximum and minimum solar wind velocity at the height of the inner boundary (8 R_⊙). Whereas the density values (n_max and n_min) are defined by the tomography data at 8 R_⊙, the range of the velocity values is unknown. This highlights the need for an optimisation process with the aim of finding optimal values of V_max and V_min. This optimisation process is described in Section 3.1. Figure 2 demonstrates the conversion of the electron density to the solar wind velocity at the inner boundary for the 2018 November (CR2210) example with V_max and V_min set to their optimal values of 480 and 220 km s⁻¹, respectively. Note that these optimal values are found in the following section.

This simple inverse linear relationship of equation (1) is likely oversimplistic compared to the true relationship between density and velocity in the corona but serves the purpose of providing a proof of concept of a direct relationship for this study. Future efforts will experiment with optimising this empirical relationship; for example, an inverse Sigmoid function or an exponential relationship may better model the likely bimodal slow and fast wind patterns in the nascent solar wind and may reveal a global model that can provide an optimised model of the solar wind over several years, or a solar cycle. We note that converting the density into an estimate of velocity is similar to that used by magnetic models, where an empirical relationship is required to transform the magnetic field distribution to velocity (e.g., Gonzi et al., 2021).

2.3 Heliospheric upwind extrapolation model

HUXt is an incompressible solar wind model that solves a reduced Burgess equation along 1D vectors of velocities in the radial domain. The “upwind” conditions only allow outward flow and thus forbids any sunward motion of the solar wind. Therefore, considerations for stream interaction are included through adjustments in solar wind speed to uphold the upwind conditions (or at a greater distance in the radial coordinate) (Riley & Lionello, 2011). For the purposes of forecasting, the results of this model compare well to full 3D MHD models for both predictions of the ambient solar wind and CME events (Owens et al., 2020; Hinterreiter et al., 2021). In this work, we use the static tomographical reconstruction and, therefore, a time-constant boundary condition.

The HUXt model includes a parameter to account for residual solar wind acceleration. Riley & Lionello (2011) used an exponential function to impose a velocity profile that approached its final value asymptotically over a distance range between the lower boundary and ≈50 R_⊙ (an acceleration parameter of 0.15, built into the HUXt model). This work uses this acceleration parameter for both slow and fast wind streams. Investigating this parameter is a central focus of our future work: it is well known that the slow wind reaches its final velocity at a greater distance from the Sun compared with that of the fast (e.g., Schwenn, 1990).

HUXt uses a five-day “spin-up” time. All longitudinal and radial model points are initialised with a value of 400 km s⁻¹. The velocity at the inner boundary condition is rotated through ≈ −66° (equal to the solar rotation over five days), and the solar wind conditions are iteratively propagated outward while the boundary condition is rotated forward with time. Therefore, after the five-day “spin-up” period, propagated solar wind velocity values have reached distances beyond the orbit of Earth. Following the spin-up period, the solar wind velocity conditions are simulated forward in time by 27.27 days (or one rotation). Results generated during the spin-up period are discarded from the model output (Owens et al., 2020). For this study, the inner boundary condition is used as a static state which is not altered during the model run.

The evolution of solar wind velocity can be plotted as a function of time at any point in space within the model. In this work, we limit the results to Earth, allowing a direct comparison with in situ measurements. We use the reduced 5-min resolution combined satellite network data provided by OMNI. Sporadic short periods of missing solar wind velocity values are linearly interpolated and then binned to form an hourly average with an associated standard deviation. This is then smoothed with a 10-h moving window average.

2.4 Dynamic time warping

Dynamic time warping (DTW) is an effective algorithm to quantify the agreement between two-time series and is used here to compare the model and in situ solar wind velocities. The DTW algorithm was initially used to aid automated speech recognition and has recently been used in a variety of fields such as economics (Franses & Wiemann, 2020), biology (Skutkova et al., 2013), and space weather (Owens & Nichols, 2021; Samara et al., 2022).

DTW requires two vectors (1D arrays, A and B) as input. Vectors A and B are not required to be the same length. The Euclidean distance is calculated between every point in set A to every point in set B and is thus assigned a cost function for every possible alignment between the two data sets. This DTW cost or “DTW distance” metric is then minimised to find the optimal alignment between data set A and B (Berndt & Clifford, 1994). Effectively, DTW is non-linearly stretching and compressing each data set by connecting like-for-like structures between the sets. More efficient versions of this algorithm have been generated for a reduced computational expense, such as the Python FastDTW algorithm package used in this study (see: Salvador & Chan, 2004).

The requirements of the DTW algorithm are as follows:

• The two sets are ordered in time.

• Both sets start and end at the same time, or A₀ is anchored to B₀, and A_n−1 is anchored to B_m−1, where n and m are the number of elements in A and B, respectively.

• Every element in data set A will be matched with at least one corresponding data point from data set B, and vice versa.

• The optimum path must not “cross” between elements (for example, if A_i is paired to B_j, then A_i+1 or later cannot be paired with B_j−1 or earlier).

A metric used in this study in order to quantify the DTW distance is the sequence similarity factor (SSF). SSF, as defined by Samara et al. (2022), is described in equation (2):

$$\mathrm{SSF}=\frac{\mathrm{DT}{\mathrm{W}}_{\mathrm{Dist}}(O,M)}{\mathrm{DT}{\mathrm{W}}_{\mathrm{Dist}}(O,\overline{O})}$$(2)

where M represents the modelled solar wind velocity, and $\overline{O}$ represents the average magnitude of the observed in situ data (O). SSF allows a direct comparison of the modelled data to that of the average in situ data and provides a surrogate score of the model. For context, if SSF is > 1, the predicted solar wind velocities are worse than that of a constant mean observed solar wind value across the full-time period of the prediction. If SSF = 0, a perfect prediction has been made, and the modelled data matches the in situ data exactly.

3 Results

3.1 Velocity optimisation using dynamic time warping

This section describes and implements a simple method to derive optimised values of both the V_min and V_max terms and shows how this approach yields a far improved agreement between the model and measured solar wind velocities at Earth, as well as the time of arrival of fast solar wind streams.

Figure 3 demonstrates the effect of using arbitrary (and inaccurate) V_max and V_min velocity terms on the solar wind model values at Earth. V_max and V_min parameters are set at 600 km s⁻¹ and 250 km s⁻¹ (shown in Fig. 3a) and 380 km s⁻¹ and 110 km s⁻¹ (shown in Fig. 3b). This figure also visualises the DTW connections between the model and OMNI measurements along the optimum DTW path (shown in red), with the overall DTW distance being visualised by the sum of the length of the red lines. Both overestimation and underestimation of the velocity parameters will cause inaccurate predictions of the solar wind at Earth and thus give a larger DTW path distance compared to a model run with optimised or “best fit” velocity parameters. For example, the magnitude of the DTW path distance metrics for the model runs seen in Figures 3a and 3b are 6.03 × 10⁴ (SSF value of 1.34) and 7.66 × 10⁴ (SSF value of 1.71), respectively. Therefore, both model runs shown in Figure 3 provide a greater DTW distance compared to a constant average observed velocity across the full Carrington rotation. In order to obtain velocity values that give an optimal fit (or lowest DTW distance) between the model and in situ data, the efficiency of the HUXt model is exploited in an exhaustive search method. The tomography/HUXt model is run repeatedly with incrementally changing V_max and V_min terms, and the total DTW distance is recorded for each run.

Fig. 3 DTW optimal path (red) between in situ data and two Tomography/HUXt model runs for CR2210 with (a) overestimated Vmax and Vmin terms of 600 and 250 km s−1 respectively, and (b) underestimated Vmax and Vmin terms of 380 and 110 km s−1, respectively.

Figure 4a shows the relationship of optimal DTW path distance between the in situ and Tomography/HUXt model solar wind velocity at 1AU with incrementally changing velocity terms in equation (1). The velocity ranges used in this study are 50–350 km s⁻¹ and 350–650 km s⁻¹ for V_min and V_max, respectively in order to account for a wide velocity range for each parameter while also constraining each parameter to values that are consistent with physical solar wind velocities. We use 30 increments between these extremes (30 × 30 = 900 model evaluations). Figure 4a shows a minimum optimal DTW path distance for the tomography-derived inner boundary condition corresponding to velocity magnitudes of 220 and 480 km s⁻¹ for V_min and V_max, respectively. This is represented by the white cross in Figure 4a. These values are lower than observed on Earth due to the HUXt heliospheric model incorporating an acceleration parameter. A comparison between smoothed OMNI satellite data and the HUXt output with the optimal V_min and V_max is shown in Figure 5a and demonstrates a very strong correlation between the OMNI data and the output of the tomography/HUXt model (for statistical analysis, see Table 1). The time of arrival of the faster solar wind streams agrees to ±1 day, with the velocities agreeing with ±20 km s⁻¹.

Fig. 4 (a) Contour plot of DTW path lengths with a varying magnitude of the Vmax and Vmin terms for the tomography-derived inner boundary condition. The minimum DTW path distance is marked by a white cross, and this point gives the optimal values of Vmax and Vmin of 480 and 220 km s−1, respectively. (b) DTW path lengths as a function of varying scale parameters for the MAS-derived inner boundary condition with the white cross showing the optimum VMmin and VMmax parameters of 280 and 580 km s−1 respectively.

Fig. 5 Further analysis of the optimal DTW path. (a) Comparison of modelled solar wind velocities at Earth (grey) with in situ data (black) and DTW optimised path (red). (b) Alignment of the modelled and in situ data along the optimised DTW path with respect to time. (c) Difference in time between observation and modelled data of aligned points along optimal DTW path (ΔT(insitu−model)). (d) Velocity difference between observation and modelled data of aligned points along optimal DTW path (ΔV(insitu−model)).

In order for a fair comparison between tomography and MAS-based inner boundary condition (see Sect. 3.2), the MAS inner boundary condition requires a similar exhaustive optimisation process to be applied. This optimisation process effectively scales the initial MAS inner boundary condition between two values (VM_min and VM_max), with the aim of minimising the DTW path distance between the modelled and in situ data. During the optimisation process and for following comparisons, the MAS inner boundary height was set at 30 R_⊙. The results of this optimisation can be seen in Figure 4b, with the minimum DTW path distance corresponding to VM_min and VM_max values of 290 and 580 km s⁻¹, respectively.

Further details of the agreement of model and observation in the context of a DTW analysis are shown in Figure 5. Figure 5a visualises the DTW optimal alignment via the red lines. Figure 5b shows that the DTW path of the optimal alignment rarely deviates over two days from the “ideal” path, which represents a near-perfect agreement between model and in situ data shown in grey. The largest disagreement is seen at 18–24 days, which is seen in both the long red lines between data points in Figure 5a and by the biggest dispersion between the DTW optimum path and the ideal path seen in Figure 5b. This is due to a model overestimation of speed during this period. The histograms demonstrate the differences in time of arrival and the magnitude of solar wind velocity between the modelled and observational data along the optimum path. Figure 5c shows that there is a bias towards negative values. This suggests the model is predicting a later time of arrival with a mean ΔT_{(insitu−model)} of −0.45 days and a standard deviation of 1.16 days. Figure 5d is also biased towards negative values, suggesting velocity overestimation by the model. The mean ΔV_{(insitu−model)} is −3.64 km s⁻¹ with a standard deviation of 19.79 km s⁻¹.

3.2 Validation of model output

Figure 6 shows a comparison of the CR2210 tomography/HUXt velocity output (Blue) with OMNI in situ measurements (Black) at 1AU and MAS/HUXt velocity output (Orange). The tomography-based inner boundary parameters V_min and V_max are set at 220 and 480 km s⁻¹, respectively. The MAS inner boundary condition has scale parameters of VM_min and VM_max of 290 and 580 km s⁻¹ as described in Section 3.1.

Fig. 6 Comparison of optimised tomography/HUXt (blue) and optimised MAS/HUXt (orange) model predictions of the solar wind velocity at Earth with insitu data (black) for CR2210.

From Figure 6, both models have similar profiles, and the main changes between slow and fast wind tend to agree. The first small peak (shown in the tomography data 2018, October 29–31) is not present in the MAS data. The solar wind velocity between the second (2018, November 2–5) and third peak (2018 November 10–14) drops to intermediate velocities (≈440 km s⁻¹) for the tomography-driven model, while the MAS model drops to slower velocity (≈350 km s⁻¹). This is significant as the in situ data, as shown in black in Figure 6, shows an intermediate velocity (≈440 km s⁻¹) of the more undisturbed solar wind between the two fast peaks (2018 November 7–10). The magnitude of solar wind velocity during this time is better represented by the tomography inner boundary condition model.

Table 1 presents statistical details of the comparison between models. The Pearson correlation coefficient between in situ measurements and models is considerably higher for the tomography/HUXt model. The mean absolute error (MAE) of velocities between measurement and MAS/HUXt is approximately 11% higher than the tomography/HUXt, while also showing a higher SSF. These results show that the general profile and magnitudes of solar wind velocity are closer to the data using the tomography boundary condition compared to the optimised MAS model. This demonstrates both the potential value of tomographical density maps as an inner boundary condition and the benefit of searching for an optimised velocity range at the inner boundary. Note that such an optimisation would not be possible without an efficient solar wind model such as HUXt.

Both models fail to give the exact time of arrival of various features. For example, the rapid rise from slow to fast wind observed on November 4 arrives approximately one day early in both models. The start of the decrease from fast to slow wind on November 13 comes late in the tomography/HUXt model, and the decrease is less rapid than observed. The main reasons for these differences are listed and discussed in Section 5. Despite the differences, the comparison of tomography/HUXt to in situ data is promising – the large-scale features of the measured velocity are present in the predicted velocities, and the timings are reasonable given the initial simple implementation of the inner boundary condition. The model fails to replicate smaller-scale structures on timescales of a day or less. One reason for this is that the tomography densities are inherently smooth – this is an unavoidable result of finding a static tomographical solution which is discussed further in Section 5. Other reasons for differences between model and measurements include the reduced physics approach of HUXt and the limited resolution of computational modelling.

3.3 Application to other dates

Here we apply the Tomography/HUXt method to two different periods during the solar cycle 24. The model is applied to 2014 May (CR2150, near solar maximum) and 2018 March (CR2202, at the start of the current solar minimum). The tomography maps for these two dates are shown in Figure 7. Figures 8a and 8c show the comparison of modelled solar wind velocity at Earth and in situ data for CR2202 and CR2150 respectfully, with optimised velocity parameters which are gained from the contour plots seen in Figures 8b and 8d.

Fig. 7 Electron density as estimated using coronal tomography at a distance of 8 R⊙, for Carrington rotations (a) 2150.5 and (b) 2202.5. The red boxes labelled A and B in (a), and the white box in (b) are relevant for the ensemble results of Section 4.

Fig. 8 (a) Comparison of model and OMNI solar wind velocities for CR 2202 for the optimal fit in terms of minimum DTW path distance (marked by a white cross in Fig. 8b), with Vmin and Vmax set to 230 and 440 km s−1 respectively. (b) Contour plot showing the DTW path distance as a function of Vmin and Vmax for CR 2202. (c) Same as (a), but for CR 2150, with Vmin and Vmax set to 240 and 400 km s−1, respectively (represented by white cross in Fig. 8d). (d) Same as (b), but for CR2150.

For CR2202, Figure 8b demonstrates a minimum optimal DTW path distance at a V_max value of 440 km s⁻¹ and a V_min value of 230 km s⁻¹. For this period, the model and observation data agree well. However, one disparity is present between 2018 March 29 – 2018 April 4, where the in situ data shows a region of higher solar wind velocity, which is not present in the modelled data. The mean velocity difference between in situ and modelled data along the optimal DTW warped path is 0.49 km s⁻¹. The time domain also shows an acceptable agreement, with a mean time difference of 1.12 days.

CR2150 spans an active period close to the height of the solar maximum. Figure 8c shows a significant fast solar wind stream in 2014, May 22–25, which is seen in both model and measurement. However, there is a large discrepancy between a peak seen in the model in 2014, May 14–16 and one seen in measurement 3 days earlier. Figure 8d demonstrates a minimum optimal DTW path distance at 400 and 240 km s⁻¹ for V_max and V_min, respectively. The velocity difference has a mean value of −0.9 km s⁻¹. These values are close to the previous case studies. However, the time difference is much greater in comparison to previous Carrington rotations, with a mean of −1.6 days and a standard deviation of 2.3 days. Such a deteriorated agreement is likely due to the increased solar activity during this time, which both disrupts the tomography process, makes the time-independent static tomography approach less valid, and increase the chance of CMEs in the in situ measurements.

Table 2 shows a statistical comparison between both optimised tomography/HUXt and MAS/HUXt models with in situ data for CR2150 and CR2202. For CR2150, the tomography/HUXt model combination yields a lower MAE (38.44 km s⁻¹) compared to the MAS/HUXt model (40.12 km s⁻¹). The MAE for CR2202 offers a significant 32% reduction for the tomography/HUXt model (39.31 km s⁻¹) compared to that of the MAS/HUXt model (57.93 km s⁻¹). For both periods, the tomography/HUXt model shows a smaller DTW path distance and a smaller SSF than MAS/HUXt model.

The model solar wind speeds for Carrington rotations 2150 and 2202 show a more gradual transition between high and low solar wind velocities than in the observational data. This is likely due to the smoothness of the density given by the tomography and the upwind dependence of HUXt. The significant velocity overestimation of ≈75–100 km s⁻¹ seen between 2014 May 27–29 (Fig. 8c) could be due to the upwind dependence of HUXt, or due to slow wind from equatorial coronal holes. Parker solar probe (PSP) has detected low density – low-velocity solar wind structures at distances close to the Sun that are thought to originate from equatorial coronal holes (see: Bale et al., 2019; Kasper et al., 2019). This specific structure of solar wind contradicts the simple inverse relationship of equation (1), which assumes a consistent low density – high-velocity relationship. This will result in overestimating solar wind velocities at the inner boundary and, therefore, at Earth. A low density – low-velocity structure could well explain the velocity overestimation of 2014 May 27–29, as this region maps back to approximately 150° longitude, a section between two coronal holes (see Fig. 7a). However, this area has obvious key implications for solar wind forecasting that demands further research.

3.4 A persistence approach

In this section, a persistence-based approach is adopted in order to attempt to predict solar wind velocities in a realistic operational context. Unlike a traditional persistence model, which predicts a near-exact repetition of the observed ambient solar wind conditions for the solar rotation prior (see Owens et al., 2013), in this case, the inner boundary condition is updated. The coronal densities (and therefore n_max and n_min terms in Eq. (1)) are extracted as described in Section 2.2. However, this inner boundary condition will not undergo the exhaustive optimisation process described in Section 3.1 but instead, use the optimised V_max and V_min terms unchanged from Carrington rotation prior. For example, Figure 9a shows the comparison between in situ data and the forecast data for CR2203 with V_max and V_min values of 440 and 230 km s⁻¹, which are the optimised parameters for CR2202. Likewise, Figure 9b compares in situ data with forecast data for CR2151 with V_max and V_min of 400 and 240 km s⁻¹ (optimised values for CR2150).

Fig. 9 Comparison of predicted solar wind conditions of (a) CR2203 (Vmax and Vmin of 440 and 230 km s−1 respectively) and (b) CR2151 (Vmax and Vmin of 400 and 240 km s−1 respectively). The Vmax and Vmin values are that of the optimised velocity terms of the Carrington rotation prior (obtained in Sect. 3.3).

Figure 9a shows a relatively good agreement with the in situ data in terms of the time of arrival of the fast solar wind streams (see Table 3 for statistical analysis). However, there are disparities in the magnitude of the peaks. For example, in the first peak (seen in the in situ data on 2018, April 21–22), the in situ data predicts a velocity of ≈600 km s⁻¹, whereas the model predicts a solar wind velocity of ≈500 km s⁻¹ . The second peak (2018 May 6–11) shows a similar difference. This suggests that the V_max term is underestimated. Likewise, the slower, more settled solar wind present in between these two peaks (2018, April 25 – 2018, May 5) is underestimated by the model, suggesting that the V_min term is also underestimated. Table 3 shows an SSF value of 0.65 for CR2203. This shows that a persistence approach to the tomography-based model can yield more accurate results compared to a single mean value.

Figure 9b shows a weaker agreement between the modelled and in situ data. The profiles are different, with the in situ data demonstrating a double peak between the dates of 2018, June 8–13, whereas the model predicts a single peak around this time. The magnitude of this peak also disagrees by ≈150 km s⁻¹. This again suggests that the V_max term is underestimated. Table 3 shows generally weaker statistics for CR2151 compared to that of CR2203. An SSF value of 1.20 infers that a mean solar wind velocity across the full-time period would yield a smaller DTW distance and potentially better predictions of solar wind velocity compared to a persistence-based approach during this time period. During periods near solar maximum (such as CR2151), we would expect a persistence-based approach to yield worse results. This is due to the coronal state changing at a significantly faster rate than during solar minimum and the tomography reconstruction failing to accurately map the true coronal density. Therefore, the inner boundary condition will differ significantly from the physical state of the solar corona. This highlights the need for a time-dependent tomography approach in an operational context.

Overall, here we show that a persistence model could potentially be used in an operational context as a worst-case scenario without updating the V_min and V_max parameters. However, the model certainly loses accuracy when deployed in this fashion, especially during periods near the solar maximum. This stresses the need for either a time-dependent boundary condition or a more complex, global relationship between coronal electron density and solar wind velocity at the height of 8 R_⊙. Both of these issues are the focus of our current efforts.

4 An ensemble approach

The high efficiency of the HUXt model allows an ensemble approach to estimating the uncertainty in model outputs based on selected uncertainties at the inner boundary. The tomography density distribution, as shown in Figures 1 and 7, shows thin elongated structures that tend to lie longitudinally and can be very narrow in latitude. Therefore, a small error in the distribution given by the tomography method, or small latitudinal deviance of the solar wind during propagation to Earth, can significantly alter the inner boundary condition and the predicted solar wind conditions at Earth. An ensemble approach is a straightforward way of investigating and quantifying the effect small variations in the latitude of the extracted tomographical data at the inner boundary can have on the resulting solar wind velocities on Earth. Another uncertainty to which we can apply the ensemble approach is the choice of V_max and V_min. Here, we use the map of the DTW cost function arising from the exhaustive search of V_max and V_min values to define a range of velocity terms at the inner boundary for 2014 May (CR2150) and 2018 March (CR2202).

4.1 Latitudinal dependence

Figure 7 shows two density maps that demonstrate two extremes of solar activity. CR2150, as shown in Figure 7a, demonstrates a complicated density distribution with many high-density streamers positioned at a wide range of latitudes and which span across the equator. Figure 7b shows a quiet solar corona where the streamer belt is longitudinally aligned near the equator and a more uniform low electron density at higher latitudes. The main density enhancements are found exclusively in the equatorial region.

A model run was conducted as described in Section 3.3 with V_max and V_min remaining fixed at the optimised values of 400 and 240 km s⁻¹ for CR2150 and 440 and 230 km s⁻¹ for CR2202. We adjust the inner boundary velocity profile by extracting densities from the tomography map with ±3° of the latitude of Earth with 1° increments. This range of latitudes was chosen with the aim of comfortably covering the latitudinal movement of Earth during one full Carrington rotation and, more importantly, the unknown drift of the solar wind in latitude between the Sun and Earth. Comparisons between the resulting model and measured solar wind conditions on Earth are shown in Figure 10. The largest variations in the ensemble velocities during CR2150 (see Fig. 10a) are present before 2014 May 13. During this time, we find that the relatively small latitude variation of ±3° leads to a wide variation in velocity. For example, the largest velocity range is ≈50 km s⁻¹, observed on 2014 May 10. For the remainder of the period, varying the latitude leads to only a small variation in velocity (≈20 km s⁻¹). This can be explained by examining the density map in Figure 7a. The solar wind streams at Earth, corresponding to the early part of the period, map back approximately to the region bounded by the red box labelled A in the tomography map. This location was calculated by considering the solar wind speed, the distance from Earth to the tomography map location, and the solar rotation rate. This region spans the boundary between a high-density streamer to the north and a low-density coronal hole to the south. Therefore, the northern latitudes in Figure 10a yield considerably slower solar wind velocities compared to southern latitudes during this time. Thus varying the latitude by small values can lead to large changes in the density profile and, therefore, the inner boundary velocity.

Fig. 10 Comparison of in situ data with the results of multiple tomography/HUXt model runs for (a) CR2150 and (b) CR2202, with the latitude of the inner boundary condition varied in one-degree increments between ±3° from the latitude of Earth.

All ensemble model runs fail to match the steep decrease in velocity during 2014 May 25. The coordinates of this decrease near the Earth map back to the region bounded by the red box labelled B in the tomography map. This is a narrow region of low density lying between longitudinally extended regions of higher density to the north and south. Either the tomography map is incorrect in this small area, and that this region should contain a higher density and thus map to slow velocities, or the simplistic inverse linear relationship mapping densities to velocities, as given by equation (1) fails in this region. If the latter, then the relationship should not be linear and should give the highest velocities only for the very lowest density features in the tomography map. Current efforts are focused on investigating an improved relationship between coronal electron density and solar wind velocity at 8 R_⊙.

All ensemble members consistently predict a higher velocity peak on a date 2014, May 14–16. This feature is not seen on these dates in the in situ data. This peak maps back to the small low-density feature at Carrington longitude 270° as seen in the equatorial region of the tomography map seen in Figure 7a. We note that the density of this feature is low but not at the minimum densities estimated for coronal holes seen in this map and others. Conversely, the peak at intermediate velocities seen in the data around 2014 May 12 is not present on this date in any of the models runs. This feature maps back to a longitude of around 305°, close to the boundary between the low-density region at 270°. It is likely that the velocity peak seen in all ensemble members (2014, May 14–16) is the same peak seen in the in situ data during 2014, May 12. The amplitude of both peaks is comparable. There are several reasons why the model predicted a later time of arrival at Earth for this specific structure, including the modelling limitations of HUXt (e.g., the acceleration parameter) or a defect in the tomography map.

Figure 11 shows how the density of the corona changes in the time of just three Carrington rotations. During CR2149, shown in Figure 11a, there is a large and very low-density coronal hole spanning the equator and reaching to high latitudes. There is a clearly defined western boundary to this coronal hole at longitude 280° near the equator. By CR2150, as shown in Figure 11b, the coronal hole has greatly reduced in size and does not reach the low densities of the previous rotation. Consequently, the western boundary is not clearly defined. The following rotation in Figure 11c shows the western streamer encroaching into the region previously occupied by the coronal hole and the coronal hole limited to southern regions. This rapid change in structure is the best explanation for the differences between the model and in situ measurements between dates 2015, May 12 and 18, and shows that a time-dependent inner boundary becomes critical during solar maximum.

Fig. 11 The rapidly changing density of the corona between (a) CR2149, (b) CR2150, and (c) CR2151.

The results for CR2202 in Figure 10b show that the variability in solar wind prediction is largest between 2018, April 1–8, where the model velocities vary by more than 100 km s⁻¹. This region maps back to the white boxed region in Figure 7b. This region spans the northern boundary of a high-density streamer, thus small variations in latitude lead to large variations in density and velocity. This result shows that even during quiet periods, large uncertainties can arise from small deviations in latitude. This is a significant problem for solar wind forecasting, particularly considering that the high-density streamer belt may actually be narrower than that reconstructed using tomography. The measured fast wind peak during 2018, on April 11 reached speeds of almost 600 km s⁻¹. The models over all latitudes consistently underestimate this peak. This feature maps back to the low-density equatorial region near Carrington longitude 145°. This is likely due to the tomography map overestimating the density at this point, a flaw in the oversimplistic relationship between density and velocity, or a combination of both.

4.2 Velocity dependence

Figure 12 shows the DTW path distance as a function of V_max and V_min, used to select the optimal parameters. In order to create a velocity-based ensemble, we wish to identify a region within this parameter space where the DTW distance is lower than a set threshold.

Fig. 12 Contour plots of the DTW path distance as a function of minimum and maximum inner boundary velocity for Carrington rotations (a) 2150 and (b) 2202. The white contour bounds the area that yields a DTW path distance under the threshold chosen (SSF ≤ 0.95), and the white pixels in the lower right corner represent invalid model parameters where Vmin > Vmax. Note that the Vmax range is altered to 300–600 km s−1 in this plot.

Figure 12 shows that there is a selection of velocity values that give an acceptable fit for in situ data, defined as SSF ≤ 0.95. This range of velocity magnitudes defines an uncertainty which is incorporated into the ensemble. The velocity terms for each ensemble member are randomly generated from a normal distribution, with a mean set at the optimal velocity parameters and a range within three standard deviations of the values that yield an SSF of 0.95 or lower. These regions are highlighted by the white contour in Figure 12. For example, CR2150 has mean (optimal) values of 400 and 240 km s⁻¹ and a standard deviation of 20 km s⁻¹ for V_max and V_min, respectively. For CR2202, mean velocity values are 440 and 230 km s⁻¹ with a standard deviation of 30 km s⁻¹. Care was taken in order to ensure the V_max term was greater than the V_min term for every ensemble member.

Figure 13 compares in situ data with seven tomography/HUXt outputs, each with a different combination of V_max and V_min magnitudes for CR2150 and CR2202. This highlights the sensitivity of the output to the velocity range. An underestimated V_max value will affect the magnitude of the solar wind peaks predicted at Earth but will also cause the fast solar wind peak to arrive later. For example, in Figure 13a, the model with V_max = 350 km s⁻¹ (red) shows a 2–3 day delay in the time of arrival of the fast solar wind on 2018 May 26, compared with a model run with V_max = 460 km s⁻¹ (green). Both these terms dictate how rapid the transition between fast and slow solar wind occurs.

Fig. 13 Comparison of in situ data with the results of multiple tomography/HUXt model runs for (a) CR2150 and (b) CR2202, with a range of Vmax and Vmin combinations.

4.3 Ensemble results

We create an ensemble with 10,000 unique inner boundary conditions generated with randomly selected pairs of V_max and V_min taken from a normal distribution around the optimal velocities and a line of latitude randomly selected between ± 3° from Earth’s latitude at the initiation of the Carrington rotation. Once generated, the 10,000 inner boundary conditions were used for the HUXt model, and the results of solar wind velocities at Earth were recorded. The results of the ensemble are shown in Figure 14. A choice of 10,000 ensemble runs was a compromise between quantifying the sensitivity of the results to the parameter uncertainties with a sensible computation time.

Fig. 14 Results of the ensemble for (a) CR2150 and (b) CR2202 with the mean predicted solar wind velocity at each time step shown as white, the light grey area representing all predicted solar wind velocities, and the darker grey representing one standard deviation from the mean. The non-ensemble model results with the optimal velocity values at the Carrington latitude of Earth is the blue line.

Figure 14a shows a strong correlation between the mean of the ensemble runs (white line) and the optimal HUXt model for the true latitude of Earth (blue line), with a maximum velocity difference of 12 km s⁻¹. The standard deviation of the ensemble around the mean is shown as the dark grey region, and the average standard deviation across the full-time series is 20.06 km s⁻¹. The most significant deviations between the in situ data and the model data lie at three intervals centred on 2014, May 12, May 15 (as discussed in Sect. 3.3), and May 27. The final disparity shows a steep transition between the fast and slow solar wind (seen around 2014 May 28 in the in situ data) that is not predicted by any of the 10,000 ensemble models runs.

CR2202, as seen in Figure 14b, also demonstrates a good agreement with the in situ data with the maximum velocity difference between the mean and the in situ data being 79 km s⁻¹. CR2202 shows a greater variability of the ensemble model runs (as shown by the height of both dark grey and light grey regions) compared to that of CR2150, with the average standard deviation across the full time-series being 41.7 km s⁻¹. This is approximately double that of CR2150. It would be expected that the value is larger than that for CR2150 due to the larger standard deviation entered into the velocity terms (20 km s⁻¹ and 30 km s⁻¹ for CR2150 and CR2202, respectively).

5 Discussion and conclusions

A new inner boundary condition for solar wind heliospheric models is derived from coronal electron density tomography maps. The tomography/HUXt model results give general good agreement with in situ data provided by the OMNI satellite network, with a significant increase in the accuracy of solar wind velocity prediction compared to a HUXt model run with a traditional MAS derived inner boundary condition for the time periods used in this study. The time periods investigated in this study were intentionally chosen to demonstrate the model at different stages of the solar activity cycle. Given further development, a future study will explore a larger dataset.

We can identify several aspects of the inner boundary condition and modelling that lead to inaccuracies and are possible to address with some further work. They are:

• The use of a static tomography reconstruction to represent a dynamic corona. We have recently developed a framework for providing time-dependent tomographical densities that will help address this, although further development is needed (Morgan, 2021).

• The overly smooth reconstruction given by tomography compared to the true density. The coronal streamer belt likely consists of very narrow high-density structures (Morgan & Habbal, 2007a, 2007b, 2010) that are highly variable on small temporal (e.g., Alzate et al., 2021) and spatial scales (Thernisien & Howard, 2006; Poirier et al., 2020). Again, developments in tomographical methods and observations can lead to reconstructions that bring us closer to these scales.

• The use of a single, fixed acceleration profile for both fast and slow wind. We are developing our own upwind model that includes both velocities and densities, and includes the acceleration profile as a search parameter. This method uses iteration to improve the fit to in situ measurements. Early results are promising and may provide a constraint on acceleration. Another approach is to use a set of tomographical maps over a range of distances (4–10 R_⊙) to constrain the early acceleration profile of the slow wind, similar to that shown by Morgan & Cook (2020).

• The simplistic inverse relationship of equation (1). We are investigating the replacement of the simple inverse relationship of equation (1) with a Sigmoid or exponential function, which gives a simple relationship, based on a small number of parameters, to model the transition from slow to fast wind as a function of the density.

• Use of a standard, possibly incorrect, coronal rotation rate. Long time series of tomographical maps can give improved estimates of the variable coronal rotation rate (Morgan, 2011b; Edwards et al. 2022). Near the equator, the coronal rotation rate may vary by a degree per day or more from the Carrington rate, leading to systematic longitudinal errors in solar wind model results.

• The omission of coronal mass ejections (CME) that may be present in the in situ data. This can be addressed to an extent using the current approach for operational forecasting: simple parameters describing a cone model of a CME can be input to HUXt, and the CME carried with the solar wind, giving an estimate of the time of arrival at Earth.

A persistence-based approach was made in order to test the operational feasibility of this model. The persistence results showed a large difference in accuracy between periods near solar minimum and maximum. The rate of evolution of the physical corona during solar maximum is such that a persistence approach will break down over the timescales of a Carrington rotation. However, during solar minimum, results were more acceptable yet less accurate than a non-persistence approach. This highlights the need for a time-dependent inner boundary condition, and a more advanced relationship between density and velocity than that given by equation (1).

The efficiency of the HUXt model allows an ensemble framework, which can quantify the uncertainties of the predicted solar wind velocities. The ensembles were based on sampling an appropriate range of latitudes and velocities at the model inner boundary. The ensemble results confirm that both these uncertainties have a significant effect on the model output velocities on Earth. Latitudinally narrow, longitudinally-aligned streamer belt structures can lead to high uncertainties in the output based on small latitudinal uncertainties at the inner boundary. The choice of velocity range at the inner boundary also has a large impact on the model results.

In the context of operational space weather forecasting, this paper shows the use of a tomography-based inner boundary condition to drive the HUXt model as part of the CORTOM module of the SWEEP project for the UK Met Office is an approach that can offer certain improvements on current systems. The use of multiple models over extended periods will enable an extensive analysis of their relative performance, and the improvements described in this paper will be implemented within SWEEP over the coming years.

Acknowledgments

We acknowledge (1) STFC grants ST/S000518/1 and ST/V00235X/1, (2) Leverhulme grant RPG-2019-361, (3) STFC PhD studentship ST/S505225/1, and (4) the excellent facilities and support of SuperComputing Wales, (5) NASA/GSFC’s Space Physics Data Facility’s OMNIWeb service and OMNI data, (6) Predictive Science Inc. for the MHDWeb service and MAS model outputs, and (7) University of Reading (UK) for the HUXt model. HUXt is available for download at https://github.com/University-of-ReadingSpace-Science/HUXt. The STEREO/SECCHI project is an international consortium of the Naval Research Laboratory (USA), Lockheed Martin Solar and Astrophysics Lab (USA), NASA Goddard Space Flight Center (USA), Rutherford Appleton Laboratory (UK), University of Birmingham (UK), Max-Planck-Institut für Sonnen-systemforschung (Germany), Centre Spatial de Liege (Belgium), Institut Optique Théorique et Appliqúee (France), and Institut d’Astrophysique Spatiale (France). The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

Fig. 1 (a) The observed coronal polarised brightness at a distance of 8 R⊙, plotted as a function of position angle (measured counter-clockwise from solar north) and time. This is used as input to the tomography method. (b) The reconstructed polarised brightness, or synthetic observations gained from the tomographical density distribution. (c) The electron density at a distance of 8 R⊙ mapped in Carrington longitude and latitude.

In the text

	Fig. 2 Coronal electron density (blue, right axis) taken from the tomographical map of Figure 1c at Earth’s latitude, and solar wind velocity at 8 R⊙ (black, left axis) gained from the density using the simple inverse relationship of equation (1). For this example, Vmax = 480 km s−1 and Vmin = 220 km s−1.
In the text

	Fig. 3 DTW optimal path (red) between in situ data and two Tomography/HUXt model runs for CR2210 with (a) overestimated Vmax and Vmin terms of 600 and 250 km s−1 respectively, and (b) underestimated Vmax and Vmin terms of 380 and 110 km s−1, respectively.
In the text

Fig. 4 (a) Contour plot of DTW path lengths with a varying magnitude of the Vmax and Vmin terms for the tomography-derived inner boundary condition. The minimum DTW path distance is marked by a white cross, and this point gives the optimal values of Vmax and Vmin of 480 and 220 km s−1, respectively. (b) DTW path lengths as a function of varying scale parameters for the MAS-derived inner boundary condition with the white cross showing the optimum VMmin and VMmax parameters of 280 and 580 km s−1 respectively.

In the text

Fig. 5 Further analysis of the optimal DTW path. (a) Comparison of modelled solar wind velocities at Earth (grey) with in situ data (black) and DTW optimised path (red). (b) Alignment of the modelled and in situ data along the optimised DTW path with respect to time. (c) Difference in time between observation and modelled data of aligned points along optimal DTW path (ΔT(insitu−model)). (d) Velocity difference between observation and modelled data of aligned points along optimal DTW path (ΔV(insitu−model)).

In the text

	Fig. 6 Comparison of optimised tomography/HUXt (blue) and optimised MAS/HUXt (orange) model predictions of the solar wind velocity at Earth with insitu data (black) for CR2210.
In the text

	Fig. 7 Electron density as estimated using coronal tomography at a distance of 8 R⊙, for Carrington rotations (a) 2150.5 and (b) 2202.5. The red boxes labelled A and B in (a), and the white box in (b) are relevant for the ensemble results of Section 4.
In the text

Fig. 8 (a) Comparison of model and OMNI solar wind velocities for CR 2202 for the optimal fit in terms of minimum DTW path distance (marked by a white cross in Fig. 8b), with Vmin and Vmax set to 230 and 440 km s−1 respectively. (b) Contour plot showing the DTW path distance as a function of Vmin and Vmax for CR 2202. (c) Same as (a), but for CR 2150, with Vmin and Vmax set to 240 and 400 km s−1, respectively (represented by white cross in Fig. 8d). (d) Same as (b), but for CR2150.

In the text

	Fig. 9 Comparison of predicted solar wind conditions of (a) CR2203 (Vmax and Vmin of 440 and 230 km s−1 respectively) and (b) CR2151 (Vmax and Vmin of 400 and 240 km s−1 respectively). The Vmax and Vmin values are that of the optimised velocity terms of the Carrington rotation prior (obtained in Sect. 3.3).
In the text

	Fig. 10 Comparison of in situ data with the results of multiple tomography/HUXt model runs for (a) CR2150 and (b) CR2202, with the latitude of the inner boundary condition varied in one-degree increments between ±3° from the latitude of Earth.
In the text

	Fig. 11 The rapidly changing density of the corona between (a) CR2149, (b) CR2150, and (c) CR2151.
In the text

Fig. 12 Contour plots of the DTW path distance as a function of minimum and maximum inner boundary velocity for Carrington rotations (a) 2150 and (b) 2202. The white contour bounds the area that yields a DTW path distance under the threshold chosen (SSF ≤ 0.95), and the white pixels in the lower right corner represent invalid model parameters where Vmin > Vmax. Note that the Vmax range is altered to 300–600 km s−1 in this plot.

In the text

	Fig. 13 Comparison of in situ data with the results of multiple tomography/HUXt model runs for (a) CR2150 and (b) CR2202, with a range of Vmax and Vmin combinations.
In the text

	Fig. 14 Results of the ensemble for (a) CR2150 and (b) CR2202 with the mean predicted solar wind velocity at each time step shown as white, the light grey area representing all predicted solar wind velocities, and the darker grey representing one standard deviation from the mean. The non-ensemble model results with the optimal velocity values at the Carrington latitude of Earth is the blue line.
In the text