© A. Bhaskar & G. Vichare, Published by EDP Sciences 2015

1 Introduction

Transient ejections from the Sun setup large-scale disturbances in the interplanetary space. These disturbances interact with the Earth’s magnetic field, resulting into the severe space weather events, such as geomagnetic storm, substorm, etc. As the present space-technology is vulnerable to the geomagnetic disturbances, predicting geomagnetic field response well in advance is an important aspect of space weather studies. Long duration southward interplanetary magnetic field injects solar wind energy into the Earth’s magnetosphere–ionosphere system mainly through reconnection (Gonzalez et al., 1994). This results in the azimuthal drift of the charged particles inside the magnetosphere, establishing ring current in the equatorial plane. Intensification (main phase) and decay (recovery) of the storm time ring current consist of different processes. The main phase is primarily controlled by the solar wind conditions, whereas decay of the ring current has a major contribution from the internal magnetospheric processes. Due to the varying nature of the storm sources, the magnetospheric dynamics and the energy budget involved in each storm differs considerably (Vichare et al., 2005). The injection of solar wind particles and transmission of solar wind electric field generate various currents in the magnetosphere–ionosphere system such as cross-tail current, field-aligned currents, partial ring current, etc. (Ohtani, 2000). Moreover, sudden variations in the dynamic pressure of the solar wind alter magnetopause current and tail current. Also, they produce transient ionospheric currents (Vichare et al., 2014; Oliveira & Samsonov, 2018). The recovery phase of the ring current during the geomagnetic storm has an influence of various nonlinear phenomenon like wave-particle interaction, charge exchange, the ionospheric outflow of O⁺ ions, particle precipitation, etc. (Daglis et al., 1999). Superposed effect of these currents and magnetospheric nonlinear processes in the magnetosphere–ionosphere system makes the prediction of storm-time temporal variations of ring current a challenging task.

Ground magnetometer measures the integrated effect of all these disturbed time and also quiet time ionospheric and magnetospheric currents. Geomagnetic indices like Disturbance storm time index (Dst) and Symmetric H-component (SYMH) index mainly represent ring current intensity during geomagnetic storms (Sugiura, 1964; Rangarajan, 1989; Wanliss & Showalter, 2006; Vichare, 2019), derived using the longitudinally distributed chain of low latitude ground-based magnetometers. SYMH is same as Dst, but it has a 1-minute temporal resolution, which is very useful to study short temporal variations during the geomagnetic disturbances. SYMH is derived by first subtracting main geomagnetic field due to internal geodynamo and external Sq induced geomagnetic field variations and then averaging residual fields. Therefore, it is a good proxy for the longitudinally symmetric component of the ring current. By removing the globally symmetric component of the magnetic field variations from geomagnetic field variations at each station, longitudinally asymmetric geomagnetic field variations are derived. The range between maximum and minimum of these subtracted fields are compiled as ASYH index. ASYH have a significant contribution from various transient currents flowing in the magnetosphere–ionosphere system such as currents associated with sudden impulses, solar flares, substorms and prompt penetration electric fields, partial ring current, field-aligned currents, magnetotail current, etc. (Clauer & McPherron, 1980; Iyemori & Rao, 1996; Singh et al., 2012, 2013). Normally, during geomagnetic storms, these asymmetric currents also get enhanced. Therefore, ASYH index is a good proxy for monitoring the currents, flowing in the magnetosphere–ionosphere system during geomagnetic storms, which are not azimuthally symmetric, but can affect mid-latitude observations. The contribution of substorms in ring current is a widely debated topic as some researchers believe to have a significant contribution and some believe it is weak. Newell & Gjerloev (2012) showed that the substorm effect in the ring current is very small. Moreover, Munsami (2000) showed that when Dst station lies under a substorm current wedge, then only they observed significant contamination of ring current due to substorms.

There are a lot of efforts to understand the relationship between the ring current (SYMH) and partial ring current (ASYH) respectively. Generally, it is observed that during the main phase of the geomagnetic storm, ring current is highly asymmetric and becomes symmetric in the late recovery phase (Siscoe et al., 2012; Jordanova et al., 2003). Liemohn et al. (2001) reported that major part of magnetic field variations during the main phase of geomagnetic storms is due to asymmetric ring current. However, there are storms which show symmetric nature of the ring current even during the main phase which remains unexplained (Newell & Gjerloev, 2012). The well known Love-Gannon relationship states that the difference between dawn and dusk disturbance-field (similar to ASYH index) at low latitudes is linearly proportional to Dst. However, (Siscoe et al., 2012) pointed out that this relationship can be explained only through field-aligned currents.

As there is a number of studies investigating the relationship between symmetric and asymmetric ring current, at the same time efforts are going on to give a more accurate prediction of these indices during geomagnetic storms. To forecast these geomagnetic indices (Dst, SYMH, AE, etc.) different approaches have been attempted (Williscroft & Poole, 1996; Wu & Lundstedt, 1996; Wu et al., 1998; Weigel et al., 1999; O’Brien & McPherron, 2000; Gholipour et al., 2004; Wei et al., 2004; Boynton et al., 2011; Rastätter et al., 2013; Revallo et al., 2014; Uwamahoro & Habarulema, 2014; Eastwood et al., 2017; Lazzús et al., 2017; Wintoft et al., 2017; Chandorkar et al., 2017; Camporeale et al., 2018; Chandorkar & Camporeale, 2018; Podladchikova et al., 2018). These methods are mainly based on physical, empirical or analytical relationships between solar wind and geomagnetic parameters, correlation, and artificial neural networks (ANNs). Rastätter et al. (2013) compared various 30 model settings with observed Dst index and found empirical models perform well in general. However, physical models given proper boundary conditions observed to be performing better than empirical models. Linear regression, statistical correlation, etc. have been proved to be useful in understanding storm time geomagnetic field variations. There are many empirical models for Dst prediction. A simple prediction algorithm for the Dst index was proposed by Burton et al. (1975), solely from a knowledge of the solar wind parameters. They assumed a constant ring current recovery time constant (e-folding time) for all the storms which may not be always a valid assumption. Iyemori & Maeda (1980) first time successfully applied linear prediction filtering method for predicting geomagnetic activity using solar wind parameters.

Artificial neural networks are being extensively used in many areas where nonlinear complexities are involved (Lippmann, 1987; Miller, 1993; Gardner & Dorling, 1998; Unnikrishnan, 2014). In last few decades artificial neural networks are used for predicting geomagnetic activity at high and low geomagnetic latitude regions (Wu & Lundstedt, 1996; Gleisner & Lundstedt, 1997; Chandorkar et al., 2017; Andriyas & Andriyas, 2017). There are many studies which attempted to predict symmetric part of the ring current and geomagnetic field variations using neural networks (Kamide & Slavin, 1986; Lundstedt & Wintoft, 1994; Wu et al., 1998; Kugblenu et al., 1999; Unnikrishnan, 2012, 2014; Wintoft et al., 2017). Lundstedt & Wintoft (1994) developed feed forward neural network to predict geomagnetic activity index, Dst, one hour in advance. They could predict the initial and main phase very well but the recovery phase was not modeled correctly. The feedforward networks do not have feedback from the output or hidden nodes (refer Sect. 2 for network architecture) which constraint them in accurately modeling time series having memory. The performance of the network is seen to be independent of whether the raw parameters are used to the derived parameters like electric field and square root of dynamic pressure. Gleisner et al. (1996) have used the time-delayed feed-forward neural network for predicting Dst index. To predict Dst recovery phase they had to use up to 20 h of solar-wind input data. They could reproduce 84% of Dst variance. Further implementation of dynamic neural networks (i.e., feedback networks) has improved the prediction accuracy for the recovery phase. Elman recurrent neural networks are implemented by Wu & Lundstedt (1997) to predict Dst during geomagnetic storms. The recovery part of the storm has modeled significantly well by Elman network as it takes feedback from the hidden layer. The predictions are generally very good for the main phase of the geomagnetic storms but are fairly good for the recovery phases. The study by Revallo et al. (2015) developed a prediction model using the neural network and analytical prediction is done for two classes of geomagnetic storms caused by Coronal Mass Ejection (CME) and Corotating Interaction Region (CIR). They observed better prediction for CME driven storms than for CIR driven storm. Recently, Lazzús et al. (2017) used swarm-optimized neural network and showed that the hybrid approach gives good predictions of Dst index. Further, Chandorkar et al. (2017) presented the probabilistic approach in accurately predicting Dst index one step ahead. The use of NARX network by Cai et al. (2009) for predicting SYMH index (equivalent to high temporal resolution Dst index) has been very successful as it is observed to be better in performance due to feedback is given from the output node.

ASYH is a very valuable index to study the asymmetric development of magnetospheric storms during the crossing of interplanetary disturbances (e.g., Huttunen et al., 2006). There is a number of studies trying to understand the physical mechanism underlying the observed asymmetry in the ring current and its origin along with the contribution of various currents in asymmetric ring current (e.g., Liemohn et al., 2001; Jordanova et al., 2003; Newell & Gjerloev, 2012). Though ASYH index is equally important during storm time dynamics, there are no reports of ANN-based model available for ASYH index till date. The ANN-based prediction of asymmetric ring current will help to understand the contribution of external/internal drivers in the observed asymmetry. Also, the early forecast of ASYH will help space weather community to have a prior knowledge of the degree of asymmetry in the ring current. Therefore, developing an ANN-based prediction model for ASYH is the main objective of the present study. The present study develops a NARX network-based model to forecast both SYMH and ASYH indices, using interplanetary parameters as inputs and feedback from the output.

