the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Probabilistic modelling of substorm occurrences with an echo state network
Ryuho Kataoka
Masahito Nosé
Jesper W. Gjerloev
The relationship between solarwind conditions and substorm activity is modelled with an approach based on an echo state network. Substorms are a fundamental physical phenomenon in the magnetosphere–ionosphere system, but the deterministic prediction of substorm onset is very difficult because the physical processes that underlie substorm occurrences are complex. To model the relationship between substorm activity and solarwind conditions, we treat substorm onset as a stochastic phenomenon and represent the stochastic occurrences of substorms with a nonstationary Poisson process. The occurrence rate of substorms is then described with an echo state network model. We apply this approach to two kinds of substorm onset proxies. One is a sequence of substorm onsets identified from auroral electrojet intensity, and the other is onset events identified from activity of Pi2 pulsations, which are irregular geomagnetic oscillations often associated with substorm onsets. We then analyse the response of substorm activity to solarwind conditions by feeding synthetic solarwind data into the echo state network. The results indicate that the effect of the solarwind speed is important, especially for Pi2 substorms. A Pi2 pulsation can often occur even if the interplanetary magnetic field (IMF) is northward, while the activity of auroral electrojets is depressed during northward IMF conditions. We also observe spiky enhancements in the occurrence rate of substorms when the solarwind density abruptly increases, which might suggest an external triggering due to a sudden impulse of solarwind dynamic pressure. It seems that northward turning of the IMF also contributes to substorm occurrences, though the effect is likely to be minor.
 Article
(2114 KB)  Fulltext XML
 BibTeX
 EndNote
