Global Earth Magnetic Field Modeling and Forecasting with Spherical Harmonics Decomposition
Panagiotis Tigas, Téo Bloch, Vishal Upendran, Banafsheh Ferdoushi, Mark C. M. Cheung, Siddha Ganju, Ryan M. McGranaghan, Yarin Gal, Asti Bhatt
GGlobal Earth Magnetic Field Modeling andForecasting with Spherical Harmonics Decomposition
Panagiotis Tigas ∗ OATMLUniversity of OxfordOxford, UK
To Bloch ∗ University of ReadingReading, UK
Vishal Upendran ∗ IUCAAPune, India
Banafsheh Ferdoushi ∗ University of New HampshireDurham, NH, USA
Yarin Gal
OATMLUniversity of OxfordOxford, UK
Siddha Ganju
NVIDIA CorporationSanta Clara, CA, USA
Ryan M. McGranaghan
ASTRA LLCLouisville, CO, USA
Mark C. M. Cheung
Lockheed MartinAdvanced Technology CenterPalo Alto, CA, USA
Asti Bhatt
SRI InternationalMenlo Park, CA, USA
Abstract
Modeling and forecasting the solar wind-driven global magnetic field perturbationsis an open challenge. Current approaches depend on simulations of computationallydemanding models like the Magnetohydrodynamics (
MHD ) model or samplingspatially and temporally through sparse ground-based stations (
SuperMAG ). Inthis paper, we develop a Deep Learning model that forecasts in Spherical Harmonicsspace , replacing reliance on MHD models and providing global coverage at one-minute cadence, improving over the current state-of-the-art which relies on featureengineering. We evaluate the performance in
SuperMAG dataset (improved by14.53%) and
MHD simulations (improved by 24.35%). Additionally, we evaluatethe extrapolation performance of the spherical harmonics reconstruction based onsparse ground-based stations (
SuperMAG ), showing that spherical harmonics canreliably reconstruct the global magnetic field as evaluated on
MHD simulation.
The space environment around Earth (geospace) does not exist in a steady state. The dynamicalchanges in geospace is termed space weather. Space weather is primarily driven by interactionsbetween the Sun’s output (in terms of solar radiation and the magnetized solar wind) with theEarth’s magnetosphere, thermosphere and ionosphere. The US National Oceanic and AtmosphericAdministration continuously monitors the solar wind upstream from Earth with satellite observatoriesorbiting Lagrangian point L1 along the Sun-Earth line. The solar wind energy is transferred to theEarth’s magnetosphere via complex mechanisms [Dungey, 1961], leading to perturbations in theEarth’s magnetic field called geomagnetic storms. ∗ equal contribution We will release all the code and models used in this work post-publication.Third Workshop on Machine Learning and the Physical Sciences (NeurIPS 2020), Vancouver, Canada. a r X i v : . [ phy s i c s . g e o - ph ] F e b olar wind parameters Predictions SampledPredictions Figure 1: Model architecture. Solar Wind data are summarized into an embedding vector usinga
Gated Recurrent Unit (GRU) which is then passed to a
Multilayer Perceptron (MLP) to outputthe predicted spherical harmonics coefficients from where we can extrapolate values for specificlocations.Geomagnetic storms drive a spectrum of potentially catastrophic disruptions to our technologically-dependent society, among the most threatening being critical disturbances to the electrical grid in theform of geomagnetically induced currents ( GICs ). Due to their proprietary nature, publicly available
GIC data are limited. However, a cohort study of insurance claims of electrical equipment providesevidence that space weather poses a continuous threat to electrical distribution grids via geomagneticstorms and
GICs [Schrijver et al., 2014, Eastwood et al., 2018].
GICs also pose threats to oilpipelines, railyways and telecommunication systems. In the case of extreme, but historically probablegeomagnetic storms, the economic impact due to prolonged power outages can exceed billions ofdollars per day [Oughton et al., 2017]. For this reason, there is urgency among public and industrystakeholders to improve monitoring and forecasting of space weather impacts like geomagneticstorms and
GICs .The geomagnetic field is continuously monitored, importantly with sparse spatial coverage, by anetwork of roughly 200 ground magnetometers [Gjerloev, 2012b] and by a constellation of 66satellites collectively known as the
Active Magnetosphere and Planetary Electrodynamics ResponseExperiment ( AMPERE ) project [AMPERE; Waters et al., 2020] in low Earth orbit. The spatially-sparse magnetometer measurements are typically synthesized into global indices (e.g.
Dst , Kp , and AE ; see appendix for more details) as measures of the geoeffectiveness of space weather perturbations.Like most indices, they are good as indicators but are too far removed from the underlying governingequations of the system. This makes interpretability and forecasting (whether by physics-based orphysics-agnostic ML models) difficult.GICs are driven by geoelectric fields through Ohms law, (cid:126)J = σ (cid:126)E (the displacement current isnegligible), where (cid:126)J is the current, σ the conductivity tensor, and (cid:126)E is the electric field. Thegeoelectric field is related to the ground magnetic field perturbations through the Faradays lawof induction −∇ × (cid:126)E = ∂ (cid:126)B∂t , where, ∂ (cid:126)B∂t is the rate of change of magnetic field on the Earth’ssurface and can be derived from the time derivative of magnetic perturbation based on ground Earthmagnetometers. Our contributions are:
1. Using simulation data from a physics-based (
MHD ) model to validate a compressed sensingtechnique to recover global maps (in spherical harmonic basis) of the geomagnetic perturba-tion from sparse measurements. This improves the temporal cadence of such maps by > We used data describing Earth’s magnetosphere and solar wind properties between and
31 Dec 2013 . We used
OMNI dataset for solar wind data, which captures the interplanetarymagnetic field components (IMF), velocity and temperature of the solar wind as well as the clock2ngle. . For Earth’s magnetic field, we used SuperMAG dataset [Gjerloev, 2012a] which consistsof geomagnetic perturbations measurements from ground earth stations located at various placesaround the globe.Finally, we used a simulation-derived dataset based on simulations conducted with Open GeospaceGeneral Circulation Model (
Open GGCM ) [Raeder et al., 2001], a magnetohydrodynamics (
MHD )model. Lacking a global and spatially-complete ground-based magnetometer array we must rely onthe first principles
MHD model to constitute a global ground truth dataset. Furthermore, this allowsus to validate our approach on finer resolution than one supported by
SuperMAG and
AMPERE dataset. However, note that the
MHD model is not a substitute for actual perturbation measurements,for many short-scale phenomena are missed by such models [Raeder et al., 2001].
In order to properly forecast GICs, we need to be able to model the relationship between the solarwind parameters and interplanetary magnetic fields (IMF) with the Earth system. Due to the scarcityof GIC data, the comparatively abundant magnetometer data and the relationship between the two, alogical step is to create models to predict ∂ (cid:126)B∂t . Such models provide complete spatial coverage andpermit the study of the nature of the connection between the solar wind and the Earth.
Figure 2: Ground truthobservations (top) andSpherical harmonic re-construction (bottom)of dbn_nez .Any scalar field over the unit sphere can be expressed as f ( θ, φ ) = (cid:80) ∞ n =0 (cid:80) nm = − n a nm Y nm ( θ, φ ) , where Y nm ( θ, φ ) := (cid:115) n + 14 π ( n − m )!( n + m )! e imθ P mn (cos( φ )) , and P mn (cos( φ )) are the associated Legendre polynomials. These functions Y nm ( θ, φ ) are solutions to Laplace equation in a 3-D spherically symmetriccoordinate system. If the sum is truncated at maximum harmonic degree N , f ( θ, φ ) is approximated as ˜ f ( θ, φ ) = (cid:80) Nn =0 (cid:80) nm = − n a nm Y nm ( θ, φ ) .Defining i = n + n + m , the double sum can be written as ˜ f ( θ, φ ) = (cid:80) ( N +1) − i =0 a i Y i ( θ, φ ) . If the 2D fields over θ, φ are unrolled as one-dimensional arrays, we have ˜ f = B (cid:126)a , where (cid:126)a = ( a i ) is the tuple ofspherical harmonic coefficients, and B = ( (cid:126)b i ) is the basis matrix whereincolumn vector (cid:126)b i represents the basis function Y nm ( θ, φ ). With
SuperMAG data, we have spatially sparse samples y which can beany component of the geomagnetic field, or its deviation from a referencefield) measured at ˜200 stations. Note that only the Northern Hemisphereis considered in this analysis, for it has an extensive coverage (in contrast to the extremely sparsecoverage in the Southern Hemisphere). The spherically symmetric reconstruction thus is ˜ f ( θ, φ ) ,and this reconstruction as evaluated at the station locations is ˜ y . Thus, we first construct a set ofcoefficients (cid:126)a from this sparse set of measurements, considered all the measurements across the globeat each time step. There are two constraints imposed for obtaining the set of coefficients:1. The coefficient set must be sparse: This is to prevent power leakage and coupling acrossmultiple modes. Solar wind measurement is performed by ACE and WIND satellites at the L1 Lagrange point of theSun-Earth system
3. As small N as possible: Since, in theory, an infinite number of modes can fit the observations,but a parsimonious model is required to limit overfitting. Furthermore, we would not wantunnatural, localized artifacts due to high number of modes, thereby motivating a constrainton n .To mitigate constraint ( ), we apply a Lasso regression technique [Tibshirani, 1996] with the sphericalharmonic functions as the basis. Figure 3: L1 error (left) and R2 (right) "knee" determination.The Lasso regression comes witha regularization term α (cid:107) (cid:126)a (cid:107) , where (cid:107) (cid:126)a (cid:107) is the L1-norm of the coeffi-cients and α > is a hyperparame-ter. This, and constraint ( ) are miti-gated by varying both α and the maxi-mum number of modes N , and search-ing for a knee in a defined error met-ric, subject to the smallest maximummodes. The sweep parameters are de-tailed in Table. 2, and two metrics are used – Maximum L1 error across all stations and time, andmaximum R2 metric (coefficient of determination) across all time steps. The maximum L1 error tellsus the worst performance across the dataset, and thus we seek the most acceptable worst possibleperformance. As shown in Fig. 3, the "knee of goodness" can be seen corresponding to α = 0 . , and20 maximum modes. These parameters are fixed in our analysis.An example reconstruction of the data from 2017, considering stations above ◦ is shown in Fig. 2.Here, we compare the north facing component of ∂d/∂t from SuperMAG and its reconstruction.Similarly, the comparison of spherical harmonic reconstruction for the
MHD simulation is shown inFig. 4. Note here, that the reconstruction is performed by sampling the
MHD simulation at locationsof
SuperMAG stations alone.
Next, we construct a forecasting model which uses solar wind data (
OMNI ) to forecast the globalmagnetic field perturbation. For this experiment we use a similar setup as Weimer [2013] model . Wefeed 25 minutes of solar wind activity into a Gated Recurrent Network (GRU) to map the sequenceinto an embedding vector. We then feed the embedding into a Multilayer Perceptron (MLP) to outputthe spherical harmonics coefficients which model the global magnetic field perturbation; specificallywe focus on north component of ∂d/∂t as a proof of concept. In contrast to Weimer [2013] model ,we use the whole sequence as input to a non-linear autoregressive model and we do not apply featureengineering to the
OMNI data. Instead we use the raw features (see appendix for detailed list offeatures). The architecture can be seen in fig. 1.We benchmark this work against the state-of-the-art empirical model by Wiemer et al. [Weimer,2013]. To evaluate the performance we first compare on the validation set on
SuperMAG but also wecompare on simulations conducted with
MHD model for two weeks worth of activity. The results aresummarized in table 1.
Table 1: Forecasting model performanceModel
SuperMAG (val) RMS (nT) ↓ MHD
RMS (nT) ↓ Ours.
Weimer [2013] model 28.35 35.72With this work, we show and evaluate the reconstruction of the global magnetic perturbation field usingspherical harmonics with LASSO regularization to promote sparsity in the coefficients. We show thatit is possible to reconstruct from sparse measurements like the ones provided by
SuperMAG stations,4igure 4: Comparison of reconstruction with the MHD simulation. Left: MHD simulation; center:MHD simulation sampled at SuperMAG station locations; right: Reconstruction.and we evaluated on
MHD data to show that the reconstruction is similar to the global field modeledby
MHD , suggesting that the results of such models can be compressed in spherical harmonics spaceby compressed sensing techniques. Additionally, we show that by using a Deep Neural Networkto forecast in the spherical harmonics space, we can improve over existing state of the art (Weimer[2013] model ) by 14.53% on
SuperMAG and 24.35% on
MHD dataset (summarized in table 1).
Broader Impact
Geomagnetic storms drive a spectrum of potentially catastrophic disruptions to our technologically-dependent society, among the most threatening being critical disturbances to the electrical grid in theform of geomagnetically induced currents ( GICs ). Due to their proprietary nature, publicly available
GIC data are limited. However, a cohort study of insurance claims of electrical equipment providesevidence that space weather poses a continuous threat to electrical distribution grids via geomagneticstorms and
GICs [Schrijver et al., 2014, Eastwood et al., 2018].
GICs also pose threats to oilpipelines, railways and telecommunication systems. In the case of extreme, but historically probablegeomagnetic storms, the economic impact due to prolonged power outages can exceed billions ofdollars per day [Oughton et al., 2017]. For this reason, there is urgency among public and industrystakeholders to improve monitoring and forecasting of space weather impacts like geomagneticstorms and
GICs . With this work, we progress towards better and faster forecasting models whichwill help us shield our infrastructure from solar-wind related hazards.
Acknowledgments and Disclosure of Funding
This project was conducted during the 2020 NASA Frontier Development Lab (FDL) program, apublic-private partnership between NASA, the SETI Institute, and commercial partners. We wishto thank, in particular, NASA, Google Cloud, NVIDIA, Intel and SRI, for supporting this project.Finally, we gratefully acknowledge the SuperMAG collaborators.
References
J. W. Dungey. Interplanetary magnetic field and the auroral zones.
Physical Review Letters , 6(2):47,1961.J. Eastwood, M. Hapgood, E. Biffis, D. Benedetti, M. Bisi, L. Green, R. Bentley, and C. Burnett.Quantifying the economic value of space weather forecasting for power grids: An exploratorystudy.
Space weather , 16(12):2052–2067, 2018.J. Gjerloev. The supermag data processing technique.
Journal of Geophysical Research: SpacePhysics , 117(A9), 2012a.J. W. Gjerloev. The SuperMAG data processing technique.
Journal of Geophysical Research: SpacePhysics , 117(A9), 2012b. ISSN 2156-2202. doi: 10.1029/2012JA017683.5. J. Oughton, A. Skelton, R. B. Horne, A. W. P. Thomson, and C. T. Gaunt. Quantifying the dailyeconomic impact of extreme space weather due to failure in electricity transmission infrastructure.
Space Weather , 15(1):65–83, Jan. 2017. doi: 10.1002/2016SW001491.J. Raeder, Y. Wang, and T. J. Fuller-Rowell. Geomagnetic storm simulation with a coupledmagnetosphere-ionosphere-thermosphere model.
Space weather , 125:377–384, 2001.C. J. Schrijver, R. Dobbins, W. Murtagh, and S. M. Petrinec. Assessing the impact of space weatheron the electric power grid based on insurance claims for industrial electrical equipment.
SpaceWeather , 12(7):487–498, July 2014. doi: 10.1002/2014SW001066.R. Tibshirani. Regression shrinkage and selection via the lasso.
Journal of the Royal StatisticalSociety: Series B (Methodological) , 58(1):267–288, 1996.C. L. Waters, B. J. Anderson, D. L. Green, H. Korth, R. J. Barnes, and H. Vanhamäki.
Science DataProducts for AMPERE , volume 17, page 141. 2020. doi: 10.1007/978-3-030-26732-2_7.D. R. Weimer. An empirical model of ground-level geomagnetic perturbations.
Space Weather-theInternational Journal of Research and Applications , 11(3):107–120, 2013.6
Reconstruction experiment sweep parameter
Table 2: Sweep parameter description for reconstructionSweep parameter DescriptionDate range 2015-06-23 to 2015-06-24Cadence 10 minMax modes sweep [5,50], steps of 5 α sweep {2,1,0.1,1e-3,1e-5}Max L1 metric max stations max time || y − ˜ y || Max R2 score max time (1 − (cid:80) stations | y − ˜ y | / variance( y )) Note that we have used a cadence of 10 min for quick computation of the metrics over the dataset –however, since the spherical harmonic generation is done at each time step, it can be performed atwhatever cadence necessary.
B Solar Wind Data -
OMNI
Dataset
Here we describe the features we used from
OMNI
Dataset.Table 3: Description of
OMNI
Dataset featuresFeature Description B T : Magnetic field magnitude ( (cid:113) ( B x + B y ) ) V SW : Solar Wind velocity magnitude ( (cid:113) ( V x + V y + V z ) ) T : Temperature of the solar wind θ c : Clock angle of the interplanetary magnetic field (IMF) F . : F10.7 measures the noise level generated by the sun at a wavelength of 10.7 cm.
C Geoeffectieness Indices
Dst: (Disturbance Storm Time Index) It is the measure of geomagnetic activity derived from nearequator ground magnetic stations providing information about the strength of ring current.
Kp:
Global geomagnetic index that is based on 3 hour measurements of mid-latitude ground magneticstations around the world.
AE: (Auroral Electrojet Index) It is the measure of auroral activity determined based on groundmagnetic stations around aurora zone.
D Weimer [2013] model
Weimer [2013] model baseline model is designed to forecast each of the three magnetic vectorcomponents. It uses solar wind data as input and it outputs the spherical harmonics coefficients,which then can be used to extract the forecasting values per CGM latitude and MLT pair.A feature vector from the solar wind data is derived7 = B T V SW t (cid:112) ( F . ) B T cos( θ c ) V SW cos( θ c ) t cos( θ c ) (cid:112) ( F . ) cos( θ c ) B T sin( θ c ) V SW sin( θ c ) t sin( θ c ) (cid:112) ( F . ) sin( θ c ) B T cos(2 θ c ) V SW cos(2 θ c ) B T sin(2 θ c ) V SW sin(2 θ c ) , (1)where B T represents the magnitude of the magnetic field, V SW the solar wind velocity, (cid:112) ( F . ) the square root of the F . feature (a measure of solar radiation), t the dipole axis angle in radiansand θ c the clock angle.Weimer [2013] model is trained on data obtained from SuperMAG stations for the year 2013 withthe loss MSE = Avg (( y − Ba ) ) , where a are the spherical harmonics coefficients, computed as a mn = ( g mn (cid:124)(cid:123)(cid:122)(cid:125) Real Part , h mn (cid:124)(cid:123)(cid:122)(cid:125) Imaginary Part ) The coefficients are computed as g mn = G mn F and h mn = H mn F , where G mn and H mn are the weightsto learn.The model is trained on 25 minute long average of solar wind data, with a lag of 20 minutes and astarget the magnetic field perturbation of 5 minutes long averages, as measured by SuperMAG stations.
E Reproducibility details of forecasting experiments
The model architecture used for the forecasting experiment is presented below. The model wastrained using Adam optimizer with learning rate lr=1e-04 and Mean Squared Error as a loss. class GeoeffectiveNet(nn.Module):def __init__(self,past_omni_length,future_length,omni_features,supermag_features,nmax,targets):super(GeoeffectiveNet, self).__init__()self.omni_past_encoder = nn.GRU(25,16,num_layers=1,bidirectional=False,batch_first=True,dropout=0.5)n_coeffs = 0 or n in range(nmax+1):for m in range(0, n+1):n_coeffs += 1n_coeffs *= 2self.encoder_mlp = nn.Sequential(nn.Linear(16, 16),nn.ELU(inplace=True),nn.Dropout(p=0.5),nn.Linear(16, n_coeffs, bias=False) )self.omni_features = omni_featuresself.supermag_features = supermag_featuresself.targets = targetsself.future_length = future_lengthdef forward(self,past_omni,past_supermag,future_supermag,dates,future_dates,**kargs):past_omni = NamedAccess(past_omni, self.omni_features)features = [] bt = (past_omni['by']**2 + past_omni['bz']**2)**.5v = (past_omni['vx']**2 + past_omni['vy']**2 + past_omni['vz']**2)**.5features.append(past_omni['bx'])features.append(past_omni['by'])features.append(past_omni['bz'])features.append(bt)features.append(v)features.append(past_omni['dipole'])features.append(torch.sqrt(past_omni['f107']))features.append(bt*torch.cos(past_omni['clock_angle']))features.append(v*torch.cos(past_omni['clock_angle']))features.append(past_omni['dipole']*torch.cos(past_omni['clock_angle']))features.append(torch.sqrt(past_omni['f107'])*torch.cos(past_omni['clock_angle']))features.append(bt*torch.sin(past_omni['clock_angle']))features.append(v*torch.sin(past_omni['clock_angle']))features.append(past_omni['dipole']*torch.sin(past_omni['clock_angle']))features.append(torch.sqrt(past_omni['f107'])*torch.sin(past_omni['clock_angle']))features.append(bt*torch.cos(2*past_omni['clock_angle']))features.append(v*torch.cos(2*past_omni['clock_angle']))features.append(past_omni['dipole']*torch.cos(2*past_omni['clock_angle']))features.append(torch.sqrt(past_omni['f107'])*torch.cos(2*past_omni['clock_angle']))features.append(bt*torch.sin(2*past_omni['clock_angle']))features.append(v*torch.sin(2*past_omni['clock_angle']))features.append(past_omni['dipole']*torch.sin(2*past_omni['clock_angle'])) eatures.append(torch.sqrt(past_omni['f107'])*torch.sin(2*past_omni['clock_angle']))features.append(past_omni['clock_angle'])features.append(past_omni['temperature'])features = torch.stack(features, -1)encoded = self.omni_past_encoder(features)[1][0]coeffs = self.encoder_mlp(encoded)predictions = torch.einsum('bij,bj->bi', future_supermag.squeeze(1), coeffs)return coeffs, predictions, None For the spherical harmonics decomposition we used scipy˙special˙sph_harm function. def basis_matrix(nmax, theta, phi):from scipy.special import sph_harmassert(len(theta) == len(phi))basis = []for n in range(nmax+1):for m in range(-n,n+1):y_mn = sph_harm(m, n, theta, phi)basis.append(y_mn.real.ravel())basis.append(y_mn.imag.ravel())basis = np.array(basis)return basis.reshape(-1, theta.shape[0], theta.shape[1]).swapaxes(0, 1).swapaxes(2, 1)