The predictions of Cai et al. (2009) for the studied storms using ACE data showed a correlation coefficient about 0.9 and RMSE of 14 nT. They showed that NARX is better in forecasting SYMH compared to Elman neural networks. The present study extends the use of the NARX network to predict ASYH index, which is the first time. Further, the present study extends the geomagnetic storm list by Cai et al. (2009) to 2015 which give a bigger dataset for the network training and testing. For this purpose, we have used interplanetary parameters, SYMH and ASYH indices during major geomagnetic storms occurred between 1998 and 2015 covering around two solar cycles. We present a prediction of SYMH and ASYH indices for the geomagnetic storm that took place on St. Patrick’s day, 2015 (intense geomagnetic storm of the current solar cycle, 24) along with few other major storms from solar cycle, 24.

The paper is arranged as follows: Section 2 describes the NARX neural network. Sections 3 and 4 introduce data and training methodology. Section 5 discusses the network performance and prediction of SYMH and ASYH indices during geomagnetic storms. Paper ends with the discussion and conclusions in Section 6.

2 NARX Neural Network

Artificial Neural Network (ANN) functions like biological neural network (Poulton, 2002). The biological neuron is composed of dendrites, the soma, and the axon. The neuron receives an input signal from other neurons which are connected to its dendrites by synapses. The soma is mainly processing unit where inputs are integrated over space and time and it activates an output depending on the total input. This output is transmitted by the axon and distributed to other neurons by the synapses at the tree structure at the end of the axon (Hérault & Jutten, 1994). The mathematical neuron functions a little simpler way since integration takes place only over space. The inputs are given at one or many nodes called input nodes. The sum of these weighted inputs is performed at the summing node which is fed to the nonlinear transform function or called an activation function to rescale the sum. The array of many nodes makes a network which can be made to learn relationships between inputs and targets, used for prediction.

For the present study, we have selected Nonlinear Auto Regressive with Exogenous inputs (NARX) model network due to its proven ability to account for the history of input and output parameters for prediction. This is feedback two-layer backpropagation network with time-delayed feedback. The basic network architecture is presented in Figure 1. As shown in the figure, inputs are shown to the network as a temporal sequence with different time lags with time delay length, d. Whereas, past outputs of the network are provided to the NARX network as feedback having history length, L, which are called as context inputs. From left to right network has an input layer, a hidden layer, and the output layer. The first layer or called as input layer receives externally provided values of input parameters. The second layer is known as the hidden layer because it does not see or act upon the external conditions. The hidden layer transforms the inputs such that the transformed inputs can be used by output layers. The output layer scales the hidden layer outputs to match the target. The dynamic behavior of the network can be formulated as:

$$O_t=F\left(O_{t-1},....O_{t-L};I_t,....I_{t-d}\right),$$(1)

where O is the output of the network, I denotes the input vector. L is the length of the feedback history and d is the length of input history. Thus, the output of the NARX network is a function of present inputs and their past values along with the history of the output. The inputs are processed by hidden nodes in the hidden layer, the output of jth hidden node is given by:

$$H_j=\tan h\left({\displaystyle \sum _{i=1}^M\mathrm{}{W}_{\mathrm{ji}}{ I}_i}+{\displaystyle \sum _{l=1}^L\mathrm{}{W}_{\mathrm{jl}}{C}_l+b_j}\right),$$(2)

where I_i is the value of input node i, M is a total number of input nodes. C denotes the context inputs. W_ji is a connecting weight of input node (i) and hidden node (j). Note that tanh (hyperbolic tangent) is the transfer function for nodes in the hidden layer. b_j is a bias of the jth neuron in the hidden layer. Complex and nonlinear relationships between inputs and output are taken care by tanh function. The output of the network (O(t)) is a linear summation over all hidden neuron outputs and output bias (b₀) which is represented by:

$$O\left(t\right)={\displaystyle \sum _{j=1}^s\mathrm{}{W}_{\mathrm{oj}}{H}_j+b_0,}$$(3)

here, W_oj is connecting the weight of hidden node to the output node and s is the number of hidden nodes.

Fig. 1 Architecture of Nonlinear Auto Regressive with eXogenous inputs (NARX) network. d is input history length and L is output feedback length.

