Global climate warming poses a significant challenge to humanity; it is associated with, e.g., rising sea level and declining Arctic sea ice. Increasing extreme events are also considered to be a result of climate warming [1, 2] and they may have widespread and diverse effects on health, agriculture, economics and political conflicts [3][4][5][6]. Still, the detection and quantification of climate change, both in observations and climate models, constitute a main focus of the scientific community. Here we develop a new approach based on network and percolation frameworks [7] to study the impacts of climate changes in the past decades using historical models and reanalysis records, and we analyze the expected upcoming impacts using various future global warming scenarios. More specifically, we classify the globe area into evolving percolation clusters and find a discontinuous phase transition, which indicates a consistent poleward expansion of the largest (tropical) cluster, as well as the weakening of the link's strength. This is found both in the reanalysis data and in the 5 (CMIP5) twenty-first century climate change simulations [8]. The analysis is based on high resolution surface (2 m) air temperature field records. We discuss the underlying mechanism for the observed expansion of the tropical cluster and associate it with changes in atmospheric circulation represented by the weakening and expansion of the Hadley cell. Our framework can also be useful for forecasting the extent of the tropical cluster in order to detect its influence on different areas in response to global warming. The series of powerful Atlantic hurricanes that hammered the Americas, the "Lucifer" heat wave that stifled Europe and the unusually dry June in Australia, all occurring in 2017, have been associated with the increased risks of extreme weather events [9]. Indeed, recent strong evidence supports the claim that this increase in extreme events is directly related to global warming [10]. Although some scientists question the significance of and the role played by human activity in global warming, the majority of the scientific community agrees that the warming is anthropogenic, due to elevated concentrations of heat-trapping (greenhouse) gases, especially triggered by the increased burning of fossil fuels and deforestation [11]. Network theory has demonstrated its potential as a useful tool for exploring the dynamical and structural properties of real-world systems from a wide variety of disciplines in physics, biology, and social science [7,[12][13][14][15][16][17]. Network approaches have been successfully implemented in climate sciences to construct "climate networks," in which the geographical locations are regarded as network nodes, and the level of similarity between the climate records of different grid points represents the network links (strength). Climate networks have been successfully used to analyze, model, and even predict climate phenomena [18][19][20][21]. Percolation theory has found to be an effective tool for understanding the resilience of connected clusters to node breakdowns through topological and structural properties [22][23][24]. The essence of the analysis is the identification of a system's different components and the connectivity between them. Percolation theory was applied to many natural and human-made systems [16,22,[25][26][27]. Here we combine climate network and percolation theory approaches to develop a framework with which to study and quantify the dynamical structure of the global climate system. Our results suggest that an abrupt first order percolation (phase) transition occurs during the evolution of the spatio-temporal climate networks. This evolution indicates the weakening and expansion of the giant component (tropical cluster) and is consistent with reported changes (expansion and weakening) of the tropical (Hadley) circulation. Similar to earlier studies [19,21], we construct a climate network based on the near-surface (monthly mean), high resolution (0.125  ), air temperature of the ERA-Interim reanalysis data [28] [Details are discussed in DATA AND METHODS]. We focus on the surface temperature field since it probably is the most commonly discussed global warming field; other variables, as well as other vertical layers, can be analyzed similarly to the way described below. Our evolving clustering process starts globally with N = 1439 x 2880 isolate nodes (the South and North Pole grid points are eliminated). We then embed the network into a two-dimensional lattice where only nearest neighbor links are considered. The links are sorted in decreasing order of strength and then added one by one according to decreasing strength W [Eq. ( 4) below]; i.e., we first choose the link with the highest weight, then the second strongest link and so on. More specifically, the nodes that are more similar (based on their temperature variations) are connected first. Existing clusters grow when a new link connects one cluster to another cluster (as small as a single node). We find that the climate network undergoes an abrupt and statistically significant phase transition, i.e., exhibiting a significant discontinuity in the order parameter G 1 , the relative size of the largest cluster. Our results indicate that links with higher similarities tend to localize into a few large components (clusters of nodes), in the tropics and in the higher latitude regions (poles) of the Northern and Southern hemispheres. [We show the dynamical evolution of the climate networks in the supplementary animation [29].] Fig. 1(a) shows the climate network component (cluster) structure in the globe map at the percolation threshold (just before the largest jump that is indicated by the orange arrow in Fig. 1(b)). We find that the network, just before this jump, is characterized by three major communities; the largest one is located in the tropical region (indicated by red color); the second and third largest are located in the high latitudes of the southern (indicated by blue color) and northern (indicated by green color) hemispheres. In the next step, a critical bond will connect and merge the tropical cluster with the Southern Hemisphere cluster, resulting in a giant component. Fig. 1(b) depicts the relative size of the largest cluster (the order parameter), G 1 , as a function of the bond/link occupied probability r in the evolution of the climate network. We find that G 1 exhibits an abrupt jump at the percolation threshold r c  0.53.
The probability density function (PDF) of the weight W i,j of links is shown in Fig. 1(d). To study the significance of these results, we constructed the PDF of (temporally or spatially) reshuffled temperature records. In this way, either the memory within each record or the crosscorrelations between the records are destroyed. We repeated the shuffling procedure 100 times and calculated their corresponding PDFs. The results, shown in Fig. 1(b) and (d), suggest that in contrast to the randomly shuffled case, where the transition is continuous, the giant component in our real climate network shows an abrupt change. The vertical line in Fig. 1(d) indicates the strength of the critical link W c at the percolation threshold r c . In addition, we find that W c is also larger than the 95% confidence level of the randomly shuffled data. During the growth of the climate network, we find, in Fig. 1(b), that there are other smaller jumps. For example, the second largest jump occurs at r  0.8, which is caused by the merging of the global cluster with the Northern Hemisphere high latitude cluster shown in blue in Fig. 1(c). S1 in [29]], we find that the tropical cluster is expanding poleward. This weakening and poleward expansion of the tropical component may be associated with global warming, as discussed below. An alternative definition of the cluster intensity is W , the average weight of links in the tropical cluster at the percolation threshold; this yields similar results (see Fig. S2 in [29]). Next, we investigate the response of the tropical component to global warming using twentyfirst century global warming experiments CMIP5 [8]. We used the Representative Concentration S1]. In Fig. 2(c) and (d S1. To study a possible climatological origin of the aforementioned results-weakening and expansion of the tropical cluster-we consider now the atmospheric circulation, especially, the Hadley cell (HC), which plays a pivotal role in the earth's climate by transporting energy and heat pole-ward. We conjecture that the tropical component of the climate network can be linked to the HC, since both of them exhibit consistent weakening and poleward expansion under global warming. In other words, the poleward expansion and weakening of the HC under global warming scenarios are linked to the expansion and weakening of the network-based tropical cluster. Below we present results that support this conjecture. An analysis of satellite observations indicates a poleward expansion of the HC by  2  latitude from 1979-2005 [30]. The main mechanisms for changes in the HC and its relation to global warming have yet to be elucidated [31]. A possible mechanism for the changes in the HC's and its relation to global warming has been proposed in [32]. Also, a possible mechanism for the changes in the HC's strength and its relation to global warming was developed in [33]. Both abovementioned observations and theories suggest a weakening and poleward expansion of the HC under global warming. To find the relationship between the evolution of the climate network and the HC, we further calculate the stream function (Eq. 8, see [34]) as a function of time, and from this, we evaluate the changes in the meridional width  H and strength  of the HC. We show the results in Fig. S3 for the ERA Interim reanalysis data and for a representative CMIP5 model under the RCP8.5 scenario. Similarly, we define   H and   as the change rates (the slope of the trend line) for the width and intensity of the HC. Fig. 3(d), (e) and (f) show the corresponding results, where, as previously reported [32,33,35], most of the CMIP5 models exhibit weakening and expansion of the HC [i.e., S4 depicts the results of a few CMIP5 models and indicates that there is a significant positive correlation between the tropical cluster intensity, W c , and the intensity of the HC, . Similar correlations have also been observed between the width of the tropical cluster, G c , and the width of the HC,  H , for each individual model. A similar significant positive relationship has been found across all models by comparing  W and   , see Fig. 3(g), (h) and (i). Fig. S5 shows the positive relationship between the width of the tropical cluster and the HC. These results indicate that the tropical component of the climate network is correlated with the HC. It should be noted that in contrast to model and theory [33,35], an increasing trend is found in the strength of HC in ERA Interim [Fig. S3(A)]. This is also confirmed by other NCEP-NCAR and ERA-40 reanalysis datasets [36]. A tentative hypothesis trying to resolve this contradiction was recently suggested [35]. It is assumed that the increasing trends are artifacts related to the fact that the tropical lapse rate in the radiosonde data is increasing rather than staying close to a moist adiabat. However, researchers still have not found a way to deal with these artifacts. Our networkpercolation approach seems to be more robust as it yields consistent decreasing trend strength for the different datasets, including models and reanalysis, as well as different time periods. This is most likely due to the fact that the network-percolation approach is based on surface temperature data only. Fig. S6 shows our network results tested on the ERA-40 reanalysis-we find the same tendency of a weakening and an expansion of the tropical component under global warming. The Intergovernmental Panel on Climate Change (IPCC) forecasts a temperature rise of 1.4 to 5.6 deg C over this century (by 2100) [1]. Here, we associate the HC circulation with the tropical component found by the network-percolation approach. The locations of the subtropical dry zones and the major tropical/subtropical deserts are associated with the subsiding branches of the HC [35,37]. Therefore, the poleward expansion of the HC may result in a drier future in some tropical/subtropical regions [32]. Another potential effect is the poleward migration of the location of tropical cyclone maximum intensity [38].
There is thus a great interest in clarifying/predicting the influenced climatological areas in response to global warming. To identify the climate change response, for simplicity but without loss of generality, we compare, in Fig. 4, the topology of the tropical component for the first and last twenty years of the twenty-first century, i.e., (2080-2100) vs. . Counting the number of models that simulate these changes , one can see that the overall pattern of stable node change is robust across most models; the patterns of adding nodes or removing nodes are shown in Fig. 4(b) and (c). We find that some regions, for example, northern India, southern Africa and western Australia have a higher probability to be influenced by the tropical component (or HC), whereas the impacts in other regions, e.g., the Northeast Pacific, will become weaker in the future. For ease of comparison, the datasets are interpolated into a 1 It has been pointed out that a random network or lattice system always undergoes a continuous percolation phase transition and shows standard scaling features during a random process [22,39]. The question whether percolation transitions could be discontinuous has attracted much attention recently in the context of interdependent networks [40] and the so-called explosive percolation models [41][42][43]. Interestingly, the dynamic evolution of our In this study, we employ the monthly 2 m near surface air temperature and 37 pressure level meridional wind velocity V of ERA-Interim [28] and ERA-40 [44] reanalysis datasets. The reso-lution is 0.125 For the climate projection, we used the RCP8.5, RCP4.5 and Historical scenarios of the CMIP5. RCP8.5 is the upper bound of the RCPs and does not include any specific climate mitigation target. The RCP8.5 assumes that greenhouse gas emissions and concentrations increase considerably over time, leading to a radiative forcing that will stabilize at about 8.5 W m 2 at the end of the twentyfirst century. All analyzed data are monthly mean surface air (2 m) temperatures and the meridional component of wind, V . The details of all 31 models we used are summarized in Table S1 [29]. The data can be downloaded in For each node i (i.e., longitude-latitude grid point), given a temperature record Ti (t), where t stands for the month, we consider the month-to-month temperature difference We take the difference since we focus on climate change and how a change in one node affects other nodes; we repeated the analysis using the original time series and obtained similar clusters. To obtain the strengths of the links between each pair of nodes i and j, we define the timedelayed cross-correlation function as, and where  is the time lag between 0 and 24 months; for the networks that span a short time of 5 years, we set   [0, 12]. We define the strength of the link as where "max", "mean" and "std" are the maximum, mean and standard deviations of the crosscorrelation [21]. Based on classical graph theory, a component is a subset of network nodes such that there exists at least one path from each node in the subset to another [7,14]. We denote S m as a series of subnetworks; specifically, S 1 indicates the largest cluster, S 2 indicates the second largest cluster, and so forth. In this study, due to the earths spherical shape, the largest component in the climate networks is defined as where  i is the latitude of node i. Since our network is finite, we use the following procedure to determine the percolation threshold. We first calculate, during the growth process, the largest change of G 1 , where M is the number of links (ordered by decreasing weight value) we added. The step with the largest jump is regarded as the phase transition point. We where  W and  G denote the increasing or decreasing trend rate, a and b are constants, and t is the time; we use a time interval of five years. The strength of the HC is computed using observed zonal-mean meridional wind in the stream function  [34], where V is the meridional velocity in pressure coordinates, R is the mean radius of the earth, and p is the pressure. The operators and [ ] stand for temporal and zonal averaging, respectively. We compute the  field, assuming  = 0 at the top of the atmosphere and integrating Eq. 8 downward to the surface. Since in the winter, the HC is stronger, we only focus on the winter HC intensity. The analyses are performed on December-January-February (DJF) for the Northern Hemisphere (NH) and June-July-August (JJA) for the Southern Hemisphere (SH) separately. We thank Ori Adam and Maaian Rotstein for helping us to retrieve the CMIP5 data and and Xiaosong Chen for helpful discussions. We acknowledge the Israel-Italian collaborative project NECST, the Israel Science Foundation, ONR, the Japan Science Foundation, BSF-NSF, and DTRA ( HDTRA-1-10-1-0014) for financial support. All authors designed the research, analyzed data, discussed results, and contributed to writing the manuscript. We thank Ori Adam and Maaian Rotstein for helping us to retrieve the CMIP5 data and and Xiaosong Chen for helpful discussions. We acknowledge the Israel-Italian collaborative project NECST, the Israel Science Foundation, ONR, the Japan Science Foundation, BSF-NSF, and DTRA ( HDTRA-1-10-1-0014) for financial support. All authors designed the research, analyzed data, discussed results, and contributed to writing the manuscript. Supplementary information is available in the online version of the paper. The authors declare no competing financial interests Supplementary information is available in the online version of the paper. The authors declare no competing financial interests