New perspectives in the study of the Earth’s magnetic field and climate connection: The use of transfer entropy

Abstract

The debated question on the possible relation between the Earth’s magnetic field and climate has been usually focused on direct correlations between different time series representing both systems. However, the physical mechanism able to potentially explain this connection is still an open issue. Finding hints about how this connection could work would suppose an important advance in the search of an adequate physical mechanism. Here, we propose an innovative information-theoretic tool, i.e. the transfer entropy, as a good candidate for this scope because is able to determine, not simply the possible existence of a connection, but even the direction in which the link is produced. We have applied this new methodology to two real time series, the South Atlantic Anomaly (SAA) area extent at the Earth’s surface (representing the geomagnetic field system) and the Global Sea Level (GSL) rise (for the climate system) for the last 300 years, to measure the possible information flow and sense between them. This connection was previously suggested considering only the long-term trend while now we study this possibility also in shorter scales. The new results seem to support this hypothesis, with more information transferred from the SAA to the GSL time series, with about 90% of confidence level. This result provides new clues on the existence of a link between the geomagnetic field and the Earth’s climate in the past and on the physical mechanism involved because, thanks to the application of the transfer entropy, we have determined that the sense of the connection seems to go from the system that produces geomagnetic field to the climate system. Of course, the connection does not mean that the geomagnetic field is fully responsible for the climate changes, rather that it is an important driving component to the variations of the climate.

About the Authors:

S. A. Campuzano

Roles Formal analysis, Investigation, Methodology, Software, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Istituto Nazionale di Geofisica e Vulcanologia (INGV), Roma, Italy

ORCID logo http://orcid.org/0000-0001-7047-5704

A. De Santis

Roles Conceptualization, Methodology, Supervision, Writing – review & editing

Affiliations Istituto Nazionale di Geofisica e Vulcanologia (INGV), Roma, Italy, Dipartimento di Ingegneria e Geologia, Università degli Studi “G. D’Annunzio”, Chieti, Italy

F. J. Pavón-Carrasco

Roles Supervision, Writing – original draft, Writing – review & editing

Affiliations Dpto. de Física de la Tierra y Astrofísica, Universidad Complutense de Madrid (UCM), Madrid, Spain, Instituto de Geociencias (IGEO) CSIC, UCM, Ciudad Universitaria, Madrid, Spain

M. L. Osete

Roles Writing – review & editing

Affiliations Dpto. de Física de la Tierra y Astrofísica, Universidad Complutense de Madrid (UCM), Madrid, Spain, Instituto de Geociencias (IGEO) CSIC, UCM, Ciudad Universitaria, Madrid, Spain

E. Qamili

Roles Writing – review & editing

Affiliation: Serco SpA, Frascati, Rome, Italy

Introduction

The possible relationship between the Earth’s climate and geomagnetic field has been highly debated in the last fifty years (e.g. [1, 2, 3, 4, 5, 6, 7, 8, 9]) but it is still an open question. The first serious proposals, that quantify this possible link, were given by Wollin et al. [1], who pointed out that low geomagnetic intensities are generally associated with warm climate periods (similar to the current situation), and by Bucha [10], who suggested that drifts of geomagnetic poles could have been responsible for displacements of a large low-pressure region of the Earth’s atmosphere associated with an increase of cyclonic activity and sudden climate changes [11].

Throughout the last few decades, other mechanisms that could explain the geomagnetic field-climate relation have been proposed (e.g. [3, 4, 6, 8]). The most plausible at long-time scale is related to the rate of galactic cosmic rays coming to the Earth’s surface. This flux of galactic cosmic rays is modulated by the intensity of both Sun and the Earth’s magnetic fields that act as a protective shield. High values of the solar (and Earth’s) magnetic field intensity reinforce the shield and then a low density of galactic cosmic rays coming to the Solar System (and in turn to Earth) is expected [12]. Entering the atmosphere, the cosmic rays could play an important role in cloud formation [13, 14] and, in this way, the geomagnetic field would be involved in climate processes. That is, a decreasing in the geomagnetic field intensity would allow a higher entrance of galactic cosmic rays to the Earth that could enhance the formation of low-lying clouds [15, 16, 17] or increase the global cloud cover leading to tropospheric cooling [3]. This mechanism was invoked to explain the possible relation between the intensity of Earth’s magnetic field and climate on glacial-interglacial timescales, since dipole moment lows (related to geomagnetic excursions) seem to occur shortly before the onset of relatively cold intervals [6, 8]. This suggests a connection between low geomagnetic intensity and climatic cooling. However, such connection could be circumstantial, as pointed out by these authors, since the variations in geomagnetic field intensity may, in fact, be linked to variations in Earth’s orbital parameters [6], which are considered the main climate-controlling factors in the Pleistocene [18]. Dergachev et al. [19] also studied the relation between short-term geomagnetic variability (jerks) and climate change, as well as the accelerated drift of the north magnetic pole and surface temperature variations. They also propose as more probable mechanism, a relation between the entrance of cosmic rays and formation of clouds.

On the other hand, Gallet et al. [4] compared the advance and retreat of the Alpine Glaciers during the last three millennia with increases and decreases of the geomagnetic field intensity in Paris estimated from archeomagnetic data (paleomagnetic data from heated archaeological artefacts). A later work with a more complete paleomagnetic intensity database corroborated a similar connection at European continental scale [20]. The results of these studies suggest a possible link between centennial-scale cooling episodes and enhanced geomagnetic intensity, the opposite to the galactic cosmic rays mechanism [3, 6, 8, 16, 17] but in agreement with the first links established in the 70’s [1, 10, 11].

Other studies point out other possible mechanisms that explain this connection, such as the experimental result of Pazur and Winklhofer [21]. They focus on the effect of the geomagnetic intensity on CO2 solubility in the ocean. They observed that low values of geomagnetic field intensity reduce the CO2 solubility in the ocean, displacing more CO2 to the atmosphere and increasing the temperature.

For shorter time scales, i.e. last 300 years, De Santis et al. [22, 23] observed a similar temporal trend between the growing South Atlantic Anomaly (SAA) area extent on the Earth’s surface and the Global Sea Level (GSL) rise. The SAA is one of the most outstanding features of the geomagnetic field. It is a large geomagnetic anomaly, presently covering a large area over the Western coast of Africa, the South Atlantic Ocean, the major part of South America and the South-eastern Pacific Ocean, which reaches lower values of intensity than expected at those geomagnetic latitudes. Several studies [24, 25, 26, 27, 28] point out that this anomaly is the response on the Earth’s surface of reversed flux patches located at the terrestrial CMB (core-mantle boundary). De Santis et al. [22] proposed three mechanisms to explain this possible link based on the entrance of charged particles from space, the possible reduction of the ozone layer in the upper stratosphere over the South Atlantic region and/or a common internal cause shared between both SAA and GSL time variations.

All these works and physical mechanisms proposed lead to the deduction that the possible link between the Earth’s climate and the geomagnetic field is far from being demonstrated and understood.

In this work, we propose to study, for the first time, the possible causal information link between two previously studied real time series by means of an innovative statistical tool for non-linear dynamic studies that measures the information flux and the sense of this flux: Transfer Entropy (TE) [29]. This measure has been used in other scientific fields for the last decades, for example in the climatic context [30, 31] or in the geomagnetic activity studies [32]. We will apply it on the SAA surface extent and GSL rise for the last 300 years following De Santis et al. [22] but on shorter scales. We choose these two time series because are important in the frame of the natural hazards. The present strong decrease of the main geomagnetic dipole field could eventually indicate a reversal (e.g. [33, 34]). As well, it plays a main role in screening most of the solar and galactic radiation from space, otherwise penetrating in a larger quantity into the atmosphere and causing possible health and environmental damages. In addition, understanding whether the present increasing trend of the GSL is continuing or not in the close future is vital because of the possible increase of new lands coverage by sea.

The present paper is structured as follows: in the first section, we expose the chosen time series to carry out this analysis. Then, we explain the details on the main methodologies applied in this work. Finally, in the discussion and conclusions we summarize the outcomes reached and their possible future implications.

Data

We analyze two time series: a) the SAA area extent at the Earth’s surface given by historical geomagnetic field models (GUFM1 model, [35]; and the later modifications [27, 36]), and b) the GSL reconstruction for the last 300 years [37]. Both time series are detailed below.

The SAA surface extent could be defined, in practice, by the area below a given intensity contour line at the Earth’s surface (here we selected the contour line of 32000 nT following De Santis et al. [22]). The SAA surface extent has been computed from the three mentioned historical geomagnetic field models covering the last 400 years. The difference between these models lies in the method used to estimate the first Gauss coefficient (g10) prior to 1840 AD, due to the lack of instrumental intensity data before that year. Jackson et al. [35] extrapolated linearly the value of this coefficient backwards from 1840 and they assumed a constant rate of temporal evolution of 15 nT/yr, which corresponds to the average time rate of g10 from 1850 to 1990. Gubbins et al. [27] modified the g10 by using the intensity paleomagnetic database [38] for the period from 1590 to 1840 to obtain a more realistic value of this coefficient. More recently, Finlay [36], using the same paleomagnetic database, applied different statistic approaches to fix again the coefficient g10 providing no rate of change for that coefficient from 1590 to 1840. Consequently, the estimations of the SAA surface extent obtained by these models differ slightly for times prior to 1840, but agree for the most recent period (see Fig 1a).

[Figure omitted. See PDF.]

Fig 1. Time series evolution.

Evolution of a) the SAA area extent (within the 32000 nT isoline of the geomagnetic field) on the Earth’s surface in km2 from three global geomagnetic field models [27, 35, 36] and b) GSL rise in mm, for the last 300 years (1700–2000). The lines represent the fits by using penalized cubic splines: (red, green, blue) SAA derived from Jackson et al. [35], Gubbins et al. [27] and Finlay [36], respectively, and (gray) GSL.

