Next Article in Journal
Trends in Nighttime Fires in South/Southeast Asian Countries
Next Article in Special Issue
Correlation of Rate of TEC Index and Spread F over European Ionosondes
Previous Article in Journal
2D and 3D Properties of Stably Stratified Turbulence
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Time-Scale Analysis of Chaos and Predictability in vTEC

by
Massimo Materassi
1,
Yenca Migoya-Orué
2,
Sandro Maria Radicella
3,
Tommaso Alberti
4,5 and
Giuseppe Consolini
5,*
1
Institute for Complex Systems of the National Research Council CNR ISC, Via Madonna del Piano 10, 50019 Sesto Fiorentino, Italy
2
STI Unit, The Abdus Salam International Centre for Theoretical Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy
3
Visiting Research Scholar, Institute of Scientific Research, Boston College, Newton, MA 02459, USA
4
Istituto Nazionale di Geofisica e Vulcanologia (INGV), Via di Vigna Murata, 605, 00143 Rome, Italy
5
INAF-Istituto di Astrofisica e Planetologia Spaziali, Via del Fosso del Cavaliere 100, 00133 Rome, Italy
*
Author to whom correspondence should be addressed.
Atmosphere 2024, 15(1), 84; https://doi.org/10.3390/atmos15010084
Submission received: 19 November 2023 / Revised: 23 December 2023 / Accepted: 3 January 2024 / Published: 9 January 2024
(This article belongs to the Special Issue Ionospheric Irregularity)

Abstract

:
Theoretical modelling of the local ionospheric medium (LIM) is made difficult by the occurrence of irregular ionospheric behaviours at many space and time scales, making prior hypotheses uncertain. Investigating the LIM from scratch with the tools of dynamical system theory may be an option, using the vertical total electron content (vTEC) as an appropriate tracer of the system variability. An embedding procedure is applied to vTEC time series to obtain the finite dimension ( m N ) of the phase space of an LIM-equivalent dynamical system, as well as its correlation dimension ( D 2 ) and Kolmogorov entropy rate ( K 2 ). In this paper, the dynamical features ( m , D 2 , K 2 ) are studied for the vTEC on the top of three GNSS stations depending on the time scale ( τ ) at which the vTEC is observed. First, the vTEC undergoes empirical mode decomposition; then ( m , D 2 , K 2 ) are calculated as functions of τ . This captures the multi-scale structure of the Earth’s ionospheric dynamics, demonstrating a net distinction between the behaviour at τ 24 h and τ 24 h . In particular, sub-diurnal-scale modes are assimilated to much more chaotic systems than over-diurnal-scale modes.

1. Introduction