Substorms are a fundamental physical phenomenon in the magnetosphere–ionosphere system. Since substorms are the main source of geomagnetic disturbances in the polar ionosphere causing geomagnetically induced currents (e.g. Viljanen et al., 2006; Wei et al., 2021; Schillings et al., 2022), prediction of substorms is an important issue. Substorms are most likely driven by the solar wind, so it is essential to understand the relationship between substorms and solarwind conditions for predicting substorms. In this context, many studies have attempted to construct predictive models of the auroral electrojet indices, AU and AL, which represent the intensities of auroral electrojets in the polar ionosphere (Davis and Sugiura, 1966). For example, Luo et al. (2013) developed a parametric model for predicting the AU and AL indices from a time history of solarwind data. There have also been a number of studies which employed machine learning approaches for predicting the AU and AL indices from given solarwind conditions (e.g. Gleisner and Lundstedy, 1997; Takalo and Timonen, 1997; Amariutei and Ganushkina, 2012). Our previous study (Nakano and Kataoka, 2022) (hereinafter referred to as NK22) also modelled the relationship between the AU and AL indices and solarwind variables using an echo state network (ESN) model, which is a kind of recurrent neural network originally introduced by Jaeger and Haas (2004).
However, a good prediction of the AU and AL indices does not necessarily guarantee successful prediction of substorm onsets. Many studies have argued that substorm onsets are triggered by internal processes of the magnetosphere–ionosphere system (e.g. Baker et al., 1996; Lui, 1996; Lyons et al., 2018; Miyashita and Ieda, 2018). Considering that substorms sometimes take place without any visible solarwind variations (e.g. Nakano and Iyemori, 2005), substorms may be caused by complex magnetospheric processes which are hard to predict. It would thus be difficult to deterministically predict individual substorm onsets from a given time history of solarwind conditions. As a matter of fact, the ESN in our previous study, NK22, did not always predict sharp decreases of the AL index due to a substorm expansion. Figure 1 compares the observed Kyoto AU and AL indices (World Data Center for Geomagnetism, Kyoto et al., 2015) and the prediction with the ESN. Although the prediction shown by the red line roughly traces the actual AU and AL indices, it did not reproduce many negative AL spikes.
In this study, we propose another approach for modelling the relationship between substorm occurrences and solarwind input. To deal with the complexity underlying substorm occurrences, we treat a substorm onset as a stochastic phenomenon. The stochastic occurrence of substorms is represented with a nonstationary Poisson process, which is a probabilistic model describing event time series (e.g. Daley and VereJones, 2003). The occurrence rate of substorms is then modelled with an echo state network (ESN) model. By introducing the ESN, we can represent the dependence of the occurrence rate on solarwind conditions. Probabilistic prediction of substorm onsets was also performed by Maimaiti et al. (2019), who predicted a substorm occurrence for next 60 min from the time history of solarwind data. In contrast, the purpose of this study is to model the response of substorm activity to given solarwind condition. Our approach thus sequentially processes the time series of solarwind data to obtain the instantaneous occurrence rate of onsets.
We applied the proposed ESNbased approach to two types of substorm onsets. One is a sequence of substorm onsets identified from the SML index (Newell and Gjerloev, 2011a, b). The SML index is a proxy of westward auroral electrojet intensity derived from the SuperMAG geomagnetic data (Gjerloev, 2012) with the same algorithm as for the AL index. However, since the events determined from the auroral electrojet intensity may contain nonsubstorm events such as DP2type convection enhancements (e.g. Nishida, 1968; Kamide and Kokubun, 1996), an increase of the auroral electrojet would not necessarily indicate a substorm onset. To identify substorm onsets without relying on the auroral electrojet intensity, we also analyse a sequence of substorm onsets identified from the Wp index (Nosé et al., 2009, 2012; World Data Center for Geomagnetism, Kyoto and Nose, 2016). The Wp index corresponds to the amplitude of lowlatitude Pi2 pulsations. An onset determined from the Wp index is thus regarded as a Pi2 event. After training the ESN to model the two types of substorms, the relationship between substorm activity and solarwind conditions is discussed by analysing the response of the ESN outputs to synthetic solarwind inputs.
We denote ν(tβ)dt as the probability of the occurrence of an event within a small time interval dt, where β is a parameter determining the shape of the function ν. The function ν is referred to as the intensity function, and it corresponds to the instantaneous occurrence rate per unit time. Given a sequence of event occurrence times ${\mathit{\tau}}_{\mathrm{1}:N}=\mathit{\{}{\mathit{\tau}}_{\mathrm{1}},{\mathit{\tau}}_{\mathrm{2}},\mathrm{\dots},{\mathit{\tau}}_{N}\mathit{\}}$, the likelihood of the parameter β is written as follows (e.g. Daley and VereJones, 2003):
where p denotes probability density, t denotes time, and N is the number of events. The loglikelihood thus becomes
Discretising Eq. (2) in time, we obtain
where ${t}_{k}={t}_{\mathrm{0}}+k\mathrm{\Delta}t$. The time interval Δt in Eq. (3) is set to be 5 min in this study. If an event occurs at t=τ_{i} during large ν, the first term on the righthand side of Eq. (2) gets larger. On the other hand, the second term on the righthand side of Eq. (2) becomes a large negative number if ν keeps large values for the whole interval. The loglikelihood, log L, thus becomes larger by choosing an intensity function which becomes large when events frequently occur and small when events do not occur. We want to obtain a good intensity function which describes the occurrence rate of events better.
We model the function ν with the ESN (Jaeger and Haas, 2004), used in NK22 and another recent study (Kataoka and Nakano, 2021). Denoting the vector consisting of the state variables of the ESN at time t_{k} as x_{k}, the ith element of x_{k}, x_{k,i} is updated at each time step according to the following equation:
where z_{k} is the vector of the input variables, w_{i} is the vector of the weights connecting the state variables, u_{i} is the vector of the weights connecting the input variables and state variables, and m denotes the dimension of the vector x_{k}. The dimension m is set to be 1000 in this study. The input vector z_{k} is given as follows:
where B_{z,k}, B_{y,k}, and B_{x,k} respectively denote the z, y, and x components of the interplanetary magnetic field in the geocentric solar magnetospheric (GSM) coordinates at time t_{k}; V_{sw,k} is the x component of the solarwind velocity in the GSM coordinates; and N_{sw,k} and Θ_{sw,k} are the solarwind density and temperature, respectively. These solarwind variables are taken from the OMNI 5 min data (King and Papitashvili, 2023). H_{k} and D_{k} respectively indicate universal time (UT) in hours and the day from the end of 2000 (i.e. D_{k}=1 on 1 January 2001) for determining the UT dependence and seasonal dependence (e.g. Cliver et al., 2000). ${S}_{{B}_{x}}$, ${S}_{{B}_{y}}$, ${S}_{{B}_{z}}$, S_{V}, S_{N}, and S_{T} are rescaling factors to adjust the value of each element of z_{k} to a similar range, and b_{V}, b_{N}, and b_{T} are for adjusting the range of each element of z_{k}. According to NK22, we assume ${S}_{{B}_{x}}={S}_{{B}_{y}}={S}_{{B}_{z}}=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{nT}$, ${S}_{V}=\mathrm{500}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, S_{N}=20 cm^{−3}, S_{T}=10^{6} K, ${b}_{V}=\mathrm{400}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, b_{N}=1 cm^{−3}, and ${b}_{T}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}\mathrm{K}$. The weights {w_{i}} and {u_{i}} are randomly given in advance and are fixed. The parameters {η_{i}} are also randomly given and are fixed. Specifically, we randomly chose 90 % of the weights {w_{i}} and {u_{i}} and set them to 0. The nonzero elements of u_{i} and η_{i} were given randomly from a standard normal distribution. The nonzero elements of w_{i} were also drawn from a normal distribution. The weights {w_{i}} were then rescaled such that the maximum singular value of the weight matrix is 0.99, where the weight matrix W is an m×m matrix defined as
By setting the maximum singular value of W to be less than unity, it is guaranteed that the ESN forgets distant past inputs and that the weights can be stably determined. Although NK22 employed a leaky echo state network, this study uses a network without a leak term, because the fitting of the substorm occurrence probability was slightly degraded when the leak term was introduced.
We represent the function ν as an exponential of the weighted sum of the state variables as
Here, we have used the parameter vector β for the weights for determining the output of the ESN. The value of β is obtained with a Bayesian approach. We take the prior distribution of β as a Gaussian distribution as
where m denotes the dimension of x_{k} and is 1000. We determine the standard deviation, σ, based on the marginal likelihood (see Appendix A) to avoid overfitting and underfitting. We set σ=0.6 for analysing the SML substorm onsets in Sect. 3 and σ=1.4 for analysing the Wp onsets in Sect. 4. We estimate β such that the following posterior probability density is maximised:
As p(τ_{1:N}) does not depend on β, we can obtain the optimal β by maximising the following objective function:
using the Newton–Raphson method.
First, we analysed the SuperMAG substorm list by Newell and Gjerloev (2011a) derived from the SuperMAG data (Gjerloev, 2012). This SuperMAG substorm list identifies substorm onset times from the SML index, which indicates westward auroral electrojet intensity. The SuperMAG list thus enumerates the events in which a westward auroral electrojet was developed. We refer to such events as auroral electrojet substorms (AE substorms). We trained the ESN with data for the 10 years from 2005 to 2014 so that the ESN output represents well the occurrence rate of the substorm onsets in the SuperMAG list. We used 5 min values of the OMNI solarwind data as the input. As in NK22, we start the comparison after spinup of the ESN for 72 steps so that the ESN satisfactorily memorises the history of the input data. If more than half of the data were missing for 1 h, we stopped the prediction and spun up the ESN again for the subsequent 72 steps. We then predicted the intensity function ν which represents instantaneous occurrence rate of substorm onsets.
Figure 2 shows an example of the prediction of the occurrence rate for the 3 d from 6 to 8 June 2015. In the top panel, the red line indicates the occurrence rate ν (h^{−1}) obtained with the ESN. The AE substorm onset times in the SuperMAG list are also plotted with blue triangles in this panel. For comparison, the SMU and SML indices (Newell and Gjerloev, 2011a, b), interplanetary magnetic field (IMF), solarwind speed, and solarwind density are shown in the second, third, fourth, and bottom panels, respectively. The IMF shown in the third panel is expressed in GSM coordinates, and the x, y, and z components are indicated with the green, blue, and red lines, respectively. The predicted occurrence rate tended to increase when substorm onsets, indicated with blue triangles, were frequently observed. This means that the ESN predicted the occurrence pattern of the AE substorms well.
To assess the performance of the ESNbased prediction, we calculate the predicted probability of the substorm occurrence for each hour, obtained as
Note that 12Δt corresponds to 1 h because Δt was taken to be 5 min. We classified the data into 10 classes by the predicted probability ($\mathrm{0}\le {P}_{k}<\mathrm{0.1}$, $\mathrm{0.1}\le {P}_{k}<\mathrm{0.2}$, …, and 0.9≤P_{k}) and calculated the actual occurrence ratio for each class. Figure 3 shows the actual occurrence ratio with respect to the predicted probability for AE substorms in the 4 years from 2015 to 2018. For reference, the gray line indicates the case where the occurrence ratio is equal to the predicted probability. The result with the red line shows that the predicted probability obtained with the ESN agrees well with the actual occurrence ratio except that the prediction slightly overestimates the occurrence ratio for $\mathrm{0.4}\le {P}_{k}<\mathrm{0.8}$.
The performance of a probabilistic prediction is often evaluated with the Brier score, defined as
where r_{k} denotes the actual result. If the substorm occurred during the period ${t}_{k}\le t<{t}_{k+\mathrm{1}}$, r_{k}=1; otherwise, r_{k}=0. For example, with a noninformative prediction when the prediction P_{k} is always 0.5, $({r}_{k}{P}_{k}{)}^{\mathrm{2}}=\mathrm{0.25}$ for the entire period, and the Brier score becomes B=0.25. The Brier score gets smaller as the prediction improves. For the data from 2015 to 2018, the Brier score of the prediction with the ESN was 0.100. This Brier score should be compared with the case where we assume a stationary Poisson process where the occurrence rate ν is constant. If a stationary Poisson process is assumed and ν is optimised to fit the same data as used for training the ESN, the probability of substorm occurrence per hour is 15.6 %. Using this value, the Brier score with the stationary Poisson process becomes 0.136, which is worse than the score with the ESN. The better score with the ESN confirms that the information on the solar wind, which is used as the input for the ESN, effectively improves the prediction about the substorm occurrence.
In the following, we conduct a prediction of Pi2 substorms with the same approach as used in the previous section. This study identifies the Pi2 substorm onset using the Wp index (Nosé et al., 2009, 2012; World Data Center for Geomagnetism, Kyoto and Nose, 2016). The Wp index corresponds to the amplitude of lowlatitude Pi2 pulsations averaged over the nightside longitudes, which is derived based on the wavelet analysis (Nosé et al., 1998). Nosé et al. (2012) proposed criteria for identifying Pi2 onsets from the Wp index. However, one of the criteria assumes quiescence of geomagnetic activity before the onset. To take into account Pi2 onsets during disturbed time, this study uses a different method for identifying Pi2 onsets. Although the original Wp index has a 1 min resolution, our identification recipe uses moving averages with window sizes of 5 min to remove the effects of noisy oscillations. We detect a time of a peak of the moving averages, T, which satisfies the following:
where $\stackrel{\mathrm{\u203e}}{\mathrm{Wp}}\left(T\right)$ denotes the moving average of the Wp index at time T. We then find the time of the minimum $\stackrel{\mathrm{\u203e}}{\mathrm{Wp}}$ within the 15 min interval before the time of the peak and define it as T−m min ($\mathrm{1}\le m\le \mathrm{15}$). If $\stackrel{\mathrm{\u203e}}{\mathrm{Wp}}\left(T\right)\stackrel{\mathrm{\u203e}}{\mathrm{Wp}}(Tm\phantom{\rule{0.125em}{0ex}}\mathrm{min})>\mathit{\theta}\phantom{\rule{0.125em}{0ex}}\mathrm{nT}$, we identify this increase of $\stackrel{\mathrm{\u203e}}{\mathrm{Wp}}$ as a Pi2 event. To determine the onset time, we calculate the difference of $\stackrel{\mathrm{\u203e}}{\mathrm{Wp}}$ between 5 min apart:
for $j=\mathrm{2},\mathrm{\dots},m$. Finding $j={j}^{\ast}$ which maximises $\mathrm{\Delta}\stackrel{\mathrm{\u203e}}{\mathrm{Wp}}(Tj\phantom{\rule{0.125em}{0ex}}\mathrm{min})$, the time $T{j}^{\ast}\phantom{\rule{0.125em}{0ex}}\mathrm{min}$ is regarded as the onset time. We hereinafter refer to the Pi2 events identified with the above criteria as Pi2 substorms, although they may contain pseudobreakups and other phenomena which are not normally classified as substorms. In our event detection method, the threshold θ is a tunable parameter which determines the sensitivity of the detection. We set θ=0.1 nT and considered two other cases: θ=0.15 nT and θ=0.2 nT for comparison, as described later.
We trained the ESN with data for the 10 years from 2005 to 2014. Again, we used 5 min values of the OMNI solarwind data as the input. We started the comparison after spinup of the ESN for 72 steps. If more than half of the data were missing for 1 h, we stopped the prediction and spun up the ESN again for the subsequent 72 steps. We then predicted the occurrence rate ν for substorm onsets. Figure 4 shows the prediction of the occurrence rate of Pi2 substorms for the 3 d from 6 to 8 June 2015, which is the same event as in Fig. 2. When identifying the substorm onsets, the threshold θ was set to be 0.1 here. In the top panel, the red line indicates the occurrence rate (h^{−1}) estimated with the ESN. The second panel shows the Wp index, and the third panel shows the SMU and SML indices. The fourth, fifth, and bottom panels show the IMF in GSM coordinates, solarwind speed, and solarwind density, respectively. Similarly to Fig. 2, the predicted occurrence rate tended to increase when Pi2 substorm onsets, indicated with blue triangles, were frequently observed.
We also calculate the predicted probability of substorm occurrence for each hour defined in Eq. (11). As in the previous section, we classified the data into 10 classes by the predicted probability ($\mathrm{0}\le {P}_{k}<\mathrm{0.1}$, …, 0.9≤P_{k}) and calculated the actual occurrence ratio for each class. Figure 5 shows the actual occurrence ratio with respect to the predicted probability for Pi2 substorms in the 4 years from 2015 to 2018. The predicted probability obtained with the ESN mostly agrees with the actual occurrence ratio although slightly underestimated. We also calculated the Brier score in Eq. (12). The Brier score values were 0.208, 0.180, and 0.150 with θ=0.1, θ=0.15, and θ=0.2, respectively. If a stationary Poisson process is assumed, the Brier score values were 0.236, 0.200, and 0.165 with θ=0.1, θ=0.15, and θ=0.2. The Brier score gets smaller as the threshold θ increases due to the number of events. If a larger threshold is taken, fewer events are identified as substorms. This situation tends to depress the predicted occurrence probability P_{k}. At the same time, as the number of events decreases, r_{k} in Eq. (12) takes a value of 0 more frequently and (r_{k}−P_{k})^{2} accordingly decreases in more cases. Thus the Brier score depends on the number of events.
Comparing with the case of AE substorm onsets, the difference in the Brier score between the nonstationary and stationary Poisson process is likely to be smaller for Pi2 substorms. This presumably indicates that predicting Pi2 substorms is more difficult than AE substorms. Figure 6 is a histogram of the frequency of P_{k} for Pi2 substorms identified with θ=0.1, Pi2 substorms identified with θ=0.2, and AE substorms. For AE substorms, P_{k} was below 0.1 in most cases. In contrast, it was relatively rare that P_{k} was lower than 0.1 for Pi2 substorms with θ=0.1. If P_{k} approaches to 1, we can confidently predict the occurrence of a Pi2 substorm. If P_{k} approaches to 0, we can confidently deny the occurrence of a Pi2 substorm. However, P_{k} for Pi2 substorms was between 0.2 and 0.6 for many cases, which means that it is usually difficult to confidently predict whether a Pi2 substorm will occur or not. If the threshold is taken to be θ=0.2, the frequency of P_{k}<0.1 increased but it was still much less than the case of AE substorms. Moreover, the frequency of P_{k}≥0.7 for Pi2 substorms with θ=0.2 is less than that for AE substorms. Overall, a Pi2 substorm onset seems to be less predictable than an AE substorm onset.
To determine what the ESN learnt from the data, we conducted an experiment analysing the response of the trained ESN to synthetic solarwind data. We used similar synthetic solarwind data to NK22 in which the solarwind parameters are fixed at constant values except that one of the parameters varying as a rectangular wave of various periods. Figure 7 demonstrates the experiment with the synthetic solarwind data over 18 d. The top panel shows the occurrence rate of AE substorms predicted with the ESN. The second panel shows the occurrence rate of Pi2 substorms predicted with the ESN. The third, fourth, and fifth panels display the IMF, solarwind speed, and solarwind density, respectively. In this experiment, the IMF B_{x} and B_{y} were set to 0 and the temperature was fixed at 5×10^{5} K. In the first 3 d, IMF B_{z} varied as a rectangular wave between −5 and 1 nT with a period of 20 min for the first day, 2 h for the second day, and 6 h for the third day, while the solarwind speed was fixed at 400 km s^{−1} and the density was fixed at 2 cm^{−3}. In the next 3 d (Day 4 to Day 6), IMF B_{z} varied with the same pattern, but the solarwind speed became 600 km s^{−1}. Next, IMF B_{z} was fixed at 1 nT, and the solarwind speed followed a similar rectangular pattern for 3 d (Day 7 to Day 9). After that, IMF B_{z} became −5 nT and the solarwind speed again followed the same rectangular pattern for 3 d (Day 10 to Day 12). The solarwind speed was then fixed at 600 km s^{−1}, and the solarwind density was perturbed with a similar rectangular pattern under a fixed IMF B_{z} of 1 nT (Day 13 to Day 15) and −5 nT (Day 16 to Day 18).
The result suggests that both AE substorms and Pi2 substorms tend to frequently occur when the IMF is southward and the solarwind speed is high. Newell et al. (2016) discussed the importance of the solarwind speed based on the analysis of substorm onsets identified from the SML index. The second panel suggests that the Pi2 substorms are also strongly dependent on the solarwind speed. Indeed, the effect of the solarwind speed seems to be more essential for Pi2 substorms. From Day 4 to Day 6 when the solarwind speed was high, the occurrence rate of Pi2 substorms tended to be high even under northward IMF situations, while the frequency of AE substorms tended to be depressed by a northward IMF. This is interpreted as showing that Pi2 pulsations can occur under geomagnetically quiet conditions without an intense auroral electrojet, which would be consistent with the result of Kwon et al. (2013). The existence of Pi2 pulsations under quiet conditions explains the poor predictability of Pi2 substorms suggested in Fig. 6. Since AE substorms are rare during northward IMF, we can confidently deny the occurrence of a substorm when the IMF is northward. In contrast, since Pi2 substorms can often occur even when the IMF is northward, it would be rare that we can confidently deny the occurrence of a substorm.
In contrast with the response to the IMF and solarwind speed, the response to the solarwind density variation is spiky, especially for AE substorms (from Day 13 to Day 18, shaded with green in the top and second panels). This may be interpreted as showing that a substorm is triggered by a sudden impulse (SI) of solarwind dynamic pressure, although it is possible that the SI effect on the geomagnetic variation is misidentified as a substorm occurrence. It is also notable that a momentary enhancement of the occurrence rate is seen at the northward turning of the IMF on Day 6 for Pi2 substorms (blue arrows in the second panel). This might suggest that Pi2 substorms tend to be triggered by a northward turning of the IMF, as suggested by Lyons et al. (1997). However, when we tried the same experiment with the ESN trained with Pi2 substorms identified with a threshold of θ=0.2, the momentary enhancement of the predicted occurrence rate at the northward turning became unclear. Thus, the response to the northward turning of the IMF is not a distinct feature. As suggested by other studies (Morley and Freeman, 2007; Wild et al., 2009), the northward turning of the IMF is not likely to be essential to substorm occurrence even though the northward turning may give favourable conditions for a substorm occurrence.
Excluding spikes of the substorm occurrence rate concurrent with jumps in the solarwind density, a higher solarwind density is likely to have increased the occurrence rate from Day 13 to Day 15 when the IMF was weak and northward. In contrast, from Day 16 to Day 18 when the IMF was southward, the occurrence rate did not show a clear change between before and after jumps of solarwind density except for spikes due to density jumps especially for AE substorms. These results suggest that the effect of solarwind density depends on the IMF. When the IMF is weak and northward, the substorm occurrence rate is higher when solarwind density is higher. Meanwhile, when the IMF is southward, the dependence of the occurrence rate on the solarwind density is not clear except that the solarwind density jumps can affect the occurrence rate. This trend is similar to the compound effect of the solarwind density and IMF on auroral electrojets as identified by NK22 and other studies (Ebihara et al., 2019; Kataoka et al., 2023). The solarwind density effect on the occurrence rate may thus contribute to the dependence of the auroral electrojet intensity on the solarwind density. However, according to NK22, a similar compound effect is also observed in the AU index, which represents an eastward electrojet. Since eastward electrojets are not considered part of the substorm current system (e.g. Kamide and Kokubun, 1996), the occurrence rate cannot fully explain the compound effect on the auroral electrojets indicated by NK22.
This study proposed an approach for analysing event time series of substorm onsets with the occurrences influenced by solarwind conditions. We treat a substorm onset as a stochastic phenomenon and represent the stochastic occurrences of substorms with a nonstationary Poisson process. The occurrence rate of substorms is then described with an echo state network model with a time series of solarwind data as an input. The echo state network allows us to diagnose the relationship between the substorm occurrence rate and solarwind conditions. We applied this proposed approach to a sequence of AE substorms and a sequence of Pi2 substorms. The ESN successfully predicted the probability of occurrences of AE and Pi2 substorms.
We also conducted experiments with a synthetic solarwind data set and obtained the following results.

