Statistical Mechanics and Information-Theoretic Perspectives on Complexity in the Earth System

Abstract

This review provides a summary of methods originated in (non-equilibrium) statistical mechanics and information theory, which have recently found successful applications to quantitatively studying complexity in various components of the complex system Earth. Specifically, we discuss two classes of methods: (i) entropies of different kinds (e.g., on the one hand classical Shannon and Rényi entropies, as well as non-extensive Tsallis entropy based on symbolic dynamics techniques and, on the other hand, approximate entropy, sample entropy and fuzzy entropy); and (ii) measures of statistical interdependence and causality (e.g., mutual information and generalizations thereof, transfer entropy, momentary information transfer). We review a number of applications and case studies utilizing the above-mentioned methodological approaches for studying contemporary problems in some exemplary fields of the Earth sciences, highlighting the potentials of different techniques.

1. Introduction

Complexity is, nowadays, a frequently used, yet (at least quantitatively speaking) poorly defined, concept [1]. It tries to embrace a great variety of scientific and technological approaches of all types of natural, artificial and social systems [1]. The concept of complexity in the Earth system pertains to processes that are not deterministic, but can be modeled with a statistical element, including fractal phenomena and deterministic chaos [2].

As in other areas of science, the last few decades of geoscientific research have been shaped by a rapid increase in the amount of existing data on relevant phenomena and processes obtained from extensive measurement campaigns (including both ground-based observations and remote sensing applications), as well as detailed simulation models making use of the meanwhile available powerful computation infrastructures. The resulting big datasets call for novel techniques, not just for data acquisition and storage, but, most prominently, also for their appropriate analysis, which is a necessary prerequisite for sophisticated process-based interpretations [3].

There is an ongoing transition from classical, yet powerful, statistical methods to such characterization of the subtle aspects of the complex dynamics displayed by the recorded spatio-temporal fluctuations of relevant variables in different components of the Earth system, such as atmosphere, biosphere, cryosphere, lithosphere, oceans, the near-Earth electromagnetic environment, etc. In fact, many novel concepts originated in dynamical systems or information theory have been developed, partly motivated by specific research questions from the geosciences, and have found a variety of different applications [4]. This continuously extending toolbox of nonlinear time series analysis highlights the paramount importance of dynamical complexity as a key to understanding the behavior of the complex system Earth and its components.

It is important to note that there is no way to find an “optimum” organization or complexity measure [5]. Moreover, a combination of complexity measures, which refer to different aspects of a system, such as structural versus dynamical properties, is the most promising way to deal with the complexity of a system. For instance, several well-known techniques have been used to extract signatures of incipient changes in the system dynamics (i.e., precursors of dynamical transitions) hidden in time series representing natural signals of the Earth system. Notable examples include applications to time series data originated in the near-Earth electromagnetic environment (magnetospheric variations associated with magnetic storms), lithosphere (preseismic electromagnetic emissions associated with earthquakes), atmosphere (meteorological or hydrological signals), as well as other areas.

As a consequence of the described scientific “trends”, there is a large variety of methodological frameworks allowing us to study, characterize and, thus, understand complexity based on observational or model data. Among the main groups of conceptual approaches, methods originating in statistical mechanics and information theory belong to the most prominent and fruitful tools of modern complex systems-based nonlinear time series analysis [6,7]. Although statistical mechanics and information theory are based on somewhat different foundations, many concepts used in both areas are conceptually closely related. Most prominently, the notion of entropy naturally arises in both disciplines and provides a versatile basis for characterizing various complementary aspects associated with structural and dynamic complexity.

In this review, we discuss some recent developments, making use of statistical mechanics and information-theoretic notions of entropy for studying complexity in the Earth system. Specifically, we do not restrict ourselves to measuring the complexity of individual variables, but the wide-spread problem of characterizing dynamical interrelationships among a set of variables. Besides reviewing the algorithmic foundations of different approaches and highlighting their conceptual similarities and differences in Section 2, we provide an overview on recent successful applications in some exemplary fields of geoscience in Section 3. A discussion and summary of the potentials of statistical mechanics and information-theoretic approaches in the Earth sciences concludes this contribution.

2. Methods

2.1. Entropy Measures from Symbolic Sequences

2.1.1. Symbolic Dynamics

Symbolic time series analysis is a useful tool for modeling and characterization of nonlinear dynamical systems. It provides a rigorous way of looking at “real” dynamics with finite precision [8,9]. Briefly, it is a way of coarse-graining or simplifying the description.

The basic idea is quite simple. One divides the phase space into a finite number of classes and labels each class with a symbol (e.g., a letter from some alphabet). Instead of representing the trajectories by sequences of (typically real-valued) numbers (iterations of a discrete map or sampled points along the trajectories of a continuous flow), one watches the alteration of symbols. Of course, in doing so, one loses a considerable amount of detailed information, but some of the invariant, robust properties of the dynamics are typically kept, such as periodicity, symmetry or the chaotic nature of an orbit [8].

In the framework of symbolic dynamics, time series are transformed into a series of symbols by using an appropriate partition, which results in relatively few symbols. After symbolization, the next step is the construction of “symbol sequences” (“words” in the language of symbolic dynamics) from the symbol series by collecting groups of symbols together in temporal order.

To be more precise, the simplest possible coarse-graining of a time series is given by choosing a single threshold (usually the mean value or median of the data considered [10]), and assigning the symbols “1” and “0” to the signal, depending on whether it is above or below the threshold (binary partition) (The generalization to obtaining a classification using more than one threshold is trivial. Moreover, in addition to this kind of “static” encoding, there are other possibilities for obtaining a “dynamic” symbolization of observed time series data, such as using order patterns (cf. Section 2.1.4) or even mixed strategies [11]). Thus, we generate a symbolic time series from a two-letter ($\lambda =2$) alphabet (0,1), e.g., $0110100110010110\dots$. We can read this symbolic sequence in terms of distinct consecutive blocks of a length of n. For example, using $n=2$, one obtains $01/10/10/01/10/01/01/10/\dots$. This reading procedure is called “lumping”.

The number of all possible kinds of words is ${\lambda }^n=2^2=4$, namely 00, 01, 10, 11. The required probabilities for the estimation of an entropy, $p_{00}$, $p_{01}$, $p_{10}$, $p_{11}$, are the fractions of the different blocks (words), 00, 01, 10, 11, in the symbolic time series, namely, 0, 4/8, 4/8 and 0, correspondingly for the specific example. Based on these probabilities, we can estimate, for example, the probabilistic entropy measure, $H_S$, introduced by Shannon [12]:

$$H_S=-\sum _{j=1}^{{\lambda }^n}{p}_j{log}_2{p}_j$$

(1)

where $p_j$ are the probabilities associated with the microscopic configurations. The binary logarithm is taken here for convenience in the case of a binary symbolization (giving rise to an interpretation of the Shannon entropy as the negative mean amount of information, measured in bits), but may be replaced by another basis (e.g., using ${log}_{\lambda }$) without loss of generality.

In the case of a real-world time series, “experimental” conditions are often rather different (e.g., in terms of the system’s energy content), which can lead to ambiguities when quantitatively comparing entropy estimates, e.g., obtained from different persons in clinical studies. In order to correct for this problem, the renormalized entropy [5,13] has been introduced as the difference between the Shannon entropies obtained for the time series under study and some reference data.

A generalization of the Shannon entropy is provided by the spectrum of Rényi entropies [14]:

$$H_q=\frac{1}{1-q}{log}_2\sum _{j=1}^{{\lambda }^n}{p}_j^q$$

(2)

with $q\in \mathbb{R}$, which give different weights to rare and frequent symbols. In this spirit, Rényi entropies allow for the studying of scaling properties associated with a non-uniform probability distribution of symbols or words as arising, e.g., in the case of multi-fractal processes. Notably, the Shannon entropy, $H_S$, is a special case of $H_q$ for $q\to 1$. In a similar way, for $q\to 0$, $H_q$ is the Hartley entropy (or Hartley information), one of the most classical information-theoretic measures of uncertainty [15]. Note that for the case of a uniform probability distribution of n-blocks, all Rényi entropies coincide with each other and, hence, the Shannon and Hartley entropies.

2.1.2. Block Entropies

In general, one can statistically analyze a symbolic time series of N symbols, $\left\{A_i\right\},i=1,2,...,N$, in terms of the frequency distribution of blocks of a length of $n\left(n<N\right)$. For an alphabet consisting of λ letters, considering a word of a length of n, there are ${\lambda }^n$ possible combinations of the symbols that may constitute a word, here, ${\lambda }^n=2^n$. The probability of the occurrence, $p_j^{\left(n\right)}$, of the j-th combination of symbols ($j={1,2,...,2}^n$) (i.e., the j-th n-word) can be approximated by the number of occurrences of this n-word in the considered symbolic sequence divided by the total number of n-words.

Based on these probabilities, one can estimate the probabilistic entropy or other related characteristics. Various tools of information theory and entropy concepts can be used to identify statistical patterns in the symbolic sequences, onto which the dynamics of the original system under analysis has been projected. For detection of an anomaly, it suffices that a detectable change in the pattern represents a deviation of the system from nominal behavior [16]. Recently published work has reported novel methods for detection of anomalies in complex dynamical systems, which rely on symbolic time series analysis, e.g., [17,18]. Entropies depending on the word-frequency distribution in symbolic sequences are of special interest, extending Shannon’s classical definition of the entropy and providing a link between dynamical systems and information theory. These entropies take a large/small value if there are many/few kinds of patterns, i.e., they decrease while the organization of patterns is increasing. In this way, these entropies can measure the complexity of a signal.

Extending Shannon’s classical definition of the entropy of a single state [12] to the entropy of a succession of states [19,20], one reaches the definition of the n-block entropy, $H\left(n\right)$, which, for a two-letter ($\lambda =2$) alphabet and word length n, is given by:

$$H\left(n\right)=-\sum _{j=1}^{2^n}{p}_j^{\left(n\right)}{log}_2{p}_j^{\left(n\right)}$$

(3)

$H\left(n\right)$ is a measure of uncertainty and gives the average amount of information necessary to predict a subsequence of length n. Consequently, $H\left(n\right)/n$ can be interpreted as the mean uncertainty per symbol and should converge for $n\to \infty$ to some stationary value if the observed dynamics is deterministic. Moreover, from a practical perspective, one is often interested in quantifying the mean information gain when increasing the word length, measured by the conditional (or differential) block entropies:

$$h\left(n\right)=H(n+1)-H\left(n\right)\mathrm{for}n\ge 1;h\left(0\right)=H\left(1\right)$$

(4)

For stationary and ergodic processes, the limit of $h\left(n\right)$ for $n\to \infty$ provides an estimator of the Kolmogorov-Sinai entropy or source entropy h of the dynamical system under study [21,22] (In a similar spirit, $H\left(n\right)/n$ provides another estimate of the source entropy for $n\to \infty$.). For a broad class of systems, $h\left(n\right)$ converges exponentially fast, with the characteristic scaling parameter being related to the third-order Rényi entropy, $H_3$. In general, exploiting the corresponding convergence properties gives rise to a widely applicable probabilistic complexity measure, the effective measure complexity [23]:

