A comparison study of extrapolation models and empirical relations in forecasting solar wind

Coronal mass ejections (CMEs) and high speed solar streams serve as perturbations to the background solar wind that have major implications in space weather dynamics. Therefore, a robust framework for accurate predictions of the background wind properties is a fundamental step towards the development of any space weather prediction toolbox. In this pilot study, we focus on the implementation and comparison of various models that are critical for a steady state, solar wind forecasting framework. Specifically, we perform case studies on Carrington rotations 2053, 2082 and 2104, and compare the performance of magnetic field extrapolation models in conjunction with velocity empirical formulations to predict solar wind properties at Lagrangian point L1. Two different models to extrapolate the solar wind from the coronal domain to the inner-heliospheric domain are presented, namely, (a) Kinematics based (Heliospheric Upwind eXtrapolation [HUX]) model and (b) Physics based model. The physics based model solves a set of conservative equations of hydrodynamics using the PLUTO code and can additionally predict the thermal properties of solar wind. The assessment in predicting solar wind parameters of the different models is quantified through statistical measures. We further extend this developed framework to also assess the polarity of inter-planetary magnetic field at L1. Our best models for the case of CR2053 gives a very high correlation coefficient ($\sim$ 0.73-0.81) and has an root mean square error of ($\sim$ 75-90 kms$^{-1}$). Additionally, the physics based model has a standard deviation comparable with that obtained from the hourly OMNI solar wind data and also produces a considerable match with observed solar wind proton temperatures measured at L1 from the same database.


INTRODUCTION
Space weather refers to the dynamic conditions on the Sun and in the intervening Sun-Earth medium that can severely influence the functioning of space-borne and ground based technical instruments thereby affecting human life. Predicting the impact of space weather thereby becomes an essential task. In particular, explosive events on the Sun that include solar flares, coronal mass ejections (CMEs) and solar energetic The various forecasting models over time have shown considerable progress in our understanding of modelling the macro-physical solar wind properties. One of the key ingredients that is critical for space weather modelling is the interplay of micro-physical turbulent and particle acceleration processes with macro-physical dynamics. Several attempts have been made to include effects of energy contained in subgrid turbulence and multi-fluid aspects in modelling of solar wind (CORHEL [8], AWSoM [37], CRONOS [47], [41], [42]). In general, handling of large separation of scales required to consistently model such an interplay of length scales is a challenging task. Such a multi-scale nature of problem is also experienced in astrophysical plasma modelling and several astrophysical codes are working in direction of developing a hybrid framework that bridges this wide gap. In particular, recent development of hybrid particle module for PLUTO code [15,43] to model particle acceleration at shocks and its subsequent non-thermal emission. Further, PLUTO code supports adaptive mesh refinement and various non-ideal MHD processes including magnetic resistivity [16] and Hall-MHD. The code also has support for anisotropic thermal conduction [44] and optical thin cooling. In last 5 years, problems pertaining to solar and magnetospheric physics have also been tackled with PLUTO code [24,32,7,27]. The goal of this work is to develop a space weather modelling framework in conjunction with PLUTO code aiming to utilize the additional functionalities of the code for modelling the micro-physical aspects of Sun-Earth environment.
In this first study, our focus is to compare the various coronal and inner heliospheric models for solar wind forecasting and the assessment of their predictive performance. The paper is arranged in the following manner -the details of the methods implemented for the forecasting framework are described in Sec 2. The results from various models adopted for forecasting of solar wind velocity for a couple of case studies along with statistical assessment of their accuracy are elaborated in Sec 3. Finally, the Sec 4 discusses the various features of the forecasting framework along with limitations and future outlook.

METHODOLOGY
The region between the solar photosphere and the Lagrangian point L1 is divided into two zones. The inner coronal zone extends from the photosphere up-to 5 R , followed by inner-heliosphere zone that extends from 5 R up-to L1 point. Data driven prediction of solar wind parameters at L1 point requires the following steps -• To calculate the magnetic fields in the coronal region through various extrapolation methods of the input observed photospheric magnetic field.
• Applying the velocity empirical relations based on the field line properties obtained from extrapolation at the outer boundary of the coronal region.
• Extending the velocity estimates from the outer boundary of the coronal region up-to Lagrangian point L1 for comparison with observations.
Detailed procedure followed for each of the above steps is described in this section.

