The frequency and magnitude of weather extreme events have increased significantly during the past few years in response to anthropogenic climate change. However, global statistical characteristics and underlying physical mechanisms are still not fully understood. Here, we adopt a statistical physics and probability theory based method to investigate the nature of extreme weather events, particularly the statistics of the day-to-day air temperature differences. These statistical measurements reveal that the distributions of the magnitudes of the extreme events satisfy a universal Gumbel distribution, while the waiting time of those extreme events is governed by a universal Gamma function. Further finite-size effects analysis indicates robust scaling behaviours. We additionally unveil that the cumulative distribution of logarithmic waiting times between the record events follows an Exponential distribution and that the evolution of this climate system is directional where the underlying dynamics are related to a decelerating release of tension. The universal scaling laws are remarkably stable and unaffected by global warming. Counterintuitively, unlike as expected for record dynamics, we find that the number of quakes of the extreme temperature variability does not decay as one over time but with deviations relevant to large-scale climate extreme events. Our theoretical framework provides a fresh perspective on the linkage of universality, scaling, and climate systems. The findings throw light on the nature of the weather variabilities and could guide us to better forecast extreme events. Natural or social systems far from equilibrium often show complex dynamic behaviours. Nevertheless, increasing empirical evidence supports that the fields as diverse as physics [1, 2], biology [3,4], ecology [5,6], human mobility [7], economics, and financial markets [8,9] may all, at least up to a certain degree of approximation, follow the principles of scale invariance and universality [10]. Scaling functions characterise the critical phenomena of phase transitions in physical systems. The following mounting evidence demonstrates these functions are also powerful tools for understanding and predicting the behaviours of various complex systems [11,12], such as the climate system [13,14] and earthquakes [15][16][17][18]. Understanding the nature and the statistical properties of climate variability across multiple spatial and temporal scales and its relationship with the ongoing climate change is a key goal of climate science [19,20]. Evaluating how climate change will impact the tail of the distribution of climatic variables is crucial to evaluate the impacts of climate change, as extreme events are often conducive to risks to human and environmental welfare [21][22][23][24]. A great deal of attention has been devoted to understanding to what extent it is possible to attribute individual extreme events to man-made climate change [25][26][27]. (WMO) reported that weather-related disasters have increased over the past 50 years, causing more damage, but early warnings save lives [28]. According to WMO, economic losses have increased sevenfold from the 1970s to the 2010s, but, thanks to the improved early warnings and disaster management, the number of deaths decreased almost three-fold over the 50-year period [28]. In particular, a growing body of literature has identified the social impacts of large temperature variability at different time scales in areas like human health [29,30], crop production [31], and economic growth [32]. Conversely, extremes reveal essential aspects of the dynamics of climate because there is a fundamental connection between suitably defined extreme events and dynamical processes of the system generating them [33][34][35][36][37][38]. Additionally, thanks to the presence of couplings acting on different temporal and spatial scales, the climate system, behaving similarly to a complex network, can feature the spread of the effects of extreme events, leading to cascades of risks [39]. This study aims at advancing our understanding of extreme events in the climate system using the viewpoint provided by statistical physics and probability theory. Specifically, we will take advantage of extreme value theory [34,[40][41][42][43], which provides a very powerful and flexible framework for studying the tails of the distribution of random variables. We will focus on near-surface temperature fields and consider different sites covering the globe. Instead of looking, as often done, at daily temperature records, we will study the extreme events of the day-to-day temperature change (both positive and negative, as well as their absolute value). We remark that extreme positive and negative events are, in general, associated with different physical processes: cold and warm weather fronts are fundamentally different at the dynamical level [44]. Looking at this climatic variable, which is the discrete approximation of the time derivative of the temperature field, amounts to studying one component of the local energy budget of near-surface air masses (the missing component being associated with moisture content). We remark that a recent study looked at the spatial statistics of daily temperature variability measured as the intramonthly standard deviation of daily surface temperature. The authors found a relevant impact of CO2 increase on this type of measure of weather variability [45]. As discussed below, our study uses a different measure of daily temperature variability, which seems better suited to representing the challenges of rapid adaptation of human and environmental systems to rapid temperature variations. We study the daily temperature variations using the differences in temperature between two adjacent discrete time periods, that is x i,d is air temperature recorded for day d at site i; D is the total number of days of the recorded data, and the analysis spans 43 years, from 1979 to 2021. N is the total number of sites (73 x 144) covering the entire Earth based on the ECMWF Reanalysis 5th Generation (ERA5) [46] dataset, see Data availability for more details. Furthermore, the absolute values of the changes provides a measure of the strength of weather variability [22]. A clear illustration of these definitions for one selected site is given in Fig. 1, which shows a 2-year time slice of daily temperature reading and day-to-day temperature variations recorded in Madrid (40.25  N, 3.75  W), Spain.
The plot shows that extreme day-to-day temperature differences are more common before the summertime, followed by high-temperature days from June to August (the summertime of the Northern Hemisphere and the sites in the South Hemisphere fol-lows similar seasonal trends, see Fig. S1 in SI). This is related to the well-known stronger weather variability in winter versus summer, associated with the presence of stronger baroclinic conversion processes [44]. To investigate the approach to the limiting behaviour in the statistics of temperature differences defined by extreme value theory, we apply the Peaks Over Threshold methods [43] to the normalised temperature positive and negative differences for all sites separately. We show the dependence of the estimates of the shape parameter and of the modified scale parameter on the threshold value for sites on land and ocean in Fig. S2 in SI. The distribution of the estimates of shape parameters is broader over the oceans, which means the asymptotic behaviours are less accurately established for ocean sites. However, most of the shape parameters of oceans and lands fall within a reasonably small range of [0.1, 0], and the asymptotic regime is established for quantiles as low as (60%, 67.5%) because the estimates of the scale parameters and of the modified scale parameter become stable as we consider higher quantiles [43,47]. This regime is considerably lower than what is usually found when looking at temperature extremes, see, e.g. [48], which is encouraging in terms of applicability of extreme value theory for day-to-day temperature variability also in conditions where data availability is limited. We consider all sites on the entire Earth together. With this consideration, we have lost the structure in space and in the magnitude scale; nonetheless, the change in the process properties with q within the extreme regime will allow us to recover some information. A positive extreme event set E + i (q) is defined when the value CV exceeds a chosen threshold associated with a quantile q. Then, we define the magnitude of the positive extreme temperature difference as the normalised value of the event, that is m + ev = CV ev +  CV ev + , ev  E + i , where  CV ev + is the standard deviation of the extreme part, CV ev + for site i. Similarly, a negative extreme event set E  i (q) is defined when the value CV is smaller than a chosen quantile q. Then, for simplicity, we define the magnitude of the negative extreme temperature difference as the minus normalised value of the event, that is After taking the minus of the values, the m  ev will be positive, which is comparable to the m + ev . Likewise, we define for the absolute values, m abs ev = |CV| ev abs  CV abs ev , ev  E abs i . The waiting time  i  T i of the extreme events in E i is the number of days between two consecutive events [49]. Therefore, we can define four types of waiting times:  + (  ,  abs ) is the waiting time between positive (negative, absolute) extreme temperature difference, i.e. the waiting time in E + (E  , E abs ). The fourth type  + ( + interchangeable) is the number of days between extreme positive and extreme negative change events (or extreme negative and extreme positive), Here, we focus on the distributions of the dimensionless magnitudes m + ev and m  ev for all sites on the Earth, that is, the positive and negative extreme parts of CV. We choose q = 90%, 92.5%, 95% and 97% as threshold quantiles within the extreme regime. The normalised PDF for m + ev and m  ev , p(m + ev ; q) and p(m  ev ; q) are shown in Figs. 2a and c, respectively. The PDF satisfies a single peak distribution and depends on our chosen threshold quantiles q. The distributions of m + ev and m  ev are expected to take the following scaling forms, and where f + and f  are two universal scaling functions, m + ev ( m  ev ) and (q) is the mean value and standard deviation associated with quantile q. We then define two rescaled parameters and  . Combining the proposed scaling forms, Eqs. ( 1) and (2), with computing values (q) of different q, yields collapse onto the respective universal scaling functions, as shown in Figs. 2b and d. In particular, we find that both f + and f  are well described by the Gumbel extreme function, and where ) is the location parameter (the mode) and  + (  ) is the scale parameter. The physical mechanism underlying our above conjecture is based on the Generalised Extreme Value (GEV) theory, see Eq. ( 9). When we fit our rescaled collapsed data with the GEV distributions (see Methods for the fitting process), we find that the shape parameter  is very close to zero (see Tab. I). The GEV functions with a small  (Eq. ( 9)) and Gumbel function (Eqs. ( 3) and ( 4)) are complemented in Figs. 2b, d and f. They are hardly distinguishable, which indicates that the f + (f  ) can be well approximated by a Gumbel distribution. We analyse the cumulative distribution function (CDF), and the data are also collapsed together. The results are consistent with the PDF's, see Fig. S4 in SI. The waiting time for extreme events is also crucial for understanding the risk management of systems. As we mentioned before, for a given site and quantile q, we introduced four types of waiting time  + ,   ,  + ( + ) and  abs . We first consider their PDFs: D( + ; q), D(  ; q), D( + ; q) and D( abs ; q). Since the multiple time scales are involved from days to many years, we work with the logarithm of  .
Figures 3a, c, e and g show the results for different q of the waiting time  + ,   ,  + and  abs , respectively. We define R as the mean value of  for a certain threshold q and rescale distributions of different thresholds with the mean rate 1/R. Figures 3b, d, f and h show the rescaling results, and this data collapse implies that (take  + as an example) where G is another universal scaling function. Similarly, we can also obtain the universal scaling functions for   ,  abs and  + . Motived by Ref. [17] for earthquakes, we find that the scaling function G can be well fitted by a generalised Gamma distribution, where , , B and C are fitting parameters. All fitting values are summarised in Tab. I. The high-quality of the fit obtained using the generalised Gamma distribution is, at first sight, surprising because the fitting parameters indicate relevant deviations from the exponential behaviour, which is instead expected in the case of chaotic flows, under the hypothesis of stationarity [34]. The weather is indeed chaotic and one expects no long time memory. But, on the other hand, the climatic data we consider are not stationary because one faces -on the time scales of relevance for the waiting times -the strong modulations due to the daily and seasonal cycle. Our data analysis procedure automatically filters out the former because we consider daily averaged data [51]. As for the latter, we observe that if we restrict the analysis to waiting times shorter than approximately 90 days, thus within the sub-seasonal range (these contribute to approximate or over 90% of the measurements), the statistics of waiting times is exponential to a high degree of accuracy. The exponential law is indicated in grey in Figs. 3 (b, d, f, h) and the statistics are given in Tab. S1 in SI. To better understand whether the scaling behaviour change with the system size and regions, we investigate the finite-size effects of the universal scaling. Our previous analysis of the dataset is based on the entire Earth with a resolution of 2.5  , yielding N = 73 x 144 = 10512 nodes. Note that the nodes are not approximately homogeneously covering the entire globe due to the Earth's ellipsoid shape. To remove the non-homogeneity, we perform the procedure developed in Ref. [52]. We define the resolution (in degree latitude) at the Equator as r 0 and calculate the and then test the goodness of fit using the 95% data through the Kolmogorov-Smirnov (KS) test. '***'means p < 0.001. The waiting time is fitted in the log scale through the least square errors method; therefore, the C = exp C . The detailed information on the fitting and test are in the Methods. The same analysis is performed on the above three regions, and identical results hold, as Fig. S7 and S8 illustrate the universal scaling of magnitudes and waiting times, respectively. We find that the universal scaling functions for magnitudes are robust and very similar to the results of the entire Earth ; The universal scaling phenomenon of the magnitudes and waiting time are very well described by the Gumbel and the Gamma functions, but with different fitting parameters for the three regions. The weather/climate system is not an equilibrium system, and 'Quakes' is a feature that can characterise complex system dynamics [53,54], which has been applied in other non-equilibrium systems, including physical systems and ecology systems [55]. The quakes of complex systems gradually change both the physical and statistical properties of the system. The record is defined as the largest value of a time signal obtained up to time t. The quake times t k is the time where the record reaches a higher value, indicating when the record jumps. It is shown that the quake event times of system variable t k undergoes record statistics satisfy the log-Poission distribution [55]. In particular, the logarithmic waiting time between quakes defined as  k = log(t k )  log(t k1 ), follows a common cumulative distribution function: Within the complex climate system, quakes are often associated with extreme events. Next, we study the evolution of the system by analysing the CDF of logarithmic waiting time between quakes. Here, we measure where the air temperature positive, negative differences (CV > 0, CV < 0) and the absolute value of all the changes (|CV|) jump as positive value, negative value and absolute value quakes in Figs. 4a and b. Figures . 4c-e show that the CDF of logarithmic waiting times follows exponential forms for positive quakes, negative quakes and quakes, respectively. Our findings are highly consistent with the results in other complex systems [55], revealing the air temperature differences of the weather system are ageing. That means it is harder to reach a new record, and the release of an intrinsic tension is decelerating even under global warming. To further investigate the potential changing properties, we apply a moving time window of 3 years to analyse whether the statistical characteristics vary with time. As the log(t k ) is anticipated to follow a Poisson distribution, the average number of events dn per unit of time decays as [56]: Figure 5 demonstrates that the average number of quakes at each site in a three-year window time decays almost as one over time; however, in the rescaled plots dn dt x t = constant, we uncover that there are some deviations throughout the evolution, as shown in Figs. 5b, d and f. These deviations are usually related to the ENSO extreme events. We highlight the duration of the El Nino and La Nina events using the light orange and light green shading, respectively. We propose the underlying mechanism is that the ENSO events seem to impact the weather variability of both the Northern and Southern Hemispheres [57], leading to stronger baroclinic disturbances. This can cause stronger day-to-day time variability [58] and make extreme weather events harder to predict.
In particular, super and extreme El Nino events may become more frequent under climate change [59,60]. Nevertheless, the overall trend is consistent with relaxational record dynamics. We further test the robustness of the PDF of daily temperature differences in the time window of 10 years, as shown in Fig. S9. We find that the distributions of magnitudes at the top 5% are very similar, indicating that the magnitudes of the defined extreme events follow a steady trend. We perform the same analysis on the CMIP6 (6) [61] model data as well. Detailed information about this dataset is summarised in the Data availability section. The robust results of the same scaling phenomena validate the universal scaling phenomenon is independent of the database, see Tab. S4 in SI. Widespread extreme weather gets a boost from anthropogenic climate change and leads to significant socioeconomic ramifications in natural and human systems. In this study, we highlight the extreme values of day-to-day temperature differences. This seldom studied climatic variable provides a measure of daily weather variability that is relevant in terms of impacts on human and environmental systems because rapid temperature variations -both positive and negative -can be extremely damaging. We apply statistical physics, extreme value and record statistics to reveal the universal scaling of magnitude and waiting time of extreme events as well as the record statistics of daily air temperature variability. We find that the distribution of the extreme magnitudes is characterised by the Gumbel distribution, while the waiting time of these extreme events is governed by the generalised Gamma distribution. This analysis is also confirmed by the extreme precipitation (see Figs. S10 and S11) and CMIP6 model data (see Tab. S4). In particular, in the CMIP6 global warming experiments, we find that our universal scaling is quite robust across different models under climate change. That is, the last 43 (2058-2100) years data can be collapsed onto the distributions of the first 43 (2015-2057) years. Further scale analysis reveals that the universal scaling phenomenon is independent of the system sizes and regions. This can be explained by the renormalisation group theory of extreme values [62][63][64]. We further investigated the evolution of extreme events of day-to-day temperature variations by analysing the 'quakes'. The cumulative distribution of logarithmic waiting time of the records in the considered data, obeys to a good approximation an exponential law, indicating the dynamics has a direction and does not deviate substantially from steady state conditions. Note that this is in stark opposition with the case of daily temperatures, where records are being beaten more and more frequently as a result of climate change. While climate change seems not to impact the statistics of extreme day-to-day temperature variations, the most relevant deviation from the exponential law for the quakes is found in concurrence with ENSO events, when records are beaten more frequently than what steady state conditions would otherwise predict. After the end of the ENSO event, the statistical properties fall back to the stationary case. We discussed the underlying mechanism for this abnormality and associated it with the impacts on the weather variability of ENSO events. The analysis of quakes is a promising tool for investigating the future impact of climate change on various climatic variables, including different measures of weather variability. As mentioned above, in [45], a clear signature of climate change was found on the intramonthly standard deviation of daily surface temperature. The same applies if one considers the temperature record itself: the ongoing climate change generates strong bias in the time statistics of quakes. The result by [45] disagrees with ours, but, as mentioned above, we are using a rather different measure of daily temperature variability. Our findings extend our understanding of the nature of extreme weather events. They could help us design schemes reducing risks associated with sudden temperature changes from one day to the next. Moreover, the universal scaling in extreme events might be applied to other complex systems, such as biological, ecological and financial systems [11]. Fitting with Maximum Likelihood Estimator. (MLE) is a general and flexible method of estimation of the unknown parameters, and each set of the parameters attaches to different probability densities of the observed data. The probability of the empirical data as a function of parameters is called the likelihood function. Values of high likelihood parameters correspond to models that give high probability to the observed data. Fitting with the least square errors. This article uses the least square errors method to fit the lines for the collapsed waiting time data. Specifically, we use the Levenberg-Marquardt algorithm as a standard tool for nonlinear fitting. It combines two minimization methods, namely gradient descent and the Gauss-Newton method. Kolmogorov-Smirnov goodness of fit test. We apply the Kolmogorov-Smirnov test (KS test) to the magnitudes, waiting times of the extreme events, and quakes to assess the goodness of fit. The Kolmogorov-Smirnov statistic is defined as D n = sup The probability density function of Generalized Extreme Value, GEV(, , ) distribution family is defined for the maximum (or minimum) values of random variables, mathematically: and the cumulative density function is: where  is the shape parameter and it can be any real number,  and  are the local and scale parameters. When  = 0, the distribution is the Gumbel distribution which is widely used in climate modelling [34]. a family distribution H(y) defined for the excesses y = X  u of a sequence of independent random variables X 1 , X 2 , X 3    with common distribution family. Theorem: (11) If convergence between the GEV and GP distributions is realised,  =  and  =  +(u). This theorem implies that the parameters of the generalized Pareto distribution of threshold excesses are uniquely determined by those of the associated GEV distribution G. As a result, once we estimate the parameters for the GP parameters, we can derive the corresponding GEV and vice versa. In particular, the  for the two distributions are the same.
The GP scale parameter depends on the chosen threshold, we usually look at the modified scale parameter which is defined as  = u u. A record within an ordered series of numerical data obtained from observations and measurements is simply a value larger (or smaller) than all preceding values. The mathematical properties of records depend on the properties of the series from which they are generated. And, the record size fluctuations of, e.g. the energy can push the system to cross the 'edge of stability, inducing irreversible configurational changes. The statistics of those changes closely follow the statistics of record fluctuations. Consider a sequence of independent random numbers drawn from the same distribution ar times 2,    , t. When this time series is sufficiently long, the number of records found between the first and the tth draw has a Poisson distribution, with an average ln t. That is: where P n (t) is the probability of having n records in t successive trials, where 1  n  t. It is a Poisson distribution when replacing usual time argument t with ln t. Therefore, the differences between two Poisson processes, the ratios ln(t k /t k1 ) are independent random numbers with an exponential distributions, where t k is the time when the kth record appears. Let n(t) be the average number of vents in time t, we have n(t) = ln t. Furthermore, the average number of events per unit of time decays as: dn dt  1 t . gram for Climate Model Diagnosis and Intercomparison' [61].The data spans from 2015 to 2100. More details of the models are described in SI. The authors would like to thank Prof. Jurgen Kurths, Dr. Xiaojie Chen and Teng Liu for many exchanges on extreme events and the access of dataset and appreciate Dr. Nanxin Wei's useful discussions on the scaling phenomena. The authors acknowledge the support of the of China ( 12135003) and Liyun Postdoctoral Program of Beijing Normal University. V.L. acknowledges the support project TiPES ( 820970) and from the EPSRC The authors would like to thank Prof. Jurgen Kurths, Dr. Xiaojie Chen and Teng Liu for many exchanges on extreme events and the access of dataset and appreciate Dr. Nanxin Wei's useful discussions on the scaling phenomena. The authors acknowledge the support of the of China ( 12135003) and Liyun Postdoctoral Program of Beijing Normal University. V.L. acknowledges the support project TiPES ( 820970) and from the EPSRC The primary dataset is the Reanalysis Air temperature ERA5 (ECMWF Reanalysis 5th Generation) 1000hPa Daily (time:00:00) from 1979-2021. We deleted the data point of 29th Feb in any leap year [46]. ECMWF is the European Centre for Medium-Range Weather Forecasts. The CMIP6 (6) model data is provided by 'The Pro-AUTHOR CONTRIBUTIONS J.F., X.C. designed the research. Q.Y. performed the analysis and prepared the manuscript, Q.Y., J.M., J.F., V.L., H.J., K.C., X.C. generated research ideas and discussed results, and contributed to writing the manuscript. Supplementary Information is available in the version of the paper. The authors declare no competing interests. The primary dataset is the Reanalysis Air temperature ERA5 (ECMWF Reanalysis 5th Generation) 1000hPa Daily (time:00:00) from 1979-2021. We deleted the data point of 29th Feb in any leap year [46]. ECMWF is the European Centre for Medium-Range Weather Forecasts. The CMIP6 (6) model data is provided by 'The Pro-AUTHOR CONTRIBUTIONS J.F., X.C. designed the research. Q.Y. performed the analysis and prepared the manuscript, Q.Y., J.M., J.F., V.L., H.J., K.C., X.C. generated research ideas and discussed results, and contributed to writing the manuscript. Supplementary Information is available in the version of the paper. The authors declare no competing interests.