Machine learning emulator for physics-based prediction of ionospheric potential response to solar wind variations

Physics-based simulations are important for elucidating the fundamental mechanisms behind the time-varying complex ionospheric conditions, such as ionospheric potential, against unprecedented solar wind variations incident on the Earth’s magnetosphere. However, carrying out an extensive parameter survey for comprehending the nonlinear solar wind density dependence of the ionospheric potential, for example, requires state-of-the-art global magnetohydrodynamic (MHD) simulations, which cannot be executed efficiently even on large-scale cluster computers. Here, we report the performance of a machine-learning based surrogate model for estimating the ionospheric potential outputs of a global MHD simulation, using the reservoir computing technique called echo state network (ESN). The trained ESN-based emulator demonstrates exceptional speed in conducting the parameter survey, which can lead to the identification of a solar wind density dependence of the ionospheric polar cap potential. Finally, we discuss future directions including the promising application for space weather forecasting.


Introduction
Machine learning (ML) techniques have been recognized as a useful tool for predicting geomagnetic activity indices, as represented by the AE, Kp, and Dst indices, using solar wind parameters as inputs, such as the solar wind speed (Vsw), proton density (Np), and interplanetary magnetic field (IMF).Many such attempts have been thoroughly reviewed by Liemohn et al. (2018).One of the latest studies used the ML technique called echo state network (ESN) to successfully predict the AE index and diagnose the nonlinearity of a magnetosphere-ionosphere coupled system using a synthetic solar wind time series (Nakano and Kataoka 2022).
Detailed patterns of ionospheric potential can be empirically modeled as a function of the solar wind parameters (e.g., Weimer 1995;2005).The empirical models, however, have limitations in that they predict the averaged observed patterns and cannot account for dynamic solar wind variations or unprecedented solar wind variations.
To overcome the limitations, a straightforward physics-based approach is to solve the idealized magnetohydrodynamic (MHD) equations for the solar wind and magnetospheric plasma flows, setting the appropriate height-integrated ionospheric conductivity layer as one of the boundary conditions (e.g., Tanaka et al. 2022a, b).
State-of-the-art global MHD simulations are also essential for understanding the mechanism behind the dynamically changing auroral oval distribution (e.g., Ebihara and Tanaka 2022).Further parameter surveys using such simulations are still necessary to understand the nonlinear effects, for example, the nonlinear density effect on the ionospheric potential or energy deposition (Khachikjan et al. 2008;Ebihara et al. 2019;Yang et al. 2020;Nakano and Kataoka 2022).However, the practical use of the global MHD simulation for such a parameter-survey purpose is still limited because it is computationally expensive.
A reasonable next approach to addressing the timeconsumption issue is to emulate the computationally expensive physics-based global MHD simulation by ML techniques, although the mechanism behind the simulation becomes a black box.Once such an ML-based emulator is developed, we can predict the dynamically changing ionospheric potential instantaneously.It is, therefore, practically possible for a small computer to make a thorough parameter survey by the ML-based emulator.
The purpose of this study is to show the very first attempt of such an ML-based emulator to imply the potential impacts.While an emulation that works with a global MHD simulation was attempted for parameter tuning (Kleiber et al. 2013;Heaton et al. 2015), here

Graphical Abstract
we develop an ESN-based emulator for predicting the dynamic ionospheric potential responses of the magnetosphere-ionosphere system to variable solar-wind time-series input.The methods of emulating global MHD simulation results by the ESN model are briefly introduced in the following Sections.In Section "Results and Discussions", we examine the obtained results from the newly developed ESN-based emulator.In Sect.4, we discuss the possible capabilities of the ESN-based emulator and some future directions.