$$EMC=\sum _{n=0}^{\infty }[h\left(n\right)-h]$$

(5)

A more extensive discussion and classification of complexity measures based on symbolic dynamics concepts has been given elsewhere [24].

2.1.3. Non-Extensive Tsallis Entropy

It has been established that physical systems characterized by long-range interactions or long-term memory or being of a multi-fractal nature are best described by a generalized statistical-mechanical formalism proposed by Tsallis [1,25,26]. More precisely, inspired by multi-fractal concepts, Tsallis introduced an entropic expression characterized by an index, q, which leads to non-extensive statistics [25,26]:

$$S_q=k\frac{1}{q-1}\left(1-\sum _{i=1}^W{p}_i^q\right)$$

(6)

where $p_i$ are probabilities associated with the microscopic configurations (i.e., the empirical frequency distribution of symbols), W is their total number, q is a real number and k is a constant, i.e., the Boltzmann’s constant from statistical thermodynamics. Notably, there is a striking conceptual similarity between Tsallis’ entropy definition and the notion of Rényi entropies.

The entropic index, q, describes the deviation of Tsallis entropy from the standard Boltzmann-Gibbs entropy. Indeed, using $p_i^{(q-1)}=e^{(q-1)ln\left(p_i\right)}\sim 1+(q-1)ln\left(p_i\right)$ in the limit $q\to 1$, we recover the usual Boltzmann-Gibbs entropy:

$$S_1=-k\sum _{i=1}^W{p}_{i}ln\left(p_i\right)$$

(7)

as the thermodynamic analog of the information-theoretic Shannon entropy.

For $q\ne 1$, the entropic index, q, characterizes the degree of non-extensivity reflected in the following pseudo-additivity rule:

$$\frac{S_q\left(A+B\right)}{k}=\frac{S_q\left(A\right)}{k}+\frac{S_q\left(B\right)}{k}+\left(q-1\right)\frac{S_q\left(A\right)}{k}\frac{S_q\left(B\right)}{k}$$

(8)

where A and B are two subsystems. If these subsystems have special probability correlations, extensivity does not hold for $q=1$ ($S_1\ne S_1\left(A\right)+S_1\left(B\right)$), but may occur for $S_q$, with a particular value of the index, $q\ne 1$. Such systems are called non-extensive [1,25]. The cases $q>1$ and $q<1$ correspond to sub-additivity or super-additivity, respectively. As in the case of Rényi entropies, we may think of q as a bias parameter: $q<1$ privileges rare events, while $q>1$ highlights prominent events [27].

We note that the parameter, q, itself is not a measure of the complexity of the system, but measures the degree of the non-extensivity of the system. In turn, the time variations of the Tsallis entropy, $S_q$, for some q quantify the dynamic changes of the complexity of the system. Specifically, lower $S_q$ values characterize the portions of the signal with lower complexity.

In terms of symbolic dynamics, the Tsallis entropy for a two-letter ($\lambda =2$) alphabet and word length n reads [18,28,29]:

$$S_q\left(n\right)=k\frac{1}{q-1}\left(1-\sum _{j=1}^{2^n}{\left[p_j^{\left(n\right)}\right]}^q\right)$$

(9)

Broad symbol-sequence frequency distributions produce high entropy values, indicating a low degree of organization. Conversely, when certain sequences exhibit high frequencies, low values of $S_q\left(n\right)$ are produced, indicating a high degree of organization.

2.1.4. Order-Pattern Based Approaches

In the symbolic dynamics framework introduced above, the transformation from a continuous to a discrete-valued time series (necessary for estimating the Shannon entropy of a discrete distribution or related characteristics) has been realized by coarse-graining the range of the respective observable. This procedure has the advantage of being algorithmically simple. However, unless block lengths $n\ge 2$ are used, the resulting entropies follow directly from the specific discretization and, thus, exclusively reflect statistical, but not dynamical, characteristics. In turn, dynamics comes into play when considering blocks of symbols.

An alternative way of introducing a symbolization in an explicitly dynamic way is making use of order patterns and related concepts from ordinal time series analysis [30]. In the simplest case (i.e., a pattern of order two), one considers a binary discretization of the underlying time series of a continuous random variable according to the sign of the difference between two subsequent values, e.g., using the symbols “0” if $s_{i+1}-s_i<0$ and “1” if $s_{i+1}-s_i>0$ (disregarding the “unlikely” case $s_{i+1}=s_i$). Notably, this kind of symbolization corresponds to the classical symbolic dynamics approach applied to the time series of increments (i.e., a first-order difference filter applied to the original data) with the increment value of zero discriminating the two classes of the resulting binary encoding. Based on this discretization, one can again define all entropies described above.

Besides using binary encodings based on order patterns of order two, the above framework can be extended to patterns of a higher order. In this case, sequences of observations are transformed into ordinal sequences (e.g., the case $s_i<s_{i+1}<s_{i+2}$ would correspond to the sequence $(1,2,3)$). Obviously, for a pattern of order q, there are $q!$ possible permutations of $(1,\cdots ,q)$ and, hence, at most $q!$ different patterns (Depending on the type of observed dynamics, there might be forbidden patterns. We do not further discuss the implications of this fact here.). Using the frequency distribution of the possible order patterns for calculating a Shannon-type entropy leads to the notion of permutation entropy [31,32,33].

Besides its direct application as a proxy for dynamical disorder based on time series data, Rosso et al. [34,35] suggested utilizing permutation entropy in combination with some associated measure of complexity for defining the so-called complexity-entropy causality plane. Here, the entropic quantity is the (Shannon-type) permutation entropy normalized by its maximum value, ${log}_{2}q!$:

$$H_S^{*}\left[P_{\pi }\right]=-\frac{1}{{log}_{2}q!}\sum _{i=1}^{q!}p\left({\pi }_i\right){log}_{2}p\left({\pi }_i\right)$$

(10)

where ${\pi }_i$ denote the different possible permutations of $(1,\cdots ,q)$. As one possible complexity measure associated with the permutation entropy, the Jensen-Shannon complexity, $C_{JS}\left[P_{\pi }\right]$, is defined as [36]:

$$C_{JS}\left[P_{\pi }\right]=Q_J[P_{\pi },P_e]H_S^{*}\left[P_{\pi }\right]$$

(11)

with the Jensen-Shannon divergence:

$$Q_J[P_{\pi },P_e]=Q_0{log}_{2}q!\left(H_S^{*}[(P_{\pi }+P_e)/2]-\frac{H_S^{*}\left[P_{\pi }\right]}{2}-\frac{1}{2}\right)$$

(12)

between the observed distribution, $P_{\pi }$, and the uniform distribution $P_e=(1/q!,\cdots ,1/q!)$, where $Q_0$ is a normalization constant, such that $0\le Q_J\le 1$.

2.2. Entropy Measures from Continuous Data

2.2.1. Approximate Entropy

Approximate entropy ($ApEn$) has been introduced by Pincus [37] as a measure for characterizing the regularity in relatively short and potentially noisy data. Specifically, $ApEn$ quantifies the degree of irregularity or randomness within a time series (of length N) and has already been widely applied to (non-stationary) biological systems for dynamically monitoring the system’s “health” status (see [38,39] and references therein). The conceptual idea is rooted in the work of Grassberger and Procaccia [40] and makes use of distances between sequences of successive observations. More specifically, $ApEn$ examines time series for detecting the presence of similar epochs; more similar and more frequent epochs lead to lower values of $ApEn$.

From a qualitative point of view, given N points, $ApEn$ is approximately equal to the negative logarithm of the conditional probability that two sequences that are similar for m points remain similar, that is, within a tolerance, r, at the next point. Smaller $ApEn$ values indicate a greater chance that a set of data will be followed by similar data (regularity). Thus, smaller values indicate greater regularity. Conversely, a larger value of $ApEn$ points to a lower chance of similar data being repeated (irregularity). Hence, larger values convey more disorder, randomness and system complexity. Consequently, a low/high value of $ApEn$ reflects a high/low degree of regularity. Notably, $ApEn$ detects changes in underlying episodic behavior not reflected in peak occurrences or amplitudes [41].

In the following, we give a brief description of the calculation of $ApEn$. A more comprehensive discussion can be found in [37,38,39].

For a time series, ${\left\{s_k\right\}}_k$ with $s_k=s\left(t_k\right)$, $t_k=kT$$,k=1,2,...,N$, and T being the sampling period, we can define $N-m+1$ vectors, each one consisting of m consecutive samples of this time series as:

$${\mathbf{X}}_i^m=\left\{s_i,s_{i+1},...,s_{i+m-1}\right\},i=1,2,...,N-m+1$$

(13)

The main idea is to consider a window of length m running through the time series and forming the corresponding vectors, ${\mathbf{X}}_i^m$. The similarity between the formed vectors is used as a measure of the degree of organization of the time series. A quantitative measure of this similarity, $C_i^m\left(r\right)$, is given by the average number of vectors, ${\mathbf{X}}_j^m$, within a distance, r, from ${\mathbf{X}}_i^m$. Here, ${\mathbf{X}}_j^m$ is considered to be within a distance r from ${\mathbf{X}}_i^m$ if $d_{ij}^m\le r$, where $d_{ij}^m$ is the maximum absolute difference of the corresponding scalar components of ${\mathbf{X}}_i^m$ and ${\mathbf{X}}_j^m$ (i.e., the two vectors have a distance smaller than r according to their supremum norm). By calculating $C_i^m\left(r\right)$ for each $i\le N-m+1$ and then taking the mean value of the corresponding natural logarithms:

$${\varphi }^m\left(r\right)={\left(N-m+1\right)}^{-1}\sum _{i=1}^{N-m+1}ln{C}_i^m\left(r\right)$$

(14)

the $ApEn$ is defined as:

$$ApEn\left(m,r\right)=\underset{N\to \infty }{lim}[{\varphi }^m\left(r\right)-{\varphi }^{m+1}\left(r\right)]$$

(15)

which, for a finite time series, can be estimated by the statistic:

$$ApEn\left(m,r,N\right)={\varphi }^m\left(r\right)-{\varphi }^{m+1}\left(r\right)$$

(16)

By tuning r, one can reach a reasonable degree of similarity for “most” vectors, ${\mathbf{X}}_i^m$, and, hence, a reliable statistics, even for relatively short time series.

In summary, the presence of repetitive patterns of fluctuation in a time series renders it more predictable than a time series in which such patterns are absent. A time series containing many repetitive patterns has a relatively small $ApEn$; a less predictable (i.e., more complex) process has a higher $ApEn$.

2.2.2. Sample Entropy

