A neural network performs a sequence of transforms by multiplying its inputs with a set of coefficients (weights) and applying nonlinear functions to provide the output. It has been shown that the neural network can approximate any continuous function (Cybenko, 1989). For a supervised network, the weights are adjusted to produce a desired output, given the inputs, known as the training phase. The training phase requires known input and target values, a cost function, and an optimization algorithm that minimizes the cost.

The Elman RNN was first developed for Dst predictions by Wu and Lundstedt (1997) and later implemented for real-time operation (Lundstedt et al., 2002). In this work, we use the term Elman network, but it is the same as the simple RNN used in the Tensorflow package that we use (Abadi et al., 2015). Using linear units at the output layer, the Elman network at time t is described by y t = b + ∑ j = 1 J v j h t j = V h t , (5) h t j = f ( a j + ∑ i = 1 I w j i x t i + ∑ k = 1 J u j k h t − 1 k ) = f ( W x t + U h t − 1 ) j , (6) with the output layer bias b, J hidden weights v j , nonlinear activation function f, J hidden layer biases a j , J × I weights w j i , and J × J recurrent weights u j k . Note that we use superscripts i , j ,  and  k as indices, not powers. The equations can be written more condensed using weight matrices W and U , where the bias terms ( a j , b ) have been collected into the matrices and increasing the lengths of x t and h t by adding a constant set to one. For example, in Lundstedt et al. (2002), there are I = 3 inputs ( B z , n , V ) and J = 4 hidden units.

A minimalistic Elman network can be constructed by using only one input unit and one linear unit in the hidden layer, thus b = 0 , v 1 = 1 , leading to y t = a 1 + w 11 x t 1 + u 11 y t − 1 , (7) which after some rearranging can be written as y t = y t − 1 + a 1 + w 11 x t 1 − ( 1 − u 11 ) y t − 1 , (8) which is identical to Eq. 2 for τ = const and Δ t = 1 , and by letting a 1 = 0 , w 11 x t 1 = Q ( t ) , and 1 − u 11 = 1 / τ . This simple network is trained using the pressure-corrected Dst index as the target. As the weights in the network are initiated with random values before training begins, there will be some variation in the final weight values if the training is repeated. We find typical values of w 11 and u 11 corresponding to a ∈ [ − 2.4 , − 2.7 ] (Eq. 3) and τ ∈ [ 14,16 ] hours, which are close to the values used by O’Brien and McPherron (2000). However, the algorithm can get stuck in local minima that results in quite different values. We provide code on Github ² for the minimalistic Elman network (see Model005.py).

The gated recurrent unit (GRU) neural network (Cho et al., 2014) has a more complex architecture than the Elman network. We implement a single GRU layer, and the output from the network is as before, given by y t = V h t (Eq. 5). The GRU layer output at unit j is h t j = z t j h t − 1 j + ( 1 − z t j ) h ˜ t j , (9) where z t j is the update gate and h ˜ t j is the candidate activation. The update gate is defined as z t j = σ ( W z x t + U z h t − 1 ) j , (10) where σ is the sigmoid function with output range 0–1. The weight matrix W z operates on the input vector x t , and the matrix U z operates on the past activation h t − 1 . The candidate activation is defined as h ˜ t j = f ( W x t + U ( r t ⊙ h t − 1 ) ) j , (11) where f is a nonlinear function with two additional weight matrices W and U . The U matrix operates on the past activation weighted by the reset gate r t j = σ ( W r x t + U r h t − 1 ) j , (12) which has a further set of weights matrices W r and U r . Clearly, the GRU network is more complex than the Elman network, and it has approximately 3 times more weights than the Elman network for the same number of units. As the update and reset gates have outputs between 0 and 1, we see that when both produce ones [ ( z t j , r t j ) = ( 1,1 ) ] , the GRU network simplifies to the Elman network. On the other hand, when z t j = 0 , no information of the current input x t is used, only the past state h t − 1 . Finally, when r t j = 0 , no information of past states goes through the candidate activation (Eq. 11); information on past states only goes through Eq. 9 and is weighted by 1 − z t j .

The long short-term memory (LSTM) neural network (Hochreiter and Schmidhuber, 1997) was introduced before GRU and has further complexity with the number of weights approximately four times that of the Elman network, given the same number of units. We will not repeat the equations here but instead refer to for example, Chung et al. (2014). The LSTM has three gating functions, instead of GRU’s two, that control the flow of information: the output gate, the forget gate, and the input gate. When they have values 1, 0, and 1, respectively, the LSTM simplifies to the Elman network.