The occurrence of space and time irregularities in terms of ionospheric variability makes it extremely difficult to build up models of the Earth’s ionosphere from first principles.
The theoretical frameworks available to describe the ionospheric medium, e.g., multi-fluid representation [1,2], give rise to physical models that are very complicated and highly demanding in terms of integration time, as well as initial and boundary data. Moreover, ionospheric space and time irregularities cannot be transparently predicted by such physical models, which may be able to explain the appearance of irregularities only if suitable initial conditions are assumed (e.g., local fluctuations from initially “seeded” plasma instabilities (Chapter 11 in [2])), despite the inability to analytically follow their full development and to consistently represent their intermittent occurrence.
The Earth’s ionosphere is highly structured throughout space and time scales such that when observed at different space or time resolutions, its behaviour is completely different [1,3,4,5,6,7,8]. At some scales and under quiet solar and geomagnetic conditions, the ionosphere resembles a fluid system that can be described by smooth local variables undergoing classical partial derivative equations. However, when smaller scales are considered or disturbed geomagnetic conditions are studied, smoothness disappears, with emerging fractal features due to high gradients of the ionisation density ( N e ) that induce irregular phenomena such as radio scintillation (Chapters 13 and 18 in [2]). Due to this wide range of behaviours of the ionospheric medium, it is advisable to imagine the system at different resolutions and under different physical conditions as represented by quite different mathematical pictures, possibly ranging from kinetic theory to fluid systems, i.e., from large sets of particles to low-dimensional dynamical systems.
A paradigm-changing approach is to generate data-driven theoretical pictures of the local ionospheric medium (LIM), that is, treating the LIM in a certain region as a basically unknown physical system to be modelled from scratch and investigating its dynamical features using real data; the set of “dynamical features” ( D ) includes whatever one needs to identify the theoretical model (this does not mean renouncing comprehensive fundamental theory; it just means working it out from experimental observations as much as possible).
For a finite-dimensional dynamical system (FDDS), the set D includes the dimension ( m N ) of its phase space ( Γ ) and the geometrical–topological entities ( g ) describing the structures in the phase space, such as the equilibrium positions, their Lyapunov exponents, etc. In the presence of chaotic regions in Γ , g must include their Hausdorff dimensions. The use of FDDSs potentially paves the way for low-dimensional chaos to describe ionospheric irregularities, without the need to assume special initial “seeding”conditions [8]. If the theoretical picture is a kinetic theory (KT) or a field theory (FT), instead, its D set is much richer, as Γ may be a functional space of measures.
The FDDS approach has been used to study the LIM around a mid-latitude location, i.e., the GNSS station in Matera (Italy) [9]: this was achieved by applying the widely known embedding procedure to the time series of the vertical total electron content (vTEC), which was assumed to represent the state of the LIM on the top of the station location. By employing two yearly vTEC time series at 30 s resolution, corresponding to solar cycle (SC) minimum and maximum conditions (2008 and 2001, respectively), the dimension m = dim Γ was found to be the same, i.e., m 2008 = m 2001 = 3 . The same embedding dimension suggests that there is a common space of states of the system where the dynamics can be confined and described using (at least) three variables. This is the lowest (minimum) dimensional dynamical system that can capture dynamical features. Then, through the Grassberger–Procaccia method [10], the correlation dimension ( D 2 ) was calculated; this D 2 is the Hausforff dimension of the trajectory of the FDDS through its m-dimensional phase space ( Γ ). The surprisingly identical result was D 2 2008 = D 2 2011 = 2.78 . The fact that D 2 2008 and D 2 2011 have identical values suggests that at large scales, the LIM can be described in a similar way, independent of the phase of the solar cycle (Although disturbed geomagnetic conditions are expected to be more frequent during solar maxima than solar minima, when considering a yearly time series, their contribution is cancelled as far as the dynamical features D = m , D 2 , K 2 are concerned). Finding D 2 close to the embedding dimension (m) means that the system is highly chaotic (unpredictable) and suggests that the trajectories are space-filling, i.e., they resemble those of a noise (for which D 2 = m would hold), but they are not exactly stochastic. Finally, the Kolmogorov entropy rate ( K 2 ) was calculated; its inverse represents the predictability horizon, i.e., the time up to which a reliable prediction can be made. If at time t 0 , the state of the FDDS in Γ is known with infinite precision, at time t 0 + K 2 1 , the state of the system is unknown, and affected by a Shannon entropy of 1 bit; the extent of this ignorance grows over time. The two yearly vTEC series have basically the same value of ( K 2 1 8 min ).
To further inspect how the elements of D = m , D 2 K 2 depend on the geomagnetic conditions, one should construct them by extracting shorter time series of geomagnetically stormy and quiet periods.
The results reported in [9] were preliminary results, paving the way for a long-term program involving the dynamical characterization ( D ) of the LIM so that important elements of ionospheric physics can be better understood. In this specific work, the question is considered as to how D = m , D 2 , K 2 depends on the time scale ( τ ) at which the vTEC is observed. As the Earth’s ionosphere is a system structured at many time scales, one may expect that D may depend on τ . This aspect is investigated starting from a 30 s resolution vTEC time series, followed by the construction of its single time-scale version via empirical mode decomposition (EMD), as illustrated in detail in Section 2. Once the τ reduction of the vTEC series is selected, the values of the elements of D are calculated using an embedding procedure applied to this τ version of the vTEC; this produces τ -dependent parameters m τ , D 2 τ and K 2 τ (see also [11,12]). When different time scales are considered with respect to the LIM, completely different dynamics appear, in terms of phase-space topology. This means that the different scale-dependent forcing processes act on the system in a very different way, changing the active number of degrees of freedom. Indeed, our analysis (Section 3) highlights that the values of m τ , D 2 τ and K 2 τ show a reduction in predictability from large to small time scales at different geographical locations; this is the central result of the present work, as discussed with physical reasoning in Section 4. Conclusions are also drawn in this section from the point of view of what can be learned about space weather physics through analysis tools and the results obtained herein, as well as possible applications in terms of ionospheric modelling. A proposal for a future line of research is also presented.

2. Data and Methods

2.1. Data

We use three vTEC yearly time series from three GNSS stations, i.e., Matera (mid-latitude, northern hemisphere; coordinates: 40.649 N , 16.704 E , Italy), Ohrid (mid-latitude, northern hemisphere; coordinates: 41.127 N , 20.794 E , North Macedonia) and Harbour (low geographic latitude, southern hemisphere; coordinates: 25.88 N , 27.70 E , South Africa) collected during the solar maximum year of 2001. The sampling time of all the time series is 30 s , as in [9].
The choice of these three GNSS stations depends on data availability, on the one hand, and on the possibility of assorting the position of the analysed LIM, expecting similar behaviour in the Matera and Ohrid vTEC time series, with some difference in the dynamics of the vTEC on the top of Harbour. The choice of the Harbour GNSS station is due to the fact that it has practically the same modified dip angle (MODIP, see [13]) and longitude as the other two stations in the northern hemisphere; any differences between Harbour and the other two stations are due to being located in the magnetically opposite hemisphere rather than the geographic latitude. Moreover, we focus only on the solar maximum year time series, since it shows similarity with those of the solar minimum [9], as long as quiet and stormy periods are not separated, despite allowing us to consider a greater fraction of time during disturbed geomagnetic conditions, which partially affects the time-scale behaviour of near-Earth space plasma (e.g., [14]).
The differences in the behaviour of the LIM on the top of the three locations is evident at a glance in Figure 1, where the raw time series are directly reported. The dynamics of Matera (red) and Ohrid (blue) are very similar, while the vTEC time series on the top of Harbour (green) exhibits considerably different behavior, even if periods of maxima and minima within the year coincide.

2.2. Empirical Mode Decomposition (EMD)