Global MHD simulation
REPPU (REProduce Plasma Universe) is an MHD simulation code developed for studying the global magnetosphere-ionosphere coupling (Tanaka 1995;2015).The REPPU code is characterized by an excellent ionospheric reproduction of fundamental auroral phenomena such as substorms (Ebihara and Tanaka 2015a;2015b), sunaligned arcs (Tanaka et al. 2017), and the theta aurora (Tanaka et al. 2018).In this study, we used an improved REPPU simulation code, including the effects of a tilted dipole axis and seasonal changes of solar zenith angles.The total number of grid cells in the magnetosphere is 30722 (horizontal) times 240 (vertical), where the unstructured grid system described by Moriguchi et al. (2008) is employed.The number of grid cells in the ionosphere is 30722.The improved REPPU simulation is essentially the same as that used for the real-time REPPU simulation of space weather forecast executed by the National Institute of Information and Communications Technology (Nakamizo and Kubota 2021).We used highresolution OMNI-2 one-min data (Bx, By, Bz, Np, Vx, Vy, Vz, and proton temperature) for input to run the REPPU simulation.Note that we used B vectors in the GSM coordinate system and we provisionally used V vectors in the GSE coordinate system.The GSE coordinates for V were originally introduced to show the simulation's capability to capture the Vy and Vz components as directly obtained from the OMNI data.Note that the OMNI data has only GSE coordinates for velocity.This study ignored the contribution of Vy and Vz, and since the dominant Vx component is common in GSE and GSM, we utilized the simulation results as they were, considering the Vy and Vz contributions as artificially added noise.
We ran the REPPU simulation to obtain several-days worth of activity outputs for training and testing the emulator.Note that the REPPU simulations have usually been done for an hour-long activity, and the several-days worth or an hundred hour-long activity output is approximately two orders of magnitude larger computations than usual.For preparing the training dataset, we selected four different long-term output datasets of moderate and intense magnetic storm events as driven by corotating interaction regions (CIR) or coronal mass ejections (CME) (see Kataoka and Miyoshi (2006) and Borovsky and Denton (2006) for the difference between CIR and CME driven storms): Intense (Dst peak = − 130 nT) CIR storm (~ 24 h from 2015/10/7 0000 UT), moderate (Dst peak = − 56 nT) CIR storm (~ 21 h from 2015/10/18 0000 UT), moderate (Dst peak = − 87 nT) CME storm (~ 34 h from 2015/11/06 1200 UT), and intense (Dst peak at -166 nT) CME storm (~ 36 h from 2015/12/19 1200 UT).Further, we utilized the testing data from non-storm time of 16.5 h data from 2015/09/06 0000 UT.In this study, we trained the ESN model using a wide range of solar wind parameters with very limited simulation data.This is why we chose the non-storm time dataset for testing.In other words, if you use storm-time data for testing, ESN cannot learn from the precious storm-time data.We discarded the first 1 h of data for each run, in which the global plasma distribution cannot yet be physically realistic.In this study, these different simulation runs are simply stitched together in time.
We used 10 min averaged output data for the ionospheric potential.After filling the data gap with linear interpolation, the solar wind data were also 10-min averaged when they were fed into the ML model.For coarsegraining purpose, we binned the ionospheric potential map for the northern polar region (> 50 deg latitude) into 15 × 32 in latitude and longitude, respectively, from the original resolution of 60 × 320.Although this binning process significantly reduces the original information especially for small-scale features, the binning results still capture the essential global patterns such as dawn-dusk convection cells.However, note that the number of grids is better to be doubled for forecasting the field-aligned current, for example, because the narrow structures of Region-1 and Region-2 patterns are better resolved by ~ 1 deg in latitude.
We then applied the principal component analysis (PCA) to the simulation results to reduce the dimensions.Note that Cousins et al. (2015) applied a similar method to evaluate variable field-aligned current.In this study we used the PCA module of Python 3 scikit-learn (https:// scikit-learn.org/ stable/ modul es/ gener ated/ sklea rn.decom posit ion.PCA.html).To reproduce the > 90% variance of the original features, we decomposed the ionospheric potential pattern into 10 PCA components.

Machine learning technique
We require an ML technique that can be trained by a relatively small training dataset because the REPPU simulation is still computationally expensive for long runs.Therefore, standard deep learning techniques are not appropriate for this study.We also need a specific ML technique suitable for the time-series prediction, because the ionosphere dynamically changes following the time history of the magnetosphere and ionosphere.Both needs can be satisfied by the reservoir computing method ESN (Jaeger 2001;Jaeger and Haas 2004), as reviewed by Tanaka et al. (2019).We used essentially the same ESN method as documented by Kataoka and Nakano (2021) and Kataoka et al. (2022).In this study, we used the ESN module of Python 3 as developed by Tanaka et al. (2022a, b) (See https:// github.com/ GTANA KA-LAB/ DTS-ESN/).
The basic ESN model used in this study is described by the reservoir sate vector x and the model output vector y (time series of the PCA components of ionospheric potential) at t = n + 1 steps as follows: Here, the weight matrices W in and W are multiplied by the input vector u (the solar wind time series) and the reservoir state vector x, respectively.In this study, we set the number of the nodes (elements of x) to be 300.These W in and W are fixed, while only W out is trained by the Ridge regularization with regularization parameter of 10 -3 .
We created the random and sparse node connections of W, where only 10% of the matrix elements are the random values between − 1.0 and 1.0, and the remaining 90% are zero.We chose the optimal spectral radius (maximum eigenvalue of W) to be 0.65, by evaluating the (1) normalized root-mean square errors using both training and testing data.As the input vectors u, the solar wind speed and density are normalized as log 10 Vsw-2.5, and log 10 Np-1.0 before training the ESN model.It is noteworthy here that both the speed and density follow lognormal distributions (Burlaga and Lazarus 2000).

Results and discussions
An example result from the ESN-based emulator compared with a REPPU simulation is shown in Fig. 1.We can see a reasonable agreement between the results from ESN-based emulator and REPPU simulation; such as two-cell convection patterns.
Figure 2  Instead of directly comparing the emulation results with the observation data, which are available only for spatially limited areas, we benchmark the emulation results against the standard empirical models (e.g., Weimer 2005, Fig. 2). Figure 3 shows an example of the IMF clock-angle dependence of ionospheric potential using the synthetic solar wind data, fixing the density at Np = 5/cc, the solar wind speed at Vsw = 400 km/s, and the IMF strength at B = 5 nT.We selected the results from quasi-steady state time steps for each IMF direction The ESN-based emulator is exceptionally fast (approximately million times faster than the REPPU simulation) to run and it is, therefore, useful for a thorough parameter survey of so-called heatmap analysis, which cannot be done by observations or simulations.As an example, we input the synthetic solar wind variations to the ESNbased emulator to examine the nonlinear density effect on the cross-polar cap potential (CPCP).Note here that we fixed the solar wind speed in the input and changed the density, so that the dynamic pressure also changes accordingly.Again, only steady state results at 150 min past, or 15 steps later, after discontinuously changing the solar wind parameters, were used when generating the pixels in the heat maps.
As shown in Fig. 4, the CPCP has a negative dependence on the density during strong southward Bz (SBZ).In fact, a similar trend has been identified by Khachikjan et al. (2008) using active-time SuperDARN observations, as they discussed that the shrinking dayside magnetopause by high dynamic pressure may be relevant to reduce the possible effect of dayside reconnection that is powering the CPCP.Recognizing that the overall IMF By dependence was fairly reproduced by the ESN-emulator as shown in Fig. 2, e.g., round dusk cells for positive By and crescent dusk cells for negative By, we further diagnose the density dependence by changing the IMF By component, and turning off the contribution of SBZ.From the heatmap analysis without SBZ component, Fig. 5 shows that the density dependence is positive when By < 4 nT, while it is negative for larger By.The By-asymmetric density effect is a new result "predicted" from the ESN-based emulator, which suggests that the nonlinear density effect on CPCP can be more complex than imagined from the SBZ reconnection hypothesis.Note that the IMF By dependence does not clearly appear if we integrate the obtained By-Np heatmap in density, which is consistent with SuperDARN observations (Mori and Kustov 2013).Related future works therefore include the observationbased identification of the By-asymmetric density effect as well as the detailed analysis of the global MHD simulation results of both magnetosphere and ionosphere to identify the exact mechanism that causes the By-asymmetric density effect.In this study, we used simulation data for CIR-and CME-driven magnetic storm events that occurred in the fall to winter seasons.The training data set is, therefore, limited to the particular situation where the auroral activity is high and the northern hemisphere is darker than the southern hemisphere.A detailed discussion of the sunlight effect associated with solar zenith angles and the general north-south asymmetry are, therefore, interesting future subjects of study when we have accumulated longer term REPPU simulation results, including different seasons.When such solar zenith angle effect is included, the future emulator will be ready for the operational space weather forecast via simply replacing the input file with the real-time solar wind data.
This paper presented the very first attempt to demonstrate the potential impact of the ESN-based emulator.Future studies using the ESN-based emulator of the REPPU simulation include an examination of the accuracy in reproducing observation data, such as SuperDARN convection maps.Any partial data or point data can also be used with cutting-edge data assimilation techniques (Nakano et al. 2020).For a data assimilation, we need to increase the ensemble number of simulation runs for integrating large probabilistic ensembles, and the ESN-based emulator will play an essential role.In this manner, the ESN-based emulator can be expected to become a basic technique for future space weather reanalysis studies.
Further, note that the OMNI data can be different from the actual incident solar wind at the magnetosphere, because it is computed using different spacecraft data, also relying on an analytic time-shift method from those spacecraft positions.The possible errors as investigated by Vokhmyanin et al. (2019) can, therefore, be used as an important probability distribution to work on the data assimilation applications in future.
shows the performance of training and testing the ESN emulator of the ionospheric potential.The input solar wind parameters are shown in the top panel.The top five PCA components and the prediction from the ESN model are also shown.The training interval is t < 664 steps.The trained ESN model reproduces the time variation of the PCA variables for both the training and the testing intervals.

Fig. 1
Fig. 1 Snapshot example of the ionospheric potential as obtained from ESN-based emulator (left) and REPPU simulation (right) at a testing time of t = 712 steps

Fig. 2
Fig. 2 Solar wind parameters (a) and top five PCA components (b-f) for ionospheric potential.Black curves show the REPPU simulation results and the red curves are from the ESN-based emulator.The testing time interval, which is not used for training, is shown by the gray hatched time interval

Fig. 3
Fig. 3 Quasi-steady state ionospheric potential patterns as obtained from the ESN-based emulator for different IMF clock angles using synthetic solar wind data of Vsw = 400 km/s, Np = 5.0/cc, and IMF strength = 5.0 nT

Fig. 4
Fig. 4 Heatmap analysis of the parameter survey to examine the nonlinear density effect on the cross polar-cap potential, changing the SBZ strength, fixing the solar wind speed at 400 km/s