Given a network with a sufficient number of weights, it can be trained to reach zero MSE; however, such a network will have poor generalizing capabilities; that is, it will have large errors on predictions on samples not included in the training data. Different strategies exist to prevent over-fitting (Goodfellow et al., 2016). We apply early stopping and dropout (Srivastava et al., 2014; Gal and Ghahramani, 2016).

In order to make a robust estimation of the performance of the networks, we apply k-fold cross-validation (Goodfellow et al., 2016). During a training session, one subset is held out for testing and the remaining k − 1 subsets are used for training and validation, and out of the k − 1 subsets, one is used for validation and the remaining k − 2 subsets for the training. During training, the validation mean squared error (MSE) is monitored, and the network with lowest validation MSE is chosen (early stopping). In practice, to know that the minimum validation MSE has been reached, the training is continued for a number of epochs after the lowest MSE has been reached. The final evaluation of the models is performed on the k different test sets (see Section 2.2 to know how the different sets are selected).

The hourly solar wind data and Dst index are obtained from the OMNI dataset (King and Papitashvili, 2005). The inputs are the solar wind magnetic field magnitude B, the y- and z-components ( B y , B z ) in the geocentric solar magnetospheric (GSM) coordinate system, the particle density n, and speed V. To provide information on diurnal and seasonal variations (O’Brien and McPherron, 2002), four additional variables are added: the day of year parameterized as sine and cosine of 2 π DOY / 365 , and local time as sine and cosine of 2 π UT / 24 . Thus, in total, nine input variables. Several previous geomagnetic prediction models also use diurnal and seasonal inputs (e.g., Temerin and Li (2006); Wintoft et al. (2017); Wintoft and Wik, (2018)). Many different coupling functions (Q) for the dayside reconnection rate have been suggested and investigated (Borovsky and Birn, 2014), but as the neural network can approximate any function, the exact function does not have to be specified as long as the relevant inputs are available.

The target variable (Dst) depends on both the solar wind and past states of the system, where past states can be described by past values of Dst itself or by past values of the solar wind. We choose to only include solar wind, thereby not relying on past observed or predicted values of Dst. For the RNN training algorithm, the data are organized so that the past T solar wind observations are presented at each time step. The input data are thus collected into a N × T × 9 tensor, and the target data have N × 1 dimension, where N is the number of samples in the set. The input history should be long enough to capture typical storm dynamics, and we found that validation errors leveled out at T = 120 hours (see also the results regarding T in Sections 2.3 and 2.4).

To implement the k-fold cross-validation (CV), the dataset must be partitioned into subsets; we perform a five-fold CV. We choose the five sets to each have similar target (Dst) mean and standard deviation so that training, validation, and testing are based on comparable data. If a more blind approach were done, then there is a high risk that training is performed on data dominated by storms, while testing is performed on more quiet conditions. Further, the samples in a subset cannot be randomly selected because there will be considerable temporal overlap between samples of different subsets due to the T = 120 hour window. Instead, we build the subsets from data covering complete years. The data we use cover the years 1995 to 2015, extending over almost two solar cycles and with few solar wind data gaps. We define five subsets based on the data for the years shown in Table 1. The datasets used for training, validation, and test are selected by cycling through the subsets. For the first CV (CV-1) subset, one is selected as test set; subsets two, four, and five for training; and subset three for validation. The process is repeated according to Table 2.

The input and target values span very different numerical ranges, whereas the training algorithm should receive input-target data that have similar numerical ranges. Therefore, the input and target data are normalized, where the normalization coefficients are found from the training set. By subtracting the mean and dividing with the standard deviation for each variable separately, the training set will have zero mean and one standard deviation on all its inputs and target variables. However, as the distributions for each variable are highly skewed, they result in several normalized values with magnitudes much larger than one. Another way to normalize is to instead rescale the minimum and maximum values to the range [-1, 1]. This guarantees that there will be no values outside this range for the training set. We found that the min–max normalization gave slightly better results, especially at the large values.

There are a number of hyperparameters (HPs) that control the model complexity and training algorithm that need to be tuned, but it is not feasible to make an exhaustive search. Initially, a number of different combinations of HP values were manually tested to provide a basic insight into reasonable choices and how the training and validation MSEs vary with epochs. In this initial exploration, we found the Tensorboard (Abadi et al., 2015) software valuable in monitoring the MSE.

The Adam learning algorithm (Kingma and Lei, 2015), which is a stochastic gradient descent method, has three parameters: learning rate ϵ and two decay rates for the moment estimates ( β 1 , β 2 ) . We fix the latter to the suggested values β 1 , β 2 = 0.9 , 0.999 and vary ε ∈ [ 5 ⋅ 10 − 4 , 1 ⋅ 10 − 3 , 2 ⋅ 10 − 3 ] .