The behavior of the vTEC series y t at different time scales is inspected via empirical mode decomposition (EMD, see for instance, [15]), through which y t is decomposed into a sum of a certain number of τ components ( υ τ t ):
y t = τ y υ τ y , t .
Note that in (1), the set of time scales ( τ ) depends on the collection of values ( y t ) of the decomposed time series, i.e., it depends on y in a functional way, as mimicked by y ; this is because EMD is a data-driven method in the sense that both the form of the functions ( υ τ ) and the characteristic scale ( τ ) of each depend functionally on y t (this is also what y , t means in Equation (1)). This is not, for instance, a feature of Fourier or wavelet decomposition, in which the period or the scale of the components is determined based on the size of the time-series domain.
Larger scales ( υ τ s) correspond to slower modes:
τ 2 > τ 1 d d t υ τ 1 > d d t υ τ 2 ,
as the larger the considered scale, the slower the mode, on average (as denoted by in (2)).
In order to get a pictorial idea of the components ( υ τ t ) into which the EMD splits the time series at hand, Figure 2 reports the different τ components for the time series on the top of the Matera GNSS station.
After τ decomposition (1), a “ τ -cumulative partial time series” ( Υ τ t ) is defined as the sum of the components ( υ τ y , t ) with a characteristic time scale smaller than or equal to τ , i.e., neglecting the modes pertaining to scales larger than τ :
Υ τ t = τ τ υ τ y , t
Of course, this definition makes Υ τ depend functionally on y t precisely as each υ τ , but this “square bracket dependence” is not highlighted in (3) in order to keep notations reasonably simple.
In the Υ τ t the time series, the original y t is considered, and all time modes “slower than” υ τ in the sense of Equation (2) are discarded. Therefore, the Υ τ t time series is the vTEC without the slow modes, and the τ scale pertains to the slowest admitted mode. If τ max is the largest time scale admitted in the modes ( υ τ ) in (1), clearly, one may write
lim τ τ max Υ τ t = y t .
As the functions ( Υ τ t ) are calculated, each of them may undergo an embedding procedure applied to the whole vTEC time series reported in [9], as well as the subsequent calculation of D 2 and K 2 . The main difference here is that this analysis of “chaos degree” and predictability of the LIM is performed at each different time scale ( τ ) for each Υ τ t (see also [11,12]). The embedding procedure à la Grassberger–Procaccia [10], applied to the partial series ( Υ τ t ) determines a trajectory ( Υ τ t ) moving through an m τ dimensional real space:
Υ τ t R m τ ,
which represents the evolution of a suitable τ FDDS and produces the τ -resolved vTEC ( Υ τ t ). This trajectory ( Υ τ t ) lives in a phase space ( Γ τ of dimension) m τ N . As underlined by the bracket in m τ , the dimension of the τ FDDS also depends on the time scale ( τ ) because observing the LIM at different time resolutions, the system may look like different dynamical systems, with many or few dynamical variables, so that each Γ τ has a different dimension.
The coarse graining of Γ τ defines the Hausdorff dimension of the curve, which is well approximated by the correlation dimension of Υ τ t (see [9,10]), namely D 2 Υ τ ( D 2 τ for simplicity). This dimension ( D 2 τ R + ) is, in practice, the fractal dimension of the attractor containing Υ τ t , where the τ FDDS dynamics live. As discussed in [9], D 2 τ is bounded from above by m τ
D 2 τ 1 , m τ :
the closer it is to 1, the more regular the dynamics of the τ FDDS; for 1 < D 2 τ < m τ , the τ FDDS becomes more chaotic as D 2 τ approaches m τ ; for D 2 τ m τ , the τ FDDS is simply stochastic.
Following the prescriptions proposed in [10], the calculation of the Kolmogorov entropy rate ( K 2 ) is applied to the single- τ time series ( Υ τ t ), i.e., to its m τ -embedded version ( Υ τ t ), yielding a value of K 2 τ , representing the rate at which the coarse-grained version of Υ τ t R m τ is located less precisely than on bit, i.e., the time ( K 2 1 τ ) after which the position ( Υ τ t + K 2 1 τ ) is known less precisely than Υ τ t of one bit. This time ( K 2 1 τ ) is understood as the predictability time horizon of the τ FDDS. For K 2 1 τ 0 , the τ FDDS is fully unpredictable, while for K 2 1 τ + , it is completely deterministic. Our τ FDDS will shows a finite value of K 2 1 τ , indicating the typical finite predictability behaviour of a chaotic system.
Once the τ -FDDS characteristics ( D τ ) are calculated as functions of the time scale ( τ ), it is possible to draw physical conclusions considering those functions ( m τ , D 2 τ , K 2 τ ). In particular, the number of dynamical degrees of freedom of the LIM as a function of the time scale, as well as degree of chaos of their dynamics, may highlight critical scales in these complex dynamics as changes in the behaviour of m τ , D 2 τ , K 2 τ take place in the correspondence with changes in the value of τ . We report on this further in our analysis presented in Section 3. In the case reported in Figure 3, the pseudospectrum
Π υ τ = τ · σ υ τ
is plotted, where σ υ τ is the standard deviation for the set of values of the time series of the single- τ components in (1) ( υ τ y , t ). This pseudospectrum is useful because its dependence on τ corresponds approximately to the value of the Fourier spectrum of the y t from which the υ τ y , t are calculated at a frequency of f τ = 1 τ . Directly using the pseudospectrum instead of the Fourier spectrum is motivated by its availability once EMD is performed and by the ease of comparison with the other results obtained via the same υ τ t functions. In Figure 3 and in all the subsequent analysis, components with approximately τ 30 min are considered, even if TEC variations with time scales of 5–20 min are linked to atmospheric waves; these components are of particular interest with respect to ionospheric dynamics. The reason to disregard those important short time scales here is related to the calculation of the correlation integrals. The latter can be affected by high-frequency contents of the time series. For this reason, modes ( υ τ ) whose time scales are below one order of magnitude the time resolution, i.e., τ 300 s , are excluded. Components ( υ τ 30 min ) can definitely be interesting in terms of ionospheric dynamics, and surely, their behaviour deserves further study, as well as additional investigation of parameters related to irregularities and radio scintillation. Moreover, the vTEC time series processed here are the result of a “calibration procedure” [16] that employs dual-frequency carrier-phase and code-delay GNSS observations related to the total electron content (TEC) along the satellite–receiver line of sight. As a consequence, the obtained vTEC is a linear combination of many TEC measurements, each including some random instrumental and processing errors. The exclusion of those high-frequency contributions ( τ 300 s ) eliminates this additive, purely stochastic process resulting from calibration.
Before moving to results, it is of use to recall the formulae due proposed by Grassberger and Procacci and used in [9] for the whole yearly time series and adapted here to treat the motions ( Υ τ t ) in Γ τ . The central quantity of all the calculations is the correlation integral ( C m τ Υ τ , r ) pertaining to the relationship between two points ( Υ τ t and Υ τ t in Γ τ ) along the trajectory of the τ FDDS, separated by an m τ -dimensional distance equal to r:
C m τ Υ τ , r = def lim N + 1 N 2 i , j = 1 N Θ r Υ τ t i Υ τ t j
In (6), the number of times considered is N; hence, the points of the trajectory ( Υ τ t ) are Υ τ t i = 1 , , N , corresponding to “an enormous number” ( N + ) of points. In practice, (6) represents the number of trajectory points closer than r to each other, which is normalized with respect to the number of point couples ( O N 2 ).
Once C m τ Υ τ , r is calculated as in (6) for the m τ -dimensional τ FDDS, D 2 τ and K 2 τ are calculated as proposed by Grassberger and Procaccia [10]:
D 2 τ = lim r 0 log C m τ Υ τ , r log r , K 2 τ = lim r 0 1 τ y log C m τ 1 Υ τ , r C m τ Υ τ , r .
The expression of D 2 τ in (7) qualifies it as the Hausdorff dimension of the r-coarse-grained trajectory ( Υ τ ) within Γ τ , while the expression of K 2 τ indicates that for each τ FDDS, one has to calculate C n Υ τ , r for various tentative embedding dimensions (n), selecting m τ as the value of n for which the growth of C n Υ τ , r reaches a plateau (see [9,10] for more details).
As a concluding remark, it is worth stressing again that all the quantities calculated by embedding the τ -limited time series ( Υ τ t ) are constrained by the limit relationship (4). As the τ max -limited time series includes all υ τ ,
D 2 τ max D 2 y , K 2 τ max K 2 y ,
where D 2 y and K 2 y are the correlation dimension and Kolmogorov entropy rate obtained by embedding the whole original time series ( y t ) into a suitable m-dimensional FDDS, respectively. Therefore, the topological characteristics of the τ FDDS, i.e., of the τ -limited Υ τ that includes all the τ -fixed components (1), are simply those of the undecomposed y t (m, D 2 and K 2 of Matera are the same as those found in [9]).