Sample entropy ($SampEn$) was proposed by Richman and Moorman [42] as an alternative that would provide an improvement of the intrinsic bias of $ApEn$ [43]. Specifically, when calculating the similarity within a time series using the vectors defined by Equation (13), $ApEn$ does not exclude self-matches, since the employed measure of similarity, $C_i^m\left(r\right)$, which is proportional to the number of vectors within a distance, r, from ${\mathbf{X}}_i^m$, does not exclude the ${\mathbf{X}}_i^m$ itself from this count. $SampEn$ is considered to be an evolution of $ApEn$ that displays relative consistency and less dependence on data length (see, e.g., [44] and the references therein).

First of all, we consider only the first $N-m$ vectors of a length of m of Equation (13), ensuring that, for $1\le i\le N-m$, ${\mathbf{X}}_i^m$ and ${\mathbf{X}}_i^{m+1}$ are defined. Then, we define $B_i^m\left(r\right)$ as ${\left(N-m-1\right)}^{-1}$ times the number of vectors, ${\mathbf{X}}_j^m$, within a distance, r, from ${\mathbf{X}}_i^m$, where $j=1,2,...,N-m$, but also $j\ne i$, excluding in this way the self-matches, and set:

$$B^m\left(r\right)={\left(N-m\right)}^{-1}\sum _{i=1}^{N-m}{B}_i^m\left(r\right)$$

(17)

Correspondingly, we define $A_i^m\left(r\right)$ as ${\left(N-m-1\right)}^{-1}$ times the number of vectors, ${\mathbf{X}}_j^{m+1}$, within a distance, r, from ${\mathbf{X}}_i^{m+1}$, where $j=1,2,...,N-m$ with $j\ne i$, and:

$$A^m\left(r\right)={\left(N-m\right)}^{-1}\sum _{i=1}^{N-m}{A}_i^m\left(r\right)$$

(18)

The $SampEn\left(m,r\right)$ is then defined as:

$$SampEn\left(m,r\right)=\underset{N\to \infty }{lim}\left\{-ln\left[A^m\left(r\right)/B^m\left(r\right)\right]\right\}$$

(19)

which, for finite time series, can be estimated by the statistic:

$$SampEn\left(m,r,N\right)=-ln\left[\frac{A^m\left(r\right)}{B^m\left(r\right)}\right]$$

(20)

2.2.3. Fuzzy Entropy

Fuzzy entropy ($FuzzyEn$) [43,44], like its ancestors, $ApEn$ and $SampleEn$ [44], is a “regularity statistic” that quantifies the (un)predictability of fluctuations in a time series. For the calculation of $FuzzyEn$, the similarity between vectors is defined based on fuzzy membership functions and the vectors’ shapes. The gradual and continuous boundaries of the fuzzy membership functions lead to a series of advantages, like the continuity, as well as the validity of $FuzzyEn$ at small values, higher accuracy, stronger relative consistency and even less dependence on the data length. $FuzzyEn$ can be considered as an upgraded alternative of $SampEn$ (and $ApEn$) for the evaluation of complexity, especially for short time series contaminated by noise.

Like $SampEn$, $FuzzyEn$ excludes self-matches. However, it follows a slightly different definition for the employed first $N-m$ vectors of a length of m, by removing a baseline, ${\overline{s}}_i$:

$${\overline{s}}_i=m^{-1}\sum _{j=0}^{m-1}{s}_{i+j}$$

(21)

i.e., for the $FuzzyEn$ calculations, we use the first $N-m$ of the vectors:

$${\mathbf{X}}_i^m=\left\{s_i,s_{i+1},...,s_{i+m-1}\right\}-{\overline{s}}_i,i=1,2,...,N-m+1$$

(22)

Then, the similarity degree, $D_{ij}^m$, between each pair of vectors, ${\mathbf{X}}_j^m$ and ${\mathbf{X}}_i^m$, being within a distance, r, from each other, is defined by a fuzzy membership function:

$$D_{ij}^m=\mu \left(d_{ij}^m,r\right)$$

(23)

where $d_{ij}^m$ is, as in the case of $ApEn$ and $SampEn$, the supremum norm difference between ${\mathbf{X}}_i^m$ and ${\mathbf{X}}_j^m$. For each vector, ${\mathbf{X}}_i^m$, we calculate the average similarity degrees with respect to all other vectors, ${\mathbf{X}}_j^m$, $j=1,2,...,N-m+1,\text{and}j\ne i$ (i.e., excluding itself):

$${\varphi }_i^m\left(r\right)={\left(N-m-1\right)}^{-1}\sum _{i=1,j\ne i}^{N-m}{D}_{ij}^m$$

(24)

Then, we evaluate:

$${\phi }^m\left(r\right)={\left(N-m\right)}^{-1}\sum _{i=1}^{N-m}{\varphi }_i^m\left(r\right)$$

(25)

and:

$${\phi }^{m+1}\left(r\right)={\left(N-m\right)}^{-1}\sum _{i=1}^{N-m}{\varphi }_i^{m+1}\left(r\right)$$

(26)

The $FuzzyEn\left(m,r\right)$ is then defined as:

$$FuzzyEn\left(m,r\right)=\underset{N\to \infty }{lim}\left[ln{\phi }^m\left(r\right)-ln{\phi }^{m+1}\left(r\right)\right]$$

(27)

which, for finite time series, can be estimated by the statistic:

$$FuzzyEn\left(m,r,N\right)=ln{\phi }^m\left(r\right)-ln{\phi }^{m+1}\left(r\right)$$

(28)

We note that by making use of vectors of different lengths, $ApEn$, $SampEn$ and $FuzzyEn$ are based on a similar rationale as the conditional block entropies (i.e., changes in the information content as the sequence length increases).

2.3. Measures of Statistical Interdependence and Causality

The above-mentioned concepts make use of the general ideas of information-theoretic or thermodynamic entropy to gain information on the statistical complexity of a single time series. A natural extension of this idea is utilizing comparable approaches for studying complex signatures of interactions between two time series, or even among $K>2$ records. In the following, we briefly review the corresponding basic concepts as used in the geoscientific context, so far.

2.3.1. Mutual Information

While Shannon entropy is a measure of the uncertainty about outcomes of one process, mutual information (MI) is a measure of the reduction thereof if another process is known (There are several other, conceptually-related symmetric measures for the strength of interdependencies between two variables, e.g., the maximal information coefficient [45]. For brevity, we do not further discuss these alternative measures here.). The Shannon-type MI can be expressed as:

$$\begin{array}{cc}\hfill I(X;Y)& =H_S\left(Y\right)-H_S\left(Y\right|X)=H_S\left(X\right)-H_S\left(X\right|Y)\hfill \end{array}$$

(29)

$$\begin{array}{cc}& =H_S\left(X\right)+H_S\left(Y\right)-H_S(X,Y)\hfill \end{array}$$

(30)

i.e., as the difference between the uncertainty in Y and the remaining uncertainty if X is already known (and vice versa). MI is symmetric in its arguments, non-negative and zero, if and only if X and Y are independent. The lagged cross-MI for two time series is given by:

$$\begin{array}{c}\hfill I_{XY}\left(\tau \right)\equiv I(X_{t-\tau };Y_t)\end{array}$$

(31)

For $\tau >0$, one measures the information in the past of X that is contained in Y, and vice versa for $\tau <0$. In analogy, the auto-MI is defined as $I(Y_{t-\tau };Y_t)$ for $\tau >0$ (for $\tau =0$, this gives the classical Shannon entropy, $H_S\left(Y\right)$).

Just like the entropy of a time series can be estimated based on symbols or patterns, MI can be estimated by simply plugging in the corresponding marginal and joint entropies into Equation (30). In contrast to pattern-based estimators, for binning estimators (i.e., estimators using frequency distributions based on a bivariate symbolization with $\lambda \gg 2$ symbols for each time series, but typically considering individual symbols, i.e., $n=1$), an adaptive strategy is preferable in which the bins are chosen, such that the marginal distributions are uniform [46,47,48]. For continuous time series, a useful class of estimators is based on nearest-neighbor statistics (see Section 2.4.1.). For example, making use of the entropy estimators developed in [49], Kraskov et al. introduced a nearest-neighbor estimator of MI [47].

In the limiting case of a bivariate Gaussian distribution, the MI is simply given by:

$$\begin{array}{c}\hfill I_{\mathrm{Gauss}}(X;Y)=-\frac{1}{2}ln\left(1-\rho {(X;Y)}^2\right)\end{array}$$

(32)

where $\rho (X;Y)$ is the linear Pearson correlation coefficient.

There are several straightforward generalizations of the mutual information concept. On the one hand, one can characterize nonlinear bivariate interrelationships between two variables by a spectrum of generalized mutual information functions obtained by replacing the Shannon entropy, $H_S$, by Rényi entropies, $H_q$, of arbitrary order q [50] (Note that this replacement does not conserve the non-negativity property of the classical MI function [50].). On the other hand, the idea behind mutual information can be directly extended to sets of $K>2$ variables, $X_1,\cdots ,X_K$, leading to the so-called (generalized) redundancies [51,52,53]:

$$I_q(X_1;\cdots ;X_K)=\sum _{k=1}^K{H}_q\left(X_k\right)-H_q(X_1,\cdots ,X_K)$$

(33)

2.3.2. Conditional Mutual Information and Transfer Entropy

An important generalization of MI is the conditional mutual information (CMI):

$$\begin{array}{cc}\hfill I(X;Y|Z)& =H_S\left(Y\right|Z)-H_S\left(Y\right|X,Z)=H_S\left(X\right|Z)-H_S\left(X\right|Y,Z)\hfill \end{array}$$

(34)

$$\begin{array}{cc}& =H_S(X,Z)+H_S(Y,Z)-H_S\left(Z\right)-H_S(X,Y,Z)\hfill \end{array}$$

(35)

which measures the mutual information between X and Y that is not contained in a third variable, Z, which can also be multivariate. CMI shares the properties of MI and is zero, if and only if X and Y are mutually independent conditionally on Z. Moreover, it is possible to use any of the existing estimators of entropic characteristics (including binning, order patterns or k-nearest neighbors) for defining CMI and generalized versions thereof. Like in Equation (32), the CMI for multivariate Gaussian processes can be expressed in terms of the partial correlation, where the correlation, $\rho (X;Y)$, is replaced by $\rho (X;Y|Z)$.

Of particular interest is the use of information-theoretic measures to quantify causal relations between variables based on time series data. A first step is to assess the directionality of information transfer between two time series. While lagged MI can be used to quantify whether information in Y has already been present in the past of X, this information could have existed in the common past of both processes and, therefore, is not a sign of a transfer of unique information in X to Y. To overcome this limitation, transfer entropy (TE) [54] measures the information in the past of X that is shared in the present of Y and not already contained in the past of Y:

$$\begin{array}{cc}\hfill I_{X\to Y}^{\mathrm{TE}}& \equiv I(X_t^{-};Y_t|Y_t^{-})\hfill \end{array}$$

(36)