The occurrence rate is enhanced under a higher solarwind speed and southward IMF for both AE and Pi2 substorms. The solarwind speed effect is more important for Pi2 substorms than for AE substorms. While a highspeed solar wind can bring Pi2 substorms even under a northward IMF, AE substorms are rarely observed under a northward IMF even if the solarwind speed is high.

The occurrence rate is enhanced at an abrupt jump of the solarwind density. A northward turning of the IMF is also likely to momentarily enhance the occurrence rate of Pi2 substorms, although the effect of a northward turning seems to be minor.

A compound effect of the solarwind density and IMF is also suggested. When the IMF is weak and northward, a higher solarwind density is likely to cause a higher occurrence rate of substorms. In contrast, when the IMF is southward, the solarwind density does not seem to make a clear effect on the occurrence rate of substorms except for spikes due to density jumps.
Parameters of a Bayesian prior distribution are often determined based on the marginal likelihood (e.g. Morris, 1983; Casella, 1985). We determine the standard deviation of the prior distribution, σ, by maximising the marginal likelihood which is defined as
where we denote the prior distribution in Eq. (8) as the conditional distribution given σ and p(βσ). Using Eq. (10), Eq. (A1) is reduced to
Denoting the optimal value which maximising J as $\widehat{\mathit{\beta}}$, it should satisfy $\mathrm{\nabla}J\left(\widehat{\mathit{\beta}}\right)=\mathrm{0}$. The Taylor expansion of J around $\widehat{\mathit{\beta}}$ is thus
where ${\mathrm{\nabla}}^{\mathrm{2}}J\left(\widehat{\mathit{\beta}}\right)$ is the Hessian matrix at $\mathit{\beta}=\widehat{\mathit{\beta}}$. This Hessian matrix is negative definite because J is maximised at $\mathit{\beta}=\widehat{\mathit{\beta}}$. This secondorder approximation of J yields an approximation of Eq. (A1) as follows:
where ${\mathrm{\nabla}}^{\mathrm{2}}J(\widehat{\mathit{\beta}}\left)\right$ is the determinant of the Hessian matrix of −J at $\mathit{\beta}=\widehat{\mathit{\beta}}$. This approximation is sometimes referred to as Laplace's approximation (e.g. Bishop, 2006). We chose the standard deviation, σ, which maximises the logarithm of the approximate marginal likelihood:
The AU and AL indices are available from the website of the WDC for Geomagnetism, Kyoto (http://wdc.kugi.kyotou.ac.jp/wdc/Sec3.html; World Data Center for Geomagnetism, Kyoto et al., 2015). The SMU and SML indices are available from the SuperMAG website (https://supermag.jhuapl.edu; Gjerloev, 2012). The SuperMAG substorm list can also be acquired through the SuperMAG website. The Wp index is available via the website of Masahito Nosé (https://www.isee.nagoyau.ac.jp/~nose.masahito/scubed/; World Data Center for Geomagnetism, Kyoto and Nosé, 2016). The OMNI solarwind data were acquired from the OMNIWeb of NASA/GSFC (https://omniweb.gsfc.nasa.gov/; King and Papitashvili, 2023).
SN conceived and conducted the analysis. RK contributed to the scientific interpretation. MN and JWG contributed to data processing and interpretation of the data.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We acknowledge the substorm timing list identified by the Newell and Gjerloev technique (Newell and Gjerloev, 2011a), the SMU and SML indices (Newell and Gjerloev, 2011b), and the SuperMAG collaboration (Gjerloev, 2012).
This research has been supported by the Japan Society for the Promotion of Science (grant no. 17H01704).
This paper was edited by Anna Milillo and reviewed by two anonymous referees.
Amariutei, O. A. and Ganushkina, N. Y.: On the prediction of the auroral westward electojet index, Ann. Geophys., 30, 841–847, https://doi.org/10.5194/angeo308412012, 2012. a
Baker, D. N., Pulkkinen, T. I., Angelopoulos, V., Baumjohann, W., and McPherron, R. L.: Neutral line model of substorms: Past results and present view, J. Geophys. Res., 101, 12975–13010, https://doi.org/10.1029/95JA03753, 1996. a
Bishop, C. M.: Pattern recognition and machine learning, Springer, New York, ISBN: 9780387310732, 2006. a
Casella, G.: An introduction to empirical Bayes data analysis, Amer. Stat., 39, 83–87, 1985. a
Cliver, E. W., Kamide, Y., and Ling, A. G.: Mountain and valleys: Semiannual variation of geomagnetic activity, J. Geophys. Res., 105, 2413–2424, 2000. a
Daley, D. J. and VereJones, D.: An introduction to the theory of point processes, Vol. I, Elementary theory and method, 2nd Edn., Chap. 7, Springer, New York, ISBN: 9780387955414, 2003. a, b
Davis, T. N. and Sugiura, M.: Auroral electrojet activity index AE and its universal time variations, J. Geophys. Res., 71, 785–801, 1966. a
Ebihara, Y., Tanaka, T., and Kamiyoshikawa, N.: New diagnosis for energy flow from solar wind to ionosphere during substorm: Global MHD simulation, J. Geophys. Res., 124, 360–378, https://doi.org/10.1029/2018JA026177, 2019. a
Gjerloev, J. W.: The SuperMAG data processing technique, J. Geophys. Res., 117, A09213, https://doi.org/10.1029/2012JA017683, 2012 (data available at https://supermag.jhuapl.edu/, last access: 19 November 2023). a, b, c
Gleisner, H. and Lundstedy, H.: Response of the auroral electrojets to the solar wind modled with neural networks, J. Geophys. Res., 102, 14269–14278, 1997. a
Jaeger, H. and Haas, H.: Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication, Science, 304, 78–80, https://doi.org/10.1126/science.1091277, 2004. a, b
Kamide, Y. and Kokubun, S.: Twocomponent auroral electrojet: Importance for substorm studies, J. Geophys. Res., 101, 13027–13046, 1996. a, b
Kataoka, R. and Nakano, S.: Reconstructing solar wind profiles associated with extreme magnetic storms: A machine learning approach, Geophys. Res. Lett., 48, e2021GL096275, https://doi.org/10.1029/2021GL096275, 2021. a
Kataoka, R., Nakano, S., and Fujita, S.: Machine learning emulator for physicsbased prediction of ionospheric potential response to solar wind variations, Earth Planets Space, 75, 139, https://doi.org/10.1186/s40623023018963, 2023. a
King, J. and Papitashvili, N.: One min and 5min solar wind data sets at the Earth's bow shock nose, NASA/GSFC [data set], https://omniweb.gsfc.nasa.gov/html/HROdocum.html (last access: 20 November 2023), 2023. a
Kwon, H.J., Kim, K.H., Jun, C.W., Takahashi, K., Lee, D.H., Jin, H., Seon, J., Park, Y.D., and Hwang, J.: Lowlatitude Pi2 pulsations during intervals of quiet geomagnetic conditions (K_{p}≤1), J. Geophys. Res., 118, 6145–6153, https://doi.org/10.1002/jgra.50582, 2013. a
Lui, A. T. Y.: Current disruption in the Earth's magnetosphere: Observations and models, J. Geophys. Res., 101, 13067–13088, https://doi.org/10.1029/96JA00079, 1996. a
Luo, B., Li, X., Temerin, M., and Liu, S.: Prediction of the AU, AL, and AE indices using solar wind parameters, J. Geophys. Res., 118, 7683–7694, https://doi.org/10.1002/2013JA019188, 2013. a
Lyons, R. R., Blanchard, G. T., Samson, J. C., Lepping, R. P., Yamamoto, T., and Moretto, T.: Coordinated observations demonstrating external substorm triggering, J. Geophys. Res., 102, 27039–27051, https://doi.org/10.1029/97JA02639, 1997. a
Lyons, R. R., Zou, Y., Nishimura, Y., GallardoLacourt, B., Angelopoulos, V., and Donovan, E. F.: Stormtime substorm onsets: occurrence and flow channel triggering, Earth Planets Space, 70, 81, https://doi.org/10.1186/s406230180857x, 2018. a
Maimaiti, M., Kunduri, B., Ruohoniemi, J. M., Baker, J. B. H., and House, L. L.: A deep learningbased approach to forecast the onset of magnetic substorms, Space Weather, 17, 1534–1552, https://doi.org/10.1029/2019SW002251, 2019. a
Miyashita, Y. and Ieda, A.: Revisiting substorm events with preonset aurora, Ann. Geophys., 36, 1419–1438, https://doi.org/10.5194/angeo3614192018, 2018. a
Morley, S. K. and Freeman, M. P.: On the association between northward turnings of the interplanetary magnetic field and substorm onsets, Geophys. Res. Lett., 34, L08104, https://doi.org/10.1029/2006GL028891, 2007. a
Morris, C. M.: Parametric empirical Bayes inference: theory and applications, J. Amer. Statist. Assoc., 78, 47–55, 1983. a
Nakano, S. and Iyemori, T.: Stormtime fieldaligned currents on the nightside inferred from groundbased magnetic data at mid latitudes: Relationships with the interplanetary magnetic field and substorms, J. Geophys. Res., 110, A07216, https://doi.org/10.1029/2004JA010737, 2005. a
Nakano, S. and Kataoka, R.: Echo state network model for analyzing solarwind effects on the AU and AL indices, Ann. Geophys., 40, 11–22, https://doi.org/10.5194/angeo40112022, 2022. a
Newell, P. T. and Gjerloev, J. W.: Evaluation of SuperMAG auroral electrojet indices as indicators of substorms and auroral power, J. Geophys. Res., 116, A12211, https://doi.org/10.1029/2011JA016779, 2011a. a, b, c
Newell, P. T. and Gjerloev, J. W.: Substorm and magnetosphere characteristic scales inferred from the SuperMAG auroral electrojet indices, J. Geophys. Res., 116, A12232, https://doi.org/10.1029/2011JA016936, 2011b. a, b
Newell, P. T., Liou, K., Gjerloev, J. W., Sotirelis, T., Wing, S., and Mitchell, E. J.: Substorm probabilities are best predicted from solar wind speed, J. Atmos. Sol.Terr. Phys., 146, 28–37, https://doi.org/10.1016/j.jastp.2016.04.019, 2016. a
Nishida, A.: Coherence of geomagnetic DP 2 fluctuations with interplanetary magnetic variations, J. Geophys. Res., 73, 5549–5559, 1968. a
Nosé, M., Iyemori, T., Takeda, T., Kamei, T., Milling, D. K., Orr, D., Singer, H. J., Worthington, E. W., and Sumitomo, N.: Automated detection of Pi2 pulsations using wavelet analysis: 1. Method and an application for substorm monitoring, Earth Planets Space, 50, 773–783, 1998. a
Nosé, M., Iyemori, T., Takeda, M., Toh, H., Ookawa, T., CifuentesNava, G., Matzka, J., Love, J. J., McCreadie, H., Tuncer, M. K., and Curto, J. J.: New substorm index derived from highresolution geomagnetic field data at low latitude and its comparison with AE and ASY indices, in: Proc. of XIIIth IAGA Workshop on Geomagnetic Observatory Instruments, Data Acquisition, and Processing, U.S. Geological Survey OpenFile Report 20091226, edited by: Love, J. J., 202–207, U.S. Geological Survey, https://doi.org/10.3133/ofr20091226, 2009. a, b
Nosé, M., Iyemori, T., Wang, L., Hitchman, A., Matzka, J., Feller, M., Egdorf, S., Gilder, S., Kumasaka, N., Koga, K., Matsumoto, H., Koshiishi, H., CifuentesNava, G., Curto, J. J., Segarra, A., and Çelik, C.: Wp index: A new substorm index derived from highresolution geomagnetic field data at low latitude, Space Weather, 10, S08002, https://doi.org/10.1029/2012SW000785, 2012. a, b, c
Schillings, A., Palin, L., Opgenoorth, H. J., Hamrin, M., Rosenqvist, L., Gjerloev, J. W., Juusola, L., and Barnes, R.: Distribution and occurrence frequency of dB/dt spikes during magnetic storms 1980–2020, Space Weather, 20, e2021SW002953, https://doi.org/10.1029/2021SW002953, 2022. a
Takalo, J. and Timonen, J.: Neural network prediction of AE data, Geophys. Res. Lett., 24, 2403–2406, 1997. a
Viljanen, A., Tanskanen, E. I., and Pulkkinen, A.: Relation between substorm characteristics and rapid temporal variations of the ground magnetic field, Ann. Geophys., 24, 725–733, https://doi.org/10.5194/angeo247252006, 2006. a
Wei, D., Dumlop, M. W., Yang, J., Dong, X., Yu, Y., and Wang, T.: Intense $\mathrm{d}B/\mathrm{d}t$ variations driven by nearearth bursty bulk flows (BBFs): A case study, Geophys. Res. Lett., 48, e2020GL091781, https://doi.org/10.1029/2020GL091781, 2021. a
Wild, J. A., Woodfield, E. E., and Morley, S. K.: On the triggering of auroral substorms by northward turnings of the interplanetary magnetic field, Ann. Geophys., 27, 3559–3570, https://doi.org/10.5194/angeo2735592009, 2009. a
World Data Center for Geomagnetism, Kyoto and Nose, M.: Geomagnetic Wp index, Data Analysis Center for Geomagnetism and Space Magnetism, Graduate School of Science, Kyoto University [data set], https://doi.org/10.17593/1343746800, available at https://www.isee.nagoyau.ac.jp/~nose.masahito/scubed/ (last access: 20 November 2023), 2016. a, b
World Data Center for Geomagnetism, Kyoto, Nosé, M., Iyemori, T., Sugiura, M., and Kamei, T.: Geomagnetic AE index, Data Analysis Center for Geomagnetism and Space Magnetism, Graduate School of Science, Kyoto University [data set], https://doi.org/10.17593/1503154800, available at http://wdc.kugi.kyotou.ac.jp/wdc/Sec3.html (last access: 19 November 2023), 2015. a
 Abstract
 Introduction
 Method
 Analysis of auroral electrojet substorms
 Analysis of Pi2 substorms
 Response to synthetic solar wind
 Concluding remarks
 Appendix A: Parameter determination with marginal likelihood
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Method
 Analysis of auroral electrojet substorms
 Analysis of Pi2 substorms
 Response to synthetic solar wind
 Concluding remarks
 Appendix A: Parameter determination with marginal likelihood
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References