3. Results

Figure 4 reports the results of the embedding procedure described in [9], applied to the time series analysed here and shown in Figure 1, where the trajectory ( y t ) associated with the whole time series ( y t ) is reported in the Γ space for Matera, Harbour and Ohri. These trajectories are simply the time vector series y t , y t + Δ , y t + 2 Δ , y t + m 1 Δ , i.e., their components are suitably delayed copies of the original observed time series y t . As explained in [9] and references therein, if the FDDS equivalent to that yielding y t is diffeomorphic to the vector of the delayed copies of the time series, it is possible to visualise the FDDS trajectories throughout the phase space. The lag Δ reduces the Δ -delayed self-mutual information of the y t to less than e 1 bit.
The three Γ s in Figure 4 all have an embedding dimension of 3. The result represented in Figure 4 is consistent with the appearance of the y t time series in Figure 1, with an apparently more elongated and flatter region of attraction for the low-geographic-latitude LIM on Harbour with respect to those on Matera and Ohrid.
The first notable result of the decomposition (1) is the behaviour with τ of the pseudospectrum ( Π υ τ ) calculated as in (5) and reported in Figure 3, where the pseudospectrum is shown for the three time series ( y Matera t , y Harbour t and y Ohrid t ). With respect to the general behaviour, the three pseudospectra are rather similar, while in the details, they differ quite a bit, with a stricter resemblance between Matera and Ohrid than relative to Harbour. Besides the differences among the three curves, one observes a considerable change in regime at τ = 1000 min , which is close to the τ day = 1441 min corresponding to the length of one terrestrial day. The slopes of Π υ τ for τ < τ day and for τ > τ day are different; in particular, the short time components ( τ < τ day ) exhibit steeper growth with τ than the slower ones. Moreover, a local maximum of Π υ τ appears at τ τ day . The pseudospectrum of the Harbour vTEC is lower than that of Matera and Ohrid vTECs, and for τ > τ day , the growth associated with τ seems “less monotonic”, with a local maximum at about τ 9000 min . Whether this maximum at a scale of about 6.25 days has a physical meaning or not and what that meaning might be deserve a much deeper study consulting more statistics of the GNSS stations. Figure 3 shows the two slope regimes separated by the time scale of the Earth’s rotation.
Once the τ -limited series ( Υ τ t ) are composed, (6) is applied to define the correlation integrals ( C m τ Υ τ , r ), as reported in Figure 5, where the curves ( C r ) correspond to C m τ Υ τ , r , and different colours correspond to different values of the time-scale τ . The value of C r varies as a function of τ for a fixed r, as do the characteristics of the τ FDDS, as they depend on C m τ Υ τ , r C r according to Formula (7). Again, the similarity between the LIM on the top of the two mid-latitude stations (Matera and Ohrid) and the difference with respect to the low-latitude station (Harbour) are evidenced in Figure 5.
Once the quantities C m τ Υ τ , r are calculated, it is possible to explore the degree of chaoticity of the different τ FDDSs by calculating the dimension ( D 2 τ ) and their predictability by calculating the entropy rate ( K 2 τ ), applying the Grassberger–Procaccia formulae (7).
The dependence ( D 2 τ ) is plotted in Figure 6 for the three GNSS stations. The 1-day scale break appears again, with (slightly) different behaviour for the three different stations. The similarities between the two mid-latitude stations and the difference between those two and the low-latitude station appear to be less important here than in the pseudospectra shown in Figure 3. It can be concluded that the representation of the LIM by τ FDDSs is limited to scales smaller than τ day , reaching values higher than 3; therefore, excluding scales larger than τ day , the effective dynamical system representing the local ionosphere needs more than three independent dynamical variables to be described. In addition to τ day , when including much slower υ τ s in Υ τ , D 2 τ rapidly decreases, reaching almost 2.8 for Matera, nearly 2.6 for Ohrid and as low as 2.3 for Harbour. This result is in agreement with what is illustrated in Figure 4, where the three whole time series ( y Matera t , y Ohrid t and y Harbour t ) appear to be embedded in a three-dimensional FDDS. Moreover, as apparent in the same Figure 4, the degree of chaos at Harbour is (slightly) lower than at the mid-latitude LIMs.
The dynamical behaviour of the vTEC on the top of the three considered GNSS stations, appearing to be richer and more chaotic at shorter than at longer time scales, with a drastic separation into τ < τ day and τ > τ day , is also reflected by the τ -dependent predictability time horizon ( K 2 1 τ ), as represented in Figure 7. For τ -limited dynamics with τ < τ day , the highest value of K 2 1 τ is about 3.8 min , corresponding to Υ 250 min t on the top of Harbour, and generally remaining below 3 min. This predictability time horizon grows abruptly by as much as 5.3 min for τ > τ day , increasing to almost 6 min for Υ 6000 min t on the top of Matera.
The dependence of D 2 and K 2 on τ reported in Figure 6 and Figure 7 not only shows a qualitative transition between τ τ day and τ τ day but also some kind of “saturation” for τ + . This may sound surprising, mainly for K 2 τ , as increasing τ values correspond to components υ τ with longer “periods”, namely increased predictability. However, here, we are calculating D 2 and K 2 for the FDDS equivalent to Υ τ t in Equation (3), which represents the cumulative contributions of υ τ τ . The series Υ τ t still encodes the chaotic properties of all the time scales ( τ τ ), and increasing τ does not exclude those highly irregular components, preventing the indefinite growth of the predictability time horizon beyond saturation.