Magnetic Field Extrapolation
Forecasting of solar wind across the domain requires accurate magnetic field solutions extrapolated from the solar surface to the outer boundaries of the coronal domain. The magnetic field extrapolations are carried out up to a distance of 5R . The inner boundary conditions for the magnetic fields at the solar surface are specified by the input magnetograms(synoptic maps) taken from the Global Oscillations Network Group(GONG) (https://gong.nso.edu/data/magmap/crmap.html). Once the inner boundary conditions are specified, the extrapolation of the magnetic fields are then carried out in a twofold manner-Just above the photospheric region, the magnetic energy density is greater than the plasma energy density and magnetic effects dominate. One can assume this region to be current free and thus, use the potential formulation for the magnetic fields [34]. This current free assumption is valid inside a sphere of radius of about 2.5R , the outer boundary of which is known as the source surface. The magnetic fields inside the source surface can be solved using the PFSS (Potential Field Source Surface) model. The outer boundary condition for this model dictates that the magnetic fields at the source surface is approximately radial [34]. For the PFSS solution, we used the python module pfsspy, which is an open source, finite difference PFSS solver. Using the observed magnetogram data as an input, potential field equations are solved in radial direction on a logarithmic grid, whereas for latitude and longitude the grid is regularly spaced in terms of cos θ (-1.0 to 1.0) and φ (0 to 2π) [38,49]. Field line tracing using the field solutions is done via Runge Kutta 4 th order method. In the region outside the source surface, the magnetic fields are extrapolated using the Schatten Current Sheet (SCS) model [33]. The SCS model extends the magnetic fields from the source surface to a distance of 5R , i.e, the outer boundary of the coronal zone. The input to the SCS model is the re-oriented output from pfsspy, i.e, if at the source surface, B r 0, the field remains unchanged, but if B r 0, B r , B θ , B φ are replaced by −B r , −B θ , −B φ . Using the same resolution in cos θ and φ plane as that of the output from pfsspy, all magnetic field components beyond the source surface can be expressed in the form of a Legendre polynomial expansion [33,25]. Following the formulations presented in the appendix of [25], we construct the matrices of coefficient of spherical harmonics g m n and h m n , where n and m are degree and order of the associated Legendre polynomial P m n (cos θ). A rather modest value of n = 15 is used to approximate the magnetic field values using the SCS model. The accuracy of SCS approximation improved only marginally with doubling the choice of n, but the computational time increased significantly. Thus, the order of Legendre polynomial used in the present work was carried out by optimising both accuracy and computational time. While tracing the field lines, one should note that the boundary conditions of PFSS and SCS are not compatible with each other. A direct combination of the PFSS solution with the SCS model at 2.5 R results in kinks and discontinuities at the model boundary which is a consequence of the different boundary requirements for the two models. To avoid this non-physical discrepancy, the input to the SCS model is given by the field values at 2.3 R rather than 2.5R . The field values in the thin radial slice between 2.3 R and 2.5 R is then overwritten by the values obtained by the SCS model. This leads to a smooth transition and a more desirable coupling between the PFSS and the SCS models [31]. The SCS model then extrapolates the fields upto 5 R . The PFSS + SCS solution together gives us a good approximation of the magnetic field structure up-to a distance of 5 R .

Velocity Empirical Relations
Based on the magnetic field structure [25] obtained by from the field extrapolation methods, we generate a velocity map for the solar wind using some empirical velocity mapping models. We employ two different empirical models for the velocity calculations described in the sections below.

Wang-Sheeley Model(WS)
The Wang-Sheeley model depends on a parameter f s , called the expansion factor of the coronal flux tubes to calculate Solar Wind velocities. The expansion factor(f s ) is given by where B P (R ss ) denotes the radial field strength at a sub-earth point P on the source surface and B P (R ) is the foot-point of the flux tube traversing P on the photosphere [45]. The expansion factor measures the amount of change in cross section of a flux tube between the photosphere and the source surface compared to a purely radial expansion. It is observed that there is a correlation between the expansion factor and the solar wind velocities. Based on this, an empirical formula for calculating the solar wind velocities can be where, v slow is the lowest expected speed as f s → ∞ and v fast is the fastest solar wind expected as f s → 1.

Wang-Sheeley-Arge Model(WSA)
The Wang-Sheeley-Arge [5] is a model that also incorporates the effect of minimum angular distance of foot point of the field line from coronal hole, along with the expansion factor. It is believed that coronal holes produce fast streams of solar wind. So the position of the field line's foot point in coronal hole plays a very important role. The empirical relation from the WSA model [30] used in the present work is - The parameters v slow and v fast corresponds to the velocity of fastest and slowest solar wind stream. θ b is the minimum angular distance for the foot point of the field line from a coronal hole boundary at the Solar surface. In the present work, we have used v slow = 250 kms −1 , v fast = 900 kms −1 , α = 1.5/9, w = 0.03, and δ = 1.5.

Extrapolation into the Heliospheric domain
The WS and WSA model gives us the solar wind maps at the outer boundary of the coronal domain i.e., at 2.5 R or 5 R based on the choice of magnetic field extrapolation method. For comparison with the observed velocities at L1 point, it is required to extrapolate these velocities into the inner heliosphere zone. This requires coupling of the coronal velocity models with heliosperic velocity extrapolation models. We have employed two such extrapolation methods, (a) Heliospheric Upwind eXtrapolation(HUX) [21] and (b) Physics based Modelling using PLUTO code. We describe these methods in the following sections.

Heliospheric Upwind eXtrapolation(HUX)
The HUX model assumes the solar-wind flow at the outer boundary of the coronal domain to be time-stationary. The extrapolation of the solar wind velocities in an r-φ grid can be then be kinetically approximated using where, ∆r = 1R and ∆φ = 1 • represent the grid spacings in r and φ directions respectively. Ω is the angular velocity of the Sun calculated assuming a rotation time period of 27.3 days. The HUX is essentially a 1D extrapolation that neglect the effects of magnetic fields, pressure gradients and gravity. The advantage being that such an extrapolation method is computationally inexpensive when compared to state-of-the-art 3D MHD models.

Physics based Modelling using PLUTO code
With an aim to incorporate the effects of some of the physical aspects in solar wind extrapolation, we describe a physics based modelling approach that involves solving a set of conservative equations using the Godunov scheme based Eulerian grid code, PLUTO [14].
For this pilot study, we assume the solar wind to be hydrodynamic and solve the following set of compressible equations in 2D polar co-ordinates (r, φ) : where, ρ is the density of the fluid, P being the isotropic thermal pressure, v is the fluid velocity, I is an identity matrix and the total energy E = 1 2 ρv 2 + ρ is sum of the kinetic energy and the internal energy. The acceleration due to gravity g = −GM /r 2 is included as a source term in the conservative momentum equation. A poly-tropic equation of state is adopted with the value of adiabatic index γ = 1.5 [17,23].
The above equations are solved in a non-inertial frame where the inner radial boundary rotates with the rate equal to the solar rotation rate. Further to simplify we neglect the Coriolis and centrifugal terms as their contribution is rather small in determining the steady flow structure of the solar wind [23]. The computational grid in polar co-ordinates ranges in the radial direction from 5R to 435R with a resolution of 512 grid cells. Same number of grid cells are used to resolve the azimuthal (φ) direction.
Initially, the computational domain is filled with a static gas with number density of 1 cm −3 and thermal sound speed of 180 km s −1 . The choice of initial conditions does not affect final steady state wind solution as the material injected from the inner boundary wash away these initial values to obtain a new steady state.
Radial velocity is prescribed at the inner boundary from the WSA mapping. The inner boundary is also made to rotate with respect to the computational domain with a rate equal to the solar rotational time period which is 27.3 days. As a general rule, boundary conditions can be specified for only those characteristics that are outgoing (away from the boundary). Since we have a supersonic inflow boundary condition, all the characteristics are pointing away from the boundary and therefore, along with the prescription of solar wind velocity, one would need to prescribe the density and pressure at the inner radial boundary as well. Following [23], we prescribe the number density (n) and pressure at the inner radial boundary in following manner : where P 0 is set to be a constant value of 2.75 nPa and n 0 = 100 cm −3 resulting in a proton temperature of 2 MK. The scaling velocity v 0 is set to be 675 km s −1 . The outer radial boundary is set to have free flowing outflow conditions. For each Carrington rotation, we carry out the simulations using velocities obtained from sub-earth points with the magnetic field parameters obtained at 5 R and employing WSA empirical relation for velocity. The results are presented in Sec 3.4.

Model Definitions
Using the combination of various magnetic field extrapolation, velocity empirical relations and method of extrapolating velocity field within the inner heliosphere, we have defined 6 different models. The combinations used for each of these models are shown in Figure 1. Models 1 and 2 involves magnetic field extrapolation using PFSS (without SCS) and velocity extrapolation in inner heliosphere using the HUX model. The difference between them is the choice of velocity empirical relation. In the other set of models 3 and 4a, the magnetic fields are extrapolated using PFSS+SCS. The empirical velocity relation employed for Model 3 is WS and that for Model 4a is WSA. Velocity field in both these models are extrapolated from 5 R to L1 using HUX. Model 4b is a variant obtained by considering an ensemble values using the same combinations as that of Model 4a (see section 3.2). Model 4c uses the PFSS+SCS field extrapolation model and WSA velocity empirical relation similar to that used by Model 4a, however, the extrapolation of velocity field is done using physics based hydrodynamic simulations.

Statistical Measures of Forecast Performance
The performance of a forecast can be determined by comparing the forecast outcome of continuous variables (e.g., velocity) to the observed values. We calculated several scalar measures of forecast accuracy which has previously been used to determine forecast performances by [26,48]. Given a set of modelled values m n and a set of corresponding observed values o n , the Mean Absolute Error (MAE) is given as the arithmetic mean of the absolute values of the differences between the model output and the observed values at each observed data point.
The Mean absolute error can be considered as a measure of the overall error in forecast of a model. Another such measure is the Mean Absolute Percentage error(MAPE) which is given by The Root Mean Square Error or RMSE is also used sometimes as a performance statistic for a model and is given by Another important measure in determining the forecast performance is the Pearson Correlation Coefficient (CC) which is a parameter used to estimate the correlation between the observed values and the model values. In addition to this, another measure is the standard deviation of the time series of each of the individual model outputs and its comparison with the estimate from observed data. These measures are estimated for all the models considered and discussed in section 3.5 for the three Carrington rotations CR2053, CR2082 and CR2104 considered in our study.

RESULTS
In this section, we present our results using the methodology described above for the three case studies spanning the declining phase of cycle 23, near the solar minimum and the rising phase of solar cycle 24. In particular, we consider CR 2053 (from 2007/02/04 to 2007/03/04), CR 2082 (from 2009/04/05 to 2009/05/03) and CR 2104 (from 2010/11/26 to 2010/12/24).We also demonstrate the assessment of performance of the different forecasting models considered for these cases.

Case Studies
CR2053 represents a relatively quiet phase of the sun during the decline of the solar cycle number 23. Six active regions were identified regions during the period of CR2053 [9] that can also be seen in our  Fig. 2). All the active region lie close to the L1 footprint and thus are expected to have pertinent effects on the solar wind velocities. As a first step, we calculated the PFSS solution for this Carrington rotation. The GONG magnetogram data which is given as input to PFSS and the output solution obtained using pfsspy are shown in Figure. 2. A 3D map of field lines joining the photosphere to the PFSS source surface at 2.5 R is shown in panel (C) of the same figure. The angle θ on the Y-axis of panels (A) and (B) is related to Carrington latitude (θ cr ) as θ = −θ cr + 90 • . It is evident from the distribution of field lines that few open field lines with both positive (red) and negative (blue) polarities have their foot points located on the photosphere that are close to solar equator. On extending these field lines from source surface up-to 5 R incorporating the SCS model, we obtain a distribution of magnetic field lines at the boundary of coronal domain. Sub-earth field lines are selected from such a distribution and the parameters like the expansion factor (f s ) and their distance to a nearby coronal hole boundary (θ b ) are estimated. These parameters obtained at 5 R are shown in panels (A) and (B) of figure 4 as a function of Carrington longitude (φ). The solar wind velocity obtained from using these parameters in the empirical velocity relations of WS (Model 3) and WSA (Model 4a) are shown in panel (C) in comparison with the observed data. As the models 1 and 2 only rely on the magnetic fields interpolated using PFSS solution, the field parameters are estimated at 2.5 R and used in the empirical relations to estimate the velocity at L1 after HUX extrapolation.
We observe that the WS model gives an inaccurate representation of the solar wind velocities and the contrast between the slow wind and the fast streams are not satisfactorily reproduced. For example, during We have also carried out the same analysis using the same set of model parameters (see Eq. 3)for the case of CR 2104 that represents the rising phase of the solar cycle 24. We find similar trend even for this case as solar wind velocity estimates from Model 2 and Model 4a that involves WSA empirical relation presents a better match with observations as compared those obtained from Model 1 and Model 3.  These values are estimated for field lines extrapolated using PFSS + SCS up-to the reference sphere of radius 5 R for CR 2082. Panel (C) shows a comparison of the observed solar wind velocities (taken from OMNI database) and the relative differences between the various outputs of the numerical models at L1 (≈ 215R ).

Ensemble Forecasting
Numerical extrapolation models like HUX are computationally less demanding than HD or 3D MHD models. Such numerical models are thus used to study a large set of initial conditions by a method known as ensemble forecasting. Ensemble forecasting has been widely used to constrain terrestrial weather and is an important tool to determine model performance and also helps to set uncertainty bounds to the model output. A solar wind forecast model is considered to be a very uncertain one if the ensemble members produce drastically different results [25]. To study the variations introduced in our model output due to an uncertainty in determining source that contribute to solar wind measured at L1, we created an ensemble of latitudes centered around our expected sub-earth latitude at a distance of ±1 • , ±2.5 • , ±5 • , ±10 • , and ±15 • . The change in sub-earth latitude changes the field lines under consideration which in turn effects the parameters associated with the field lines (e.g, θ b , f s ) that are used as input to the WSA model. This reflects as a change in the final output velocities at L1. Our ensemble of sub-earth points provides us with a set of 11 different velocities at each Carrington longitude. We consider the ensemble median to be the preferred average as an ensemble of sub-earth points are known to produce highly skewed solutions and thus the ensemble mean may not be an accurate representation of the results [25]. One should note that the goal of ensemble forecasting in this case is to provide a measure of the uncertainty in the obtained model velocities rather than the ensemble results leading towards a better forecast [20,10,25]. Figure: 7 shows the forecast obtained from the ensemble of sub-earth latitudes by using Model 4a for all three cases CR 2053, CR 2082 and CR 2104.The dashed line in the figures represent the ensemble median for the respective cases. We have also plotted the velocity obtained by considering zero uncertainties in the the sub-earth latitudes and the same is plotted as a black solid line. The dark and light shades around the median represent the 1σ and 2σ variations respectively around the ensemble average. A wider spread in the 2σ values indicates a larger uncertainly in the velocities forecast by the model. In the case of CR2053 and CR2082, a large portion of the observed velocity profile lies within the 2σ error bounds indicating an excellent forecast performance. For the case of CR 2104, the high speed stream is not predicted well with any of our models and some portion of the observed values lie outside of the 2σ error bounds. This is mainly due to the underestimation of standard deviation from the models as compared to observations. However, the structure is well correlated as evident from high CC (section 3.5).

Assessing Interplanetary Magnetic field Polarity
In addition to predicting the solar wind velocity, we also assess the polarity of the magnetic field lines at L1 assuming a Parker spiral magnetic field configration between Sun and Earth (see [11]) We evaluate the corresponding longitude at L1 for each CR longitude at 5 R using the standard streamline equation. The solar wind velocity is assumed to be constant and corresponds to the value obtained using WSA formulation. The polarity of radial magnetic field obtained at 5 R from PFSS + SCS extrapolation is then coverted to GSE/GSM co-ordinate frame to compare with the observed polarity of B x component of interplanetary magnetic field at L1 from OMNI database for all three Carrington rotations. The convention of polarity followed at L1 in this case is -1 (outward) and +1 (inward). Figure: 6 shows the observed polarity (daily averaged) of inter-planetary magnetic field B x in GSE/GSM co-ordinates and its comparison with that obtained for all three cases CR 2053, CR 2082 and CR 2104 using the field extrapolations and velocity from Model 4a. The inter-planetary magnetic field polarity obtained for all three cases show a good agreement with observed daily averaged values. In particular for the case of CR 2104 shown in panel (C) of figure 6, significant agreement exists between model prediction and observation.

Forecasting from Physics based Modelling
Two dimensional hydrodynamic simulation runs with initial conditions described in section 2.3.2 provide a time dependent solar wind forecasting framework. Plasma is injected in the computational domain from One can distinctively observe a spiral pattern whereby high velocity streams tend to have higher proton temperature and lower density values than its surrounding low or moderate velocity streams. The velocity measured at L1 from simulations during the considered Carrington rotation period is shown as blue solid line in panel D. In comparison, hourly averaged observed values are shown in black dashed lines, while velocity estimates from models 4a and 4b are shown in magneta solid and green dashed lines respectively for comparison. The velocity estimate from the simulation runs have more variations as compared to its counterparts from model 4a and 4b. This comparison of models for velocity prediction demonstrates that addition of pressure gradient term and incorporating the energy equation captures the interaction of streams leading to intermittent variable patterns. Such patterns do not appear for velocities obtained from model 4a and 4b as the kinematic extrapolation disregards physical effects due to presence of velocity shear and pressure gradient terms. Statistical analysis of velocity pattern at L1 obtained from HD simulations indicate a good match with observed values (see Table: 1). Similar variable velocity pattern is also obtained after using the empirical velocity values from CR2082 in the inner boundary (figure not shown). Statistical comparison with other models and observed values is presented in Table:

Quantifying Predictive Performance using Statistical Analysis
Statistical measures (see section 2.5) quantifying forecast performance for CR 2053, CR 2082 and CR 2104 are presented in Table: 1, Table: 2 and Table: 3 respectively. The standard deviation of the observed time series for CR2053 is 119.4 kms −1 . Model 4c was able to reproduce this standard deviation quite satisfactorily (123.6 kms −1 ). All models for CR 2053 involving the WS model failed to produce reasonable agreement with the observations, which can be inferred from the extremely low value of the correlation coefficient (∼ 0). We believe that this is due to the proximity of the coronal holes to the equator, whereby open field lines emerging from these holes drive high speed streams. The WS empirical model disregards the size of coronal holes in its empirical relation and fails to capture this effect. The presence of a significant number of coronal holes on the magnetogram that are close to the sub-earth latitudes thus necessitates the use of the WSA model for velocity mapping on the inner boundary. CR2082 on the other hand has a significantly lesser number of coronal holes [22]. And thus even the models involving the WS velocity mapping show acceptable CC values (0.42-0.54). The highest CC for CR 2053 is produced by model 4b (CC= 0.81) which represents the ensemble median of the sub-earth latitude variations, however the same cannot be said for CR 2082 as all the other models produce a greater CC and lower error values than that of the ensemble median. This reinforces the statement that the ensemble average does not always produce a more accurate representation of the velocity profile.
A graphical representation and analysis of various forecasting models implemented in this work are represented in a Taylor diagram [39]. The Taylor plots for CR 2053, CR 2082 and CR 2104 are given in Figure: 9 The radial distance from the origin represents the standard deviation, the azimuthal coordinate represent the correlation coefficient (CC). In addition to this, one more set of arcs are drawn, centered at the reference (observed) standard deviation on the X axis that represent the RMS error for the models. In general a solar wind model is considered to be a good model if it has a high correlation coefficient(CC), low RMS error, and having a similar standard deviation to the observed time series of velocities [26] and in a Taylor plot it would lie close to the smallest circle on the X axis. For the case of CR 2053, Model 1 and Model 3 lie very close to the vertical axis indicating very low CC value and thus are considered to have poor forecast performance. Model 4b has the highest CC value and comparatively low errors, but it lies relatively far from the observed standard deviation curve. The physics based models, Model 4c is seen to have standard deviation profiles well in agreement with the observed values along with a good CC. We thus assert that Model 4b and 4c are the best representations of the velocity profiles for CR 2053. For the case of CR 2082, Model 4c and Model 2 has a standard deviation closest to that of the observed value of ∼ 60 kms −1 and Model 2 has the lowest error values (in terms of RMSE, MAE and MAPE). The CC of model 3 is one of the highest in the models considered followed closely by Model 2. Model 3 however has larger values of all the errors considered showing that even though its correlated well with the observed values, the model is prone to errors which would result in inaccurate output solar wind velocity profiles.

DISCUSSIONS AND SUMMARY
In this pilot study, we have developed a python module towards constructing a robust framework for accurate predictions of a steady state background solar wind model using various empirical and extrapolation formulations. This framework is also integrated with physics based modelling using PLUTO code. The entire workflow is a combination of (a) extrapolating the magnetic fields to the outer coronal domain using magnetic models such as PFSS and SCS, (b) mapping the velocities in the outer coronal domain using models such as WS and WSA and (c) extrapolating the velocities to L1 using extrapolation techniques such as HUX as well as hydrodynamic propagation of the velocities using PLUTO. We have studied various combinations of the coronal magnetic models as well as velocity extrapolation models in view of three different Carrington rotations, CR 2053, CR 2082 and CR 2104. We were successfully able to generate velocity profiles well in agreement with the observed values.
Even though the the models involving the WS velocity mapping performed relatively well for the case of CR 2082 when compared to the equivalent models of the other CRs in this study, it performed very poorly when applied to CR 2053 and CR 2104. The WSA model on the other hand has shown consistently superior performance for all the CRs with lesser error measures in these cases than the WS model. The WS model also failed to capture the contrast between the slow and fast solar wind streams and thus produced velocity profiles with significantly lesser standard deviation than the observed profile. We infer from this that the WSA model is superior to the WS model for the cases and parameters considered in this work.
In general, the models with a combination of PFSS + SCS for magnetic field extrapolation paired with the WSA model for velocity mapping (models 4a, 4b and 4c) had the best performance in all the cases, CR 2053, CR 2082 as well as CR 2104. This can be seen as near identical standard deviations of the output velocity pattern in Model 4a and Model 4c paired with a very high correlation coefficient (∼ 0.73-0.81) and relatively low errors for the case of CR 2053 (Table: 1) and CR 2104 (Table: 3). For CR 2082, Model 4a and Model 4c had standard deviations which are close to that of the observed profile. We note that Model 3 performed quite well in case of CR 2082 with a high CC (∼ 0.53) but with higher errors than the other models. However, it showed a lack of consistency in performance with very poor results in case of CR 2053. We observe that Model 2 also performed well for all three CRs which can be seen in its high CC values and low errors.
Model 4c that employs PFSS+SCS and WSA velocity empirical relation with physics based extrapolation, gives us an additional advantage of being able to determine the thermal properties of the solar wind which also showed good correlation with the observed values (see Figure: 10). This model also naturally takes into account the acceleration of the solar wind during its propagation. The standard deviation of the Model 4c velocity profile also had good agreement with the observations with high CC values (0.73 for CR 2053, 0.46 for CR 2082 and 0.85 for CR 2104). The Model 4c thus has its own pros and cons. The pros being the ability to determine additional solar wind properties than what HUX or HUXt can offer, and the cons being that it takes significantly more computational resources.
Ensemble forecasting helps provide a rather clear quantification of the uncertainty associated with the predictions. We have chosen a total of 11 realisations covering a latitude range of ±15 • about the expected sub-earth latitude. We can assert that the variation helped us accurately capture the uncertainty in the predictions as a bulk of the observed velocity curve was entirely within the 2σ bounds of the model output for CR 2053 as well as CR 2082 (Figure: 7). Underestimation in the solar wind velocity is seen in case of CR 2104 particularly for the high speeds streams. We also note that the ensemble median did not necessarily provide a better estimate of the model velocity and in general, had sub-par performance when compared to some of the other models (Model 4a and Model 4c). The radial magnetic field polarity obtained at 5 R has also been extrapolated to L1 assuming a Parker spiral field profile for Model 4a. A good agreement between the observed and predicted field polarity is seen in all three cases with CR 2082 showing misses only for three data points obtained from observations. Even though ensemble forecasting can provide a rather clear statistical uncertainty prediction of the model, it has been seen that various other uncertainties can creep in that are solely dependent on the choice of the input magnetograms. Riley et.al [29] highlighted that magnetograms from various observatories show significant differences in magnetic field measurements. In-spite of the fact that one can quantify the conversion factors between various independent datasets, there is no particular observatory that can produce a ground truth dataset for the input magnetograms. Additionally, one can improve the WSA model performance using time-dependent flux transport model based magnetograms [e.g., 35,3,2]. Quantifying the effects of input magnetograms on model performance is beyond the scope of this paper and shall be addressed in a future study. The statistical parameters represented in the Taylor diagrams indicate the accuracy of our current forecasting framework. Typically, for an accurate forecasting model, low RMSE, high CC and similar standard deviation with observed values are expected. The forecasting presented in this paper uses the same set of parameters (except θ b and f s ) for empirical velocity estimate of all the Carrington rotations. We find that these parameters give a more accurate forecast of the bulk solar wind speeds at L1 for the case of CR 2053. Even though the same parameters result in good forecast for CR 2082 and CR 2104 as well, we believe that the performance may be improved with a changes in the free parameters and radial distance of reference sphere [48].
We would like to point out that the free parameters used for the WS and WSA empirical relations cannot be universally applied for all CRs and must be individually tuned for a particular case scenario. This is evident from our use of the same set of free parameters for all the CRs which results in significantly different performances. Additionally, the statistical estimates have indicated that in case of CR 2082, the model with PFSS+WSA have performed slightly better than its counterpart including SCS. In this particular case, this could perhaps be related to over-estimation of magnetic field lines contributing to solar wind speeds at L1 by SCS. However, through the ensemble modelling we have demonstrated the estimation from Model 4a is within the 2σ bounds. Improving the accuracy and performance of SCS extrapolation will be considered in our further studies. The caveat of hydrodynamic models (Model 4c) is that it solves the equations on a plane and does not contain the magnetic information of the sun-earth system on which state-of-the-art 3D MHD model are based. However, even without the magnetic field information and having a 2D geometry, the HD models capture many important aspects of the solar wind. Model 4c facilitate better physical insight in the behaviour of solar wind than the HUX models (4a and 4b). Model 4c (HD) can thus be seen as a halfway point between HUX and full MHD models.
In summary, we have successfully produced velocity maps for the cases considered and also matched our results with additional observable e.g, proton temperature. This pilot study is the first step in developing an indigenous space weather framework which is an absolute need of hour as Indian Space Research Organisation (ISRO) is coming up with Aditya-L1 to study the properties of Sun at L1 (https://www. isro.gov.in/aditya-l1-first-indian-mission-to-study-sun). A full-fledged 3D-MHD run of solar wind using PLUTO code will be carried out in our subsequent paper for predicting other observable quantities like the number density and IMF magnetic field. Further, including the propagation of CMEs in such a realistic solar wind background would be carried out with an aim to study its impact on space weather.