https://doi.org/10.1371/journal.pone.0207270.g001

For the Global mean Sea Level (GSL), we use a reconstruction since 1700 based on the longest available tide-gauge records [37] (http://www.psmsl.org/products/reconstructions/jevrejevaetal2008.php), where the effects of vertical land movement induced by the glacial isostatic adjustment of the solid Earth have been removed. Jevrejeva et al. [37] extended the record backwards from 1850 using three of the longest (though discontinuous) tide-gauge records available, being the error of the reconstruction higher in this epoch (Fig 1b).

We have smoothed both the SAA and the GSL series by using penalized cubic splines in order to avoid future mathematical artefacts resulting from the differences in the reconstruction prior and after 1850. For both records, the fitting was carried out using knot points every 5 years from 1700 to 2000 and a spline damping parameter of 10 yr4/km4 and 10 yr4/mm2 for the SAA and GSL time series, respectively. These optimal values were estimated according to the root mean square (rms) error (see Fig A in S1 File).

In general, the Transfer Entropy (TE) is applied on stationary time series [39]. However, as evident from Fig 1, both SAA and GSL series cannot reasonably be assumed as stationary, being both curves almost monotonically increasing. For this reason, we will apply the TE to the anomaly time series after removing the best-fit long-period trend (see Fig 2). In our case, we choose the simplest polynomial function that accounts for the time evolution of the series: a second order polynomial, which seems the best compromise to remove a reasonable trend and not to completely destroy some similar short-period fluctuations in both series. A positive/negative anomaly would mean that the SAA area extent or GSL rise grow more/lesser than expected.

[Figure omitted. See PDF.]

Fig 2. Evolution of the time series anomalies.

Red, green, blue lines correspond to SAA anomalies derived from Jackson et al. [35], Gubbins et al. [27] and Finlay [36], respectively. Grey line represents the GSL anomalies. See text for further details. Both time series have been normalized to zero mean and unit variance.

https://doi.org/10.1371/journal.pone.0207270.g002

Methods

TE is an information theoretic measure introduced by Schreiber [29] as a generalization of the mutual information [40]. While the mutual information contains neither dynamics nor directional information, the TE takes into account the dynamics of information transport between two systems. This allows quantifying both the exchange of information to the predominant sense of this flow.

The foundations of the TE are to be found in the basic works of the theory of information [41]. The Shannon entropy is given by:(1)where i represents the states that the process I can assume and p(i) the probability distribution which they follow. This quantity measures the average amount of information needed to encode a process optimally.

From finite-order Markov processes, Schreiber [29] introduced a measure to quantify information transfer between two different time series, based on appropriately conditioned transition probabilities instead of static probabilities. Assuming that the system under study can be approximated by a stationary Markov process of order k, the transition probabilities describing the evolution of the system are p(in+1|in, …, in−k+1). If two processes I and J are independent, then the generalized Markov property(2)holds, where and l indicates the number of conditioning states for process J.

Schreiber [29] proposed, using the Kullback entropy for conditional probabilities [42, 43], to measure the incorrectness of assuming the generalized Markov property (Eq [2]), i.e. I and J are independent, which results in:(3)denoted as transfer entropy (a schematic representation of the TE can be found in Fig B in S1 File). The TE can be understood as the excess amount of information that must be used to encode the state of a process by erroneously assuming that the actual transition probability distribution function is , instead of .

The TE computation on real time series has some shortcomings and limitations that must be addressed in the best possible way: 1) the choice of the strategy followed to calculate the TE: discretization method and optimal parameters. The results depend on the different parameters used and it is important to check that we find approximately invariant results with different sets of them. 2) The finite sample size of the real time series: it is always necessary check that the number of data is enough to apply the TE. By examining the log posterior probability for the optimal number of bins S used to discretize the time series, it is possible to verify whether one possesses sufficient data, and it is the method used in this work. 3) The interpretation of the TE results: Smirnov [44] pointed out the inability of the TE to differentiate indirect influences from direct influences. In general, the most widely used interpretation of the TE is to consider that, if it exists, this means that there is an information flow or transfer between the two time series analyzed (I, J). James et al. [45] found that transfer-like entropies could both overestimate information flow and underestimate influence. They proposed a new interpretation of the transfer entropy as a measure of the reduction in uncertainty about one time series given another, instead of as information flow or transfer, which is understood as the existence of information that is currently in I is caused solely by J’s past.

There are different strategies to calculate the TE from the analysis of real data. Here, we use the method based on the discretization of the time series, which was explained in detail by Sandoval Jr [46]. This method consists in dividing the data in a number of bins S, by assigning a numeric symbol to each bin from 1 to S. Each symbol corresponds to a range of values of data series, which are replaced by the symbols assigned (from 1 to S).

Obviously, the calculation of TE will depend on the specific partition chosen S. In order to obtain the optimal number of bins S, we consider the approach proposed by Knut [47], where S is given by the maximization of the posterior probability p(S|N, nk). Given a uniform bin-width histogram for a statistical data set of N samples, the posterior probability p(S ∣ N, nk) is given by:(4)where nk is the number of samples in the kth bin, V is the data range length, and Γ is the Gamma function. In optimization problems, it is common to maximize the logarithm of the Eq [4] [47], also because from the behaviour of the logarithm one can study if the chosen time series are sufficiently long to be analyzed with a tool like the TE [48]. For this reason, we maximize the logarithm of the posterior probability to, firstly, determine if the chosen time series are long enough and then, estimate the optimal number of bins.

Once we have checked that the number of data is enough and estimated the optimal number of bins S, we discretize the time series as we explained above, and compute directly the TE from the Eq [3] given by Schreiber [29], with and representing both involved series. The choice of the embedding dimension k and l is a key point in the computation of the TE. If the dimension is too low, the information contained in the past time (or memory) of the series I might be assigned to come from J. In order to avoid this, we must get that the series I is independent from itself with a delay k. Therefore, we base the selection of this parameter on the determination of the mutual information between the time series I and itself with a delay k [49]:(5)being Ik the time series I with delay k. The value of k associated with the first local minimum reported in the Eq [5] is considered the optimal embedding dimension.

For the dimension of the embedding l of the J series, it is usually considered l = 1 or l = k [29, 39]. In a conservative approach we consider l = 1. To calculate the different probabilities of the Eq [3] we simply count the number of times that a symbol or sequence of symbols appears in our time series.

Due to finite size of the time series and the reduced data number, the establishment of a threshold at which the result can be considered significant is essential. In order to establish the statistical significance of our results we calculate the TE with the data points of the J series, which represents the source of the presumed information flow, shuffled randomly [39, 50]. The objective of this procedure is to destroy all potential relations between the two series, I and J, and hence the observed TE should be zero. In finite time series this value rarely is zero due to the finite sample effects, and we obtain the threshold value of TE above which is significant. Practically, we create 1000 surrogate time series of J by using the Iterated Amplitude Adjusted Fourier Transform technique (IAAFT) [51, 52, 53]. This procedure assures that the surrogate time series have the same mean, variance, autocorrelation function and therefore, power spectrum as the original series but destroys the non-linear relations and, therefore, the information actually significant transferred from J to I series. To consider the original TE significant we consider the 5% null hypothesis being the null hypothesis that the transfer entropy between the two original time series is not significant. Whether the 95% of the new TEs values, calculated from surrogate series J, are lesser than the original one, then we consider the original TE significant.

Results and discussion

The analysis of the logarithm of Eq [4] (log posterior) in function on the number of bins provides useful information: a) both time series are long enough to apply the TE and b) the selection of the optimal number of bins S according to the maximum in the log posterior function (see Fig 3a and 3b). The log posterior of SAA anomalies (Fig 3a) increases sharply according to the number of bins considered, reaching a peak (corresponding to the optimal number of bins S = 5) and then decreasing. Respect to the GSL anomalies series (Fig 3b), the log posterior also decreases gradually but the maximum is not so clear. These behaviors indicate a sufficient amount of data to develop this analysis with the TE, but finite sample effects could be important. Due to the lack of an obvious peak in the GSL anomalies series, we establish an agreement between the log posterior curve and the main characteristics of the histogram of the time series. In view of Fig 3d, we consider that with S = 4 we have captured the main information of this series (see also Fig C in S1 File). Finally, in order to avoid a future bias in the computation of the TE, we choose the same number of bins S for both time series i.e., equal to 4 (see Table 1 and Fig C in S1 File) due to larger bin sizes (smaller S) are usually favored in the literature because show the differences more sharply [46].

[Figure omitted. See PDF.]

Fig 3. Evaluation of the length of time series and optimal number of bins.

Log posterior curves in function on the number of bins S: a) for the SAA anomalies computed from Jackson et al. [35] and b) for the GSL anomalies. The subplots c) and d) represent, in orange and cyan respectively, the chosen discretization (S = 4) taking into account the results given in a) and b), as well as the main characteristics of the probability density of both systems (see red and blue bars in c) and d) plots). The error bars indicate the standard deviation of the bin heights.

https://doi.org/10.1371/journal.pone.0207270.g003

[Figure omitted. See PDF.]

Table 1. Optimal parameters.

https://doi.org/10.1371/journal.pone.0207270.t001

As indicated in the methodology, the selection of the embedding dimension k for both series was carried out using the mutual information given by Eq [5]. Results are plotted in Figs Da and Db and contained in the Table Ab in the S1 File. For the GSL anomaly series the optimal dimension was obtained for kGSL = 13, while different values were obtained for the 3 SAA anomalies series (24 for the SAA anomalies series of Jackson et al. [35]; and 26 for the other two series). Nevertheless, since different embedding dimensions can generate TE bias [54], we have fixed the dimension kSAA in 26 for all the 3 SAA anomalies series because a slight over-embedding does not compromise the detection of significant TE [55]. To corroborate the different value of dimension k for GSL and SAA series, we have also calculated the autocorrelation function since the simplest estimate of an optimal k is the first zero of the autocorrelation function [56, 57]. The problem is that these estimates generally yield too large k values for stochastic dynamical systems [58]. In fact, the first minimum reported for the SAA anomalies is given in kSAA = 29 and for the GSL anomalies in kGSL = 17 (Figs Dc and Dd in S1 File). As provided by the mutual information (kSAA = 26 and kGSL = 13), the autocorrelation functions also indicate a lower memory for the GSL series than the three SAA series.

In order to evaluate how the selection of these parameters (S, k) affects the results, we have performed several tests using different sets of them. The results are detailed in the S1 File along with the Tables A and B. In addition, we have performed different tests to study the impact of a) the use of a different detrending approach to define the anomalies (Fig E and Table C in S1 File) and b) the use of an unsmooth GSL time series (Fig F and Table D in S1 File). Detailed information about these tests could be also found in the S1 File. We find that these changes can slightly affect the statistical significance of our results but not the sense of the information flow between the two time series.

For the chosen parameters, the TE results (Eq [3]) are given in Table 2 and Figs 4 and 5. As it can be observed, there is a significant information flow from SAA to GSL anomalies by considering the 5% null hypothesis when the most recent geomagnetic field models given by Gubbins et al. [27] and Finlay [36] are used. Anyway, the significant levels calculated following the IAATF approach are clarifying, with percentages around the 90% in all cases for the TE from SAA to GSL anomalies. This outcome indicates that the SAA anomalies add great predictability to the GSL anomalies by suggesting interactions between the two time series of anomalies at a time scale lower or equal to two consecutive data, i.e. one year. However, more investigations must be carried out about the time delay that needs the influence to propagate between both series (e.g. [32, 30, 59]).

[Figure omitted. See PDF.]

Table 2. Results of transfer entropy analysis.

https://doi.org/10.1371/journal.pone.0207270.t002

[Figure omitted. See PDF.]

Fig 4. Results of transfer entropy analysis.

Transfer entropy by measuring the information flow from SAA to GSL anomalies and from GSL to SAA anomalies, by using the three historical models for the geomagnetic field to compute the SAA surface extent. In brackets, the significant level indicates the percentage of TEs calculated from surrogate series that are lesser than the original TE.

https://doi.org/10.1371/journal.pone.0207270.g004

[Figure omitted. See PDF.]

Fig 5. Statistical significance of the transfer entropy results.

Transfer entropy calculated from surrogate series a), c) and e) of SAA anomalies from Jackson et al. [35] (SAAJ), Gubbins et al. [27] (SAAG) and Finlay [36] (SAAF) respectively and b), d) and f) GSL anomalies. The results show that the statistical significance is higher when the sense of the information goes from SAA to GSL anomalies, also registering greater values of the TE.

https://doi.org/10.1371/journal.pone.0207270.g005

In view of these results, it would be expected that a future SAA anomaly taking into account our selected trend generates a GSL anomaly with a time lag of one year or less. Several physical mechanisms are proposed to explain this possible coupling [22]. The first of them is that an increase of the SAA area facilitates the entrance of charged particles from space. If the SAA area extent grows more than it is expected (positive anomaly), then this entrance is favored. As a result we have a warmer atmosphere, which, in turn, implies a consequent melting of major ice caps (Antarctica and Greenland) that finally would cause a greater increasing of the global sea level (positive anomaly). Recent works (e.g. [31, 60, 61, 62, 63]) have found interesting correlations between solar and galactic cosmic rays periodic variations and climatic (such as temperature and rainfall) variations in the region where the SAA is located. The entrance of galactic cosmic rays at the atmosphere depends both solar and Earth’s magnetic fields, hence these correlations could also be influenced by a factor depending on the low geomagnetic intensity due to the SAA presence in the region and its continuous increasing for the last centuries.

Another mechanism proposed is that a possible reduction of the ozone layer in the upper stratosphere over the South Atlantic region can modify the radiative flux at the top of the atmosphere and hence can cause changes in the weather and climate patterns, including cloud coverage. Solanki et al. [64] propose a similar mechanism to explain relation between solar activity and climate based on the fact that the variations in solar activity during an 11-year cycle are more intense at shorter wavelengths, which include UV radiation. The variations in UV radiation modify the concentrations of ozone and lead to changes in the atmospheric circulation dynamics.

As we can observe, these two mechanisms relate the solar activity, the galactic cosmic rays production and the geomagnetic field with the Earth’s climate, by suggesting that all of them can work together and be needed to completely explain the found outcomes.

Finally, an internal mechanism was presented by which a convective dynamism in the outer core could cause a variation of the magnetic field and an elastic deformation at the Earth’s surface [65].

In the analyzed case study, we have shown that the sense of the information goes from SAA to GSL time series (Fig 5). This would discard any physical mechanism in which the climate controls the geomagnetic field and support the mechanisms caused by the presence of the SAA.

Conclusions

We have applied for the first time a recent statistical tool, transfer entropy, to shed light on the question of a possible link between the Earth’s magnetic field and climate and provide new perspectives in its future analysis. In this work, we have analyzed two real time series with an analogous evolution for the last 300 years, the South Atlantic Anomaly area extent on the Earth’s surface and the Global Sea Level rise. We have analyzed the anomalies of both time series, after removing the long term trend. The results seem to support the existence of an information flow between SAA and GSL anomalies, with larger information transferred from SAA to GSL and a confidence level about 90%. The found connection does not mean that the geomagnetic field is fully responsible of the climate changes, rather that it is an important driving component to the variations of the climate. This result is especially relevant because could help to find a physical mechanism able to explain this connection by discarding those in which the climate controls the geomagnetic field and supporting the mechanisms associated to the geomagnetic field.

Although this work seems to provide a favorable argument to this link, future investigations are needed to completely exploit this issue, for example to check other time series at longer timescales.

Supporting information

[Figure omitted. See PDF.]

S1 File. Detailed description of the supplementary material.

Complementary figures and tables. Results from different tests that support the main outcomes described in the main text and some case studies that help to better understand the TE results.

https://doi.org/10.1371/journal.pone.0207270.s001

(DOCX)

Acknowledgments

S.A.C., F.J.P.C. and M.L.O. are grateful to the Spanish research project CGL2017-87015-P of the Spanish Ministerio de Economía y Competitividad and the FPI grant BES-2012-052991, which has allowed S.A.C. two 3-month stays at INGV in Rome in 2014 and 2015 (EEBB-I-14-09023 and EEBB-I-15-10151). A.D.S. and S.A.C. also thank the ESA-funded Projects TEMPO and LIMADOU, for providing partial financial support to this research, and the INGV for providing computational and service facilities. S.A.C. thanks J. Jarillo for the help and advices in the development of this work. Authors are very grateful to two anonymous reviewers for their careful and detailed reviews, which have helped to improve the quality of the manuscript substantially. All algorithms have been developed in Matlab codec (Matlab 7.11.0, R2010b) and R-project (R 2.12.2) along with the figures. The data used are listed in the references, tables and main manuscript.

Citation: Campuzano SA, De Santis A, Pavón-Carrasco FJ, Osete ML, Qamili E (2018) New perspectives in the study of the Earth’s magnetic field and climate connection: The use of transfer entropy. PLoS ONE 13(11): e0207270. https://doi.org/10.1371/journal.pone.0207270

References

1. Wollin G, Ericson DB, Ryan WBF. Magnetism of the Earth and climatic changes. Earth Planet. Sci. Lett. 1971; 12: 175–183.

2. Wagner G, Livingstone DM, Masarik J, Muscheler R, Beer J. Some results relevant to the discussion of a possible link between cosmic rays and the Earth’s climate. J. Geophys. Res. 2001; 106: 3381–3387.

3. Christl M, Mangini A, Holzkämper S, Spötl C. Evidence for a link between the flux of galactic cosmic rays and Earth’s climate during the past 200,000 years. J. Atmos. Sol.-Terr. Phys. 2004; 66 (3–4): 313–322.

4. Gallet Y, Genevet A, Fluteau F. Does Earth’s magnetic field secular variation control centennial climate change? Earth Planet. Sci. Lett. 2005; 236: 339–347.

5. Courtillot V, Gallet Y, Le Mouël J-L, Fluteau F, Genevey A. Are there connections between the Earth’s magnetic field and climate? Earth Planet. Sci. Lett. 2007; 253: 328–339.

6. Thouveny N, Bourlès DL, Saracco G, Carcaillet JT, Bassinot F. Paleoclimatic context of geomagnetic dipole lows and excursions in the Brunhes, clue for an orbital influence on the geodynamo? Earth Planet. Sci. Lett. 2008; 275 (3–4): 269–284.

7. Knudsen MF, Riisager P. Is there a link between Earth’s magnetic field and low-latitude precipitation? Geology. 2009; 37: 71–74.

8. Kitaba I, Hyodo M, Katoh S, Dettman DL, Sato H. Midlatitude cooling caused by geomagnetic field minimum during polarity reversal. PNAS. 2013; 110 (4): 1215–1220. pmid:23297205

9. Rossi C, Mertz-Kraus R, Osete ML. Paleoclimate variability during the Blake geomagnetic excursion (MIS 5d) deduced from a speleothem record. Quaternary Science Reviews. 2014; 102: 166–180.

10. Bucha V. Variations of the geomagnetic field, the climate and weather. Stud. Geophys. Geod. 1976; 20:149–167.

11. Bucha V. Correlations of the changes in the Earth’s magnetic field—the climate and weather, electromagnetic field of the Earth (A79-2601309-46). Veda-Vydavatelstvo Slovenskej Akademie Vied, Bratislava. 1978; 59–67.

12. Snowball I, Muscheler R. Palaeomagnetic intensity data: an Achilles heel of solar activity reconstructions. Holocene. 2007; 17, 6: 851–859.

13. Duplissy J, Enghoff MB, Aplin KL, Arnold F, Aufmhoff H, Avngaard M et al. Results from the CERN pilot CLOUD experiment. Atmos. Chem. Phys. 2010; 10: 1635–1647.

14. Kirkby J, Curtius J, Almeida J, Dunne E, Duplissy J, Ehrhart S et al. Role of sulphuric acid, ammonia and galactic cosmic rays in atmospheric aerosol nucleation. Nature. 2011; 476: 429–433. pmid:21866156

15. Svensmark H, Friis-Christensen E. Variation of cosmic ray flux and global cloud coverage—a missing link in solar-climate relationships. J. Atm. Solar Terr. Physics. 1997; 59 (11): 1225–1232.

16. Svensmark H. Influence of Cosmic Rays on Earth’s Climate. Phys. Rev. Lett. 1998; 81 (22): 5027–5030.

17. Usoskin IG, Kovaltsov GA. Cosmic rays and climate of the Earth: Possible connection. C. R. Geosci. 2008; 340 (7): 441–450.

18. Hays JD, Imbrie J, Shackelton NJ. Variations in the Earth’s Orbit: Pacemaker of the Ice Ages. Science. 1976; 194 (4270): 1121–1132. pmid:17790893

19. Dergachev V., Vasiliev S., Raspopov O., Jungner H. Impact of the geomagnetic field and solar radiation on climate change. Geomagnetism and Aeronomy. 2012; 52(8): 959–976.

20. Pavón-Carrasco FJ, Osete ML, Torta JM, Gaya-Piqué LR. A Regional Archaeomagnetic Model for the Palaeointensity in Europe for the last 2000 Years and its Implications for Climatic Change. Pure Appl. Geophys. 2008; 165: 1209–1225.

21. Pazur A, Winklhofer M. Magnetic effect on CO2 solubility in seawater: A possible link between geomagnetic field variations and climate. Geophys. Res. Lett. 2008; 35 (16).

22. De Santis A, Qamili E, Spada G, Gasperini P. Geomagnetic South Atlantic Anomaly and global sea level rise: A direct connection? J. Atmos. Sol.-Terr. Phys. 2012; 74: 129–135.

23. De Santis A, Qamili E, Wu L. Toward a possible next geomagnetic transition? Nat. Hazards Earth Syst. Sci. 2013; 13: 3395–3403.

24. Gubbins D. Mechanism for geomagnetic polarity reversals. Nature. 1987; 326: 167–169.

25. Hulot G, Eymin C, Langlais B, Mandea M, Olsen N. Small-scale structure of the geodynamo inferred from Ørsted and Magsat satellite data. Nature; 2002; 416: 620–623. pmid:11948347

26. Olson P, Amit H. Changes in Earth’s dipole. Naturwissenschaften. 2006; 93: 519–542. pmid:16915369

27. Gubbins D, Jones AL, Finlay CC. Fall in Earth’s Magnetic Field is erratic. Science. 2006; 312 (5775): 900–902. pmid:16690863

28. Finlay CC, Jackson A, Gillet N, Olsen N. Core surface magnetic field evolution 2000–2010. Geophys. J. Int. 2012; 189 (2): 761–781.

29. Schreiber T. Measuring Information Transfer. Phys. Rev. Lett. 2000; 85, 2: 461–464. pmid:10991308

30. Das Sharma S., Ramesh D., Bapanayya C., Raju P. Sea surface temperatures in cooler climate stages bear more similarity with atmospheric co2 forcing. Journal of Geophysical Research: Atmospheres (1984–2012). 2012; 117(D13).

31. Bhaskar A., Ramesh D., Vichare G., Koganti T., Gurubaran S. Quantitative assessment of drivers of recent global temperature variability: an information theoretic approach. Climate Dynamics. 2017; 49 11–12: 3877–3886.

32. De Michelis P., Consolini G., Materassi M., Tozzi R. An information theory approach to the storm-substorm relationship. Journal of Geophysical Research: Space Physics. 2011; 116(A8).

33. Constable C, Korte M. Is the Earth’s magnetic field reversing? Earth Planet. Sci. Lett. 2006; 246: 1–16.

34. Tarduno JA, Watkeys MK, Huffman TN, Cottrell RD, Blackman EG, Wendt A et al. Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo. Nat Commun. 2015; 6: 7865. pmid:26218786

35. Jackson A, Jonkers ART, Walker MR. Four centuries of geomagnetic secular variation from historical records. Philos. Trans. R. Soc. Lond. A. 2000; 358 (1768): 957–990.

36. Finlay CC. Historical variation of the geomagnetic axial dipole. Phys. Earth Planet. Interiors. 2008; 170: 1–14.

37. Jevrejeva S, Moore JC, Grinsted A, Woodworth PL. Recent global sea level acceleration started over 200 years ago? Geophys. Res. Lett. 2008; 35: L08715.

38. Korte M, Genevey A, Constable CG, Frank U, Schnepp E. Continuous geomagnetic field models for the past 7 millennia: 1. A new global data compilation. Geochem. Geophys. Geosyst. 2005; 6: Q02H15.

39. Marschinski R, Kantz H. Analysing the information flow between financial time series: An improved estimator for transfer entropy. Eur. Phys. J. B. 2002; 30: 275–281.

40. Shannon CE. A mathematical theory of communication. Bell System Tech. J. 1948; 27: 379–423.

41. Shannon CE, Weaver W. The Mathematical Theory of Information. University of Illinois Press, Urbana, IL. 1949.

42. Kullback S, Leibler RA. On information and sufficiency. Ann. Math. Statist. 1951; 22: 79–86.

43. Kullback S. Information Theory and Statistics. Wiley, NewYork. 1959.

44. Smirnov DA. Spurious causalities with transfer entropy. Phys. Rev. E. 2013; 87(4): 042917.

45. James RG., Barnett N., Crutchfield JP. Information flows? A critique of transfer entropies. Phys. Rev. Lett. 2016; 116(23): 238701. pmid:27341264

46. Sandoval L Jr. Structure of a Global Network of Financial Companies Based on Transfer Entropy. Entropy. 2014; 16(8): 4443–4482.

47. Knut KH. Optimal Data-Based Binning for Histograms. arXiv:physics/0605197v2 [physics.data-an]. 2013.

48. Knut KH, Gotera A, Curry CT, Huyser KA, Wheeler KR., Rossow WB. Revealing Relationships among Relevant Climate Variables with Information Theory. In The Earth-Sun System Technology Conference, NASA, Adelphi, Md. 2005.

49. Dimpfl T, Peter FJ. Using transfer entropy to measure information flows between financial markets. Stud. Nonlinear Dyn. E. 2013; 17(1): 85–102.

50. Sensoy A, Sobaci C, Sensoy S, Alali F. Effective transfer entropy approach to information flow between exchange rates and stock markets. Chaos Soliton. Frac. 2014; 68: 180–185.

51. Theiler J, Eubank S, Longtin A, Galdrikian B, Farmer JD. Testing for nonlinearity in time series: the method of surrogate data. Physica D. 1992; 58: 77–94.

52. Kugiumtzis D. Surrogate data test for nonlinearity including nonmonotonic transforms. Phys. Rev. E. 2000; 62 (1).

53. Schreiber T, Schmitz A. Surrogate time series. Physica D. 2000; 142: 346–382.

54. Kraskov A, Stogbauer H, Grassberger P. Estimating mutual information. Phys Rev E Stat Nonlin Soft Matter Phys. 2004; 69(6 Pt 2): 066138. pmid:15244698

55. Lindner M, Vicente R, Priesemann V, Wibral M. TRENTOOL: A Matlab open source toolbox to analyse information flow in time series data with transfer entropy. BMC Neurosci. 2011; 12 (119): http://www.biomedcentral.com/1471-2202/12/119.

56. Abarbanel HI. Analysis of Observed Chaotic Data. Springer, New York. 1996.

57. Kantz K, Schreiber T.: Nonlinear Time Series Analysis. Cambridge University Press, Cambridge. 1997.

58. Ragwitz M, Kantz H. Markov models from data by simple nonlinear time series predictors in delay embedding spaces. Phys. Rev. E. 2002; 65: 056201. pmid:12059674

59. Wibral M., et al. Measuring information-transfer delays. PloS One. 2013; 8(2): e55809. pmid:23468850

60. Souza Echer MP., Echer E., Nordemann DJR., Rigozo NR., Prestes A. Wavelet analysis of a centennial (1895–1994) southern Brazil rainfall series (Pelotas, 31°46’19”S 52°20’33”W). Clim. Change. 2008; 87: 489–497.

61. Rampelotto PH., Rigozo NR., da Rosa MB., Prestes A., Frigo E., Souza Echer MP., Nordemann DJR. Variability of Rainfall and Temperature (1912–2008) from Santa Maria (29°41’S, 53°48’W) and its Connection with Natural Influences. J. Atmos. Solar-Terr. Phys. 2012; 77: 152–160.

62. Frigo E., Pacca IG., Pereira-Filho AJ., Rampelloto PH., Rigozo NR. Evidence for cosmic ray modulation in temperature records from the South Atlantic Magnetic Anomaly region. Ann. Geophys. 2013; 31: 1833–1841. https://doi.org/10.5194/angeo31-1833-2013.

63. Frigo E., Antonelli F., Silva DSS., Rampelloto PH., Lima PCM., Pacca IIG., Bageston JV. Effects of solar activity and galactic cosmic ray cycles on the modulation of the annual average temperature at two sites in southern Brazil. Ann. Geophys. 2018; 36: 555–564. https://doi.org/10.5194/angeo36-555-2018.

64. Solanki SK., Krivova NA., Haigh JD. Solar Irradiance Variability and Climate. Annu. Rev. Astron. Astr. 2013; 51: 311–351. https://doi.org/10.1146/annurev-astro-082812-141007.

65. Greff-Lefftz M, Pais MA, Le Mouel J-L. Surface gravitational field and topography changes induced by the Earth’s fluid core motions. J. Geodesy. 2004; 78: 386–392.