4. Discussion and Conclusions

Herein, we present here a step forward in describing the local ionospheric medium (LIM) in a brand-new and deeply dynamic theoretical manner.
Data-conditioned representations of the LIM are traditionally empirical models constructed through the collection of huge datasets referring to suitably assorted ionospheric locations and conditions [17,18]. Despite their success in describing the average climatological behaviour of systems, they are neither able to take into account suddenly developing irregularities nor focus on theoretical explanations of LIM dynamics. In this study, we interrogated ionospheric data in order to evaluate quantities pointing to describe the LIM as a dynamical system from scratch. This attitude, on the one hand, allows results to be obtained based on reality rather than (even physically reasonable) assumptions and, on the other hand, the LIM to invoke representations with the goal of tracing detailed physical mechanisms.
The physical data representing LIM evolution are time series of the vertical total electron content (vTEC) on the tops of three considered stations. Our aim was to describe the LIM dynamics as those of a finite dimensional dynamical system (FDDS). These time series are embedded into a suitable phase space according to the Grassberger–Procaccia procedures [10], and the Hausdorff dimensions of their strange attractors, as well as the Kolmogorov entropy rate, are calculated.
In [9], raw yearly time series of a vTEC ( y t ) were analysed so to ensure an equivalent FDDS of trajectory y t R m , m N ; here, we decompose time series y t into empirical modes ( υ τ t ), each with a characteristic time scale ( τ ), then use them to define τ -limited time series ( Y τ t ) as in (3) so that Y τ t includes all the modes ( υ τ t ) from y t faster than υ τ t . Each of the τ -limited vTEC series ( Y τ t ) is ultimately embedded as Y τ t R m τ , and the dynamical proxies ( m τ , D 2 τ , K 2 τ ) of the scale-dependent embedding dimension, correlation dimension and Kolmogorov entropy rate are calculated. Excluding from Y τ all the slower components ( υ τ > τ ) reveals the faster motion details of the LIM dynamics. In fact, slower modes tend to have larger amplitudes, masking the faster, weak-amplitudevariability.
The interest in assigning a distinct equivalent FDDS to each Y τ is twofold. On the one hand, the LIM may look very different at different time scales; identifying the dependence of the dynamical quantities ( m τ , D 2 τ and K 2 τ ) on time scales rephrases this statement in a precise and quantitative way. On the other hand, values of τ discriminating different behaviours of m τ , D 2 τ or K 2 τ should represent characteristic time scales of the local ionospheric physics that are useful in singling out important processes influencing the LIM.
Our results suggest that (1) there is an expected similarity between the vTEC dynamics above Matera and Ohrid, while the vTEC above Harbour behaves a slightly differently, indicating hemispheric variation and (2) a time scale of about 1 day represents a threshold at which the characteristics of the τ FDDS change.
The take-home message of these results is that multi-scale analysis calculating the quantities of m τ , D 2 τ and K 2 τ unveils very important new details about LIM dynamics as represented by vTEC time series. According to the findings presented here, the LIM is a complex dynamical system that has finite dimensional representations depending on the time scale. The smaller the time scale ( τ ), the higher the number of dynamical variables necessary to describe the LIM at that scale. Moreover, all τ FDDSs representing the LIM evolve along chaotic trajectories ( Y τ t ), filling D 2 τ -dimensional attractors with D 2 τ N and increasing with smaller scales; therefore, the smaller the τ , the more chaotic the LIM-equivalent τ FDDS. Lastly, the predictability time horizon ( K 2 1 τ ) decreases with smaller scales; therefore, the smaller the τ , the more unpredictable the LIM-equivalent τ FDDS. The picture presented herein is that of a chaotic and unpredictable system at all time scales, although it is more chaotic and less predictable in its fast modes than in its slow modes. When examining the LIM at a more detailed time scale, one is able to distinguish the action of independent dynamical variables that are undetectable at large time scales, i.e.,
d m τ d τ < 0 ,
and one finds that such “fast modes” are invisible at large τ s, resulting in a higher level of chaos, i.e.,
d D 2 τ d τ < 0 ,
and unpredictability, i.e.,
d K 2 1 τ d τ > 0 .
Equations (9)–(11) refer to the “trends” of m τ , D 2 τ and K 2 1 τ , while oscillations and local peaks and valleys with τ may appear. These tendencies abruptly increase at about 1440 min . In accordance with Earth’s rotation period, m τ and D 2 τ decrease abruptly with increasing τ , while K 2 1 τ exhibits a sudden increase. This has a precise physical meaning: sunlight is the natural time scale that separates fast modes from slow modes in the vTEC dynamics, which sounds very reasonable due to the role played by insolation in ionospheric physics.
All the conclusions just reported apply slightly differently to the two northern hemisphere, mid-latitude locations than to the southern hemisphere location examined herein, with the caveat that such a difference is an interhemispheric difference. However, geomagnetic latitude should be expected to influence the detailed behaviour of m τ , D 2 τ and K 2 1 τ , probably more than other local conditions. This calls for a systematic extension of the present study in order to confirm our findings.
Based on the results reported in [9] and those presented here, our attempt to construct an FDDS equivalent to the local ionospheric medium portrays the system as chaotic and unpredictable to a certain extent, appearing richer and less predictable in its fast modes than in its slow modes, corresponding to the period of the Earth’s daily rotation, i.e., the separation time scale. The details of these findings depend on the location under consideration. The emergence of this complex behaviour is witnessed here for yearly time series ( y t ) of vTEC above the three investigated stations. Therefore, we suggest (1) an extension of this analysis to many more locations in order to more consistently how the m τ , D 2 τ and K 2 τ vary with geomagnetic coordinates; (2) time-limited studies concentrating on geomagnetically homogeneous periods in order to determine how the multi-scale complexity of vTEC series depends on present geomagnetic activity, as well as techniques compatible with shorter time series and the expression of τ FDDS quantities as functions of the present (e.g., m = m τ , t ; see [19]), which can help to monitor LIM complexity above a given location throughout the various phases of a geomagnetic storm [20].
The very difficult challenge of this line of research is to pass from knowledge of m , D 2 , K 2 for an FDDS to a possible set of m-coupled ordinary differential equations:
d Y d t = f Y ,
where Y = Y 1 , , Y m represent independent dynamical variables moving along a D 2 -dimensional attractor within R m , with a K 2 Kolmogorov entropy rate [21]. Mapping m , D 2 , K 2 into the suitable system (12) is all but a trivial challenge in general, and in the case of the LIM, one has a different FDDSs, i.e., a different (12) at each scale ( τ ) and probably also differing under the various geomagnetic conditions. The theoretical sense of arriving at a system like (12) for the LIM is to produce, for the ionosphere, what the Lorenz system equations [22] represent for the physics of the neutral atmosphere.
In principle, within a geomagnetically homogeneous period, after collecting all the dynamics ( Y τ t s) and finding what (12) looks like for Y τ ,
d Y τ d t = f τ Y τ
(that is, the Lorenz equivalent of the LIM system at scale (τ)), keeping track of the function Y τ Y τ , (13) is an analytical/numerical predictor of the τ -limited vTEC within the time interval under consideration. Despite the loftiness of the goal, this program, when fully pursued, may represent a striking breakthrough in ionospheric modelling.