The learning algorithm updates the weights in batches of samples from the training set, where the number of samples in each batch N B is much smaller than the total number of training samples ( N B ≪ N ) . We test batch sizes of N B ∈ [ 32,64,128 ] . One training epoch includes approximately N / N B training iterations in which the weights are updated at each iteration.

The model capacity is determined by the number of weights and the network architecture. In this work, we have one input layer, a recurrent layer (hidden layer), and a single output. Thus, the capacity is determined by the network type ( M ∈ [ Elman , GRU , LSTM ] ) and the number of hidden units ( N H ) .

The current state ( h t ) in the RNN depends on both its inputs ( x t ) and the past state ( h t − 1 ) . For computational performance reasons, past states are not kept indefinitely; instead, there is a limit T on the length of the memory. We explored T ∈ [ 48,72,96,120 ] hours and found that the validation MSE decreased with increasing T, but that it leveled out for large T. We therefore set T = 120 hours. This also means that any dynamical processes extending past 120 h cannot be modeled internally by the RNN. The choice of T for the Elman and GRU networks is studied on simulated Dst data in the next section.

The dropout is controlled by parameters that specify the fraction of network units in a layer that are randomly selected per epoch and temporarily disregarded. The dropout can be applied to all layers: the input layer ( d i ) , the recurrent layer ( d r ) , and the hidden layer ( d h ) . The dropout is a number between 0 and 1, where 0 means all units are included and 1 that all units are unused.

For each combination of HP that we explore, we train three networks initiated with different random weight values as there is no guarantee that the training algorithm will find a good local minimum. The network with the lowest validation error is selected. Note that here validation refers to the split into training and validation sets used during training, which is different from the cross-validation sets that make up the independent test set.

It is interesting to study the RNNs on data generated from a known function relating solar wind to Dst, and for this purpose, we use the AK1 model (O’Brien and McPherron, 2000). Using the datasets defined in the previous section, we apply the AK1 model to the solar wind inputs and create the target data. Thus, there exists an exact relation between input and output, and the learning process of the RNN will only be limited by the amount of data, network structure (type of RNNs), and network capacity (size of RNNs). We showed that the minimalistic Elman network (Eq. 7) can model the pressure-corrected Dst. The AK1 model also includes the pressure term, and its inputs are B z , n, and V. The five-fold CV is applied to Elman and GRU networks, and we vary the time window T and the network size N H .

In Figure 1, the validation errors as function of T are shown for the Elman and GRU networks. At each T, the optimal networks with respect to N H are used. We see that for small T, the RMSE is large, but it is similar for the two network types. At small T, only part of the storm recovery phase can be modeled. But as T is increased, the RMSE becomes much smaller for the GRU network than for the Elman network. It is likely that the Elman network suffers from the vanishing gradient problem (Bengio et al., 1994): the reason for introducing GRU and LSTM networks. We also see that the GRU network reaches an RMSE of lower than 0.6 nT, which could be further decreased by increasing T. Thus, the GRU network can learn the AK1 model using the observed solar wind data.

As described in Section 2.2, the inputs to the Dst model are solar wind magnetic field ( B , B y , B z ) , density (n), and speed (V); the day of year parameterized as sine and cosine of 2 π DOY / 365 ; and local time as sine and cosine of 2 π UT / 24 . The DOY and UT are added to model the seasonal and diurnal variations in Dst (O’Brien and McPherron, 2002).

We perform a search in the hyperparameter space as described above and conclude that training is not very sensitive on the learning rate (ϵ) or batch size ( N B ) , and therefore fix them at ( ε , N B ) = ( 10 − 3 , 128 ) .

For each of the five splits, we select the corresponding training set (Tables 1,2), and RNNs are trained with different number of hidden units ( N H ) and different dropout rates ( d i , d r , d h ) . For each combination of ( N H , d i , d r , d h ) , three networks are trained starting from different random initial weights. During training, the validation error is monitored, and the network with lowest validation error is selected. The training is stopped 20 epochs after the minimum validation error has been reached, and the network at minimum validation error is saved. Typically, the minimum validation RMSE is found after 40 to 80 epochs. This results in five different networks for each HP combination that can be tested using the CV approach. See Appendix for software used and typical training times.