with the upper script “${}^{-}$” denoting the (infinite) past vector, e.g., $X_t^{-}=(X_{t-1},X_{t-2},\dots )$. It should be noted that the notion of causality employed here is in the predictive sense of Granger causality, which was introduced for linear processes in [55] and later generalized to nonlinear frameworks [56,57]. TE and Granger causality can actually be shown to be equivalent for Gaussian variables [58].

2.3.3. Graphical Models

The definition of TE leads to the problem that infinite-dimensional densities have to be estimated, which is commonly called the “curse of dimensionality”. In the usual naive estimation of TE, the infinite vectors are simply truncated at some ${\tau }_{max}$, which has to be chosen at least as large as the supposed maximal coupling delay between X and Y. This can lead to very large estimation dimensions, since the delay is not known a priori. As shown in [59,60], a high estimation dimension severely affects the reliability of the inference of directionality. In [59], the problem of high dimensionality is overcome by utilizing the concept of graphical models.

While the TE discussed before is a bivariate measure of directionality, a more comprehensive causality assessment requires the inclusion of multiple variables to be able to exclude common drivers and indirect causalities. In the graphical model approach [61,62,63], the conditional independence properties of a multivariate process are visualized in a graph; in our case, a time series graph. As depicted in Figure 1a, each node in that graph represents a subprocess of a multivariate discrete-time process, $\mathbf{X}$, at a certain time, t. Nodes $X_{t-\tau }$ and $Y_t$ are connected by a directed link “$X_{t-\tau }\to Y_t$”, if and only if $\tau >0$ and:

$$\begin{array}{c}\hfill I(X_{t-\tau };Y_t|{\mathbf{X}}_t^{-}\setminus \left\{X_{t-\tau }\right\})>0\end{array}$$

(37)

i.e., if they are not independent conditionally on the past of the whole process, which implies a lag-specific causality with respect to $\mathbf{X}$. If $Y\ne X$, we say that the link “$X_{t-\tau }\to Y_t$” represents a coupling at lag τ, while for $Y=X$, it represents an auto-dependency at lag τ. Nodes $X_t$ and $Y_t$ are connected by an undirected contemporaneous link “$X_t-Y_t$” (visualized by a line, cf. in Figure 1a) [63], if and only if:

$$\begin{array}{c}\hfill I(X_t;Y_t|{\mathbf{X}}_{t+1}^{-}\setminus \{X_t,Y_t\})>0\end{array}$$

(38)

where also the contemporaneous present ${\mathbf{X}}_t\setminus \{X_t,Y_t\}$ is included in the condition. Note that stationarity implies that “$X_{t-\tau }\to Y_t$” whenever “$X_{t^{\prime }-\tau }\to Y_{t^{\prime }}$” for any $t^{\prime }$. In the resulting graph, the parents and neighbors of a node $Y_t$ are defined as:

$$\begin{array}{cc}\hfill {\mathcal{P}}_{Y_t}& \equiv \{Z_{t-\tau }:Z\in \mathbf{X},\tau >0,Z_{t-\tau }\to Y_t\}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\mathcal{N}}_{Y_t}& \equiv \{X_t:X\in \mathbf{X},X_t-Y_t\}\hfill \end{array}$$

(40)

In [59], an algorithm for the estimation of time series graphs by iteratively inferring the parents is introduced. The supplementary material of [59] also provides a suitable shuffle test and a numerical study on the detection and false positive rates of the algorithm.

Figure 1. (Color online) Causal interactions in a multivariate process ${\mathbf{X}}_t={(X,Y,Z)}_t$. (a) The time series graph (see definition in text). The set of parents, ${\mathcal{P}}_{Y_t}$ (blue box, gray nodes), separates $Y_t$ from the past of the whole process, ${\mathbf{X}}_t^{-}\setminus {\mathcal{P}}_{Y_t}$, which is used in the algorithm to estimate the graph. (b) The process graph, which aggregates the information in the time series graph for better visualization (edge labels denote the associated lags).

Figure 1. (Color online) Causal interactions in a multivariate process ${\mathbf{X}}_t={(X,Y,Z)}_t$. (a) The time series graph (see definition in text). The set of parents, ${\mathcal{P}}_{Y_t}$ (blue box, gray nodes), separates $Y_t$ from the past of the whole process, ${\mathbf{X}}_t^{-}\setminus {\mathcal{P}}_{Y_t}$, which is used in the algorithm to estimate the graph. (b) The process graph, which aggregates the information in the time series graph for better visualization (edge labels denote the associated lags).

Time series graphs encode the existence of a lag-specific causality, but they are not meant to assess the causal strength in a meaningful way. This topic is discussed in [64], where the concept of momentary information transfer (MIT) [65] is utilized for attributing a well-interpretable coupling strength solely to the inferred links of the time series graph. MIT between X at some lagged time, $t-\tau$, in the past and Y at time t is the CMI that measures the part of the source entropy of Y that is shared with the source entropy of X:

$$\begin{array}{cc}\hfill I_{X\to Y}^{\mathrm{MIT}}\left(\tau \right)& \equiv I(X_{t-\tau };Y_t|{\mathcal{P}}_{Y_t}\setminus \left\{X_{t-\tau }\right\},{\mathcal{P}}_{X_{t-\tau }})\hfill \end{array}$$

(41)

The attribute momentary is used, because MIT measures the information of the “moment” $t-\tau$ in X that is transferred to $Y_t$. Similar to the definition of contemporaneous links in Equation (38), we can also define a contemporaneous MIT:

$$\begin{array}{cc}\hfill I_{X-Y}^{\mathrm{MIT}}& \equiv I(X_t;Y_t|{\mathcal{P}}_{Y_t},{\mathcal{P}}_{X_t},{\mathcal{N}}_{X_t}\setminus \left\{Y_t\right\},{\mathcal{N}}_{Y_t}\setminus \left\{X_t\right\})\hfill \end{array}$$

(42)

MIT is well interpretable, because it excludes past information, not only from Y, but also from X, i.e., information transmitted by serial correlations. Often X and Y are entangled by multiple links, and only the exclusion of both pasts truly yields a coupling measure that depends solely on the mutual coupling strength. This is demonstrated analytically and numerically in [64].

2.4. Linkages with Phase Space Methods

2.4.1. Reinterpreting Entropies As Phase Space Characteristics

In their most common formulation, many of the aforementioned characteristics are explicitly based on a discretization of the system under study (either via coarse-graining and subsequent symbolization or using order patterns or related approaches), giving rise to an information-theoretic interpretation of entropies and related complexity measures. Consequently, these measures are often considered as providing one major field of nonlinear time series analysis methods, which is opposed to another class of well-developed theoretical frameworks making use of the notion of the phase space of a dynamical system. In the latter case, one commonly takes sets of variables (or time-lagged replications of one variable) into account, which are assumed to span a metric space in which the system’s behavior can be fully described by means of distances between sets of “state vectors” sampled from the observed trajectory.

Reconsidering the entropic measures described above in more detail, one recognizes that there are also multiple quantifiers that make explicit use of distances (i.e., approximate, sample and fuzzy entropy), which implies that the latter methods can also be considered as phase space-based approaches if one identifies subsequent observations as individual “coordinates” of some abstract space (i.e., a “naive” embedding space of variable dimension). In fact, as already stated above, the mentioned class of entropy characteristics is based on ideas formulated by Grassberger and Procaccia about 30 years ago, which are among the pioneering works utilizing the concept of phase space for nonlinear time series analysis.

In a similar spirit, state-of-the-art estimators of mutual information and causality concepts based on this framework make use of the phase space concept by replacing the traditional histogram-based symbolizations or order patterns by nearest-neighbor relationships between sampled state vectors. Specifically, k-nearest neighbor estimators of entropies [49], MI [47] and CMI [66] have been found practically useful and have desirable numerical properties not shared by other estimators. Here, we just briefly illustrate the basic idea of k-nearest neighbor estimates for CMI.

As the first step, the number of neighbors, k, to be considered around each sampled state vector in the joint space of $(X,Y,Z)$ is chosen as a free parameter. For each sample corresponding to observations at time t, the supremum norm distance, ${\epsilon }_t$, to the k-th nearest neighbor defines a cube of a length of $2{\epsilon }_t$ around the corresponding state vector. Next, the numbers of points in the respective subspaces with a distance less than ${ϵ}_t$ are counted, yielding $k_{z,t}$, $k_{xz,t}$ and $k_{yz,t}$. Then, the CMI is estimated as:

$$\begin{array}{c}\hfill \widehat{I}(X;Y|Z)=\psi \left(k\right)+\frac{1}{T}\sum _{t=1}^T\left[\psi \left(k_{z,t}\right)-\psi \left(k_{xz,t}\right)-\psi \left(k_{yz,t}\right)\right]\end{array}$$

(43)

where ψ is the Digamma function. Smaller k result in smaller cubes, and since (as assumed in the estimator’s derivation) the density is approximately constant in these cubes, the estimator has a low bias. Conversely, for large k, the bias is stronger, but the variance gets smaller. Note, however, that for independent processes, the bias is approximately zero, i.e., $\widehat{I}(X;Y|Z)\approx 0$; large k are therefore better suited for independence tests, as needed in the algorithm to infer the time series graph.

2.4.2. Other Entropy Measures from Phase Space Methods

Besides the above-mentioned examples, entropic quantities frequently occur in various phase space-based methods of nonlinear time series analysis. One prominent example is the Grassberger-Procaccia algorithm for estimating the correlation dimension of chaotic attractors [67,68]. As in the calculation of $ApEn$, $SampEn$ and $FuzzyEn$, let us consider the m-dimensional delay vectors, ${\mathbf{X}}_i^m$, where the embedding dimension, m, and the mutual temporal offset between the components of ${\mathbf{X}}_i^m$ have been optimized, so as to yield coordinates, which are as independent as possible and represent the “true” dimensionality of the recorded system appropriately. Computing the pairwise distance of the resulting delay vectors, we may distinguish between mutually “close” and “distant” vectors by comparing the distance with some predefined global threshold value, r, which is summarized in the recurrence matrix [69]:

$${\mathbf{R}}_{ij}\left(r\right)=\Theta (r-\parallel {\mathbf{X}}_i^m-{\mathbf{X}}_j^m\parallel )$$

(44)

where $\Theta (·)$ is the Heaviside “function”. The fraction of mutually close pairs of vectors:

$$C^m\left(r\right)=\frac{2}{N^{\prime }(N^{\prime }-1)}\sum _{i<j}{\mathbf{R}}_{ij}\left(r\right)$$

(45)

(with $N^{\prime }=N-m+1$) is called the correlation sum or recurrence rate of the observed trajectory. For a dissipative chaotic system with attractor dimension $<m$, the latter obeys a characteristic scaling behavior as:

$$C^m\left(r\right)\sim r^{D_2}exp(-m{H}_2\Delta t)$$

(46)

with $\Delta t$ denoting the sampling interval of the time series data and $D_2$ the attractor’s correlation dimension. By varying r and m, it is thus possible to estimate the second-order Rényi entropy, $H_2$, together with the correlation dimension from the correlation sum.

Alternatively, other properties of the recurrence matrix may be exploited for estimating different entropic quantities. Among others, the length distribution, $p\left(l\right)$, of “diagonal lines” formed by non-zero entries in $\mathbf{R}$ provides a versatile tool for estimating $H_2$ from time series data [70]. A similar approach can be used for estimating the generalized MI of second order. The Shannon entropy, $ENTR$, of the discrete distribution, $p\left(l\right)$, has been used as a heuristic complexity measure in the framework of the so-called recurrence quantification analysis [69]. However, it has been shown that this definition does not allow for properly distinguishing between different “degrees” of chaotic dynamics in the desired way. In turn, the Shannon entropy associated with the length distribution, $p\left(\overline{l}\right)$, of diagonal lines in $\mathbf{R}$ formed exclusively by zero entries provides a measure of dynamical disorder that increases with the complexity of chaotic dynamics as measured by the system’s largest Lyapunov exponent [71]. Even more, the latter quantity shows an excellent agreement with the classical Shannon entropy estimate based on symbolic dynamics, but does not share some of the disadvantages, due to the coarse-graining, which arise in the latter case.

It is worth noting that, recently, there have been the first successful attempts to transfer the idea of “recurrences” in phase space to symbolic data [11,72,73]. A special case is order pattern recurrence plots [74,75], which have been successfully applied in the analysis of neurophysiological time series.

2.4.3. Causality from Phase Space Methods

For the sake of completeness, we briefly mention that similar to the characterization of dynamical complexity with entropic characteristics, also the complementary bi- and multi-variate problem of detecting causal interrelationships can be complementarily tackled by information-theoretic and phase space methods. Distance-based characteristics allowing one to reliably detect directional couplings between different variables have been developed and applied by various authors, e.g., [76,77,78,79,80,81,82]. A specific class of approaches is based on the recurrence matrix, making use of conditional recurrence probabilities [83,84] or asymmetries in coupled network statistics [85]. Since these methods are not within the focus of this review, we refer to the aforementioned references for further details.

3. Applications

3.1. Space Weather and Magnetosphere

During the past two decades, it has become increasingly evident that the magnetosphere cannot be viewed as an autonomous system, but should rather be examined in terms of driven nonlinear dynamical systems [86]. As early as 1990, Tsurutani et al. [87] provided one of the first indications of a nonlinear process in the magnetosphere, claiming that there is a nonlinear response of the AE index to the southward component of the interplanetary magnetic field (IMF). At the same time, Baker et al. were able to provide a mechanical analogue to model the driving of the magnetosphere, concluding that using tools of chaos theory, we may gain considerable insight into the storm-time behavior of the magnetosphere [88]. The discussion on low-dimensional chaos in magnetospheric activity was continued in many important works [89,90,91,92]. On the other hand, Chang and co-workers [93,94,95,96] suggested that the magnetospheric dynamics is that of an infinite-dimensional nonlinear system near criticality, a hypothesis that was later supported by statistical analyses of the AE index [97,98] and auroral observations obtained by the Ultraviolet Imager (UVI) on the Polar spacecraft [99,100].

In situ observations revealed that scale-invariance, turbulence, non-Gaussian fluctuations and intermittency are essential elements of plasma sheet dynamics [101,102,103], while both in situ and global observations pointed out the possible emergence of dynamical phase transitions and symmetry breaking phenomena in coincidence with the occurrence of local and global relaxation processes [104,105,106,107].

Due to the above, it becomes obvious that the classical magnetohydrodynamic (MHD) approach is not sufficient to describe the more dynamical aspects of many magnetospheric phenomena, since the emergence of long-range correlations, as well as the coupling between various scales demand the utilization of novel concepts that were developed with respect to the multi-scale dynamics of out-of-equilibrium systems [28,29,108,109,110].

In this spirit, multi-fractal models have been employed for the description of space weather [111]. In addition, the emergence of self-organized criticality (SOC) has been examined in both AE index and solar wind data [112,113], as a means of understanding the global energy storage and release in the coupled solar wind–magnetosphere system [114,115] and as a paradigm for the magnetotail dynamics [116]. Moreover, the probability distributions of physical variables in turbulent space plasmas are scale-dependent and perfectly well modeled by non-extensive one-parameter (Tsallis) distributions, where q measures the degree of correlations (long-range interactions) in the system [117].

Magnetic storms are the most prominent global phenomenon of geospace dynamics, interlinking the solar wind, magnetosphere, ionosphere, atmosphere and, occasionally, the Earth’s surface [118,119,120]. Magnetic storms produce a number of distinct physical effects in the near-Earth space environment: acceleration of charged particles in space, intensification of electric currents in space and on the ground, impressive aurora displays and global magnetic disturbances on the Earth’s surface [121]. The latter serve as the basis for storm monitoring via the hourly disturbance storm time (Dst) index, which is computed from an average over four mid-latitude magnetic observatories [122].

The Dst data used here for illustrative purposes (Figure 2) include two intense magnetic storms, which occurred on March 31, 2001, and November 6, 2001, with minimum Dst values of −387 nT (nanoTesla) and −292 nT, respectively, as well as a number of smaller events (e.g., May and August, 2001, with Dst <−100 nT in both cases).

All entropy calculations were performed on 256 sample long windows sliding along the time series with a step of one sample. For the case of symbolic Tsallis entropy, the value $q=1.84$ [123] was found for the non-extensivity parameter, q. For the $FuzzyEn$ calculations, the exponential function $\mu \left(d_{ij}^m,r\right)=exp\left(-{\left(d_{ij}^m/r\right)}^n\right)$ with $n=2$ has been used as the fuzzy membership function, using $m=2$ and $r=0.65·STD$, where $STD$ is the standard deviation of the analyzed time series fragment, allowing fragments with different amplitudes to be compared quantitatively. The same values for m and r were also used in the case of $ApEn$ and $SampEn$ calculations. The results of the corresponding analyses are shown in Figure 2.

Approximate entropy, sample entropy, Fuzzy entropy and non-extensive Tsallis entropy sensitively trace the complexity dissimilarity among different “physiological” (quiet times) and “pathological” (intense magnetic storms) states of the magnetosphere. They imply the emergence of two distinct patterns: (1) a pattern associated with the intense magnetic storms, which is characterized by a higher degree of organization; and (2) a pattern associated with normal periods, which is characterized by a lower degree of organization. In general, all four entropy measures clearly distinguish between the different complexity regimes in the Dst time series (see the red part of the corresponding plots in Figure 2).

Figure 2. From top to bottom: Dst time series from January 1, 2001, to December 31, 2001, associated with two intense magnetic storms (marked in red) and the corresponding $ApEn$, $SampEn$, $FuzzyEn$ and non-extensive Tsallis entropy, $S_q$, for sliding windows in time. The common horizontal axis is the time (in days), denoting the relative time position from the beginning of the Dst data. The triangles denote five time intervals in which the first, third and fifth intervals correspond to higher entropies, whereas the second and fourth time windows exhibit lower entropies (cf. the parts of the entropy plots shown in red).

Figure 2. From top to bottom: Dst time series from January 1, 2001, to December 31, 2001, associated with two intense magnetic storms (marked in red) and the corresponding $ApEn$, $SampEn$, $FuzzyEn$ and non-extensive Tsallis entropy, $S_q$, for sliding windows in time. The common horizontal axis is the time (in days), denoting the relative time position from the beginning of the Dst data. The triangles denote five time intervals in which the first, third and fifth intervals correspond to higher entropies, whereas the second and fourth time windows exhibit lower entropies (cf. the parts of the entropy plots shown in red).

Figure 2 indicates that $ApEn$ and $S_q$ yield superior results in comparison with the other entropy measures regarding the detection of dynamical complexity in the Earth’s magnetosphere (i.e., they offer a clearer picture of the transition). A possible explanation for this is that, on the one hand, $S_q$ is an entropy obeying a non-extensive statistical theory, which is different from the classical Boltzmann-Gibbs statistical mechanics. Therefore, it is expected to better describe the dynamics of the magnetosphere, which is a non-equilibrium physical system with large variability. On the other hand, $ApEn$ is more stable when dealing with non-stationary signals from dynamical systems (such as the magnetospheric signal) than the other entropy measures presented here.

Figure 2 could also serve for placing boundary values or thresholds for $ApEn$ and Tsallis entropy in order to distinguish the different magnetospheric states. Along these lines, we suggest the limits of 0.2 and 0.7 for $ApEn$ and Tsallis entropy, respectively. Similarly, the values of 0.15 and 0.1 could play the role of indicative limits for the transition between the quiet-time and stormy magnetosphere in the cases of $SampEn$ and $FuzzyEn$, respectively.

3.2. Preseismic Electromagnetic Emissions

Under specific conditions, the Earth’s heterogeneous crust experiences sudden releases of energy following large-scale mechanical failures. These phenomena of a potentially catastrophic nature are known as earthquakes (EQs). In a simplified view, a significant EQ is what happens when the two surfaces of a major fault slip, one over the other, under the stresses rooted in the motion of tectonic plates. Due to the catastrophic nature of some significant EQs and their direct impact on human life, there is a great interest in all related phenomena, especially those that are systematically observed prior to EQ occurrence and could therefore be regarded as potential EQ precursors.

Electromagnetic (EM) phenomena associated with EQs have been repeatedly recognized as soon as human history began to be recorded. The EQ light phenomenon has been reported all around the world already since ancient times [124], given that no specific instruments, but only eye-witnesses were necessary for the observations. However, since the existence of appropriate instrumentation for their measurement, a wide variety of EM phenomena have been observed [125]. For example, the existence of magnetic effects associated with earthquakes or volcanic activity was first suggested by Kalashnikov in 1954 [126].

Focused research in the field, which now is usually referred to as “Seismo-Electromagnetics”, practically started in the 1980s, extending from direct current (DC) to very high frequency (VHF) bands, while it concerns the theoretical and experimental study of either EM signals supposedly emitted from the focal zone or anomalous transmission of EM waves over the epicentral region [124]. From that time onwards, a large number of telluric current anomalies (e.g., [127,128]), ultra-low frequency (ULF) anomalies (e.g., [129,130,131,132,133]), very low frequency (VLF) to VHF anomalies (e.g., [133,134,135,136,137,138] and references therein), as well as a variety of atmospheric-ionospheric anomalies attributed to the phenomenon known as lithosphere-atmosphere-ionosphere (LAI) coupling (e.g., [124,139,140,141,142] and the references therein) have been reported, analyzed and interpreted. More information, historical data and reviews can be found in [124,125,143,144,145] and their references.

The identification of an observed EM anomaly, such as an EQ-related one, is a challenging research topic. The time-lead of the EM anomaly relative to the EQ occurrence per se is not credible evidence. Due to the multi-facet nature of EQ preparation processes, a multidisciplinary analysis of the EM anomalies is needed before its classification as EQ-related for the identification of features that characterize the EQ preparation processes.

Since the early work of Mandelbrot [146], the self-affine nature of faulting and fracture is widely documented based on the analysis of data from both field observations and laboratory experiments; cf. [147] and the references therein. Moreover, several characteristics of seismicity and fracture (like criticality, complexity, frequency-magnitude laws, etc.) have been investigated both in the laboratory and at the geophysical scale (e.g., [148,149,150,151,152,153] and the references therein). These findings led to the hypothesis that any EM emissions related to EQ preparation should present specific characteristics, like fractal scaling, criticality, phase transition characteristics, symmetry breaking, frequency-size laws, universal structural patterns of fracture and faulting, complexity characteristics, long-range correlations, memory, etc., depending on the phase of the EQ preparation process to which they are related. A number of analyses in field-recorded preseismic EM signals have verified this hypothesis (e.g., [18,136,137,138,148,153,154,155,156,157,158,159] and references therein).

As previously stated, for systems exhibiting long-range correlations, memory or fractal properties, non-extensive Tsallis entropy is considered an appropriate mathematical tool [1,25,26,160]. A central property of the EQ preparation process is the occurrence of coherent large-scale collective behavior with a very rich structure, resulting from repeated nonlinear interactions among the constituents of the system [147,161]. Consequently, Tsallis entropy is also a physically meaningful tool for investigating the launch of a fracture-induced EM precursor.

A way to examine transient phenomena like preseismic EM emissions is to analyze the corresponding time series into a sequence of distinct time windows. The aim is to discover a clear difference in the dynamical characteristics as the catastrophic event is approaching. As already pointed out, $ApEn$, $SampEn$ and $FuzzyEn$ are especially useful for the analysis of short time series contaminated by noise. Since the recorded EM time-series are usually strongly contaminated by the background noise at the location of the field measurement station, $ApEn$, $SampEn$ and $FuzzyEn$ are particularly appropriate for their analysis in the form of distinct (short) time windows.

Seismicity is a critical phenomenon [147,148,155]. Thus, it is expected that a significant change in the statistical pattern, namely, the appearance of entropy “drops”, represents a deviation from the normal behavior of the EM recordings, revealing the presence of an EM anomaly. Specifically, we aim to reveal any increase of the organization level (entropy drop) in the observed EM time series as approaching the time of the EQ occurrence, by applying the different entropic metrics of symbolic Tsallis entropy, $ApEn$, $SampEn$ and $FuzzyEn$.

As an illustrative example, we consider here the preseismic kilohertz (kHz) EM activities associated with the Athens (Greece) earthquake ($M=5.9$) that occurred on 7 September 1999 [137,138,156,157,162,163,164,165,166]. All entropy calculations were performed on 1,024 sample long windows sliding along the time series with 25% mutual overlap, i.e., with a step of 768 samples. For the symbolic Tsallis entropy, the value $q=1.8$ [137] was found for the non-extensivity parameter, q. Here, the results for words of a length of $n=3$ are presented. $ApEn$, $SampEn$ and $FuzzyEn$ have been estimated using the same setting as for the Dst time series in Section 3.1.

Figure 3 shows the same analyses as those presented for the Dst index (Figure 2) for the 10 kHz preseismic EM time series associated with the Athens earthquake. It can be observed that lower entropy values, compared to that of the background noise, can be identi?ed between the time points ∼ 292,600 s and ∼729,000 s (green-colored part in Figure 3), although sparsely distributed in time. In turn, the entropy values suddenly drop during the two strong EM pulses at the tail of the analyzed excerpt of the Athens kHz signal, signifying a different behavior. This has been attributed to a new distinct phase in the tail of the EQ preparation process, which is characterized by a significantly higher degree of organization and lower complexity in comparison to that of the preceding phase (e.g., [153,157]).

Figure 3. From top to bottom: Part of the recorded time series of the 10 kHz (east-west) magnetic field strength (in arbitrary units) covering an 11-days period from August 28, 1999, 00:00:00 UT (Universal Time) , to September 7. 1999, 23:59:59 UT, associated with the Athens EQ, and the corresponding $ApEn$, $SampEn$, $FuzzyEn$ and non-extensive Tsallis entropy for sliding windows in time (see text). The common horizontal axis is the time (in days), denoting the relative time position from the beginning of the analyzed part of the EM recording. The preseismic EM emissions exhibit lower entropy values, which are observed in two distinct phases: the first phase characterized by lower entropy values sparsely distributed in time (cf. the part of the time series and entropy plots displayed in green), followed by a phase characterized by continuous significantly lower entropy values (time series and entropy plotted in red).

Figure 3. From top to bottom: Part of the recorded time series of the 10 kHz (east-west) magnetic field strength (in arbitrary units) covering an 11-days period from August 28, 1999, 00:00:00 UT (Universal Time) , to September 7. 1999, 23:59:59 UT, associated with the Athens EQ, and the corresponding $ApEn$, $SampEn$, $FuzzyEn$ and non-extensive Tsallis entropy for sliding windows in time (see text). The common horizontal axis is the time (in days), denoting the relative time position from the beginning of the analyzed part of the EM recording. The preseismic EM emissions exhibit lower entropy values, which are observed in two distinct phases: the first phase characterized by lower entropy values sparsely distributed in time (cf. the part of the time series and entropy plots displayed in green), followed by a phase characterized by continuous significantly lower entropy values (time series and entropy plotted in red).

From the obtained results, it seems that all employed entropic metrics display a comparable performance. However, $ApEn$, $SampEn$ and $FuzzyEn$ generally present stronger fluctuations with time, providing a higher resolution and a clearer difference of the EM anomalies, especially the two strong EM pulses, from the background noise baseline. It should be noted that the Tsallis entropy calculations presented in Figure 3 were performed without an explicit check of stationarity, considering that the window length of 1,024 points is small enough to provide pseudo-stationarity conditions.

3.3. Climate and Related Fields

Due to the relatively large amount of data on climate variability and the closely related aspects of Earth system dynamics (for example, hydrology and biogeochemistry), there have been numerous applications of statistical mechanics and information-theoretic concepts, including various entropic characteristics and complexity measures, in this field. Due to the large number of potential applications in this area, we restrict ourselves in the following to summarizing some examples of recent efforts related to certain specific climate-related studies, which can be considered representative for a great variety of potential research questions.

3.3.1. Complexity of Present-Day Climate Variability

The quantitative characterization of dynamical complexity in today’s Earth climate has been the subject of many studies. Notably, entropic characteristics and associated measures of complexity are particularly well suited for this purpose. For example, Paluš and Novotná [167] demonstrated the use of redundancies (i.e., multivariate generalizations of MI), in combination with surrogate data preserving the linear correlation properties of the original time series, for testing for the nonlinearity of different climatic variables (surface air pressure and temperature, geopotential height data), a necessary prerequisite for the emergence of chaotic dynamics.

On the global scale, von Bloh et al. [168] studied surface air temperature observations, as well as the corresponding outcomes of a global circulation model. In order to reveal the spatial characteristics associated with the nonlinear dynamics of the climate system, they considered the spatial patterns of the second-order Rényi entropy, $H_2$, estimated from recurrence plots (A similar analysis has been recently presented by Paluš et al. [169] in terms of a Gaussian process entropy rate computed for globally distributed surface air temperature anomalies). Their results suggested that state-of-the-art climate models (at the time of this study) have not been able to trace the fundamental spatial signatures of nonlinearity in observed climate dynamics, such as marked latitudinal, as well as land-ocean entropy gradients. The latter finding can be considered a well-established fact, which has been described for other nonlinear characteristics, such as the Hurst exponent describing the degree of long-range dependence in temperature fluctuations (see, e.g., [170,171,172,173]).

One reason for the observed complexity of temperature fluctuations is that the climate exhibits relevant variability on a multitude of different time-scales, from diurnal to seasonal and inter-annual scales. The longer time-scales are particularly highlighted by nonlinear large-scale oscillatory (dipole) patterns, such as the El Niño/Southern Oscillation (ENSO). Investigating the spatial profile and temporal variability of such patterns and characterizing their degree of dynamical disorder can provide a vast amount of complementary information not provided by classical linear analysis techniques. Among other entropic characteristics, Tsallis’ non-extensive entropy has been found particularly useful for this purpose, once more making use of the fact that the climate system is typically far from equilibrium.

Ausloos and Petroni [174] studied the normalized variability of the southern oscillation index (SOI) on different time-scales. Their results indicated a qualitative change in the non-extensivity as the considered scales vary, leading to a classical Boltzmann-Gibbs scaling ($q\to 1$) at larger scales. Beyond exclusive consideration of $S_q$, in the non-extensive statistical mechanics framework, characteristic scaling exponents can also be found for other quantities, which is formalized by Tsallis’ q-triplet [175], which has been obtained for daily ENSO fluctuations [176], but also other geoscientific variables, such as the depth of the stratospheric ozone layer [177], various meteorological variables, earthquake inter-occurrence times, magnetospheric dynamics, etc. [178].

In addition to non-extensive entropy, also other information-theoretic entropy concepts have great potentials for discriminating between climatological records regarding their respective dynamical complexity. As a recent example, we mention here the consideration of the complexity-entropy causality plane for the characterization of daily streamflow time series [36], vegetation dynamics based on remote sensing data [179] or the degree distribution of a complex network obtained from correlations among surface air temperature records in the ENSO region in the equatorial Pacific [180].

3.3.2. Dynamical Transitions in Paleoclimate Variability

In order to obtain reliable estimates of the dynamical complexity associated with the variability of large-scale climate modes, such as ENSO, on long time-scales, instrumental records may provide insufficient time coverage. Instead, it is feasible to utilize paleoclimate proxy data for obtaining information on nonlinear variability during the last few centuries, millennia or even the deeper past. In such situations, entropic characteristics can again provide interesting insights into qualitative changes in historical climate variability.

Saco et al. [181] studied the histogram-based Shannon entropy, as well as permutation entropy for a proxy record obtained from a sedimentary sequence from Laguna Pallcacocha, Ecuador. Their study revealed a distinct time interval with markedly reduced permutation entropy corresponding to a known period of large-scale rapid climate change (RCC) between 9,000 and 8,000 years before present (BP). The latter period includes the globally observed 8.2 kyr (kiloyear) event, which is believed to have originated from a strong meltwater pulse into the North Atlantic [182]. In terms of dynamical entropy, the obtained results indicate a considerable degree of regularity in the variability of environmental conditions. We hypothesize that this regularity indeed reflects the marked variability profile during the 8.2 kyr event, which exceeds the typical magnitude of “background fluctuations”. Notably, in the case of other, weaker RCC episodes during the last 8,000 years [182], there is no such clear decrease in permutation entropy. Instead, the corresponding entropy values are found to display a marked low-frequency variability, where even some of the temporary maxima of entropy coincide with known RCC episodes.

The special characteristics of the time interval between 9,000 and 8,000 years BP are supported by complementary analyses of Ferri et al. [183] using Tsallis’ q-triplet [175], who found that the three characteristic scaling exponents from the triplet change their relative ordering during this time interval. The aforementioned results once more underline the relevance of entropic characteristics (especially in the non-extensive framework) as indicators of approaching critical transitions, which is additionally supported by the results of another study on some Antarctic ice core records [184].

The previous results highlight the great potentials of nonlinear methods of time series analysis for the investigation of paleoclimate variability, which are not restricted to statistical mechanics and information-theoretic approaches, but have also been demonstrated for phase space-based characteristics, such as different quantifiers based on recurrence plots [85,185,186,187,188,189]. As an important cautionary note, we need to mention the presence of severe methodological problems, due to the non-uniform sampling and time-scale uncertainty of paleoclimate proxy data [190,191,192], which often do not permit a straightforward use of existing (linear and nonlinear) methods of time series analysis, but call for the development of specifically tailored estimators. Notable exceptions are tree rings, annually laminated (varved) sediments and other layer-forming geological archives (e.g., some ice cores, corals, etc.), where a reliable age-depth model can eventually be inferred based on layer counting.

3.3.3. Hydro-Meteorology and Land-Atmosphere Exchanges

Another climate-related recent field of application of information-theoretic concepts is the investigation of cross-scale processes in biogeochemistry and ecohydrology, i.e., the investigation of interdependencies between atmospheric dynamics and hydrology, on the one hand, and the biosphere (especially vegetation), on the other hand. In general, the biosphere is considered an important part of the Earth system, which exhibits mutual interactions with other components on, or even between, various spatial and temporal scales. The identification of such interdependencies, along with the associated relevant scales, is a challenging task for which bi- and multi-variate characteristics, based on information theory, appear particularly well suited.

Realizing the latter fact, the application of methods, such as mutual information or transfer entropy, has recently gained increasing importance. Brunsell and co-workers [193,194] applied entropy concepts for quantifying the information gained by surface vegetation as a function of the time-scale covered by the input precipitation field, revealing a clear relationship between information content and data resolution. In addition, they were able to characterize the sensitivity of evapotranspiration estimates to the spatial scales of vegetation and soil moisture dynamics obtained from satellite measurements that have been used for deriving these estimates. Stoy et al. [195] utilized the Shannon entropy, as well as relative entropy (Kullback-Leibler divergence) for assessing the amount of information contributed by the aggregation of scales in remote sensing data. Based on the inferred loss of information, they were able to define an “optimum” pixel size for working with such data. Subsequently, Brunsell [196] applied Shannon entropy, relative entropy and mutual information content for investigating the spatial variation in the temporal scaling of daily precipitation data, revealing the existence of distinct scaling regions in precipitation not discriminated before by traditional analysis techniques.

While in the aforementioned examples, information content and information transfer have been considered for different datasets with the same spatial and/or temporal resolution, a more general framework was put forward recently [197]. In this context, Shannon entropy and relative entropy have been used for identifying the dominant spatial scales in land-atmosphere H${}_2$O fluxes between surface and atmosphere obtained from remote sensing data and investigating the interactions of corresponding processes for different stages of the phenological cycle (i.e., the developmental stages of vegetation). The obtained results provided a deep and unique look into the multi-scale spatial structure of evapotranspiration, one of the most important climate variables determining vegetation growth. In a similar spirit, Cochran et al. [198] used relative entropy together with a wavelet-based multi-resolution analysis of ground-based eddy covariance measurements from so-called flux towers for studying the CO${}_2$ exchange between the atmospheric boundary layer and the free troposphere. The aforementioned studies exemplify the relevance of entropic characteristics for studying various aspects of climate-vegetation feedback on local to global scales, thereby providing previously unrecognized information on the terrestrial carbon cycle and, hence, the Earth’s biogeochemistry.

As a complementary approach to studying interdependencies between the biosphere and climate system, Ruddell and Kumar put forward the idea of information process networks in ecohydrology, characterizing different types of coupling corresponding to synchronization, feedback and forcings [199,200]. Here, the quantification of information flows between simultaneously measured variables at various time lags by means of transfer entropy provides weighted network representations of the corresponding interaction structures, which can be further quantified using proper statistical measures. Kumar and Ruddell [201] successfully applied their formalism for inferring feedback loops from flux tower measurements in seven ecosystems in different climatic regions. As a particularly relevant result, they provided evidence that feedbacks tend to maximize information production, since variables participating in feedback loops exhibit “moderated” variability. The amount of information production was found to increase with gross ecosystem productivity and revealed a clear dependence of the corresponding processes on the specific ecological and climate setting, seasonal patterns, etc. Specifically, drier and colder ecosystems were found to respond more intensely and immediately to small changes in water and energy supply, which is in accordance with the general experience in ecology.

3.3.4. Interdependencies and Causality between Atmospheric Variability Patterns

One important recent field of application of mutual information is the construction of climate networks from spatially distributed climate records. Here, individual time series corresponding to data from distinct measurement stations or grid points in a climate model are considered as variables, the strength of statistical associations between which can then be characterized by a variety of linear, as well as nonlinear measures. Selecting the strongest and/or statistically most significant pairwise associations provides a discrete set of spatial interdependencies, which can be further analyzed by network-theoretic concepts. Besides classical linear correlations and concepts referring to different notions of synchronization [202,203,204,205] (which have proven their potentials for bivariate analyses of climate data [206,207] before applications to network construction), application of mutual information based on symbolic dynamics [208,209] or order patterns [210,211] has attracted considerable interest and leads to a multi-scale description of the backbone of relevant spatial interdependencies in the climate system. Recently, the corresponding framework has been successfully generalized to studying directed interrelationships by making use of conditional mutual information [60]. We need to mention that the climate network approach based on information-theoretic characteristics is a subject of ongoing methodological discussions, referring to effects, such as the intrinsic dynamic complexity (which can be characterized itself by concepts, such as Gaussian process entropy rates [169]) or data artifacts resulting in spurious nonlinearities [212].

In addition to climate network construction and analysis, there have been first successful case studies applying information-theoretic concepts for coupling and causality analysis associated with some more specific research questions from climatology. Pompe and Runge [65] used a technique based on permutation entropy for studying the mutual interdependency between sea-surface temperature (SST) anomalies at two specific sites in the North Atlantic, motivated by the path of the North Atlantic Current. Their analysis revealed two distinct, relevant time-scales of oceanic coupling in addition to an instantaneous interrelationship; the latter is most likely mediated by atmospheric processes, namely the North Atlantic Oscillation (NAO) dipole pattern. These non-zero coupling delays were not as clearly detectable using the traditional mutual information or linear cross-correlation. In a similar study, graphical models extending the established transfer entropy framework have proven their skills in obtaining climatologically meaningful coupling patterns from the station data of daily sea-level pressure over Central and Eastern Europe during winter [59]. Another recent study making use of the same formalism investigated the vertical heat transport between the sea surface and upper troposphere based on monthly air temperature anomalies at different altitudes, having distinct implications for our mechanistic picture of the functioning of the atmospheric Walker circulation [64].

As an illustrative case study, we recall here the example originally provided by [213], who demonstrated the estimation algorithm of MIT using the linear partial correlation and its application to studying the Pacific Walker circulation based on reanalysis data with monthly resolution. In the following, we re-examine this example using the nonlinear conditional mutual information estimated with the k-nearest neighbor approach described in Section 2.4.1.

The basic mechanism of the Walker circulation [214,215,216,217] suggests that heating of the western equatorial Pacific, which is accompanied by low surface air pressure, promotes upwelling of waters in the eastern Pacific. This leads to easterly trade winds along the surface of the equatorial Pacific. The moist air then rises above the western Pacific, and the circulation is closed by the dry air masses descending along the central and eastern Pacific.

To test whether the feedback between the eastern Pacific (EPAC, represented here by the average SST over ${80}^{\mathrm{o}}$–${100}^{\mathrm{o}}$ W, $5^{\mathrm{o}}$S–$5^{\mathrm{o}}$N) and western Pacific (WPAC, average sea-level pressure anomalies over the region ${130}^{\mathrm{o}}$–${150}^{\mathrm{o}}$ E, $5^{\mathrm{o}}$S–$5^{\mathrm{o}}$N) was mediated via the surface of the central equatorial Pacific (CPAC, average SST over ${150}^{\mathrm{o}}$–${120}^{\mathrm{o}}$ W, $5^{\mathrm{o}}$S–$5^{\mathrm{o}}$N), we study the time series graph and MIT of the three-variable process (EPAC, CPAC, WPAC).

Figure 4 shows the lagged MI and the significant values of MIT after the time series graph was estimated using a significance level of $\alpha =95$%. The MI lag function is significant between all pairs of variables, making a precise assessment of the causal structure impossible. For example, we find broad peaks between EPAC → CPAC and CPAC → WPAC, while the MIT values are significant only at lags 1 and 0 (indicating a contemporaneous dependency). Very interestingly, while MI features a broad peak in EPAC → WPAC, the time series graph shows no link (even for a lower significance level of 90%). Indeed, this link is obviously mediated via the surface of the central equatorial Pacific. In turn, the MIT of the link WPAC → EPAC at lag 1 is still significant, which demonstrates that the link back takes a different path, not via the surface of the equatorial central Pacific, which is in line with the climatological theory that the air masses return into the high troposphere.

In summary, the time series graph approach confirms the Walker circulation pattern, and the MIT further gives an estimate of the strengths of the underlying causal dependencies. We find generally very strong auto-MIT values, especially for the temperature variables, CPAC and EPAC, which indicates their strong inertia due to the large specific heat capacity of the ocean. Furthermore, the MIT values of the couplings EPAC → CPAC and CPAC → WPAC are comparable, whereas the MI peak of the former is significantly larger. This effect can be explained by the strong auto-dependencies within EPAC and CPAC, which increase MI, but are `filtered out’ in MIT. Another interesting point is that the link WPAC → CPAC that describes the descending mechanism is recovered here, while it was missing in the linear analysis in [213]. This finding suggests a nonlinear character of this dependency. However, in the simple Walker circulation picture, one would not expect the link CPAC → EPAC, which could be explained by the fact that the above described Walker circulation pattern is different during El Niño events, where the region of atmospheric updrafts shifts towards the central Pacific and also broadens markedly. Here, we used the whole time sample to test the hypothesis that the “average” influence is mediated via the central Pacific; a time-resolved analysis will be the subject of future research.

Figure 4. (Left) Lag functions of mutual information (gray) and significant momentary information transfer (MIT) values (black) for all pairs of variables (western Pacific (WPAC), central equatorial Pacific (CPAC), eastern Pacific (EPAC)). Note that the (conditional) mutual information ((C)MI) values have been normalized to the scale of (partial) correlations according to Equation (32). The gray horizontal line denotes the significance threshold for MI only. (Right) Process graph as in Figure 1b, where the node color represents the auto-MIT strength at lag 1 and the link color the cross-MIT strength, with the label denoting the lag. Straight lines mark contemporaneous links. The most important finding is the vanishing link EPAC → WPAC, which shows that the influence of the east on the west Pacific is mediated via the surface of the equatorial central Pacific. In turn, the link WPAC → EPAC is conserved, implying that this influence was not mediated via this region.

Figure 4. (Left) Lag functions of mutual information (gray) and significant momentary information transfer (MIT) values (black) for all pairs of variables (western Pacific (WPAC), central equatorial Pacific (CPAC), eastern Pacific (EPAC)). Note that the (conditional) mutual information ((C)MI) values have been normalized to the scale of (partial) correlations according to Equation (32). The gray horizontal line denotes the significance threshold for MI only. (Right) Process graph as in Figure 1b, where the node color represents the auto-MIT strength at lag 1 and the link color the cross-MIT strength, with the label denoting the lag. Straight lines mark contemporaneous links. The most important finding is the vanishing link EPAC → WPAC, which shows that the influence of the east on the west Pacific is mediated via the surface of the equatorial central Pacific. In turn, the link WPAC → EPAC is conserved, implying that this influence was not mediated via this region.

4. Discussion

The examples discussed in Section 3 demonstrate that entropic characteristics provide versatile means to detect not only spatial patterns, but also temporal changes of dynamical complexity, which can be particularly relevant for anticipating approaching qualitative changes in the dynamical regime of the system under study. From a statistical point of view, such critical phenomena (extreme events) should be characterized (among other features, [218,219,220,221,222]) by the following crucial symptoms: (i) a high degree of organization or high information content; (ii) strong persistency, indicating the presence of a positive feedback mechanism in the underlying phenomenon that drives the systems out of equilibrium; (iii) the existence of a clear preferred direction of the underlying activities; and (iv) the absence of any footprint of a second-order transition in equilibrium.

In this review, we have focused on the first of the aforementioned symptoms in terms of entropies of different kinds: on the one hand, non-extensive Tsallis entropy using symbolic dynamics techniques and, on the other hand, approximate entropy, sample entropy and fuzzy entropy. We note that it has been established that physical systems exhibiting critical phenomena, such as magnetic storms, earthquakes or climate shifts, which are characterized by long-range interactions or long-term memories, or are of a multi-fractal nature, are best described by the generalized statistical mechanics formalism proposed by Tsallis [1,25]. Therefore, although the various introduced measures of dynamical organization or information content can reveal crucial features associated with the launch of an extreme event, namely, low complexity, high degree of organization or high information content, the use of Tsallis entropy seems to be a particularly appropriate measure from a physical point of view.

We clarify that the use of Tsallis entropy leads to double information. The entropic index, q, describes the deviation of Tsallis entropy from the standard Boltzmann-Gibbs entropy, namely, the appearance of long-range interactions, long-term memory and/or multi-fractal behavior. However, the index, q, itself is not a measure of the complexity or organization of the system. In turn, the time variation of the Tsallis entropy, $S_q$, for a given q quantifies the dynamic changes of the complexity or organization of the system. Lower $S_q$ values characterize the portions of the signal with lower complexity or higher organization.

In contrast to the non-extensive statistical mechanics framework, $ApEn$, $SampEn$ and $FuzzyEn$ are especially useful for the analysis of short time series contaminated by noise. For example, since the recorded preseismic EM emissions time series are usually strongly contaminated by such background noise, $ApEn$, $SampEn$ and $FuzzyEn$ are particularly appropriate for their analysis in the form of distinct (short) time windows. We expect similar benefits for the analysis of certain time series in a climatic context, as well.

Mounting empirical evidence has been supporting the possibility that a number of systems arising in disciplines as diverse as physics, biology, engineering and economics may have certain quantitative features that are intriguingly similar. These properties can be conveniently grouped under the headings of scale invariance and universality [223,224,225,226,227,228]. For instance, de Arcangelis et al. [228] have shown that the stochastic processes underlying apparently different phenomena, such as solar flares and earthquakes, have universal properties.

Balasis et al. [229] have recently established the principle of universality in solar flares, magnetic storms and earthquake dynamics. The aforementioned similarity is quantitatively supported by the observation of power-laws in the distributions of solar flare and magnetic storm energy related to a non-extensive Tsallis formalism that gives the Gutenberg-Richter law for the earthquake magnitude distribution as a special case. In a later study [230], Balasis et al. extended their previous work [229] to include preseismic EM emissions. They, thus, provide evidence for universal behavior in the Earth’s magnetosphere, solar corona and in the coupled Earth’s ionosphere, atmosphere and lithosphere system, where magnetic storms, solar flares and preseismic kHz EM anomalies occur, respectively.

Another interesting observation that witnesses universality among different systems is the emergence of discrete scale invariance. Discrete scale invariance has been documented so far in ruptures, seismicity, financial systems (for a review see [147]) and, also, quite recently, for Dst time series that represent magnetospheric activity [231].

The evidence for universal statistical behavior suggests the possibility of a common approach to the forecasting of space weather and earthquakes. In any case, the transfer of ideas and methods of seismic forecasting to the prediction of solar flares and magnetic storms could improve space weather forecasting. Similar concepts could be applied to the treatment of space weather and earthquakes along with extreme events arising from meteorological hazards. For instance, weather forecasting techniques could be applied and tested to predict the evolution of space weather or earthquakes and vice versa.

Besides the quantitative characterization of dynamical complexity, various recent examples summarized in this review point to a considerable importance of entropy-based characteristics for identifying and quantifying (possibly nonlinear) interdependencies between different (geophysical or other) variables, variability at different scales, etc. Besides obtaining estimates of the strengths of bivariate interrelationships between two variables and assessing the corresponding significance by means of suitable statistical tests, unveiling and characterizing asymmetric interdependencies is a crucial point for developing a better “mechanistic” understanding of the governing processes, which can exhibit symmetric interactions between different quantities, but also distinct “driver-response” relationships. The latter are of particular relevance when thinking about causality issues, e.g., trying to understand how certain dynamical patterns in one variable affect the behavior of another one. As a specific geoscientific example for this kind of research question, we mention the not yet systematically understood effect of extreme climatic conditions on vegetation dynamics in different types of ecosystems [232].

Among other complications for appropriate analysis and modeling of geoscientific data, cross-scale variability is a common problem in Earth sciences [3] when it comes to the detection, quantification and attribution of physical processes interlinking such “components”. Notably, such dynamical interdependencies foster complex nonlinear dynamics in terms of feedbacks or various kinds of synchronization. Beyond certain climatic examples mentioned in this review, we envision the potential of modern approaches related to transfer entropies, momentary information transfer and graphical models for tackling contemporary research questions in various fields of Earth system science. Among others, these information-theoretic functionals have great potentials for applications to problems, such as the time-dependent coupling between solar wind and the Earth’s magnetosphere and the geospace (storm–substorm relationship and interdependency [233]).

Finally, we envision the systematic application of information-theoretic approaches for the explicit study of spatio-temporal dynamics. Specifically, combinations with complex network analysis, where concepts like ordinary or conditional mutual information have already been successfully utilized for inferring the most relevant couplings among sets of variables in a climate context [60,208,209,210,211], have great potentials for unveiling and quantifying possibly unknown dynamical interrelationships, responses and feedbacks. Besides their application in the context of climate networks, the latter concepts have also been successfully utilized for studying the spatio-temporal backbone of seismicity [234], suggesting a much wider range of potential research questions to be addressed with similar approaches.

5. Conclusion

Entropy characteristics based on statistical mechanics and information-theoretic considerations exhibit a large degree of conceptual similarity. For example, the classical Boltzmann-Gibbs entropy is known to have its analog in the Shannon entropy originally introduced in a telecommunications science context. This duality, which is further elaborated on in terms of the recently developed non-extensive statistical mechanics framework, provides the basic foundation and rationale for the systematic application of information-theoretic notions of entropies and associated complexity measures and, also, complementary phase space-based entropy characteristics to investigating dynamical complexity in “abstract” (mathematical), as well as observational (physical) systems. More specifically, the idea of quantifying the disorder of a configuration in a microscopic thermodynamic framework is transferred to “configurations” representing observable dynamics, i.e., patterns arising due to the evolution of the system under study, which allows quantifying dynamical complexity by various types of entropies and associated complexity measures.

In this review, we have identified four main fields of application of such characteristics in the Earth sciences:

• the measurement of complexity of stationary records (e.g., for identifying spatial patterns of complexity in climate data);
• the identification of dynamical transitions and their preparatory phases from non-stationary time series (i.e., critical phenomena associated with approaching “singular” extreme events, like earthquakes, magnetic storms or climatic regime shifts known from paleoclimatology, but expected to be possible in the future climate, associated with the presence of climatic tipping points [235,236]);
• the characterization of complexity and information transfer between variables, subsystems or different spatial and/or temporal scales (e.g., couplings between solar wind and the magnetosphere or the atmosphere and vegetation);
• the identification of directed interdependencies between variables related to causal relationships between certain geoscientific processes, which are necessary for an improved process-based understanding of the coupling between different variables or even systems, a necessary prerequisite for the development of appropriate numerical simulation models.

For tackling the aforementioned, as well as conceptually similar research questions, different classes of entropic characteristics have been discussed, illustrating their respective strengths for different kinds of data and problems. Notably, methodological potentials and limitations of the considered approaches are closely associated with the specific properties of the available time series:

• Symbolic dynamics approaches (e.g., block entropies, permutation entropy, Tsallis’ non-extensive entropy) typically require a considerable amount of data for their accurate estimation, which is mainly due to the systematic loss of information detail in the discretization process. Consequently, we suggest that such approaches are particularly well suited for studying long time series of approximately stationary processes (e.g., climate or hydrological time series) and high-resolution data from non-stationary and/or non-equilibrium systems (e.g., electromagnetic recordings associated with seismic activity).
• Tsallis’ non-extensive entropy is specifically tailored for describing non-equilibrium phenomena associated with changes in the dynamical complexity of recorded fluctuations, as arising in the context of critical phenomena, certain extreme events or dynamical regime shifts.
• Distance-based entropies ($ApEn$, $SampEn$ and $FuzzyEn$) provide a higher degree of robustness to short and possibly noisy data than other entropy characteristics based on any kind of symbolic discretization. This suggests their specific usefulness for studying changes of dynamical complexity across dynamical transitions in a sliding windows framework.
• Directional bivariate measures and graphical models allow for a statistical evaluation of causality between variables. In this spirit, this class of approaches provides a versatile and widely applicable toolbox, where the graphical model idea gives rise to a multitude of bivariate characteristics (including conditional mutual information and transfer entropy as special cases), but may exhibit a considerable degree of algorithmic complexity in practical estimation.

We emphasize that the research questions and data properties mentioned above are relatively generic and can be found in various fields of the Earth sciences, but also other scientific disciplines. Learning from successful applications of entropy concepts originating in statistical mechanics or information theory in one field can therefore provide important insights or conceptual ideas for other areas (in the Earth sciences or beyond) or even stimulate new research questions and approaches. Some corresponding examples for prospective problems to be studied by such means have been mentioned in Section 4. In this spirit, this review is intended to provide such a stimulation, without aiming to provide a complete summary of all available approaches or recent applications.