Reconstructing Solar Wind Profiles Associated With Extreme Magnetic Storms: A Machine Learning Approach

The lack of data on solar wind have prevented a detailed understanding of extreme magnetic storms. To address this issue, we apply a machine learning technique in the form of an Echo State Network (ESN) to reconstruct solar wind data for several extreme magnetic storms for which little or no solar wind data were previously available. Multiple geomagnetic activity indices are used as the input data for the ESN, which produces a continuous time series of solar wind parameters as output. As a result, the solar wind parameters for the largest storm event in March 1989 are obtained, and the minimum Bz is estimated to be −95 nT ±10 nT. Two different types of solar wind profiles are discussed for the extreme magnetic storms, a sheath‐driven profile and a magnetic cloud‐driven profile. The results reported here will be highly useful as input data for future simulation studies modeling extreme magnetic storms.

used for evaluation are presented in Section 4, along with the results of a superposed epoch analysis of two different types of extreme magnetic storms. In Section 5, we discuss our findings and their implications. Concluding remarks are presented in Section 6.

Data Set
The solar wind parameters were obtained from the OMNI2 data set, which is available at https://omniweb.gsfc. nasa.gov/ow.html. OMNI2 gives the long-term time series of the solar wind parameters at the Earth's position. The time cadence is 1 hr, which was calibrated by time shifting and averaging from available spacecraft data. The available OMNI2 solar wind data covers the period from November 28, 1963, although there are frequent data gaps before 1994. In this study, we use the north-south component of the interplanetary magnetic field (Bz) in the geocentric solar magnetosphere (GSM) coordinate system, the number density of protons (Np), and the solar wind speed (Vsw). These three parameters are the most important parameters driving the geomagnetic activities via exciting the magnetospheric and ionospheric current systems (Akasofu, 1981;Borovsky, 2008;Newell et al., 2007).
The most widely used indices of geomagnetic activity are AE, Kp, and Dst, each of which responds differently to the solar wind parameters because of their respective physical meanings: auroral electrojet (AE), subauroral disturbances (Kp), and ring current development (Dst). For example, the primary parameter affecting the Dst index is Bz, while for the AE and Kp indices the primary parameter is Vsw. It is expected that the slight difference in the geomagnetic indices includes information related to the solar wind parameters, which, we propose, can be extracted by the machine learning technique developed in this study (and described in detail below). There are many examples of attempts to reproduce the time series of AE, Kp, and/or Dst from solar wind parameters Bz, Np, and Vsw using various methods, including machine learning techniques (Liemohn et al., 2018).
Geomagnetic activity indices have been available since before the space age. For example, values of the Kp index are available at a 3-hr time cadence for all past years since 1932, at GFZ Potsdam (Matzka et al., 2021). In this study, we converted the 3-hr Kp index to consecutive hourly values in the Kp*10 index, following the OMNI2 format, e.g., Kp = [0+, 3−, 1+] becomes Kp*10 = [0, 0, 0, 26.7, 26.7, 26.7, 13.3, 13.3, 13.3]. The Dst index is available at a 1-hr time cadence for the years 1957 through 2014 at WDC Kyoto. A provisional Dst index is available for three subsequent years, 2015, 2016, and 2017. The hourly data for the AE index are also available for the years 1957 through 2014 at WDC Kyoto. However, there are data gaps in the AE index for the period from January 1, 1976 to December 31, 1977 and the period from July 1, 1988to December 31, 1989(except for March 1989. A provisional AE index is available for the years 1990 through 2017. After 2018, only quick-look data are available for both the Dst and AE indices. Given the data availability described above, we have limited our analysis to the time period from January 1, 1957 to December 31, 2017.

Method of Analysis: Echo State Network
The Echo State Network (ESN) was developed by Jaeger (2001) as a type of recurrent neural network. The architecture has been investigated and used by many authors (Hart et al., 2020;Tanaka et al., 2019). ESN is now known as a novel method that works with a relatively small set of training data, which is especially suitable to the analysis of times series. The first application of ESN to space physics appears in Nakano and Kataoka (2021). In ESN, the input vector u, reservoir state vector x, and output vector y are defined by an N-point time series where n = 1, 2, 3, … N.
The reservoir state vector x consists of a large number of nodes, which are updated in time with the input vector u and the previous state of the nodes as follows: We use the hyperbolic tangent function as function f and fix the weight matrices W in and W. To create the random and sparse node connections of W, we set the number of nodes, Nx, equal to 10 3 , where only 10% of the matrix elements are random values between −1.0 and 1.0, and the remaining 90% are zero, that is, the network connection density = 0.1. We chose 0.95 as the spectral radius (maximum eigenvalue) of W in order to satisfy the echo state property guaranteeing the independence of the reservoir state from the initial values when the spectral radius is below unity (Jaeger, 2001).
Thus, there are several parameters to change the model structure. In the present problem, selecting different number of nodes and/or the connection density are not sensitive to change the structure, and the selection of spectral radius is the most sensitive. Figure S2 in Supporting Information S1 shows the dependence on the spectral radius from 0.5 to 1.0. In this paper, we select 0.95 as the best value of the spectral radius to minimize the errors of both test data and training data for Bz and Vsw. To minimize the error of Np, however, smaller spectral radius gives better scores for the test data. Such a separate optimization for each parameter can therefore be a future work to improve the model, although it is beyond the scope of this paper.
The output vector y is calculated using the linear combination of the output weight matrix and the reservoir state vector as follows: where out E W is the output weight matrix. We train only the output weight matrix by the set of T-point time series of input vectors X and desired output vectors D: The least squares method for minimizing the difference between the outputs y and d can be represented by a standard linear regression as follows: We use the AE, Kp, and Dst indices as the three-dimensional input vector (Nx = 3), and Bz, Np, and Vsw as the three-dimensional output vector (Ny = 3). In this study, we take the logarithms of the density Np and speed Vsw to avoid negative outputs. It is noteworthy here that both the density and speed follow log-normal distributions (Burlaga & Lazarus, 2000). Note, too, that it is the solar wind that drives geomagnetic activities and not the other way around. Considering such a cause-and-result relationship, the time stamp is reversed when we apply ESN, that is, all of the N-point time series are sorted by n = N, N − 1, N − 2, … , 3, 2, 1. Here, we assumed that OMNI2 data was perfectly time-shifted already to match the geomagnetic indices and solar wind drivers in time.
To prepare the training data set for X and D in Equation 4, we constructed a special time series by selecting and connecting the storm events, which was intended to ensure that the model was reasonable, especially during disturbed periods. In so doing, we first detected the storm days using the 1-day minimum amplitude of the Dst index below −100 nT, with the previous 3 days of Dst > −100 nT. We then discarded the storm days if there were any data gaps in the solar wind parameters during the 7-day time period spanning the storm day, from 3 days before to 3 days after. Finally, the time series of selected storm events were connected in time. The training data set was constructed from the 28 years of data from 1990 to 2017, a period during which the OMNI2 solar wind data has only infrequent data gaps. We identified a total of 74 storm events during the selected period (1990-2017) to create a time series of 74 × 7 × 24 = 12,432 data points. The percentage of how much of storm events had complete data is ∼51% as we discarded the other 70 data-gap storm events. The training data set and the event list of storms are provided in Data Set S1.

Reconstruction of the Solar Wind Parameters in March 1989
The largest observed magnetic storm of minimum Dst = −589 nT occurred in March 1989. Figure 1 shows the corresponding geomagnetic activity indices and the solar wind parameters available in the OMNI2 data. As can 4 of 10 be seen here, there are large data gaps in the solar wind parameters, especially around the magnetic storm on March 13 and 14. The red lines show the reconstructed solar wind parameters using the ESN model. The most notable result is that the modeled minimum Bz reaches −95 nT during the magnetic storm. It is also noteworthy that the ESN model fails to predict the northward Bz (∼10 nT was predicted to be ∼5 nT on March 27) and underestimates density (∼40/cc was predicted to be ∼20/cc on March 27) and speed (>700 km/s was predicted to be ∼600 km/s on March 19). Although these errors are not small, as described in the next subsection, the example illustrates the first step in reconstructing the entire time series for such large magnetic storms. Futher examples of storm-time comparisons are shown in Figures S3-S6. To evaluate the errors in our analysis, the standard metrics were calculated using the storm-time data (a total of 70 storm events) for the time period from 1990 to 2017, in which the gap-less training data were excluded. The cross-correlation coefficients were 0.72, 0.42, and 0.61 for Bz, Np, and Vsw, respectively (as documented in Figure S1 in Supporting Information S1). The root mean square errors (RMSE) for Bz, Np, and Vsw were 4.0 nT, 7.2/cc, and 104 km/s, respectively. Normalizing the RMSE values by the respective standard deviations produced normalized root mean square errors (NRMSE) of 0.71, 1.0, and 0.79. The error is large for the large values, for example, ∼10 nT for |Bz| >10 nT events, as shown in Figure S2 in Supporting Information S1.

Occurrence Distribution of Extreme Events
To investigate the occurrence distribution of extreme events, the solar wind parameters of all of the storm events from 1957 to 2017 were reconstructed. The absolute values of the minimum Bz (SBZ) and the density at the minimum Bz timing were extracted for each storm event to create the scatter plot in Figure 3. Figure 3 shows the cumulative distribution function (CDF) of the SBZ values. After extending the CDF curve by fitting the power-law distribution (Kataoka, 2013(Kataoka, , 2020, the one-in-100-year-event SBZ is estimated to be 137 nT. It is noteworthy that such a large amplitude of SBZ ∼ 100 nT was observed at the STEREO-A spacecraft during the July 2012 event (Russell et al., 2013) in which an extreme CME missed the Earth by only a few days.

Superposed Epoch Analysis for Extreme Magnetic Storms
Based on the density values scattered around ∼10/cc as shown in the top panel Figure 3, possible sheath-like (high density > 15/cc) and cloud-like (low density < 5/cc) large-amplitude (Bz < −25 nT) events were selected for further superposed epoch analysis. Note that this is a biased selection, not based on a cluster analysis, just to discuss possibly different populations. We omitted the medium-density events (5-15/cc) because they can belong either high or low events considering the large error in density estimation of our model.
The left panels in Figure 4 show the result for the sheath-like events. The basic characteristic of the rapid evolution of these storms main phase with spiky AE enhancement (less than ∼6 hr) is clearly indicated by the superposed epoch analysis for sheath-type selection.
The right panels in Figure 4 show the result for the cloud-like events. The main phase and duration of strong SBZ are several times longer than those of the sheath-like storms. Moreover, the peak amplitude of SBZ is >50 nT, 7 of 10 which is also a few times larger than that of the sheath-like storms. The resultant magnetic storms of this population are, therefore, largest. The latest example of this type of storm occurred in November 2003 (Kataoka, Fairfield, et al., 2005). The most prominent difference in these extreme solar wind-driven storms arises from their low-Mach number state, as discussed in the next section.

Discussions
The AE and Kp indices in sheath-like events are comparable to those of cloud-like events, as shown in Figure 4, despite of the smaller amplitude of SBZ. This is consistent with the fact that the high-beta sheath regions have higher solar wind-magnetosphere coupling efficiencies than magnetic clouds (Myllys et al., 2016). In fact, the majority of the largest Kp values are related to the sheath regions (Kilpua, Isavnin, et al., 2013). The density effect on the AE index was also recently discussed by Ebihara et al. (2019) and Nakano and Kataoka (2021).  Tables S1-S2 in Supporting Information S1.
We turn now to a discussion of the most extreme situation of cloud-like events, which becomes especially important for discussing the 100-year event shown in Figure 3. Assuming the fraction of α particles to be a constant value of 5% for simplicity, we can calculate the Alfven Mach number as follows: The Alfven Mach number ranges from 7 to 11 in a typical high-speed solar wind, but is less than half that during CMEs or magnetic clouds (Lavraud & Borovsky, 2008). The extreme magnetic storms, therefore, tend to occur under a low Alfven Mach number.
In the special situation of a low Mach number, the bow shock has been observed far from the magnetopause (Fairfield et al., 2001). The magnetospheric interaction with the low-Mach solar wind can also be different from the usual high-Mach state because of the weak bow shock and resultant low-beta flow in the magnetosheath. Global MHD simulations predict that due to the relatively strong magnetic force, the cross-section of the magnetopause shape and the surrounding magnetosheath flow becomes asymmetric (Lavraud & Borovsky, 2008). In fact, Lavraud et al. (2007) reported in-situ observation data showing that the magnetosheath flow was accelerated beyond 1,000 km/s from the upstream solar wind speed of 650 km/s when the Alfven Mach number was between 2.0 and 4.4. Other examples of low-Mach events include the inner magnetospheric response (Lugaz et al., 2016). Recent MHD simulation results indicate that the saturation of cross polar cap potential is dependent on the Alfven Mach number (Lakka et al., 2018).
Even more extreme cases of sub-Alfvenic (M A < 1) solar wind exist. It is suggested that the magnetosphere has a different shape, called Alfven wings (Chane et al., 2012), rather than the usual comet tail shape. Ridley (2007) showed global MHD simulation results indicating that the northern and southern lobes were separated. The formation of the Alfven wings is one of the proposed mechanisms to explain the saturation of cross polar cap potential (Kivelson & Ridley, 2008;Ridley, 2007).
Future MHD simulations based on the average solar wind parameters obtained in this study should provide essential information for understanding the energy flow from the solar wind to the magnetosphere and ionosphere during extreme magnetic storms. However, for a detailed modern analysis, solar wind data with a time resolution higher than hourly (for example, 5 min values) will be necessary as input parameters for the simulations. Reproducing such high time resolution data remains a future work; however, there are several physical constraints regarding cloud-like and sheath-like events. For cloud-like events, it is possible to simply interpolate the hourly data into five-minute data, assuming that the magnetic field inside the magnetic cloud rotates smoothly. On the other hand, the fluctuation level of the magnetic field is strong in the sheath region (Kilpua, Hietala, et al., 2013) due to the existence of Alfvenic discontinuities (Kataoka, Watari, et al., 2005). The rapid changes of the disturbed condition embedded in the sheath region will also be important for reproducing the large geomagnetically induced currents.
It should be noted that the estimated minimum value of Bz = −95 nT in the March 1989 event (Figure 1) is significantly larger than the value of Bz = −50 nT previously estimated by Nagatsuma et al. (2015). They assumed the high-Mach number regime of M >> 1 in which the compression ratio at the bow shock was approximated as 4. The observed strength of ∼200 nT in the magnetosheath can, therefore, be converted to ∼50 nT in the solar wind. However, the compression ratio becomes ∼2, if M ∼ 2. Our estimate is thus consistent with the assumption of a low Mach number regime of M ∼ 2.
Our results provide a continuous data set for the ∼60-year period since 1957, well before the OMNI2 data set. The data should prove extremely useful for interdisciplinary studies of the space climate, including efforts such as modeling the possible long-term impact of radiation belts and energetic electron precipitation events to the middle atmosphere via atmospheric chemistry (Seppala et al., 2015). While there are several other space climate studies focused on the reconstruction of long-term solar wind data using the geomagnetic activity indices (e.g., Lockwood, 2013;Svalgaard & Cliver, 2007), the main difference between those studies and ours is the time resolution. Whereas the other studies discuss a time variation of 27-day averaged values, our study focuses on hourly time variation and represents the first step in high-resolution (hourly time cadence) space climate study.
Filling the data gaps is another important aspect of the ESN method applied here. We recognize that there have been a number of excellent prior studies on this topic. For example, using a gap-filling method based on Singular Spectrum Analysis (Kondrashov et al., 2010, Kondrashov et al., 2014 report NRMSE values of 0.77, 0.48, and 0.10 for Bz, Np, and Vsw, respectively, which are better scores than our results, except for Bz. To facilitate further comparisons with other models, the hourly data produced in our study for Bz, Np, and Vsw covering the period from 1957 to 2017 are provided in Data Set S2. The output data from our study gives a starting point for improving the metrics. There are several possible pathways forward. For example, using other machine learning techniques appears promising. Using other datasets, such as the SME index of the SuperMAG project, also holds promise. Finally, we note that the methodology developed in this study can be generally used to complement datasets containing data gaps in various fields of science where there are limited observation data available.

Conclusions
With the help of an Echo State Network, a type of recurrent neural network, we discussed sheath-like and cloudlike average patterns in solar wind profiles associated with extreme magnetic storms. We also constructed a continuous solar wind data set for the years 1957 through 2017. Using our approach, the SBZ amplitude during the extreme March 1989 event was estimated to be 95 nT ± 10 nT, while for the one-in-100-year-event, the SBZ amplitude was estimated to be 137 nT.