The coupling function from solar wind to observed Dst is subject to a number of uncertainties, and to provide a few examples: the solar wind data have been measured at different locations upstream of Earth, mostly from orbit around the L1 location, and then shifted to a common location closer to Earth (King and Papitashvili, 2005); We rely on a point measurement; there may be both systematic and random errors in the derived Dst index. The uncertainties introduce errors in the input–output mapping, and to reduce their effect and improve generalization, we apply dropout. From a search of different combinations of ( N H , d i , d r , d h ) , it was found that dropout on the inputs ( d i ) always led to poor performance, which can be understood as several inputs individually are critical, for example B z . Therefore, we set d i = 0 . The performance improved when dropout was applied on the recurrent and hidden layers. Figure 2 shows the training and validation RMSE as a function of N H for different dropouts. In the case with no dropout (left panel), it is seen that the GRU and LSTM validation errors are similar and significantly below the Elman validation errors. There is also a large gap between the training and validation errors, indicating over-fitting on the training set. When dropout is introduced (middle and right panels), the network size must be increased to reach similar validation RMSE as when no dropout is applied, which is expected as only a fraction of units are active at any one time. But we also see that the gap between training and validation errors decreases. We also applied dropout on the Elman network, but the validation errors became large when d r > 0 ; therefore, the results are not included in the middle and right panels. When d r = 0.5 and d h = 0.5 (right panel), the optimal GRU and LSTM networks have N H = 50 and N H = 40 , respectively. In terms of the number of weights, they are of similar sizes, 9,051 and 8,041, respectively. When dropout is applied, the number of active weights drops to 2,651 and 2,421

For each CV set, we select the GRU and LSTM networks with minimum validation RMSE with and without dropout, and run them and collect the 5 CV sets into one set. We thereby get an estimate of the generalization performance for the whole 1995 to 2015 period. Figure 3 shows scatterplots of predicted Dst vs. observed Dst on the test sets for different networks. Table 3 summarizes the performance on the training, validation, and test sets. The 95% confidence intervals have been estimated by both assuming independent data points and taking into account the autocorrelation (Zwiers and von Storch, 1995). It is clear that using dropout significantly improves the generalization capability. We also see that there is no significant difference between the GRU and LSTM networks. The bias (mean of errors) and linear correlation coefficient are computed on the test set and shown in Table 4.

The performance of the networks varies with the level of Dst; the errors have a tendency to increase with the magnitude of Dst. Figure 4 shows the RMSE binned by observed Dst. Down to D s t = − 300  nT , the networks with dropout have the lowest RMSE. The bin at D s t = − 400  nT has too few samples to be interpretable. The main reason the errors increase with the magnitude of Dst is that there are very few samples around the extremes; thus, the uncertainty of the function estimation will be large. In Wintoft et al. (2017), this problem was addressed by using an ensemble of networks; the predictions from several networks with different weights were averaged. In this work, we study the use of dropout not only in the training phase but also in the prediction phase. The algorithm that temporarily cancels units at random during training can also be applied during prediction. This means that there is practically an indefinitely number of weight combinations that can be used to produce an arbitrary number of predictions at each time step. For the GRU network with ( N H , d r , d h ) = ( 50,0.5 , 0.5 ) , there are more than 10 28 possible combinations. There is a Bayesian interpretation of dropout as the estimation of model uncertainty (Gal and Ghahramani, 2016b). The idea is that the weights are random variables leading to a distribution of predictions for fixed inputs. For each sample, a large number of predictions can be generated, which randomly use different combinations of network units. For each sample, we generate 100 predictions and compute the mean and standard deviation. Figure 5 shows two examples, the first a severe geomagnetic storm and the second a major storm. The mean predictions with dropout come close to the predictions without dropout. The prediction uncertainty is small during quiet times (Dst close to zero) and increases with storm magnitude. Again, this is a result of the greater uncertainty in parameter estimates in regions which are poorly sampled.

As time of day and season are included in the inputs, the network can model diurnal and seasonal variations in Dst. These variations are not strong, and the left panel in Figure 6 shows Dst for all years averaged over month and UT hour. Running the GRU networks on the test data from the five CV sets reveals a very similar pattern (second panel from left). Thus, the network shows similar long-term statistics considering that it is driven only by solar wind and time information. The two left panels contain contributions from all levels of Dst from quiet conditions to storm conditions. But we may now simulate solar wind conditions that we can define as quiet conditions. The two right panels show predicted Dst, assuming solar wind flowing out from the Sun (GSEQ system) along the Parker spiral with a 45 ∘ angle at L1 at two different speeds, 350 km/s and 400 km/s, respectively. In this configuration, B z = 0 in the solar coordinate system, but via geometric effects (Sun’s and Earth’s tilts with respect to the ecliptic and Earth’s dipole tilt), B z will be nonzero in the GSM system showing diurnal and seasonal variations (Lockwood et al., 2020).