Author Contributions

Conceptualization, M.M. and S.M.R.; methodology, T.A. and G.C.; software, T.A.; validation, M.M., Y.M.-O. and S.M.R.; formal analysis, M.M. and T.A; data curation, Y.M.-O.; writing—original draft preparation, M.M.; writing—review and editing, Y.M.-O., S.M.R., T.A. and G.C. All authors have read and agreed to the published version of the manuscript.

Funding

T.A. acknowledges funding from the “Bando per il finanziamento di progetti di Ricerca Fondamentale 2022” of the Italian National Institute for Astrophysics (INAF)—Mini Grant: “The predictable chaos of Space Weather events”.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The vTEC data series were obtained from the ICTP Calibrated GNSS TEC data service and are available at https://arplsrv.ictp.it, accessed on 10 January 2022.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Kelley, M.C. The Earth’s Ionosphere—Plasma Physics and Electrodynamics; Elsevier Inc.: Amsterdam, The Netherlands, 1989; ISBN 978-0-12-404013-7. [Google Scholar]
  2. Materassi, M.; Forte, B.; Coster, A.; Skone, S. (Eds.) The Dynamical Ionosphere—A Systems Approach to Ionospheric Irregularity; Elsevier Inc.: Amsterdam, The Netherlands, 2020; ISBN 978-0-12-814782-5. [Google Scholar]
  3. Goodman, J.M. The Ionosphere. In Space Weather & Telecommunications; The International Series in Engineering and Computer Science; Springer: Boston, MA, USA, 2005; Volume 782. [Google Scholar] [CrossRef]
  4. Hunsucker, R.D. Radio Techniques for Probing the Terrestrial Ionosphere; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013; Volume 22. [Google Scholar]
  5. Lilensten, J.; Blelly, P.L. The TEC and F2 parameters as tracers of the ionosphere and thermosphere. J. Atmos. Sol.-Terr. Phys. 2002, 64, 775–793. [Google Scholar] [CrossRef]
  6. Bolzan, M.J.A.; Tardelli, A.; Pillat, V.G.; Fagundes, P.R.; Rosa, R.R. Multifractal analysis of vertical total electron content (VTEC) at equatorial region and low latitude, during low solar activity. Ann. Geophys. 2013, 31, 127–133. [Google Scholar] [CrossRef]
  7. Chandrasekhar, E.; Prabhudesai, S.S.; Seemala, G.K.; Shenvi, N. Multifractal detrended fluctuation analysis of ionospheric total electron content data during solar minimum and maximum. J. Atmos. Sol.-Terr. Phys. 2016, 149, 31–39. [Google Scholar] [CrossRef]
  8. Kumar, K.S.; Kumar, C.V.A.; George, B.; Renuka, G.; Venugopal, C. Analysis of the fluctuations of the total electron content (TEC) measured at Goose Bay using tools of nonlinear methods. J. Geophys. Res. 2004, 109, A02308. [Google Scholar] [CrossRef]
  9. Materassi, M.; Alberti, T.; Migoya-Orué, Y.; Radicella, S.M.; Consolini, G. Chaos and Predictability in Ionospheric Time Series. Entropy 2023, 25, 368. [Google Scholar] [CrossRef] [PubMed]
  10. Grassberger, P.; Procaccia, I. Characterization of strange attractors. Phys. Rev. 1983, 50, 346. [Google Scholar] [CrossRef]
  11. Alberti, T.; Consolini, G.; Ditlevsen, P.D.; Donner, R.V.; Quattrociocchi, V. Multiscale measures of phase-space trajectories. Chaos Interdiscip. J. Nonlinear Sci. 2020, 30, 123116. [Google Scholar] [CrossRef] [PubMed]
  12. Alberti, T.; Donner, R.V.; Vannitsem, S. Multiscale fractal dimension analysis of a reduced order model of coupled ocean-atmosphere dynamics. Earth Syst. Dyn. 2021, 12, 837–855. [Google Scholar] [CrossRef]
  13. Rawer, K. Meteorological and Astronomical Influences on Radio Wave Propagation; Landmark, B., Ed.; Academic Press: New York, NY, USA, 1963; pp. 221–250. [Google Scholar]
  14. Consolini, G.; Alberti, T.; De Michelis, P. On the Forecast Horizon of Magnetospheric Dynamics: A Scale-to-Scale Approach. J. Geophys. Res. 2018, 123, 9065–9077. [Google Scholar] [CrossRef]
  15. Huang, N.E.; Shen, Z.; Long, S.R.; Wu, M.C.; Shih, H.H.; Zheng, Q.; Yen, N.-C.; Tung, C.C.; Liu, H.H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. Ser. A 1998, 454, 903. [Google Scholar] [CrossRef]
  16. Ciraolo, L.; Azpilicueta, F.; Brunini, C.; Meza, A.; Radicella, S.M. Calibration errors on experimental slant total electron content (TEC) determined with GPS. J. Geod. 2007, 81, 111–120. [Google Scholar] [CrossRef]
  17. Bilitza, D. IRI the International Standard for the Ionosphere. Adv. Radio Sci. 2018, 16, 1–11. [Google Scholar] [CrossRef]
  18. Radicella, S.M. The NeQuick model, genesis, uses and evolution. Ann. Geophys. 2009, 52, 417–422. [Google Scholar]
  19. Alberti, T.; Faranda, D.; Lucarini, V.; Donner, R.V.; Dubrulle, B.; Daviaud, F. Scale dependence of fractal dimension in deterministic and stochastic Lorenz-63 systems. Chaos Interdiscip. J. Nonlinear Sci. 2023, 33, 023144. [Google Scholar] [CrossRef] [PubMed]
  20. Ogunsua, B.O.; Laoye, J.A.; Fuwape, I.A.; Rabiu, A.B. The comparative study of chaoticity and dynamical complexity of the low-latitude ionosphere, over Nigeria, during quiet and disturbed days. Nonlinear Process. Geophys. 2014, 21, 127–142. [Google Scholar] [CrossRef]
  21. Brunton, S.L.; Proctor, J.L.; Kutz, J.N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. USA 2016, 113, 3932–3937. [Google Scholar] [CrossRef] [PubMed]
  22. Lorenz, E.N. Deterministic non periodic flow. J. Atmos. Sci. 1963, 20, 130–141. [Google Scholar] [CrossRef]
