- Email: [email protected]
The nonlinear dynamics of space weather
Pergamon www.elsevier.nl/locate/asr
Adv. SpaceRes. Vol. 26, No. 1, pp. 197-207,200O 0 2000 COSPAR. Published bv Elsevier Science Ltd. All rightsreserved Printedin &eat Britain 0273-1177/00 $20.00 + 0.00 Pll: SO273-1177(99)01050-9
THE NONLINEAR DYNAMICS OF SPACE WEATHER
D. Vassiliadisl,
A. J. Klimas2, J. A. Valdivia3, D. N. Baker4
’ Universities Space Research Association, at: NASA/GSFC, Code 692, Greenbelt, MD 20771, USA 2Luboratory for Extraterrestrial Physics, NASA/GSFC, Greenbelt, MD 20771, USA 3National Research Council, at: NASA/GSFC, Greenbelt, MD 20771, USA ‘Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80309-0590, USA
ABSTRACT Studies of the nonlinear magnetospheric dynamics have led to several directions useful in understanding space physics processes, in particular those related to magnetospheric currents, and making space weather forecasts possible. Four such directions are identified: (a) empirical time series prediction with nonlinear autoregressive moving-average (ARMA) models, for which an example is given in terms of a geosynchronous electron flux index model. (b) Measurement of physical properties of the currents from the coefficients of the ARMA models, such as characteristic time scales and coupling strengths. Using a Dst index model we measure the ring current decay time as a function of storm phase and activity, and identify new oscillation time scales correlating with the substorm injection activity. (c) Spatiotemporal nonlinear modeling predicts the amplitude and location of the disturbance as a function of space, as well as its time evolution. An example is given in terms of predicting the longitudinal and temporal profile of midlatitude geomagnetic disturbances. (d) Nonlinear dynamical models can be coupled to other physical or empirical approaches to build comprehensive space weather models. SWIFT, an ionospheric space weather model, is discussed as an example. Other applications of nonlinear dynamics are briefly discussed. 0 2000 COSPAR. Published by Elsevier Science Ltd. 1. INTRODUCTION Nonlinear dynamics in space physics has grown steadily in scope and depth. We discuss this growth in terms of four directions whose results have useful implications for theoretical and applied physics, in particular space weather. While earlier research focuses on classification, empirical prediction, and first-principles modeling of geomagnetic time series (review: Klimas et al., 1996), recently three more directions have been developed: closed-form timeseries models which are data-derived, but in addition have a physical interpretation; spatial, in addition to temporal, prediction of space-related disturbances, and modeling from multiple measurement probes in the same region, or in several regions; and combinations of the nonlinear dynamics to other, physical or empirical, models for more comprehensive space weather forecasting. We discuss these directions in the above order. The starting point is a brief review of basic nonlinear dynamics concepts and some of the early results in Section 2. The following section discusses empirical time series modeling and illustrates it with the prediction of a newly proposed geosynchronous electron flux index. The next two sections (4 and 5) are devoted to physical time series modeling and spatiotemporal modeling. Following that, Sec. 6 discusses the synergy between nonlinear dynamical and other models for space weather applications, and finally ways of making predictions and other results available on-line.
197
D. Vassiliadis et al.
198
2. WHENCE
NONLINEAR
DYNAMICS?
A basic goal of nonlinear dynamics as applied to space physics is to model current systems in the magnetosphere and ionosphere and the geomagnetic disturbances they cause. The configuration and intensity of the basic current systems determine the global magnetospheric structure and activity in the predominantly low-p plasma (Potemra, 1984). The justification for applying nonlinear dynamics methods is the conspicuous coherence of these currents, and other aspects of magnetospheric activity. However, a theoretical understanding of how self-organization arises in plasmas for realistic ranges of plasma parameters is missing. Regarding the currents, their coherent nature is characterized by their (a) being permanent and stable features, and (b) by their having large interaction lengths and fast interaction timescales. The ring, electrojet, magnetopause, and tail currents are permanent features of the magnetosphere. Following a small-scale disturbance they rapidly return to their equilibrium values, except if a stability threshold is‘exceeded (e.g. in the case of the tail current). Because of their long-range magnetic effect, and the rapid signal speed (VA 2 L/T where L is a scale of the magnetospheric subsystem and T is a time scale for current growth), the magnetosphere is more organized spatially and responds faster than, e.g., the neutral atmosphere, a fluid of much smaller volume. An example of coherent activity in various magnetospheric locations is shown in Figure 1. Solar wind driving (from a high-speed stream at the end of solar cycle 22) excites simultaneous geomagnetic disturbances in polar, auroral, and midlatitude locations (whose activity is normally characterized by quite different time scales) while substorm injections often’ produce high geosynchronous particle fluxes. (c) Additional observational evidence about magnetospheric coherence has been discussed by Siscoe (1993). (d) Finally the predictability of magnetospheric-current processes, as will be reviewed below, is another body of evidence for the coherence. Hence the observed magnetospheric collective response is a direct argument for using nonlinear dynamics methods. The underlying assumption for most methods is that the spatial characteristics of the system (here: the geomagnetic activity) are due to a small number of spatial modes (corresponding to the geomagnetic effect of distinct current systems), and the complex activity can be described by the interaction, probably nonlinear, of these modes. The need for nonlinear data analysis became obvious after the linear moving-average (MA) filters were applied (Clauer, 1986). In examining the effect of solar wind driving, represented by the VB, electric field, on the energization of the magnetosphere, represented by the AL index, linear prediction filters were used, and they showed the first predictive capability
T G(t)
=
H(r)
* [VB,](t
- T)dT
J 0
Z(t Fig. 1. Coherent magnetospheric response during the month of 1995 to a solar wind high-speed stream (not shown). From top to bottom: Solar wind input; Geosynchronous orbit: index for the energetic electron flux (GEF); logro(GEF) showing electron depletions due to growth phase dipolarizations; Polar latitudes: index measuring the ground geomagnetic disturbance near the cap center; Mid-latitudes: Dst geomagnetic index.
(1)
t At) = 2
H(jAt)[VB,](t
- jAt)
j=O
January
where the first equation is for continuous ond for discrete-time data of resolution denotes the filter output.
and the secAt, and A^L
Now the filter’s impulse response function, H(r), where T is the time lag relative to the onset of solar
199
Nonlinear Dynamics of Space
wind input, is calculated from the VB, and AL data by inverting (1). Then the filter can be convolved with the same solar wind data to reproduce AL (in-sample), or to different solar wind data to predict the corresponding AL (outof-sample). The method assumes linearity, namely that the impulse response is the same function independently of the input. However, soon after the first applications of the method it was found that H(r) varies quite systematically with the magnetospheric (and therefore the solar wind’s) activity level (Bargatze et al., 1985). Based on the linear filter results, a nonlinear model was suggested (Baker et al., 1990) whose response was irregular because it was dependent on the activity level. The global magnetosphere was likened to a dripping faucet (Baker et al., 1990), a strongly nonlinear oscillator where sensitive dependence on the present state allows oscillations with a large variety of amplitudes and phases to appear, so that this aperiodic “dripping” would represent global reconfigurations (substorms). At the same time there was some evidence that the AE index was self-similar in a way consistent with a chaotic system (Vassiliadis et al., 1990), although it was later shown that the selfsimilarity was biased by the autocorrelation of the high-time-resolution data (Prichard and Price, 1992). rather than, e.g., correlation between successive substorms. In these early efforts the magnetosphere was represented as an autonomous system, where the solar wind input was neglected or set to a constant, but after Prichard and Price (1992) pointed out this shortcoming, an input was included explicitly and used for prediction with nonlinear time series models (see next section). The input can be either the upstream solar wind (Vassiliadis et al., 1995; Valdivia et al., 1996) or a magnetospheric precursor like the PC index (Vassiliadis et al., 1996b). Similar studies were carried out with neural networks (e.g. Wu and Lundstedt, 1996). At the same time the nonlinear “dripping faucet” model (Baker et al., 1990) was developed into a plasmaphysics/electrodynamics’model consisting of a small number of first-principles equations which are reviewed by Klimas et al. (1996).
3. TIME
SERIES
PREDICTION:
A GEOSYNCHRONOUS
ELECTRON
FLUX
INDEX
MODEL
Time series models of geomagnetic activity mentioned above have been the nonlinear dynamics methods with the most predictive power to date (e.g. Vassiliadis et al., 1995; 1996; Valdivia et al., 1996), and they continue to be useful in space physics and space weather studies. We illustrate them by introducing an empirical geosynchronous electron flux model. Geosynchronous electron injections constitute a characteristic type of substorm signature. The available energy for them is the lobe (open-field-line) magnetic energy. Open lines, driven by the solar wind, participate in the large-scale convection of the polar cap; the amount of open flux is directly related to polar cap size and less directly to convection speed. Polar Hall currents and the convection-related electrojets and associated closure currents, can be used for approximates measure of the available energy and dissipation rate in the magnetosphere. Such currents are measured by the polar cap (PC) geomagnetic index. The PC index is constructed from the horizontal magnetic variations close to the center of the northern polar cap. Its aim is to measure Hall currents over the polar cap from which one can infer the convection speed and electric field (Troshichev et al., 1988), but the index is also affected by other polar and aurora1 current systems. In winter, in particular, Hall currents wane and the dominant geomagnetic effect is due to closure-providing field-aligned currents. Also, under Northward BZ conditions, when the polar cap contracts, a fainter effect due to the aurora1 electrojets is measured. The end result is a good correlation between PC and AL or AU which is useful in prediction of the latter indices. The geosynchronous electron flux (GEF) index is compiled from >50 keV particle measurements by five LANL geosynchronous spacecraft at a time resolution of 1 min. The index value is the highest flux measured among those spacecraft within 2 hours in Ml_T of the midnight meridian. If no spacecraft is in that region, the index value is filled in by linear interpolation in time. The construction of GEF resembles that of the AL/AU envelopes of aurora1 geomagnetic activity. In dispersionless injections GEF provides estimates of the minimum duration of the preceding current sheet thinning (by measuring the particle depletion) and estimates for injection’s onset time and intensity. These signatures are significantly reduced if the index is derived from dispersed-injection observations. The correspondence between PC and GEF is clear in intervals like those of Figure 1. For many events the relation is more quantitative: weak polar currents during the growth phase coincide with particle depletion and subsequent
200
D. Vassiliadis et al.
rapid intensification coincides with injections during expansion (Figure 2a; Pi2 pulsations provide an additional time check for substorm onset). Therefore we use the PC index as the input to a nonlinear dynamical model for geosynchronous injections. Each time series is run through a moving-average filter of 30 min. The correlation coefficient between PC and GEF over the l-month interval is 48.9% peaking at a lag time of 1 min. This time is the average time difference between the PC expansion peak and the injection peak, since these are the dominant parts of the respective time series. The correlation coefficient measures dynamical equation it is necessary to times. That can be accomplished by More accurate than MA is an ARMA
the relation between GEF(t) and PC(t) at the same time step. To obtain a express GEF(t) as a function of PC(t-jAt) at the present as well as previous a “linear prediction filter” or moving-average (MA) model similar to Eq. 1. model which includes both previous GEF and previous PC measurements: m-l
GZ(t
l-l
oi * GEF(t
+ At) = C
- iAt) + Cbj
i=O
* PC(t
- jAt)
(2)
j=O
The number of delays (terms in the summations) is m = 2 = 1, while the coefficients ai, bj can either be constants (linear ARMA model) or allowed to change with activity level (nonlinear ARMA). The ARMA filter can be mapped to an ordinary differential equation (Klimas et al., 1998), such as that for a (nonlinear) oscillator. Such an oscillator,
swinging
in discrete-time
according
representing the dissipative character of injections and of ionospheric current decay. In fact the linear ARMA model for GEF turns out to be overdamped, i.e. the oscillation period is longer than the damping time, because an impulsive PC input results in a single observed injection rather than a sequence of injections; (b) its output has a phase difference relative to the input, because of a response lag in the tail region to energy made available in the magnetosphere, plus a delay for the duration of the injection; and (c) it is in general nonlinear, i.e. the frequency, phase, and output amplitude depend on the input amplitude. The ai coefficients give the characteristic frequencies and damping rates of the system and the bj give the coupling coefficients and phase lags relative to the input (the mapping between the discrete-time ARMA and the continuous-time oscillator is described in: Klimas et al., 1997, 1998). The numerical calculation of ai, bj involves a deconvolution, similar in principle to the linear prediction filters (Clauer, 1986), but using the singular value decomposition method (Press et al., 1993). Single-step prediction uses the observation at step t+l to correct the prediction before predicting t+2; that type of prediction with a linear ARMA model reproduces the training GEF data with a correlation of 99%. Iterative prediction, however, (which would be the one to use in a realistic setting ivhere observations would not be available) gives a modest 51% correlation with observations. Nonlinear ARMA filters are expected to be more accurate.
to (2), has the following
0 a.
features:
5
IO Time
(a) it is damped,
15 [hours]
20
Fig. 2. a. Variation of GEF (in logarithmic scale) and PC indices during substorm activity in a day of the interval shown in Fig. 1. The growth (G) and expansion (E) phases are indicated as well as the Pi2 pulsation index. b. Observed and predicted GEF index (upper panel) from PC index (lower) for the first 7 days of January 1995. Note the relation of the peak values in the PC index to peaks in the injection.
NonlinearDynamicsof Space
201
The nonlinear ARMA model represents the dependence of injection dynamics with activity level, and therefore it should be calculated for each level separately, rather than averaging over all of them as it happens for a linear model. The dependence of ai, bj on activity level provides effectively a nonlinear filter. One way is to subdivide the dataset in intervals of increasing activity (as done by Bargatze et al., 1985). A frequently-used geometric method is to construct a state-input space (Vassiliadis et al., 1995) and calculate the coefficients at every point, by selecting a group of geomagnetic/solar wind conditions similar to the instantaneous level of activity with no explicit dependence on the time it occurred (the so-called “nearest neighbors”). A similar alternative is to embed the time-series data in a phase space (similar to the state-input space), divide it into regions using a grid, and calculate the coefficients in each region separately (see next section). Using the geometric representation of a state space or a phase space is advantageous, because (a) “activity level” is defined for an instant, rather than a long time interval, and the filter changes continuously (although not necessarily smoothly) with activity; (b) visualization of the data is in terms of a phase-space-like “portrait”, rather than a time series; (c) under certain conditions there is a well-defined mapping between the dynamics of the physical system and those reconstructed in the phase space (Sauer et al., 1991, and references therein). (d) The filters obtained from the state space, or state-input space for input-output systems such as PC-GEF data, can be mapped to nonlinear ODES (Klimas et al., 1997, 1998). The state-input space technique has been applied to a variety of data pairs in the solar wind-magnetosphere coupling: VBs to AL or AU (Vassiliadis et al., 1995;1996a); VBs to PC (Vassiliadis et al., 1996b); VBs to D,r (Valdivia et al., 1996). Nonlinear relations have been found between VB, and AL or D,t, and weakly nonlinear between VBs and AU.
4. FROM EMPIRICAL MODELING .TO A PHYSICAL INTERPRETATION: RING CURRENT TIME CONSTANTS FROM A HIGH-RESOLUTION DST INDEX While empirical prediction is a primary application of nonlinear dynamics methods, their prediction accuracy suggests that their coefficients can be related to the physical parameters. We will illustrate this by means of the mid-latitude geomagnetic index, D,t, which is a usual proxy for the strength of the ring current (McPherron, 1997; but see also criticism of the index construction by Campbell, 1996). The l-hour index is a standard diagnostic of the long-term average level of activity in the magnetosphere; however as will be shown below the higher-time-resolution (below: 5-min) data gives significantly more information about the storm phases and the current’s interaction with substorms. Most of the space weather damages on the ground occur during great storms, but the large dB/dt variations that cause them are localized under the substorm current system. Nonlinear dynamical methods applied to D,t have shown the capability for prediction (Valdivia et al., 1996). They can be considered extensions of a linear model (Burton et al., 1975; hereafter “B75”). That model is a linear first-order differential equation and here we extend it to a second-order ARMA model, which is the discrete-time equivalent of a second-order differential equation, a nonlinear oscillator, under certain conditions (Klimas et al., 1997). 1-I 773-l Est(t + At) = C ci * D,t(t - iAt) + C bj . VB,(t - jAt) (3)
i=o
j=o
where m = 2 = 1. The ring current dynamics changes for each storm phase, so in order to construct a nonlinear (activity-dependent) model, we divide the data in regions and calculate the coefficients for each region separately. The data are also smoothed with a 25-min running average. We construct the 3-dimensional space (Dst, dD,t(t)/dt, VB,(t)) which is the phase space for a D,t oscillator driven by an input VBs (in principle there is also a fourth dimension for dVB,/dt, but the running-average smoothing renders this derivative negligible). A spatial grid divides the phase space in rectangular regions (2D projections of the grid are shown in Figures 3a and 4a). The grid is nonuniform and along each direction (Dsr, dD,t/dt, and VBS) the spacing of the grid is such that it contains the same number of data. The coefficients ai, bj are calculated in each cell separately.
202
D. Vassiliadis et al.
The Laplace transform
of Eq. (3) gives the transfer function E,&Ll)
rl(w)
=
I-1 C
Q-j
II
m-l
p -
domain:
= I?(w)vFls(“)
l-l Jgo
in the frequency
C a;pmi i=o
bjpm-tj+l)
(4)
j=o m-l
go
(P-
PiI
where p E ezp(wAt), and p; E ezp(w;At) is one of the m roots of the polynomial in the denominator, or poles of the transfer function, while wi is one of the m characteristic frequencies of the system (Klimas et al., 1997, 1998; Vassiliadis et al., “The D,t Geomagnetic Response”. submitted to JGR. December 1998). c
.5
E s
tom 5
A
C
Y
r
0
5
-5
5
-10
% -15
L -200-150-100
-50
0
50
Main phase
‘J,t [nTl
Fig. 3. a. Projection of the (D,,, dD,/dt, VB,) phase space in its first two dimensions. The subspace is divided in 2x2 = 4 cells, indicating early SYSC phase (Compression, C), late SYSC (return, r), main phase (M), and recovery (R).
While H(w) gives the characteristic modes of oscillation/decay for D,t, its inverse Laplace transform gives the impulse response function H(kAt), k = 0, 1, 2, _.. in the time domain (cf. Eq. (1)):
.
As before H(kAt) determines the coupling wind input, and the length of the storm.
O)
to the solar
To each phase space region, or “cell”, there corresponds a set of poles and characteristic frequencies of the equation, as well as a response function, i.e. a model for that region. One can time-integrate such a model from an initial condition that starts in the cell. When the trajectory crosses over to another cell, the corresponding model is selected, etc. The combination of these local-linear models is, of course, a nonlinear model. For 3a) the the
a “coarse” nonlinear model (2x2~2 grid, Figure the impulse response functions (5) are similar to average response function, corresponding to that of B75 model (Figure 3b). The average-model decay
Fig. 3. b. Impulse response functions H(t) for the four regions of Figure 3a. The thick curve is added for comparison and gives the linear (average) response, qualitatively similar to the B75 model. Deviations from the linear model are shown for VB, < 1 mV/m (dotted lines) and VB, > 1 mV/m (solid lines). Note the presence of a slowly growing response function in the early SVSC phase and two weakly oscillatory functions in the late ST/SC phase.
time is 4.87 hours, which can be interpreted as the decay time of the ring current (for pressure-corrected Dst the time becomes 7.26 hours, comparable to B75’s 7.7 hours); the rise time is much shorter, 0.86 hours, justifying the B75 model where the rise time is taken to be 0. Most of the individual response functions are overdamped, but during the early SYSC there a growing solution and in the late SI/SC the responses become weakly oscillatory. The growing and oscillatory responses are new D,r modes that do not appear in the linear ARMA model.
203
Nonlinear Dynamics of Space
A higher phase space resolution (4x4~4 grid, Figure 4a) shows that the oscillations are present in the main phase. Figure 4b shows the phase space cells where oscillations appear; these are in the compression and main phase, but virtually absent from the recovery phase. The average frequency increases with solar wind activity, from l/2.78 to l/l .43 hours-‘. Figure 4c gives an example of the amplitude and time profile of the response for the cases marked with an ‘x’ in Figure 4b. The oscillation is also found in the time series of the _IH perturbation at all D,t magnetometers. The oscillations are identifiable in longitudinal profiles of AH (such as constructed by the method of Valdivia et al., 1999; see also Section 5). Spectrograms i.e. in the WSC and main phases. 0.61< wss c 1.21 d/m 0c WIc0.6,“V/m ..-,
?
10
E
5
p 2
0 -5
g $
-10 -15 -200-150-100
-50
0
50
o,i [nTI Fig. 4. a. Same as Figure 3a, but with the (D,,, dD,,/dt) subspace divtded in 4x4 = 16 cells. The two rightmost cells indicate early St/SC commencement (EC) and late SYSC (LC). The other fourteen are divided between main phase (M) and recovery (R).
show that the oscillations
are generally
before the Dst minimum,
*..
[nTlmV/m]
Fig. 4. b. D,, oscillations appear in certain regions (“cells”) of the phase space in Figure 4a. The oscillation frequency, Im(w), is the z axis. The four panels correspond to four increasing solar wind VB, levels. The average Im(w) increases with solar wind activity. From each graph we arbitrarily select a cell marked with ‘Ix” whose impulse response function we plot in Fig. 4c.
Time [ho”,.]
Time [hours]
Fig. 4. c. Absolute value of the impulse response functions for the cell denoted by *x’ in Figures 4a, b. The oscillations are evident for the linear (left) and logarithmic (right) scale.
There are two explanations put forward for the oscillatory response. Because it was initially found in D,t data that were not corrected for the positive geomagnetic effect of solar wind pressure increases, one of the candidates is a magnetospheric “ringing” response to solar wind pressure pulses. However, this explains only a small part of oscillations, those with periods of 5-10 min. However, most of the oscillations remain even after the data are pressure-corrected. Also the division in cells shows that the oscillations do not usually appear in the early compression (SVSC) phase, but rather at its end, as well as in the early main phase. The second explanation is that these ring current enhancements are due to substorm injections. Indeed examination of the AL index shows that its intensifications coincide with the times where the oscillations occur in dD,t/dt (Vassiliadis et al., “The D,t Geomagnetic Response”, submitted to JGR, December 1998.).
D. Vassiliadis et al.
204
The above modeling method has been based on selecting similar-activity data (nearest neighbors) using a phase space. Alternatively cne can construct a model from temporally close data (time-nearest neighbors) (Klimas et al., 1997, 1998). By selecting an interval which contains a storm or other period of interest, one can obtain model coefficients from Eq. (3) for all these times, i.e. a time-local model for a sequence of observations around the current point. In fact the coefficients not only evolve smoothly in time, but consist of nontrivial functions of the input and output. The constitutive expressions, of the form ai = f(Dst(t), dD,t(t)/dt, VB,(t)), or closure relations (Klimas et al., 1997), can be combined with Eq. (3) to produce a truly nonlinear model. Input from other intervals can be run through the model, and the predicted and observed output can be compared. 5. SPATIOTEMPORAL MODELING: PREDICTING THE LOCATION 0~ GEOMAGNETIC DISTURBANCES The magnetospheric currents that we model depend on several solar wind parameters, and simultaneously affect many locations on the ground. It is therefore natural to integrate multiple inputs and multiple outputs related to the same current system in a single nonlinear model. If the outputs are simultaneous measurements of the same physical variable at different locations, the resulting dynamical model is called a spafiotemporulmodel. A spatiotemporal model leas been developed for the mid-latitude disturbances (Valdivia et al., 1999) which are in the largest part due to ring current variations. The energy density distribution of current carriers (ions) is nonuniform in the azimuthal direction because of injections and drifting particles. Through the Dessler-Parker-Sckopke relation (McPhetron, 1997; also extension by Valdivia lza rm et al., 1999), where the geomagnetic disturbance is proportional to the totaltietgy of the particles, the energy nonuniformity produces a longitudinally-varying profile of the North-South fleld component. We characterize the mid-latitude field via its geomagnetic N-S component as measured at six magnetometers (Boulder and Tashkent added to the four standard Dpt-index magnetometers). Together with the solar wind input they are iterated in a spatiotemporal model for the r-th magnetometer measurement:
Tj$)(t + &]
=
2
.......... ..... j_iid ............. I;$-@................ ..
fz
y!i; .“.....
a(‘”
-B
....
..
..-
....
.
.
.
.....
.......
“..-“..”
_.
.....
.
.
...
..
I
s=l (6)
l-l
+
~~j%,(t
-jAt)
j=O
where m = 6 is the ntilhber of magnetometers, and 1 is the number of lags in the input VBs. The a@) coefficients are due to the correlation between magnetometer stations (increasing with the intensity and speed of propagation of a current density disturbance) and the bO)i represents the couphng to the injection function, represented here simply & the solar wind electric field (see Eq. 3). Again the coefficients can be obtained by deconvolving G,. (6j. Using all of the data set results in a linear model, and dividing the data according to geomagnetic&~ wind activity level results in a nonlinear model.
Fii. 5. Spatiotemporalprediction of midlatitude stations for the April 4, 1979 storm (DOY 114). First two panels: solar wind VB, aud pressure measured by BEE-3 and propagatedto subsolar magnetopause. Subsequent four panels: observed and predicted AH at the four standardmagnetometer stations (linear model (6): blue; nonlinear: red). Last panel: observed Dd index compared to the value obtained from the individual magnetometer predictions (linear: green; nonlinear: red).
205
Nonlinear Dynamics of Space
Given a sequence of solar wind input values, Eq. (6) can be iterated from an initial condition and can be compared to the observed AH(r) . Such comparisons show that the nonlinear model is reproducing AH(‘) than the linear spatiotemporal model (Figure 5). The predicted perturbations to form a predicted E’st which is closer to observations than that predicted from a nonlinear (Valdivia et al., 1999).
-(r) its output AH more accurate in can be compiled temporal model
Spatiotemporal prediction is currently being extended to 2D high-latitude geomagnetic predictions. Based on the predictive accuracy of nonlinear models for the AL/AU indices, similar in structure to Eqs. (2, 3), one can extend them to spatiotemporal models similar to Eq. (6) for local perturbations measured by one or more A first approach analyzed Canopus’s east magnetometer arrays as a function of MLT and invariant latitude. meridional magnetometer array in conjunction with solar wind input from WIND at the L1 point (Valdivia et al., 1998). Here, too, the nonlinear spatial model is more accurate than the linear one. Spatiotemporal modeling will be extended to more than one magnetometer arrays. Yellow Knife, November
Novamkr 1-3.1078
E
Fig. 6. Coupling of nonlinear dynamical models and a ballistic propagation scheme as part of the ionospheric SWIFT model. First three panels: propagation of BEE-3 data to the subsolar magnetopause and comparison to IMP-8 near-Earth measurements. The last two panels are AL predictions (grey lines) from the BEE-3 input (IMP-8 gives similar results) and from the PC index.
6. COMBINATION OF NONLINEAR DYNAMICS TOWARDS GEOMAGNETIC AND IONOSPHERIC
5.1979
from both AL and Ati
Fig. 7. Effect of AU parameter on the SWIFT model output. At around 3 am on November 5, 1978, Yellowknife is west of the Harang discontinuity so the effect of the eastward_electrojet is significant. Using only the nonlinear-predicted AL in the geomagnetic part of SWIFT gives the wrong sign for the ground disturbance (dash-dotted line). When the predicted E is included, the prediction is improved significantly (dotted line).
WITH OTHER FORECASTING
MODELS
Nonlinear dynamical models are usually combined with other physical or empirical models. For example, the solar wind input for geomagnetic time series models is the measurement at L1 processed by a propagation program which simulates the transport to the subsolar magnetopause. Similarly when a nonlinear model is sufficiently accurate, its output can be used by other models as a proxy for the modeled magnetospheric disturbance. In this way more comprehensive models can be built. An example for such a combination of models for solar wind-magnetosphere -ionosphere coupling is the Space Weather Ionospheric Forecasting Technologies (SWIFT) model. The solar wind drives nonlinear dynamics and neural networks models, and in addition the solar wind input and predicted indices drive empirical ionospheric
206
D. Vassiliadis et nl.
models which give the electric-field patterns, currents, and Joule dissipation, and geomagnetic disturbance as a function of invariant latitude and MLT. An early test of SWIFT shows the propagation and index-prediction parts (Figure 6). The solar wind input is measured at L1 by ISEE- and propagated to the subsolar magnetopause where it is compared with IMP-g. There is fair agreement between the model and observed solar wind. Subsequently both the propagated solar wind and the PC index are used in turn as inputs to predict AL. Other tests have been carried out to measure the individual models’ sensitivity to perturbations to or uncertainties in the initial state. After the indices have been calculated, they are used in conjunction with the solar wind input to obtain the ionospheric electrodynamic variables. Since different current systems are parametrized through indices and other model variables, it is possible to optimize the model further by weighting the relative contribution of each current system to increase agreement with training data. Figure 7 shows the effect of adding the parameter of the eastward electrojet, the AU index, to the model optimization. A sequence of intensifications, starting at O4:OOUT, November 4 (when the measuring station of Yellowknife is in the dusk section), is predicted with the wrong sign. Inclusion of the AU parameter reduces the error by at least 50 and in some cases more than 150 nT. Index predictions are produced daily based on the solar wind input (originally from WIND and since May 1998 from ACE) and the PC index. A homepage has been set up for AL and AU, http:/Aepgst.gsfc.nasa.gov/ people/vassi/htmls/alprediction.html. For data-derived models it is important to maintain a database of intervals representative of as many activity levels and profiles as possible. Recently we have created an on-line database for AL/AU predictions during magnetic cloud passages in 1997 at h~~p://epgs~.gsfc..nasa.gov/people/vassilhhnls/97/even~s97.html. Remarkably the predictions are less accurate compared to the standard, non-cloud predictions. The explanation is that during the cloud events, the amplitude is larger and the variability is smaller than in the training intervals. Since the training dataset does not contain cloud conditions, it is difficult to predict AL/AU correctly. We plan to test this hypothesis by adding the cloud intervals to the database and repeating the runs. 7. SUMMARY Nonlinear dynamics techniques &e an essential set of data analysis methods which cover the middle ground between single-event and statistical analysis. We discussed a sequence of methods to quantitatively predict and physically The complexity of the model building can increase in a commensurate model a space physical environment. way with the variability of the data sets: From single-time-series prediction for scalar indices, we have gone to methods of physical interpretation of the models (e.g. Klimas et al., 1997), where we measure characteristic time scales (growth/decay/oscillation) and coupling strengths. The instantaneous state of the current system can be more accurately defined if measurements at several locations are available, from which one can construct spatiotemporal models and predictions. Finally if the accuracy of the prediction is sufficiently high, the model can be used in conjunction with other empirical or physical models. The evaluation of nonlinear dynamics models in terms of error statistics provides a good guide in two areas not covered in this paper: model verification and metric development. In developing space weather models, it is essential to verify the predictive capability by using out-of-sample (cross-validation) testing, and comparison between the candidate model and reference models or competing models to extract the prediction skill, i.e. the prediction capability relative to other models (issues discussed by Doggett, 1996). These are standard topics for nonlinear model development (e.g. the reference model is the linear model). Similarly the error analysis in nonlinear dynamics models (Vassiliadis and Klimas, 1996) can also be used as a guide to develop metrics necessary for progress evaluation in space weather forecasting (Wolf et al., Metrics for the National Space Weather Program, January 1999). 8. ACKNOWLEDGMENTS The following individuals and/or organizations were essential in providing the data we used in this study: Hightime-resolution D,t index: Y. Kamide, R. Nakamura, S. Sharma. GEF index: G. D. Reeves. Pi2 pulsation index: T. Iyemori, M. Nose’. ACE and WIND magnetic field and plasma data: Spacecraft PIs and NOAA/SEC. PC index: E. Friis-Christensen and S. Vennerstroem. Support from NASA/SR&T and NSF/GEM grants is gratefully acknowledged.
Nonlinear Dynamics of Space
REFERENCES Baker, D. N., A. J. Klimas, R. L. McPherron, and J. Buchner, The Evolution from Weak to Strong Geomagnetic Activity: an Interpretation in Terms of Deterministic Chaos, Geophys. Res. Len., 17,41 (1990). Bargatze, L. F.. D. N. Baker, R. L. McPherron, and E. W. Hones, Jr., Magnetospheric Impulse Response for Many Levels of Geomagnetic Activity, J. Geophys. Res., 90,6387 (1985). Burton, R. K., R. L. McPherron, and C. T. Russell, An Empirical Relationship Between Interplanetary Conditions and Dst, J. Geophys. Res., 80,4204 (1975). Campbell, W. H., Geomagnetic Storms, the Dst Ring-Current Myth, and Lognormal Distributions, f. Arm. Terr. Phys., 58, 1171 (1996). Clauer, C. R., The Technique of Linear Prediction Filters Applied to Studies of Solar Wind-Magnetosphere Coupling, in Solar Wind-Magnetosphere Coupling, edited by Y. Kamide and J. Slavin, pp. 39-57, Terra Scientific Publ. Co., Tokyo (1986). Doggett, K. A., (ed.), The Evaluation of Space Weather Forecasts. Proceedings of a Workshop in Boulder, Colorado, June 19-21, 1996, NOAA/SEC, Boulder, CO (1996). Horton W., and I. Doxas, A Low-Dimensional Dynamical Model for the Solar Wind-Driven Geotail-Ionosphere System J. Geophys. Res., 103,456l (1998). Klimas, A. J., D. Vassiliadis, D. N. Baker, and D. A. Roberts, The Organized Nonlinear Dynamics of the Magnetosphere, J. Geophys. Res., 101, 13089 (1996). Klimas, A. J., D. Vassiliadis, and D. N. Baker, Data-Derived Analogues of the Magnetospheric Dynamics, J. Geophys. Res., 102,26993 (1997). Klimas, A. J., D. Vassiliadis, and D. N. Baker, Dst Index Prediction Using Data-Derived Analogues of the Magnetospheric Dynamics, J. Geophys. Res., 103,20435 (1998). McPherron, R. L., The Role of Substorms in the Generation of Magnetic Storms, in Magnetic Storms, edited by B. T. Tsurutani, W. D. Gonzalez, Y. Kamide, and J. K. Arballo. pp. 131-147, American Geophysical Union, Washington, DC (1997). Potemra, T. A. (ed.), Magnetospheric Currents, American Geophysical Union, Washington, DC (1984). Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. V. Vetterling, Numerical Recipes: the Art of Scientific Computing (second ed.), Cambridge University Press, Cambridge, UK (1993). Prichard, D., and C. P. Price, Spurious Dimension Estimates from Time Series of Geomagnetic Indices, Geophys. Res. L&t., 19, 1623 (1992). Sauer, T., J. A. Yorke, and M. Casdagli, Embedology, J. Star. Phys., 65,579 (1991). Siscoe, G., Manifestations of global coherence, Eos Transactions AGU, Fall Meeting Suppl., 501 (1993). Troshichev, 0. A., V. G. Andrezen, S. Vennerstrom, and E. Friis-Christensen, Magnetic Activity in the Polar-Cap A New Index, Planet. Space Sci., 33, 1095 (1988). Valdivia, J. A., A. S. Sharma, and K. Papadopoulos, Prediction of Magnetic Storms by Nonlinear Models, Geophys. Res. L-e%.,23, 2899 (1996). Valdivia, J. A., D. Vassiliadis, A. J. Klimas, and A. S. Sharma, The Electrojet Currents: Understanding Their SpatioTemporal Multivariate Properties, in Proceedings of the Fourth International Conference on Substorms, edited by Y. Kamide and S. Kokubun, pp. 669-672, Terra Scientific Publ. Co., Tokyo (1998). Vassiliadis, D., and A. J. Klimas, Evaluation of Nonlinear Systems for Geomagnetic Index Forecasting and Modeling, in The Evaluation of Space Weather Forecasts. Proceedings of a Workshop in Boulder, Colorado, June 19-21, 1996, edited by K. A. Doggett, pp. 29-33, NOAA/SEC, Boulder, CO (1996). Vassiliadis, D., A. S. Sharma, T. E. Eastman, and K. Papadopoulos, Low-Dimensional Chaos in Magnetospheric Activity from AE Time Series, Geophys. Res. Zen., 17, 1841 (1990). Vassiliadis, D., A. J. Klimas, D. N. Baker, and D. A. Roberts, A Description of the Solar Wind-Magnetosphere Coupling Based on Nonlinear Filters, J. Geophys. Res., 100,3495 (1995). Vassiliadis, D., A. J. Klimas, D. N. Baker, and D. A. Roberts, The Nonlinearity of Models of the vBsouth - AL Coupling, J. Geophys. Res., 101, 19779 (1996a). Vassiliadis, D., V. Angelopoulos, D. N. Baker, and A. J. Klimas, The Relation Between Northern Polar Cap and Aurora1 Electrojet Geomagnetic Indices in the Wintertime, Geophys. Res. Lett., 23,278 1 (1996b). Wu, J.-G., and H. Lundstedt, Prediction of Geomagnetic Storms from Solar Wind Data Using Elman Recurrent Neural Networks, Geophys. Res. Zett., 23,3 19 (1996).
207
Adv. SpaceRes. Vol. 26, No. 1, pp. 197-207,200O 0 2000 COSPAR. Published bv Elsevier Science Ltd. All rightsreserved Printedin &eat Britain 0273-1177/00 $20.00 + 0.00 Pll: SO273-1177(99)01050-9
THE NONLINEAR DYNAMICS OF SPACE WEATHER
D. Vassiliadisl,
A. J. Klimas2, J. A. Valdivia3, D. N. Baker4
’ Universities Space Research Association, at: NASA/GSFC, Code 692, Greenbelt, MD 20771, USA 2Luboratory for Extraterrestrial Physics, NASA/GSFC, Greenbelt, MD 20771, USA 3National Research Council, at: NASA/GSFC, Greenbelt, MD 20771, USA ‘Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80309-0590, USA
ABSTRACT Studies of the nonlinear magnetospheric dynamics have led to several directions useful in understanding space physics processes, in particular those related to magnetospheric currents, and making space weather forecasts possible. Four such directions are identified: (a) empirical time series prediction with nonlinear autoregressive moving-average (ARMA) models, for which an example is given in terms of a geosynchronous electron flux index model. (b) Measurement of physical properties of the currents from the coefficients of the ARMA models, such as characteristic time scales and coupling strengths. Using a Dst index model we measure the ring current decay time as a function of storm phase and activity, and identify new oscillation time scales correlating with the substorm injection activity. (c) Spatiotemporal nonlinear modeling predicts the amplitude and location of the disturbance as a function of space, as well as its time evolution. An example is given in terms of predicting the longitudinal and temporal profile of midlatitude geomagnetic disturbances. (d) Nonlinear dynamical models can be coupled to other physical or empirical approaches to build comprehensive space weather models. SWIFT, an ionospheric space weather model, is discussed as an example. Other applications of nonlinear dynamics are briefly discussed. 0 2000 COSPAR. Published by Elsevier Science Ltd. 1. INTRODUCTION Nonlinear dynamics in space physics has grown steadily in scope and depth. We discuss this growth in terms of four directions whose results have useful implications for theoretical and applied physics, in particular space weather. While earlier research focuses on classification, empirical prediction, and first-principles modeling of geomagnetic time series (review: Klimas et al., 1996), recently three more directions have been developed: closed-form timeseries models which are data-derived, but in addition have a physical interpretation; spatial, in addition to temporal, prediction of space-related disturbances, and modeling from multiple measurement probes in the same region, or in several regions; and combinations of the nonlinear dynamics to other, physical or empirical, models for more comprehensive space weather forecasting. We discuss these directions in the above order. The starting point is a brief review of basic nonlinear dynamics concepts and some of the early results in Section 2. The following section discusses empirical time series modeling and illustrates it with the prediction of a newly proposed geosynchronous electron flux index. The next two sections (4 and 5) are devoted to physical time series modeling and spatiotemporal modeling. Following that, Sec. 6 discusses the synergy between nonlinear dynamical and other models for space weather applications, and finally ways of making predictions and other results available on-line.
197
D. Vassiliadis et al.
198
2. WHENCE
NONLINEAR
DYNAMICS?
A basic goal of nonlinear dynamics as applied to space physics is to model current systems in the magnetosphere and ionosphere and the geomagnetic disturbances they cause. The configuration and intensity of the basic current systems determine the global magnetospheric structure and activity in the predominantly low-p plasma (Potemra, 1984). The justification for applying nonlinear dynamics methods is the conspicuous coherence of these currents, and other aspects of magnetospheric activity. However, a theoretical understanding of how self-organization arises in plasmas for realistic ranges of plasma parameters is missing. Regarding the currents, their coherent nature is characterized by their (a) being permanent and stable features, and (b) by their having large interaction lengths and fast interaction timescales. The ring, electrojet, magnetopause, and tail currents are permanent features of the magnetosphere. Following a small-scale disturbance they rapidly return to their equilibrium values, except if a stability threshold is‘exceeded (e.g. in the case of the tail current). Because of their long-range magnetic effect, and the rapid signal speed (VA 2 L/T where L is a scale of the magnetospheric subsystem and T is a time scale for current growth), the magnetosphere is more organized spatially and responds faster than, e.g., the neutral atmosphere, a fluid of much smaller volume. An example of coherent activity in various magnetospheric locations is shown in Figure 1. Solar wind driving (from a high-speed stream at the end of solar cycle 22) excites simultaneous geomagnetic disturbances in polar, auroral, and midlatitude locations (whose activity is normally characterized by quite different time scales) while substorm injections often’ produce high geosynchronous particle fluxes. (c) Additional observational evidence about magnetospheric coherence has been discussed by Siscoe (1993). (d) Finally the predictability of magnetospheric-current processes, as will be reviewed below, is another body of evidence for the coherence. Hence the observed magnetospheric collective response is a direct argument for using nonlinear dynamics methods. The underlying assumption for most methods is that the spatial characteristics of the system (here: the geomagnetic activity) are due to a small number of spatial modes (corresponding to the geomagnetic effect of distinct current systems), and the complex activity can be described by the interaction, probably nonlinear, of these modes. The need for nonlinear data analysis became obvious after the linear moving-average (MA) filters were applied (Clauer, 1986). In examining the effect of solar wind driving, represented by the VB, electric field, on the energization of the magnetosphere, represented by the AL index, linear prediction filters were used, and they showed the first predictive capability
T G(t)
=
H(r)
* [VB,](t
- T)dT
J 0
Z(t Fig. 1. Coherent magnetospheric response during the month of 1995 to a solar wind high-speed stream (not shown). From top to bottom: Solar wind input; Geosynchronous orbit: index for the energetic electron flux (GEF); logro(GEF) showing electron depletions due to growth phase dipolarizations; Polar latitudes: index measuring the ground geomagnetic disturbance near the cap center; Mid-latitudes: Dst geomagnetic index.
(1)
t At) = 2
H(jAt)[VB,](t
- jAt)
j=O
January
where the first equation is for continuous ond for discrete-time data of resolution denotes the filter output.
and the secAt, and A^L
Now the filter’s impulse response function, H(r), where T is the time lag relative to the onset of solar
199
Nonlinear Dynamics of Space
wind input, is calculated from the VB, and AL data by inverting (1). Then the filter can be convolved with the same solar wind data to reproduce AL (in-sample), or to different solar wind data to predict the corresponding AL (outof-sample). The method assumes linearity, namely that the impulse response is the same function independently of the input. However, soon after the first applications of the method it was found that H(r) varies quite systematically with the magnetospheric (and therefore the solar wind’s) activity level (Bargatze et al., 1985). Based on the linear filter results, a nonlinear model was suggested (Baker et al., 1990) whose response was irregular because it was dependent on the activity level. The global magnetosphere was likened to a dripping faucet (Baker et al., 1990), a strongly nonlinear oscillator where sensitive dependence on the present state allows oscillations with a large variety of amplitudes and phases to appear, so that this aperiodic “dripping” would represent global reconfigurations (substorms). At the same time there was some evidence that the AE index was self-similar in a way consistent with a chaotic system (Vassiliadis et al., 1990), although it was later shown that the selfsimilarity was biased by the autocorrelation of the high-time-resolution data (Prichard and Price, 1992). rather than, e.g., correlation between successive substorms. In these early efforts the magnetosphere was represented as an autonomous system, where the solar wind input was neglected or set to a constant, but after Prichard and Price (1992) pointed out this shortcoming, an input was included explicitly and used for prediction with nonlinear time series models (see next section). The input can be either the upstream solar wind (Vassiliadis et al., 1995; Valdivia et al., 1996) or a magnetospheric precursor like the PC index (Vassiliadis et al., 1996b). Similar studies were carried out with neural networks (e.g. Wu and Lundstedt, 1996). At the same time the nonlinear “dripping faucet” model (Baker et al., 1990) was developed into a plasmaphysics/electrodynamics’model consisting of a small number of first-principles equations which are reviewed by Klimas et al. (1996).
3. TIME
SERIES
PREDICTION:
A GEOSYNCHRONOUS
ELECTRON
FLUX
INDEX
MODEL
Time series models of geomagnetic activity mentioned above have been the nonlinear dynamics methods with the most predictive power to date (e.g. Vassiliadis et al., 1995; 1996; Valdivia et al., 1996), and they continue to be useful in space physics and space weather studies. We illustrate them by introducing an empirical geosynchronous electron flux model. Geosynchronous electron injections constitute a characteristic type of substorm signature. The available energy for them is the lobe (open-field-line) magnetic energy. Open lines, driven by the solar wind, participate in the large-scale convection of the polar cap; the amount of open flux is directly related to polar cap size and less directly to convection speed. Polar Hall currents and the convection-related electrojets and associated closure currents, can be used for approximates measure of the available energy and dissipation rate in the magnetosphere. Such currents are measured by the polar cap (PC) geomagnetic index. The PC index is constructed from the horizontal magnetic variations close to the center of the northern polar cap. Its aim is to measure Hall currents over the polar cap from which one can infer the convection speed and electric field (Troshichev et al., 1988), but the index is also affected by other polar and aurora1 current systems. In winter, in particular, Hall currents wane and the dominant geomagnetic effect is due to closure-providing field-aligned currents. Also, under Northward BZ conditions, when the polar cap contracts, a fainter effect due to the aurora1 electrojets is measured. The end result is a good correlation between PC and AL or AU which is useful in prediction of the latter indices. The geosynchronous electron flux (GEF) index is compiled from >50 keV particle measurements by five LANL geosynchronous spacecraft at a time resolution of 1 min. The index value is the highest flux measured among those spacecraft within 2 hours in Ml_T of the midnight meridian. If no spacecraft is in that region, the index value is filled in by linear interpolation in time. The construction of GEF resembles that of the AL/AU envelopes of aurora1 geomagnetic activity. In dispersionless injections GEF provides estimates of the minimum duration of the preceding current sheet thinning (by measuring the particle depletion) and estimates for injection’s onset time and intensity. These signatures are significantly reduced if the index is derived from dispersed-injection observations. The correspondence between PC and GEF is clear in intervals like those of Figure 1. For many events the relation is more quantitative: weak polar currents during the growth phase coincide with particle depletion and subsequent
200
D. Vassiliadis et al.
rapid intensification coincides with injections during expansion (Figure 2a; Pi2 pulsations provide an additional time check for substorm onset). Therefore we use the PC index as the input to a nonlinear dynamical model for geosynchronous injections. Each time series is run through a moving-average filter of 30 min. The correlation coefficient between PC and GEF over the l-month interval is 48.9% peaking at a lag time of 1 min. This time is the average time difference between the PC expansion peak and the injection peak, since these are the dominant parts of the respective time series. The correlation coefficient measures dynamical equation it is necessary to times. That can be accomplished by More accurate than MA is an ARMA
the relation between GEF(t) and PC(t) at the same time step. To obtain a express GEF(t) as a function of PC(t-jAt) at the present as well as previous a “linear prediction filter” or moving-average (MA) model similar to Eq. 1. model which includes both previous GEF and previous PC measurements: m-l
GZ(t
l-l
oi * GEF(t
+ At) = C
- iAt) + Cbj
i=O
* PC(t
- jAt)
(2)
j=O
The number of delays (terms in the summations) is m = 2 = 1, while the coefficients ai, bj can either be constants (linear ARMA model) or allowed to change with activity level (nonlinear ARMA). The ARMA filter can be mapped to an ordinary differential equation (Klimas et al., 1998), such as that for a (nonlinear) oscillator. Such an oscillator,
swinging
in discrete-time
according
representing the dissipative character of injections and of ionospheric current decay. In fact the linear ARMA model for GEF turns out to be overdamped, i.e. the oscillation period is longer than the damping time, because an impulsive PC input results in a single observed injection rather than a sequence of injections; (b) its output has a phase difference relative to the input, because of a response lag in the tail region to energy made available in the magnetosphere, plus a delay for the duration of the injection; and (c) it is in general nonlinear, i.e. the frequency, phase, and output amplitude depend on the input amplitude. The ai coefficients give the characteristic frequencies and damping rates of the system and the bj give the coupling coefficients and phase lags relative to the input (the mapping between the discrete-time ARMA and the continuous-time oscillator is described in: Klimas et al., 1997, 1998). The numerical calculation of ai, bj involves a deconvolution, similar in principle to the linear prediction filters (Clauer, 1986), but using the singular value decomposition method (Press et al., 1993). Single-step prediction uses the observation at step t+l to correct the prediction before predicting t+2; that type of prediction with a linear ARMA model reproduces the training GEF data with a correlation of 99%. Iterative prediction, however, (which would be the one to use in a realistic setting ivhere observations would not be available) gives a modest 51% correlation with observations. Nonlinear ARMA filters are expected to be more accurate.
to (2), has the following
0 a.
features:
5
IO Time
(a) it is damped,
15 [hours]
20
Fig. 2. a. Variation of GEF (in logarithmic scale) and PC indices during substorm activity in a day of the interval shown in Fig. 1. The growth (G) and expansion (E) phases are indicated as well as the Pi2 pulsation index. b. Observed and predicted GEF index (upper panel) from PC index (lower) for the first 7 days of January 1995. Note the relation of the peak values in the PC index to peaks in the injection.
NonlinearDynamicsof Space
201
The nonlinear ARMA model represents the dependence of injection dynamics with activity level, and therefore it should be calculated for each level separately, rather than averaging over all of them as it happens for a linear model. The dependence of ai, bj on activity level provides effectively a nonlinear filter. One way is to subdivide the dataset in intervals of increasing activity (as done by Bargatze et al., 1985). A frequently-used geometric method is to construct a state-input space (Vassiliadis et al., 1995) and calculate the coefficients at every point, by selecting a group of geomagnetic/solar wind conditions similar to the instantaneous level of activity with no explicit dependence on the time it occurred (the so-called “nearest neighbors”). A similar alternative is to embed the time-series data in a phase space (similar to the state-input space), divide it into regions using a grid, and calculate the coefficients in each region separately (see next section). Using the geometric representation of a state space or a phase space is advantageous, because (a) “activity level” is defined for an instant, rather than a long time interval, and the filter changes continuously (although not necessarily smoothly) with activity; (b) visualization of the data is in terms of a phase-space-like “portrait”, rather than a time series; (c) under certain conditions there is a well-defined mapping between the dynamics of the physical system and those reconstructed in the phase space (Sauer et al., 1991, and references therein). (d) The filters obtained from the state space, or state-input space for input-output systems such as PC-GEF data, can be mapped to nonlinear ODES (Klimas et al., 1997, 1998). The state-input space technique has been applied to a variety of data pairs in the solar wind-magnetosphere coupling: VBs to AL or AU (Vassiliadis et al., 1995;1996a); VBs to PC (Vassiliadis et al., 1996b); VBs to D,r (Valdivia et al., 1996). Nonlinear relations have been found between VB, and AL or D,t, and weakly nonlinear between VBs and AU.
4. FROM EMPIRICAL MODELING .TO A PHYSICAL INTERPRETATION: RING CURRENT TIME CONSTANTS FROM A HIGH-RESOLUTION DST INDEX While empirical prediction is a primary application of nonlinear dynamics methods, their prediction accuracy suggests that their coefficients can be related to the physical parameters. We will illustrate this by means of the mid-latitude geomagnetic index, D,t, which is a usual proxy for the strength of the ring current (McPherron, 1997; but see also criticism of the index construction by Campbell, 1996). The l-hour index is a standard diagnostic of the long-term average level of activity in the magnetosphere; however as will be shown below the higher-time-resolution (below: 5-min) data gives significantly more information about the storm phases and the current’s interaction with substorms. Most of the space weather damages on the ground occur during great storms, but the large dB/dt variations that cause them are localized under the substorm current system. Nonlinear dynamical methods applied to D,t have shown the capability for prediction (Valdivia et al., 1996). They can be considered extensions of a linear model (Burton et al., 1975; hereafter “B75”). That model is a linear first-order differential equation and here we extend it to a second-order ARMA model, which is the discrete-time equivalent of a second-order differential equation, a nonlinear oscillator, under certain conditions (Klimas et al., 1997). 1-I 773-l Est(t + At) = C ci * D,t(t - iAt) + C bj . VB,(t - jAt) (3)
i=o
j=o
where m = 2 = 1. The ring current dynamics changes for each storm phase, so in order to construct a nonlinear (activity-dependent) model, we divide the data in regions and calculate the coefficients for each region separately. The data are also smoothed with a 25-min running average. We construct the 3-dimensional space (Dst, dD,t(t)/dt, VB,(t)) which is the phase space for a D,t oscillator driven by an input VBs (in principle there is also a fourth dimension for dVB,/dt, but the running-average smoothing renders this derivative negligible). A spatial grid divides the phase space in rectangular regions (2D projections of the grid are shown in Figures 3a and 4a). The grid is nonuniform and along each direction (Dsr, dD,t/dt, and VBS) the spacing of the grid is such that it contains the same number of data. The coefficients ai, bj are calculated in each cell separately.
202
D. Vassiliadis et al.
The Laplace transform
of Eq. (3) gives the transfer function E,&Ll)
rl(w)
=
I-1 C
Q-j
II
m-l
p -
domain:
= I?(w)vFls(“)
l-l Jgo
in the frequency
C a;pmi i=o
bjpm-tj+l)
(4)
j=o m-l
go
(P-
PiI
where p E ezp(wAt), and p; E ezp(w;At) is one of the m roots of the polynomial in the denominator, or poles of the transfer function, while wi is one of the m characteristic frequencies of the system (Klimas et al., 1997, 1998; Vassiliadis et al., “The D,t Geomagnetic Response”. submitted to JGR. December 1998). c
.5
E s
tom 5
A
C
Y
r
0
5
-5
5
-10
% -15
L -200-150-100
-50
0
50
Main phase
‘J,t [nTl
Fig. 3. a. Projection of the (D,,, dD,/dt, VB,) phase space in its first two dimensions. The subspace is divided in 2x2 = 4 cells, indicating early SYSC phase (Compression, C), late SYSC (return, r), main phase (M), and recovery (R).
While H(w) gives the characteristic modes of oscillation/decay for D,t, its inverse Laplace transform gives the impulse response function H(kAt), k = 0, 1, 2, _.. in the time domain (cf. Eq. (1)):
.
As before H(kAt) determines the coupling wind input, and the length of the storm.
O)
to the solar
To each phase space region, or “cell”, there corresponds a set of poles and characteristic frequencies of the equation, as well as a response function, i.e. a model for that region. One can time-integrate such a model from an initial condition that starts in the cell. When the trajectory crosses over to another cell, the corresponding model is selected, etc. The combination of these local-linear models is, of course, a nonlinear model. For 3a) the the
a “coarse” nonlinear model (2x2~2 grid, Figure the impulse response functions (5) are similar to average response function, corresponding to that of B75 model (Figure 3b). The average-model decay
Fig. 3. b. Impulse response functions H(t) for the four regions of Figure 3a. The thick curve is added for comparison and gives the linear (average) response, qualitatively similar to the B75 model. Deviations from the linear model are shown for VB, < 1 mV/m (dotted lines) and VB, > 1 mV/m (solid lines). Note the presence of a slowly growing response function in the early SVSC phase and two weakly oscillatory functions in the late ST/SC phase.
time is 4.87 hours, which can be interpreted as the decay time of the ring current (for pressure-corrected Dst the time becomes 7.26 hours, comparable to B75’s 7.7 hours); the rise time is much shorter, 0.86 hours, justifying the B75 model where the rise time is taken to be 0. Most of the individual response functions are overdamped, but during the early SYSC there a growing solution and in the late SI/SC the responses become weakly oscillatory. The growing and oscillatory responses are new D,r modes that do not appear in the linear ARMA model.
203
Nonlinear Dynamics of Space
A higher phase space resolution (4x4~4 grid, Figure 4a) shows that the oscillations are present in the main phase. Figure 4b shows the phase space cells where oscillations appear; these are in the compression and main phase, but virtually absent from the recovery phase. The average frequency increases with solar wind activity, from l/2.78 to l/l .43 hours-‘. Figure 4c gives an example of the amplitude and time profile of the response for the cases marked with an ‘x’ in Figure 4b. The oscillation is also found in the time series of the _IH perturbation at all D,t magnetometers. The oscillations are identifiable in longitudinal profiles of AH (such as constructed by the method of Valdivia et al., 1999; see also Section 5). Spectrograms i.e. in the WSC and main phases. 0.61< wss c 1.21 d/m 0c WIc0.6,“V/m ..-,
?
10
E
5
p 2
0 -5
g $
-10 -15 -200-150-100
-50
0
50
o,i [nTI Fig. 4. a. Same as Figure 3a, but with the (D,,, dD,,/dt) subspace divtded in 4x4 = 16 cells. The two rightmost cells indicate early St/SC commencement (EC) and late SYSC (LC). The other fourteen are divided between main phase (M) and recovery (R).
show that the oscillations
are generally
before the Dst minimum,
*..
[nTlmV/m]
Fig. 4. b. D,, oscillations appear in certain regions (“cells”) of the phase space in Figure 4a. The oscillation frequency, Im(w), is the z axis. The four panels correspond to four increasing solar wind VB, levels. The average Im(w) increases with solar wind activity. From each graph we arbitrarily select a cell marked with ‘Ix” whose impulse response function we plot in Fig. 4c.
Time [ho”,.]
Time [hours]
Fig. 4. c. Absolute value of the impulse response functions for the cell denoted by *x’ in Figures 4a, b. The oscillations are evident for the linear (left) and logarithmic (right) scale.
There are two explanations put forward for the oscillatory response. Because it was initially found in D,t data that were not corrected for the positive geomagnetic effect of solar wind pressure increases, one of the candidates is a magnetospheric “ringing” response to solar wind pressure pulses. However, this explains only a small part of oscillations, those with periods of 5-10 min. However, most of the oscillations remain even after the data are pressure-corrected. Also the division in cells shows that the oscillations do not usually appear in the early compression (SVSC) phase, but rather at its end, as well as in the early main phase. The second explanation is that these ring current enhancements are due to substorm injections. Indeed examination of the AL index shows that its intensifications coincide with the times where the oscillations occur in dD,t/dt (Vassiliadis et al., “The D,t Geomagnetic Response”, submitted to JGR, December 1998.).
D. Vassiliadis et al.
204
The above modeling method has been based on selecting similar-activity data (nearest neighbors) using a phase space. Alternatively cne can construct a model from temporally close data (time-nearest neighbors) (Klimas et al., 1997, 1998). By selecting an interval which contains a storm or other period of interest, one can obtain model coefficients from Eq. (3) for all these times, i.e. a time-local model for a sequence of observations around the current point. In fact the coefficients not only evolve smoothly in time, but consist of nontrivial functions of the input and output. The constitutive expressions, of the form ai = f(Dst(t), dD,t(t)/dt, VB,(t)), or closure relations (Klimas et al., 1997), can be combined with Eq. (3) to produce a truly nonlinear model. Input from other intervals can be run through the model, and the predicted and observed output can be compared. 5. SPATIOTEMPORAL MODELING: PREDICTING THE LOCATION 0~ GEOMAGNETIC DISTURBANCES The magnetospheric currents that we model depend on several solar wind parameters, and simultaneously affect many locations on the ground. It is therefore natural to integrate multiple inputs and multiple outputs related to the same current system in a single nonlinear model. If the outputs are simultaneous measurements of the same physical variable at different locations, the resulting dynamical model is called a spafiotemporulmodel. A spatiotemporal model leas been developed for the mid-latitude disturbances (Valdivia et al., 1999) which are in the largest part due to ring current variations. The energy density distribution of current carriers (ions) is nonuniform in the azimuthal direction because of injections and drifting particles. Through the Dessler-Parker-Sckopke relation (McPhetron, 1997; also extension by Valdivia lza rm et al., 1999), where the geomagnetic disturbance is proportional to the totaltietgy of the particles, the energy nonuniformity produces a longitudinally-varying profile of the North-South fleld component. We characterize the mid-latitude field via its geomagnetic N-S component as measured at six magnetometers (Boulder and Tashkent added to the four standard Dpt-index magnetometers). Together with the solar wind input they are iterated in a spatiotemporal model for the r-th magnetometer measurement:
Tj$)(t + &]
=
2
.......... ..... j_iid ............. I;$-@................ ..
fz
y!i; .“.....
a(‘”
-B
....
..
..-
....
.
.
.
.....
.......
“..-“..”
_.
.....
.
.
...
..
I
s=l (6)
l-l
+
~~j%,(t
-jAt)
j=O
where m = 6 is the ntilhber of magnetometers, and 1 is the number of lags in the input VBs. The a@) coefficients are due to the correlation between magnetometer stations (increasing with the intensity and speed of propagation of a current density disturbance) and the bO)i represents the couphng to the injection function, represented here simply & the solar wind electric field (see Eq. 3). Again the coefficients can be obtained by deconvolving G,. (6j. Using all of the data set results in a linear model, and dividing the data according to geomagnetic&~ wind activity level results in a nonlinear model.
Fii. 5. Spatiotemporalprediction of midlatitude stations for the April 4, 1979 storm (DOY 114). First two panels: solar wind VB, aud pressure measured by BEE-3 and propagatedto subsolar magnetopause. Subsequent four panels: observed and predicted AH at the four standardmagnetometer stations (linear model (6): blue; nonlinear: red). Last panel: observed Dd index compared to the value obtained from the individual magnetometer predictions (linear: green; nonlinear: red).
205
Nonlinear Dynamics of Space
Given a sequence of solar wind input values, Eq. (6) can be iterated from an initial condition and can be compared to the observed AH(r) . Such comparisons show that the nonlinear model is reproducing AH(‘) than the linear spatiotemporal model (Figure 5). The predicted perturbations to form a predicted E’st which is closer to observations than that predicted from a nonlinear (Valdivia et al., 1999).
-(r) its output AH more accurate in can be compiled temporal model
Spatiotemporal prediction is currently being extended to 2D high-latitude geomagnetic predictions. Based on the predictive accuracy of nonlinear models for the AL/AU indices, similar in structure to Eqs. (2, 3), one can extend them to spatiotemporal models similar to Eq. (6) for local perturbations measured by one or more A first approach analyzed Canopus’s east magnetometer arrays as a function of MLT and invariant latitude. meridional magnetometer array in conjunction with solar wind input from WIND at the L1 point (Valdivia et al., 1998). Here, too, the nonlinear spatial model is more accurate than the linear one. Spatiotemporal modeling will be extended to more than one magnetometer arrays. Yellow Knife, November
Novamkr 1-3.1078
E
Fig. 6. Coupling of nonlinear dynamical models and a ballistic propagation scheme as part of the ionospheric SWIFT model. First three panels: propagation of BEE-3 data to the subsolar magnetopause and comparison to IMP-8 near-Earth measurements. The last two panels are AL predictions (grey lines) from the BEE-3 input (IMP-8 gives similar results) and from the PC index.
6. COMBINATION OF NONLINEAR DYNAMICS TOWARDS GEOMAGNETIC AND IONOSPHERIC
5.1979
from both AL and Ati
Fig. 7. Effect of AU parameter on the SWIFT model output. At around 3 am on November 5, 1978, Yellowknife is west of the Harang discontinuity so the effect of the eastward_electrojet is significant. Using only the nonlinear-predicted AL in the geomagnetic part of SWIFT gives the wrong sign for the ground disturbance (dash-dotted line). When the predicted E is included, the prediction is improved significantly (dotted line).
WITH OTHER FORECASTING
MODELS
Nonlinear dynamical models are usually combined with other physical or empirical models. For example, the solar wind input for geomagnetic time series models is the measurement at L1 processed by a propagation program which simulates the transport to the subsolar magnetopause. Similarly when a nonlinear model is sufficiently accurate, its output can be used by other models as a proxy for the modeled magnetospheric disturbance. In this way more comprehensive models can be built. An example for such a combination of models for solar wind-magnetosphere -ionosphere coupling is the Space Weather Ionospheric Forecasting Technologies (SWIFT) model. The solar wind drives nonlinear dynamics and neural networks models, and in addition the solar wind input and predicted indices drive empirical ionospheric
206
D. Vassiliadis et nl.
models which give the electric-field patterns, currents, and Joule dissipation, and geomagnetic disturbance as a function of invariant latitude and MLT. An early test of SWIFT shows the propagation and index-prediction parts (Figure 6). The solar wind input is measured at L1 by ISEE- and propagated to the subsolar magnetopause where it is compared with IMP-g. There is fair agreement between the model and observed solar wind. Subsequently both the propagated solar wind and the PC index are used in turn as inputs to predict AL. Other tests have been carried out to measure the individual models’ sensitivity to perturbations to or uncertainties in the initial state. After the indices have been calculated, they are used in conjunction with the solar wind input to obtain the ionospheric electrodynamic variables. Since different current systems are parametrized through indices and other model variables, it is possible to optimize the model further by weighting the relative contribution of each current system to increase agreement with training data. Figure 7 shows the effect of adding the parameter of the eastward electrojet, the AU index, to the model optimization. A sequence of intensifications, starting at O4:OOUT, November 4 (when the measuring station of Yellowknife is in the dusk section), is predicted with the wrong sign. Inclusion of the AU parameter reduces the error by at least 50 and in some cases more than 150 nT. Index predictions are produced daily based on the solar wind input (originally from WIND and since May 1998 from ACE) and the PC index. A homepage has been set up for AL and AU, http:/Aepgst.gsfc.nasa.gov/ people/vassi/htmls/alprediction.html. For data-derived models it is important to maintain a database of intervals representative of as many activity levels and profiles as possible. Recently we have created an on-line database for AL/AU predictions during magnetic cloud passages in 1997 at h~~p://epgs~.gsfc..nasa.gov/people/vassilhhnls/97/even~s97.html. Remarkably the predictions are less accurate compared to the standard, non-cloud predictions. The explanation is that during the cloud events, the amplitude is larger and the variability is smaller than in the training intervals. Since the training dataset does not contain cloud conditions, it is difficult to predict AL/AU correctly. We plan to test this hypothesis by adding the cloud intervals to the database and repeating the runs. 7. SUMMARY Nonlinear dynamics techniques &e an essential set of data analysis methods which cover the middle ground between single-event and statistical analysis. We discussed a sequence of methods to quantitatively predict and physically The complexity of the model building can increase in a commensurate model a space physical environment. way with the variability of the data sets: From single-time-series prediction for scalar indices, we have gone to methods of physical interpretation of the models (e.g. Klimas et al., 1997), where we measure characteristic time scales (growth/decay/oscillation) and coupling strengths. The instantaneous state of the current system can be more accurately defined if measurements at several locations are available, from which one can construct spatiotemporal models and predictions. Finally if the accuracy of the prediction is sufficiently high, the model can be used in conjunction with other empirical or physical models. The evaluation of nonlinear dynamics models in terms of error statistics provides a good guide in two areas not covered in this paper: model verification and metric development. In developing space weather models, it is essential to verify the predictive capability by using out-of-sample (cross-validation) testing, and comparison between the candidate model and reference models or competing models to extract the prediction skill, i.e. the prediction capability relative to other models (issues discussed by Doggett, 1996). These are standard topics for nonlinear model development (e.g. the reference model is the linear model). Similarly the error analysis in nonlinear dynamics models (Vassiliadis and Klimas, 1996) can also be used as a guide to develop metrics necessary for progress evaluation in space weather forecasting (Wolf et al., Metrics for the National Space Weather Program, January 1999). 8. ACKNOWLEDGMENTS The following individuals and/or organizations were essential in providing the data we used in this study: Hightime-resolution D,t index: Y. Kamide, R. Nakamura, S. Sharma. GEF index: G. D. Reeves. Pi2 pulsation index: T. Iyemori, M. Nose’. ACE and WIND magnetic field and plasma data: Spacecraft PIs and NOAA/SEC. PC index: E. Friis-Christensen and S. Vennerstroem. Support from NASA/SR&T and NSF/GEM grants is gratefully acknowledged.
Nonlinear Dynamics of Space
REFERENCES Baker, D. N., A. J. Klimas, R. L. McPherron, and J. Buchner, The Evolution from Weak to Strong Geomagnetic Activity: an Interpretation in Terms of Deterministic Chaos, Geophys. Res. Len., 17,41 (1990). Bargatze, L. F.. D. N. Baker, R. L. McPherron, and E. W. Hones, Jr., Magnetospheric Impulse Response for Many Levels of Geomagnetic Activity, J. Geophys. Res., 90,6387 (1985). Burton, R. K., R. L. McPherron, and C. T. Russell, An Empirical Relationship Between Interplanetary Conditions and Dst, J. Geophys. Res., 80,4204 (1975). Campbell, W. H., Geomagnetic Storms, the Dst Ring-Current Myth, and Lognormal Distributions, f. Arm. Terr. Phys., 58, 1171 (1996). Clauer, C. R., The Technique of Linear Prediction Filters Applied to Studies of Solar Wind-Magnetosphere Coupling, in Solar Wind-Magnetosphere Coupling, edited by Y. Kamide and J. Slavin, pp. 39-57, Terra Scientific Publ. Co., Tokyo (1986). Doggett, K. A., (ed.), The Evaluation of Space Weather Forecasts. Proceedings of a Workshop in Boulder, Colorado, June 19-21, 1996, NOAA/SEC, Boulder, CO (1996). Horton W., and I. Doxas, A Low-Dimensional Dynamical Model for the Solar Wind-Driven Geotail-Ionosphere System J. Geophys. Res., 103,456l (1998). Klimas, A. J., D. Vassiliadis, D. N. Baker, and D. A. Roberts, The Organized Nonlinear Dynamics of the Magnetosphere, J. Geophys. Res., 101, 13089 (1996). Klimas, A. J., D. Vassiliadis, and D. N. Baker, Data-Derived Analogues of the Magnetospheric Dynamics, J. Geophys. Res., 102,26993 (1997). Klimas, A. J., D. Vassiliadis, and D. N. Baker, Dst Index Prediction Using Data-Derived Analogues of the Magnetospheric Dynamics, J. Geophys. Res., 103,20435 (1998). McPherron, R. L., The Role of Substorms in the Generation of Magnetic Storms, in Magnetic Storms, edited by B. T. Tsurutani, W. D. Gonzalez, Y. Kamide, and J. K. Arballo. pp. 131-147, American Geophysical Union, Washington, DC (1997). Potemra, T. A. (ed.), Magnetospheric Currents, American Geophysical Union, Washington, DC (1984). Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. V. Vetterling, Numerical Recipes: the Art of Scientific Computing (second ed.), Cambridge University Press, Cambridge, UK (1993). Prichard, D., and C. P. Price, Spurious Dimension Estimates from Time Series of Geomagnetic Indices, Geophys. Res. L&t., 19, 1623 (1992). Sauer, T., J. A. Yorke, and M. Casdagli, Embedology, J. Star. Phys., 65,579 (1991). Siscoe, G., Manifestations of global coherence, Eos Transactions AGU, Fall Meeting Suppl., 501 (1993). Troshichev, 0. A., V. G. Andrezen, S. Vennerstrom, and E. Friis-Christensen, Magnetic Activity in the Polar-Cap A New Index, Planet. Space Sci., 33, 1095 (1988). Valdivia, J. A., A. S. Sharma, and K. Papadopoulos, Prediction of Magnetic Storms by Nonlinear Models, Geophys. Res. L-e%.,23, 2899 (1996). Valdivia, J. A., D. Vassiliadis, A. J. Klimas, and A. S. Sharma, The Electrojet Currents: Understanding Their SpatioTemporal Multivariate Properties, in Proceedings of the Fourth International Conference on Substorms, edited by Y. Kamide and S. Kokubun, pp. 669-672, Terra Scientific Publ. Co., Tokyo (1998). Vassiliadis, D., and A. J. Klimas, Evaluation of Nonlinear Systems for Geomagnetic Index Forecasting and Modeling, in The Evaluation of Space Weather Forecasts. Proceedings of a Workshop in Boulder, Colorado, June 19-21, 1996, edited by K. A. Doggett, pp. 29-33, NOAA/SEC, Boulder, CO (1996). Vassiliadis, D., A. S. Sharma, T. E. Eastman, and K. Papadopoulos, Low-Dimensional Chaos in Magnetospheric Activity from AE Time Series, Geophys. Res. Zen., 17, 1841 (1990). Vassiliadis, D., A. J. Klimas, D. N. Baker, and D. A. Roberts, A Description of the Solar Wind-Magnetosphere Coupling Based on Nonlinear Filters, J. Geophys. Res., 100,3495 (1995). Vassiliadis, D., A. J. Klimas, D. N. Baker, and D. A. Roberts, The Nonlinearity of Models of the vBsouth - AL Coupling, J. Geophys. Res., 101, 19779 (1996a). Vassiliadis, D., V. Angelopoulos, D. N. Baker, and A. J. Klimas, The Relation Between Northern Polar Cap and Aurora1 Electrojet Geomagnetic Indices in the Wintertime, Geophys. Res. Lett., 23,278 1 (1996b). Wu, J.-G., and H. Lundstedt, Prediction of Geomagnetic Storms from Solar Wind Data Using Elman Recurrent Neural Networks, Geophys. Res. Zett., 23,3 19 (1996).
207