3 Database

Different types of forecast models were studied prior to deciding the input and output database for NARX network. Interplanetary magnetic field, solar wind density, and velocity are most crucial parameters controlling the storm profile. Also, the history has a significant influence on the prediction accuracy. Therefore, we considered the total interplanetary magnetic field (B) and its components (B_y and B_z), solar wind density (N_sw) and solar wind speed (V_sw) as input parameters. SYMH and ASYH indices are considered as the target for the training two independent networks.

The study is carried out considering 67 major geomagnetic storms (minimum SYMH < −85 nT) during 1998–2005 (Storm list is adapted from: Cai et al., 2009) and 34 geomagnetic storms (minimum SYMH < −85 nT) identified between 2006 and 2015 (listed in Table 1). This period (1998–2015) covers 23rd and ongoing 24th solar cycles having total 101 geomagnetic storms of minimum SYMH < −85 nT. It also includes 17th March geomagnetic storm which is a major storm of the 24th solar cycle to date. Total 92 storms (1998–2013) are used for training and nine storms occurred during 2014–2015 are used to predict SYMH and ASYH indices. The utilized Solar wind parameters and indices were acquired from CDAWEB database (http://cdaweb.gsfc.nasa.gov/). One minute time resolution data was converted to 5-minute resolution for reducing the computation time. The missing data were interpolated using piecewise cubic Hermite polynomial. The total data length of ∼685 days having 5-minute resolution was used for developing the network.

The data of 92 geomagnetic storms between 1998 and 2013 is used for learning the network, which is divided into the three parts: (1) training (75%), (2) validation (15%) and (3) test (10%). As stated earlier the training data is used to learn the relationship between inputs and output. The validation of the network is determined through the identification of minimum error using 15% of the data. Validation data is used to stop the network from over-fitting the target. The test data was used to evaluate the performance of the network. Further, to check the prediction performance of the networks nine geomagnetic storms are used, which occurred during January, 2014–July, 2015.

4 Method and analysis

ANN-based prediction model consists of mainly three steps: training, validation and testing (Haykin & Network, 2004). The present network consists of one input layer with 30 external input nodes and 24 context inputs from the output, one hidden layer with 16 neurons and one output node. We have performed the trial and error on the number of hidden nodes and arrived at the optimum number of hidden nodes which was used in our network. The network performs better as the complexity of the time series is learned better with more number of hidden layers. However, it is to be noted that with an increased number of hidden layers, the computation time also increases largely. Therefore, to reduce the computation time we have used only one hidden layer. We have selected the input parameters (B, B_y, B_z, N_sw, and V_sw) based on the physical understanding of the solar wind-magnetosphere coupling. IMF, solar wind speed and density are primary parameters which affect reconnection and transfer of solar wind energy to the magnetosphere. It is reported by earlier researchers that ring current history of about 2 h is adequate for predicting SYMH index (Cai et al., 2009). Also, the communication time of the interplanetary electric field from the Bow-shock nose to the equatorial ionosphere is observed to be ∼20 min (Bhaskar & Vichare, 2013). Therefore, the input history of 30 min and output feedback length of 120 min having a 5-minute temporal resolution is used in the network. There are total 54 input nodes for each network. Input history of 30 min implies the need of total 30 external input nodes for five input parameters (B, B_y, B_z, N_sw, and V_sw) each parameter having six nodes. Similarly, 120-minute feedback length from output makes 24 context inputs.

The network is presented with the inputs to produce the desired output. For training the network, we have used a most popular back-propagation algorithm (Rumelhart et al., 1985) learning. In this algorithm the weights are updated by using delta rule which is given by:

$$\Delta w\left(i+1\right)=-\eta \frac{\mathrm{d}E}{\mathrm{d}w}+\alpha .\Delta w\left(i\right),$$(4)

here, w represents the weight of the nodes, i is epoch, a and n denote the momentum parameter and learning rate respectively. Momentum parameter is used to avoid local minimum, whereas learning rate controls the learning speed of the network. They range between 0 and 1. For optimization of the speed of learning the n is adjusted in each iteration according to the performance of the network. For initialization, small random values are assigned to the network weights. Initially a = 0.9 and n = 0.01 were considered for the networks training. E is network error which is estimated by using the following equation which is also known as a cost function,

$$E=\frac{1}{2}{\displaystyle \sum _{k=1}^N\mathrm{}{\left(O^k-T^k\right)}^2}.$$(5)

O is the output of the network, T is the target value and N is the total number of training samples. The error is minimized during epoch to epoch of the training for obtaining final trained network. To avoid the over-fitting the validation dataset is used. During the training, error on the validation data is continually monitored. The training is terminated when the validation error reaches a minimum and then increases for next six epochs consecutively. It is known and observed that initialization of network weights and number of nodes in the hidden layer affect the performance of neural networks. Hence, we trained the network multiple times by changing the initial weights and number of hidden layer nodes and selected the one which gave the best results for prediction.

Further, Root Mean Square Error (RMSE) was estimated to evaluate the performance of the network on test data consisting nine geomagnetic storms occurred during January, 2014–July, 2015 including the geomagnetic storm of March 17, 2015 (see Table 1). The root mean square error can be computed as,

$$\mathrm{RMSE}={\left[\frac{1}{N}{\displaystyle \sum _{i=1}^N\mathrm{}(O^i-T^i{)}^2}\right]}^{1/2}.$$(6)

Also, the cross-correlation coefficient (R) was estimated using equation (7) to quantify the similarities between the time series of the observed and predicted SYMH/ASYH index.

$$R=\frac{{\displaystyle \sum _{i=1}^N}\mathrm{}(T^i-\overline{T})(O^i-\overline{O})}{\sqrt{{\displaystyle \sum _{i=1}^N\mathrm{}(T^i-\overline{T}{)}^2}}\sqrt{{\displaystyle \sum _{i=1}^N\mathrm{}(O^i-\overline{O}{)}^2}}}$$(7)

5 Results

5.1 Network performance

Figure 2 shows the performance of the trained networks. Figure 2a shows the performance for SYMH whereas, Figure 2b shows performance for ASYH index. The figure presents the performance of all the steps i.e. training, validation and the test. As a part of the learning of the network, after each iteration, the mean squared error of both the networks initially decreases. This characteristic is observed in both the panels of Figure 2 i.e. initially the error in estimated SYMH and ASYH decreases with increasing epochs in a similar fashion and then remain steady during training, validation and testing. One more common feature observed during training, validation, and testing of SYMH and ASYH is that the errors of testing and validation converge to a smaller value compared to the training. The best test performance of the SYMH network is achieved at Mean Squared Error (MSE), ∼6 and that for ASYH network is ∼18. This implies the training of SYMH network is better compared to ASYH. Note that, to prevent the network from over-fitting, the training was stopped when validation error increased continuously for the next six iterations. This is achieved at epoch = 20 and 40 for SYMH and ASYH networks respectively.

Fig. 2 The training performance of the neural network is shown for both (a) SYMH and (b) ASYH indices.

Figure 3a, b shows the linear regression of targets (SYMH/ASYH) and predicted outputs of the networks for best training epochs. The correlation values are almost the same R ∼ .99 for both the networks output-target pairs. However, it is evident that the scatter is better for SYMH as compared to ASYH. The slope value close to unity and low value of intercept of the linear fit between target and output indicate that training is impressive for both the networks.

Fig. 3 The regression of the target and the modeled output by the networks are shown for (a) SYMH and (b) ASYH for training data.

5.2 Prediction

5.2.1 Geomagnetic storms −85 > SYMH > −210 nT

To test the prediction capability of the networks developed here, we have used the geomagnetic storms which were not considered in the training process. The details of these major storms used for prediction are presented in Table 2. Figure 4 shows the predicted (blue line) and observed (dashed red line) values of SYMH for first 8 geomagnetic storms listed in the table. The storm on March 17, 2015, is discussed in detail in the next subsection. In general, all the storms show a very good match between the predicted and observed SYMH profiles. One can notice that even finer features of timescales 10–30 min are reflected in the predicted profiles. The prediction of the minimum SYMH matches very well with the observed strength of the storms occurred on Feb 27, 2014; Jan 07, 2015; Jun 21, 2015, and Jul 04, 2015. However, the minimum SYMH of storms occurred on Apr 11, 2014, and Jun 07, 2014 are underestimated by the network. Transient variations like storm sudden commencement (SSC) are reproduced well by the network during storms on Jan 07, 2015 and Jul 04, 2015. Nevertheless, SSC occurred on June 21, 2015, and Feb 27, 2014, is not predicted by the model. Also, a two-step decrease observed during the main phase of storms (Jun 07, 2015, Jun 21, 2015, and Jul 04, 2015) is reproduced by the network to the good extent.

Fig. 4 Predicted and observed SYMH index for test storms listed in Table 2. The number indicates the geomagnetic storm number listed in the table. Dotted red curve is observed SYMH and solid blue curve is the predicted SYMH by the network.

Figure 5 shows the predicted (blue line) and observed (dashed red line) values of ASYH index during the storms. The predicted profiles of ASYH match well with the observed profiles. In general, the prediction model underestimates the amplitude of ASYH index. However, the overall temporal profile of ASYH index is well predicted by the network. Note that, finer structures of smaller timescales ∼10–30 min are also well mimicked by the model predictions except for Apr 11, 2014 storm.

Fig. 5 Predicted and observed ASYH index for test storms listed in Table 2. The number indicates the geomagnetic storm number listed in the table. Dotted red curve is observed ASYH and solid blue curve is the predicted ASYH by the network.

Tsurutani & Gonzalez (2013) noted a variability of around −30 nT in Dst index, which they considered as a threshold/noise level for geomagnetic storms (Munsami, 2000). Therefore, here we have considered ±30 nT as the threshold even for SYMH/ASYH, although it is possible that for higher time resolution indices the noise level might be larger. The residuals of SYMH and ASYH are estimated by subtracting model values from the observed, which are presented in Figure 6 for the selected storms. The noise levels (±30 nT) are marked by the red dashed lines and the main phase bounded by two vertical dashed lines, in each panel. The figure shows that for most of the storms the residuals lie within the noise level, in general. However, for the storm occurred on June 21, 2015, the residuals are well above noise level (>100 nT), in particular during the main phase. Note that the model estimates are based on the interplanetary (external) inputs. Therefore, the higher residual values could be ascribed to the influence of the magnetospheric origin (internal), which is not modeled by the present network. Also, one can notice that compared to the residuals in SYMH index, the ASYH residuals are higher in magnitude. This may imply the larger contribution in the ASYH index due to magnetospheric sources, compared to that in SYMH index.

Fig. 6 The residual of predicted and observed values for SYMH and ASYH indices during test storms listed in Table 2. Dotted horizontal dashed lines mark ±30 nT noise level (i.e., quiet time background). The vertical dashed lines from left represent onset and the end of main phase of the storm respectively.

Further, to quantify how good are the networks in predicting these indices, the correlation coefficient (R) and RMSE are estimated for the predicted geomagnetic storms which are listed in Table 2. A good performance of the networks is more evident from the observed high mean correlation coefficients, R ∼ 0.9 and R ∼ 0.7 for SYMH and ASYH indices respectively (see Table 2). The table clearly shows the cross-correlation is high between predicted and observed SYMH as compared to the correlation between predicted and observed ASYH. This smaller value of the correlation coefficient of ASYH index could be due to the presence of very high-frequency fluctuations in ASYH index compared to SYMH index (refer Figs. 4 and 5). Also, as discussed earlier, ASYH is more complex in nature due to various currents, unlike SYMH. The developed forecast model is compared with the base model. Here the base model is assumed as persistence model one hour ahead (i.e., the index at time t is predicted to be equal to itself at time t−1 h. The estimated correlation coefficient and RMSE for the predicted geomagnetic storms are also shown in the table. The comparison of base model and NARX mdoel shows that later performs better.

5.2.2 St. Patrick’s day 2015 geomagnetic storm (SYMH ∼ −234 nT)

Recent storm of 17th March, 2015 (known as St. Patrick’s day 2015 geomagnetic storm) is a great storm of the ongoing 24th Solar cycle (see Table 2). The parameters of interplanetary disturbance during this storm are presented in Figure 7. Interplanetary magnetic field (B, B_y, B_z, nT), Solar wind density (N_sw) and Velocity (V_sw) show the clear enhancement at the onset of the storm on March 17, 2015, time ∼ 05: 00 UT. It is evident that long duration (∼24 h) southward B_z with amplitude ∼20 nT has given rise to an intense storm having minimum SYMH ∼ −234 nT and maximum ASYH of ∼250 nT. The storm shows two distinct steps during the main phase, which might be associated with the sudden north-south turnings of the IMF occurred on March 17, 2015, at ∼11 UT. ASYH index is enhanced during the main and early recovery phase of the storm. The ASYH index is almost of the same magnitude that of SYMH during the main phase. The prediction results obtained from the present models along with the observed indices are displayed in Figure 8. The trained NARX network for SYMH predicts the observed SYMH very well including the sudden storm commencement, two steps in the main phase and small transient fluctuations during the recovery phase of the storm. The predicted storm time minimum SYMH is close to the observed one. The correlation between predicted and observed SYMH is very high (R ∼ 0.93) and RMSE value is low (∼21 nT). The match between predicted and observed ASYH is excellent during the main phase. For ASYH, correlation between predicted and observed ASYH is R ∼ 0.7 and RMSE value is ∼20 nT.

Fig. 7 Interplanetary parameters of intense geomagnetic storm of March 17, 2015: Interplanetary magnetic field (B) and its components By, Bz; solar wind speed (Vsw) and density (Nsw). Geomagnetic indices: ASYH and SYMH. The dashed vertical line marks the onset of the storm.

Fig. 8 Predicted and observed (a) SYMH and (b) ASYH indices for March 17, 2015 geomagnetic storm.

6 Discussion and conclusions

As ASYH index is of paramount importance to unravel the information about the asymmetric response of the magnetosphere especially during geomagnetic storms, the present study attempts to predict ASYH index for the first time. We have applied NARX neural network to 92 geomagnetic storms occurred between 1998 and 2013. The developed networks successfully predict SYMH and ASYH indices about an hour prior to the onset of storm provided the real-time upstream solar wind data is available. The temporal horizon of the forecast is the propagation time of the solar wind from the spacecraft ACE/WIND to the Earths bow shock which is about 1 h. The propagation time is not constant and varies based on the solar wind speed and location of the spacecraft. However, time estimates can be done based on the know information of the spacecraft location and solar wind speed. To implement the developed network the real-time solar wind data will be time-shifted to the bow-shock nose following the procedure adopted for generating OMNI data (https://omniweb.gsfc.nasa.gov/html/ow_data.html). To demonstrate the feasibility of this approach we have time-shifted ACE solar wind and magnetic field data to bow shock for the forecasted events. The time-shifted and OMNI interplanetary magnetic field and solar wind data comparison is shown in Figure 9. There are various methods being used to time-shift spacecraft data to the Earth and each have their own merits and limitations (see for details Cameron & Jackel, 2016 and references therein). However, to keep it consistent with training data of the network, we have followed the original method used by Coordinated Data Analysis Web (CDAWeb) to time-shift solar wind data to the Earths bow shock. The time-shift is calculated as follows,

$$\Delta t=\frac{X}{V_{\mathrm{sw}}}\times \frac{[1+(Y\times W)/X]}{[1-V_{\mathrm{e}}\times W/V_{\mathrm{sw}}]},$$(8)

where, Δt is the time shift in seconds, X and Y are GSE X and Y components of the spacecraft position vector, in km, V_sw is the observed solar wind speed in km/s, V_e is the speed of the Earth’s orbital motion (30 km/s), W = tan [0.5 × atan (V_sw/428)] is parameter related to the assumed orientation of the phase front relative to the Earth-sun line. Figure 9 clearly shows that almost same OMNI interplanetary magnetic field and solar wind data could be generated applying this time-shift. Note that solar wind density (N_sw) data is missing for some initial period in ACE but present in OMNI. This is due to the fact that OMNI is comprised of ACE and WIND data. The parameters estimated at bow shock at (time = t + Δt) are used as inputs to ANN and predict the SymH and AsyH indices at time t + Δt. Thus, using the real time data at L1 point, we predict the SymH/AsyH indices after time shift of Δt. This assures that the trained network using OMNI data can be implemented for real-time forecasting provided the data is time-shifted to the bow shock using the same method as CDAWEB used to generate OMNI data.

Fig. 9 Time-shifted ACE data and OMNI data comparison for July 04, 2015 geomagnetic storm.

The need of 30 min input history and 120-minute feedback for better predictions imply the role of preconditioning of magnetosphere i.e. the future state of contributing currents in indices depends on their present and past values. We have examined the prediction for nine geomagnetic storms from solar cycle 24, occurred during January, 2014–July, 2015. These storms include the major storm occurred on St. Patrick’s day, 2015, which is the most intense storm occurred so far in solar cycle 24.

The ability of NARX having feedback from output enabled us to model these indices quite accurately. The temporal variations of the order of 10–30 min are well predicted by both the networks. The network trained for SYMH index predicts SYMH very well and observed average correlation between predicted and observe SYMH is high (R ∼ 0.88) i.e. almost ∼77% variations of SYMH are modeled by the network. The average RMSE is about 13.98 nT and matches with the observations by Cai et al. (2009). Therefore, the prediction performance of the network is almost the same as that of ANN constructed by Cai et al. (2009). However, as noted earlier the prediction accuracy varies from storm to storm. Munsami (2000) also observed mismatch between predicted and observed Dst index, which they thought to be due to other than external drivers such as substorms. However, they did not observe noticeable improvement in Dst prediction even by considering inputs from substorm activity. This could be due to the contribution to Dst from other processes such as wave-particle interaction, charge exchange, the ionospheric outflow of O⁺ ions, particle loss to the atmosphere and magnetopause (Daglis et al., 1999; Liemohn et al., 2001). For the present study, the residuals between observed and predicted values of SYMH generally lie within the noise level of ±30 nT. However, note that sometimes high residual values are observed above the noise level, especially during the main phase of geomagnetic storms. For SYMH the prediction results are almost the same as that Cai et al. (2009). However, more dataset is used in the present study which helped in the generalization of the network.

As an extension of the earlier studies we have extended NARX network to ASYH index. In general, the prediction of ASYH index is very good, within the noise level of ±30 nT. More than 50% variations of ASYH are explained by the present network. This implies that the variation in asymmetric ring current could be explained by solar wind parameters. However, during the main and early recovery phase of storms, the residuals (observed-modeled) are above the noise level, which could be ascribed to the internal magnetospheric processes such as field-aligned currents, particle loss. The study shows that the prediction is better for SYMH as compared to ASYH. This implies that there is something more than external solar wind driver i.e. internal dynamics which is contributing in ASYH index. If both the variations are caused by the same current system then prediction of ASYH and SYMH would have been very much similar at least in the main phase of the storms. However, this is not the case and hence one may conclude that SYMH and ASYH variations during the main phase are not necessarily related to the same current system.

The present study demonstrates that developed networks are capable of predicting SYMH and ASYH indices and hence can be implemented for the real-time forecasting. The network performs better as compared to the assumed base model. This is the first attempt to forecast ASYH index using the neural network with the use of a long dataset of geomagnetic storms spanning more than one solar cycle. Interestingly, even though ASYH is a good proxy for the internal variability of the asymmetric ring current ANN could model large part of the variations using external (solar wind) parameters. The reliable forecast of SYMH and ASYH indices will help space weather community and space programs to get early information on the strength of geomagnetic disturbances and their asymmetric geomagnetic response. The ASYH index represents the magnitude of the asymmetric geomagnetic response during geomagnetic disturbances. The main contribution comes from the dawn-dusk asymmetry. The forecast of ASYM will enable to understand how strong asymmetry expected in the magnetosphere prior to disturbance hits the magnetosphere. SYMH forecast will only enable to understand the global average response of the magnetosphere, addition of this index will help to get a complete picture of the geomagnetic response during geomagnetic storms. In future, the predictions may be improved by considering inputs representative of internal magnetospheric dynamics.

Acknowledgments

The solar wind parameters, interplanetary magnetic field, and geomagnetic indices used in this paper are obtained from CDAWEB (http://cdaweb.gsfc.nasa.gov/). We thank the ACE Science Center for making data available in the public domain. Authors extend thanks to Indian Institute of Geomagnetism which supported the work. A.B. is presently supported by the NASA Living with a Star Jack Eddy Postdoctoral Fellowship Program, administered by UCARs Cooperative Programs for the Advancement of Earth System Science (CPAESS). The editor and authors thank K. Unnikrishnan and an anonymous referee for their assistance in evaluating this paper. Authors also thank Jean Lilensten, Enrico Camporeale, and Richard Boynton for their valuable comments which helped authors to improve the manuscript.

References

All Tables

All Figures

	Fig. 1 Architecture of Nonlinear Auto Regressive with eXogenous inputs (NARX) network. d is input history length and L is output feedback length.
In the text

	Fig. 2 The training performance of the neural network is shown for both (a) SYMH and (b) ASYH indices.
In the text

	Fig. 3 The regression of the target and the modeled output by the networks are shown for (a) SYMH and (b) ASYH for training data.
In the text

	Fig. 4 Predicted and observed SYMH index for test storms listed in Table 2. The number indicates the geomagnetic storm number listed in the table. Dotted red curve is observed SYMH and solid blue curve is the predicted SYMH by the network.
In the text

	Fig. 5 Predicted and observed ASYH index for test storms listed in Table 2. The number indicates the geomagnetic storm number listed in the table. Dotted red curve is observed ASYH and solid blue curve is the predicted ASYH by the network.
In the text

	Fig. 6 The residual of predicted and observed values for SYMH and ASYH indices during test storms listed in Table 2. Dotted horizontal dashed lines mark ±30 nT noise level (i.e., quiet time background). The vertical dashed lines from left represent onset and the end of main phase of the storm respectively.
In the text

	Fig. 7 Interplanetary parameters of intense geomagnetic storm of March 17, 2015: Interplanetary magnetic field (B) and its components By, Bz; solar wind speed (Vsw) and density (Nsw). Geomagnetic indices: ASYH and SYMH. The dashed vertical line marks the onset of the storm.
In the text

	Fig. 8 Predicted and observed (a) SYMH and (b) ASYH indices for March 17, 2015 geomagnetic storm.
In the text

	Fig. 9 Time-shifted ACE data and OMNI data comparison for July 04, 2015 geomagnetic storm.
In the text