Figure 1. The three yearly time series of vTEC during solar maximum year 2001: from the tops of the Matera, Harbour and Ohrid GNSS stations.
Figure 1. The three yearly time series of vTEC during solar maximum year 2001: from the tops of the Matera, Harbour and Ohrid GNSS stations.
Atmosphere 15 00084 g001
Figure 2. The different components ( υ τ 1 , , 14 t ) into which the EMD splits the time series ( y Matera t ) plotted in Figure 1. The convention here is to count the τ s from the shortest value to the longest one so that τ i < τ i + 1 .
Figure 2. The different components ( υ τ 1 , , 14 t ) into which the EMD splits the time series ( y Matera t ) plotted in Figure 1. The convention here is to count the τ s from the shortest value to the longest one so that τ i < τ i + 1 .
Atmosphere 15 00084 g002
Figure 3. Plot of the pseudospectrum as defined in (5) for the three time series reported in Figure 1. The remarkable feature of these plots is the break at τ = 1440 min 1 day , corresponding to the rotation period of the planet, resulting in the insolation tempo (see text for details).
Figure 3. Plot of the pseudospectrum as defined in (5) for the three time series reported in Figure 1. The remarkable feature of these plots is the break at τ = 1440 min 1 day , corresponding to the rotation period of the planet, resulting in the insolation tempo (see text for details).
Atmosphere 15 00084 g003
Figure 4. The three-dimensional trajectories resulting from the embedding procedure as described in [9] applied to the whole yearly time series plotted in Figure 1. Note the strong similarity between the trajectories ( y t ) for Matera and Ohrid, while the LIM on the top of Harbour appears to correspond to a differently shaped attractor. See the text for details.
Figure 4. The three-dimensional trajectories resulting from the embedding procedure as described in [9] applied to the whole yearly time series plotted in Figure 1. Note the strong similarity between the trajectories ( y t ) for Matera and Ohrid, while the LIM on the top of Harbour appears to correspond to a differently shaped attractor. See the text for details.
Atmosphere 15 00084 g004
Figure 5. The correlation integrals ( C m τ Υ τ , r C r ) depending on the time scale ( τ ) for the three considered GNSS stations.
Figure 5. The correlation integrals ( C m τ Υ τ , r C r ) depending on the time scale ( τ ) for the three considered GNSS stations.
Atmosphere 15 00084 g005
Figure 6. The dependence ( D 2 τ ) for the τ FDDS representing the LIM on the top of the three GNSS stations (Matera, Ohrid and Harbour). The horizontal dashed line represents the value of D 2 τ for a large τ value for the Matera LIM, taken as a reference; and the vertical dashed line marks the τ day scale. See the text for details.
Figure 6. The dependence ( D 2 τ ) for the τ FDDS representing the LIM on the top of the three GNSS stations (Matera, Ohrid and Harbour). The horizontal dashed line represents the value of D 2 τ for a large τ value for the Matera LIM, taken as a reference; and the vertical dashed line marks the τ day scale. See the text for details.
Atmosphere 15 00084 g006
Figure 7. The dependence ( K 2 τ ) for the a τ FDDS representing the LIM on the top of the three GNSS stations (Matera, Ohrid and Harbour). The horizontal dashed line represents the value of K 2 τ for a large τ value for the Matera LIM, taken as a reference; and the vertical dashed line marks the τ day scale. See the text for details.
Figure 7. The dependence ( K 2 τ ) for the a τ FDDS representing the LIM on the top of the three GNSS stations (Matera, Ohrid and Harbour). The horizontal dashed line represents the value of K 2 τ for a large τ value for the Matera LIM, taken as a reference; and the vertical dashed line marks the τ day scale. See the text for details.
Atmosphere 15 00084 g007
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Materassi, M.; Migoya-Orué, Y.; Radicella, S.M.; Alberti, T.; Consolini, G. Multi-Time-Scale Analysis of Chaos and Predictability in vTEC. Atmosphere 2024, 15, 84. https://doi.org/10.3390/atmos15010084

AMA Style

Materassi M, Migoya-Orué Y, Radicella SM, Alberti T, Consolini G. Multi-Time-Scale Analysis of Chaos and Predictability in vTEC. Atmosphere. 2024; 15(1):84. https://doi.org/10.3390/atmos15010084

Chicago/Turabian Style

Materassi, Massimo, Yenca Migoya-Orué, Sandro Maria Radicella, Tommaso Alberti, and Giuseppe Consolini. 2024. "Multi-Time-Scale Analysis of Chaos and Predictability in vTEC" Atmosphere 15, no. 1: 84. https://doi.org/10.3390/atmos15010084

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop