# Emergence of universal scaling in weather extreme events

Qing Yao  School of Systems Science/Institute of Nonequilibrium Systems,
Beijing Normal University, 100875 Beijing, China  Blackett Laboratory and
Centre for Complexity Science, Imperial College London, South Kensington
Campus, London SW7 2AZ, United Kingdom  Jingfang Fan<EMAIL_ADDRESS>School of Systems Science/Institute of Nonequilibrium Systems, Beijing Normal
University, 100875 Beijing, China  Potsdam Institute for Climate Impact
Research, 14412 Potsdam, Germany  Jun Meng  School of Science, Beijing
University of Posts and Telecommunications, 100876 Beijing, China  Potsdam
Institute for Climate Impact Research, 14412 Potsdam, Germany  Valerio
Lucarini  Department of Mathematics and Statistics, University of Reading,
Reading, RG6 6AX, United Kingdom  Centre for the Mathematics of Planet Earth,
University of Reading, Reading, RG6 6AX, United Kingdom  School of Systems
Science/Institute of Nonequilibrium Systems, Beijing Normal University, 100875
Beijing, China  Henrik Jeldtoft Jensen  Department of Mathematics and Centre
for Complexity Science, Imperial College London, South Kensington Campus,
London SW7 2AZ, United Kingdom  Institute of Innovative Research, Tokyo
Institute of Technology, 4259, Nagatsuta-cho, Yokohama, 226-8502, Japan  Kim
Christensen  Blackett Laboratory and Centre for Complexity Science, Imperial
College London, South Kensington Campus, London SW7 2AZ, United Kingdom
Xiaosong Chen  School of Systems Science/Institute of Nonequilibrium Systems,
Beijing Normal University, 100875 Beijing, China

(August 29, 2024)

###### Abstract

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.

## Introduction

Natural or social systems far from equilibrium often show complex dynamic
behaviours. Nevertheless, increasing empirical evidence supports that the
fields as diverse as physics odor2004universality ; fan2020universal , biology
munoz2018colloquium ; di2018landau , ecology sole2007scaling ;
ross2021universal , human mobility song2010modelling , economics, and
financial markets bouchaud1990anomalous ; mantegna1995scaling  may all, at
least up to a certain degree of approximation, follow the principles of scale
invariance and universality stanley2000scale . 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
sornette2006critical ; west2017scale , such as the climate system
ashkenazy2003nonlinearity ; fan2021statistical  and earthquakes bak2002unified
; saichev2006universal ; corral2004long ; livina_memory_2005 .

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 ghil2020physics ;
franzke2020structure . 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 katz1992extreme ; karl1995trends ; Easterling2000 ;
Ipcc2012 . 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 Allen2003 ; trenberth2015attribution ; Otto2017 . The World
Meteorological Organization (WMO) reported that weather-related disasters have
increased over the past 50 years, causing more damage, but early warnings save
lives WMO . 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 WMO . 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 Poumadere2005 ; shi2015impacts , crop
production wheeler2000temperature , and economic growth kotz2021day .
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
albeverio2006extreme ; lucarini2016extremes ; Galfi2017 ; Faranda2017 ;
Hochman2019 ; Galfi2021 . 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 helbing2013globally .

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 Leadbetter1983 ; Katz2002 ; Ghil2011 ; coles2001introduction ;
lucarini2016extremes , 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 Holton2004 . 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
Kotz2021b . 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.

### Extreme Events

We study the daily temperature variations using the differences in temperature
between two adjacent discrete time periods, that is
$\mathrm{CV}_{i,d}=x_{i,d}-x_{i,d-1},\leavevmode\nobreak\
i=1,\leavevmode\nobreak\ \dots,N,d=2\leavevmode\nobreak\ \dots D$, where
$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\times 144$) covering
the entire Earth based on the ECMWF Reanalysis 5th Generation (ERA5)
hersbach2020era5  dataset, see Data availability for more details.
Furthermore, the absolute values of the changes
$\lvert\mathrm{CV}_{i,d}\rvert=\lvert x_{i,d}-x_{i,d-1}\rvert$ provides a
measure of the strength of weather variability karl1995trends .

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^{\circ}$N,
$3.75^{\circ}$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 follows 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 Holton2004 .

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 coles2001introduction  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 coles2001introduction ; galfi2017convergence .
This regime is considerably lower than what is usually found when looking at
temperature extremes, see, e.g. Zahid2017 , 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.

Figure 1: The illustration of weather extreme events defined for day-to-day
air temperature variability. a and c are air temperatures of single days
(April 1, 2020 and August 14, 2021) visualised through nearside perspective
projection; b and d are the air temperature differences between the next day
and the chosen days. We choose a single site in Madrid ($40.25^{\circ}$N,
$3.75^{\circ}$W), Spain, to illustrate the definitions of extreme events for
day-to-day air temperature variability. e are the daily air temperatures (AT)
from January 2020 to Dec 2021. f is the day-to-day temperature differences,
$\mathrm{CV}$, and the red bars and navy bars represent the positive and
negative differences between the two adjacent days, respectively. The
horizontal yellow (sea green) lines are the threshold of $97.5\%$ positive
(negative) changes. The day-to-day temperature differences are greater (lower)
the thresholds are marked by vertical yellow (sea green) lines, which indicate
the happening of the extreme events defined for the given thresholds. The
width of the yellow (sea green) rectangular is the waiting time between
positive (negative) extremes.

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
$\mathcal{E}^{+}_{i}(q)$ is defined when the value $\mathrm{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}=\frac{\mathrm{CV}_{ev^{+}}}{\sigma_{CV_{ev^{+}}}},\leavevmode\nobreak\
ev\in\mathcal{E}^{+}_{i}$, where $\sigma_{CV_{ev^{+}}}$ is the standard
deviation of the extreme part, $\mathrm{CV}_{ev^{+}}$ for site $i$. Similarly,
a negative extreme event set $\mathcal{E}^{-}_{i}(q)$ is defined when the
value $\mathrm{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
$m^{-}_{ev}=-\frac{\mathrm{CV}_{ev^{-}}}{\sigma_{CV_{ev^{-}}}},\leavevmode\nobreak\
ev\in\mathcal{E}^{-}_{i}$. 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}=\frac{\lvert\mathrm{CV}\rvert_{ev^{abs}}}{\sigma_{CV_{ev}^{abs}}},\leavevmode\nobreak\
ev\in\mathcal{E}^{abs}_{i}$.

The waiting time $\tau_{i}\in\mathrm{T}_{i}$ of the extreme events in
$\mathcal{E}_{i}$ is the number of days between two consecutive events
boers2019complex . Therefore, we can define four types of waiting times:
$\tau^{+}$ ($\tau^{-}$, $\tau^{abs}$) is the waiting time between positive
(negative, absolute) extreme temperature difference, i.e. the waiting time in
$\mathcal{E}^{+}$ ($\mathcal{E}^{-}$, $\mathcal{E}^{abs}$). The fourth type
$\tau^{+-}$ ($\tau^{-+}$ interchangeable) is the number of days between
extreme positive and extreme negative change events (or extreme negative and
extreme positive), i.e. the waiting time in
$(\mathcal{E}^{+}\cup\mathcal{E}^{-})$ and the detailed description can be
found in Fig. S5 in SI.

## Results

### Universal scaling for magnitudes of extreme temperature differences

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 $\mathrm{CV}$. We choose
$q=90\%,92.5\%,95\%$ and $97\%$ as threshold quantiles within the extreme
regime.

Figure 2: The magnitudes of extreme day-to-day temperature, a, c and e are
probability density distribution for different chosen quantiles (thresholds)
for positive, negative and absolute extreme events respectively; b, d and f
are universal scaling functions. The fitting method is Maximum Likelihood
Estimation (MLE), and we use the Kolmogorov–Smirnov(KS) test for the
cumulative distribution of $95\%$. The fitting parameters are shown in Tab. 1.
The fitted lines by GEV are plotted in grey, and the lines which presented
Gumbel distribution are in dashed black lines. These two distributions are
hardly distinguishable in b, d and f for PDF. The associated cumulative
distribution functions are in SI, see Fig. S4.

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,

$p(m^{+}_{ev};q)=\frac{1}{\sigma(q)}f^{+}\left(\frac{m^{+}_{ev}-\langle
m^{+}_{ev}\rangle}{\sigma(q)}\right),$ (1)

and

$p(m^{-}_{ev};q)=\frac{1}{\sigma(q)}f^{-}\left(\frac{m^{-}_{ev}-\langle
m^{-}_{ev}\rangle}{\sigma(q)}\right),$ (2)

where $f^{+}$ and $f^{-}$ are two universal scaling functions, $\langle
m^{+}_{ev}\rangle$ ($\langle m^{-}_{ev}\rangle$) and $\sigma(q)$ is the mean
value and standard deviation associated with quantile $q$. We then define two
rescaled parameters $\omega^{+}=\frac{m^{+}_{ev}-\langle
m^{+}_{ev}\rangle}{\sigma(q)}$ and $\omega^{-}=\frac{m^{-}_{ev}-\langle
m^{-}_{ev}\rangle}{\sigma(q)}$. Combining the proposed scaling forms, Eqs. (1)
and (2), with computing values $\sigma(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,

$f^{+}(z^{+})=\frac{1}{\beta}e^{-e^{-z^{+}}-z^{+}},$ (3)

and

$f^{-}(z^{-})=\frac{1}{\beta}e^{-e^{-z^{-}}-z},$ (4)

where $z^{+}=(\omega^{+}-\mu^{+})/\beta^{+}$ and
$z^{-}=(\omega^{-}-\mu^{-})/\beta^{-}$. $\mu^{+}$ ($\mu^{-}$) is the location
parameter (the mode) and $\beta^{+}$ ($\beta^{-}$) 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 $\xi$ is very close to zero (see
Tab. 1). The GEV functions with a small $\xi$ (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.

### Universal scaling for waiting time between extreme events

Figure 3: Waiting time of extreme day-to-day temperature differences. The four
columns are $\tau^{+}$, $\tau^{-}$, $\tau^{+-}$ and $\tau^{abs}$,
respectively. The fitting parameters are in Tab. 1. a, c, e and g are log-
binned probability density distribution for different chosen quantiles with
the base of $1.3$ (pruessner2012self, , Ch. 7, p. 233). The dashed grey
vertical lines indicate $365$ days, while the solid grey vertical lines mark
$90$ days. b, d, f and h are the universal scaling functions. The fitting
parameters for the generalised Gamma function, Eq. (6), are summarised in Tab.
1. The statistics of the mean value $R$ of waiting times and fitting
parameters for exponential function can be found in Tab. S1 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 $\tau^{+}$, $\tau^{-}$,
$\tau^{+-}$ ($\tau^{-+}$) and $\tau^{abs}$. We first consider their PDFs:
$D(\tau^{+};q)$, $D(\tau^{-};q)$, $D(\tau^{+-};q)$ and $D(\tau^{abs};q)$.
Since the multiple time scales are involved from days to many years, we work
with the logarithm of $\tau$. Figures 3a, c, e and g show the results for
different $q$ of the waiting time $\tau^{+}$, $\tau^{-}$, $\tau^{+-}$ and
$\tau^{abs}$, respectively. We define $R$ as the mean value of $\tau$ 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 $\tau^{+}$ as an example)

$D(\tau^{+};q)=\frac{1}{R}\mathcal{G}\left(\tau^{+}/R\right),$ (5)

where $\mathcal{G}$ is another universal scaling function. Similarly, we can
also obtain the universal scaling functions for $\tau^{-}$, $\tau^{abs}$ and
$\tau^{+-}$. Motived by Ref. corral2004long  for earthquakes, we find that the
scaling function $\mathcal{G}$ can be well fitted by a generalised Gamma
distribution,

$\mathcal{G}(u;\gamma,\eta,B,C)=C\frac{1}{u^{1-\gamma}}\exp\left(-\frac{u^{\eta}}{B}\right),$
(6)

where $\gamma,\eta,B$ and $C$ are fitting parameters. All fitting values are
summarised in Tab. 1.

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
lucarini2016extremes . 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 111Note that in continental and arid areas, day-to-night
excursions of more than 30 ${}^{\circ}C$ are far from atypical. 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.

| Magnitude
---|---
| $\xi$ | $\mu$ | $\beta$ | $R^{2}$ | KS-Dn($95\%$)
$m^{+}_{ev}$ |  $-0.0155$ $(0.0001)$ |  $-0.4487$ $(0.0001)$ |  $0.7962$ $(0.0001)$ | $0.997$ |  $0.56$***
$m^{-}_{ev}$ |  $0.0410$ $(0.0001)$ |  $-0.4647$ $(0.0001)$ |  $0.7472$ $(0.0001)$ | $0.992$ |  $0.65$***
$m^{abs}_{ev}$ |  $-0.0004$ $(0.0001)$ |  $-0.4538$ $(0.0001)$ |  $0.7828$ $(0.0001)$ | $0.994$ |  $0.64$***
| Waiting time
| $\gamma$ | $\eta$ | $B$ | $C^{\prime}$ | $R^{2}$
$\tau^{+}$ |  $1.04$ $(0.36)$ |  $0.56$ $(0.11)$ |  $0.44$ $(0.23)$ |  $1.03$ $(1.24)$ | $0.966$
$\tau^{-}$ |  $1.36$ $(0.50)$ |  $0.47$ $(0.11)$ |  $0.29$ $(0.17)$ |  $2.22$ $(2.03)$ | $0.964$
$\tau^{+-}$ |  $0.73$ $(0.54)$ |  $0.39$ $(0.11)$ |  $0.33$ $(0.25)$ |  $1.59$ $(2.37)$ | $0.982$
$\tau^{abs}$ |  $0.40$ $(0.18)$ |  $0.85$ $(0.13)$ |  $1.49$ $(0.72)$ |  $-0.75$ $(0.41)$ | $0.964$
Table 1: Table of fitting parameters. The numbers in the brackets are the
standard errors of the nonlinear fitting methods. We fit the magnitudes with
the Maximum Likelihood Estimator (MLE) 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^{\prime}$. The detailed
information on the fitting and test are in the Methods.

### Universal scaling of finite-size analysis

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^{\circ}$, yielding $N=73\times 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. meng_percolation_2017 . We define the resolution
(in degree latitude) at the Equator as $r_{0}$ and calculate the number of
nodes $n_{0}=360/r_{0}$. Then the number of nodes in latitude $mr_{0}$ is
$n_{m}$ $=n_{0}\cos\left(mr_{0}\right)$, where
$m\in\left[-90/r_{0},90/r_{0}\right]$. The total number of nodes is
$N^{\prime}=\sum_{m=0}^{m=90/r_{0}}2n_{m}-n_{0}$. Here, we choose $r_{0}$ to
be $2.5^{\circ}$, which yields $N^{\prime}=6570$. Furthermore, we divide the
globe (6570 nodes) into three ‘Tropics’ ($25^{\circ}$S to $25^{\circ}$N),
‘Mid-latitude’ ($60^{\circ}$S to $27.5^{\circ}$S and $27.5^{\circ}$N to
$60^{\circ}$S), as well as ‘High latitude’ regions ($90^{\circ}$S to
$62.5^{\circ}$S and $62.5^{\circ}$N to $90^{\circ}$N), see Fig. S6 in SI.

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
(Fig. 2); 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.

### Evolution of the climate system

The weather/climate system is not an equilibrium system, and ‘Quakes’ is a
feature that can characterise complex system dynamics sibani1993slow ;
jensen2013stochastic , which has been applied in other non-equilibrium
systems, including physical systems and ecology systems anderson2004evolution
. 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
anderson2004evolution . In particular, the logarithmic waiting time between
quakes defined as $\tau_{k}=\log(t_{k})-\log(t_{k-1})$, follows a common
cumulative distribution function:

$P(\tau_{k})=1-e^{-\tau_{k}}.$ (7)

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 ($\mathrm{CV}>0$, $\mathrm{CV}<0$)
and the absolute value of all the changes ($\lvert\mathrm{CV}\rvert$) jump as
positive value, negative value and absolute value quakes in Figs. 4a and b.

Figure 4: The quake times of the temperature differences and CDF of their
logarithmic waiting time for a site around Madrid, which is from the same time
signal as Fig. 1e. a, The red bars represent the positive changes,
$\mathrm{CV}>0$, while the blue bars represent negative changes,
$\mathrm{CV}<0$. The red line represents the largest positive value of
$\mathrm{CV}>0$ up to time $t$, and the blue line represents the largest
negative value (i.e. the smallest values) of $\mathrm{CV}<0$ up to time $t$.
The upper plot is the zoom-in part of the year 1979. b, The purple bars
represent the absolute values of a, the purple line indicates the largest
value of $\rvert\mathrm{CV}\rvert$ up to time $t$, and the b.1 is the zoom in
part of the year 1979. The red (dark blue) circles in c, blue circles in d and
purple circles in d represent the one minus cumulative distribution of $\log$
waiting time of all the site records of air temperature differences for the
positive, negative and absolute value quakes, respectively. The three plots
follow an exponential, represented by the black lines.

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 anderson2004evolution , 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 $\mathrm{d}\bar{n}$ per unit of
time decays as sibani2009record :

$\frac{\mathrm{d}\bar{n}}{\mathrm{d}t}\propto\frac{1}{t}.$ (8)

Figure 5: The average number of quakes at a site in a unit time against the
time in a linear-linear plot. a, c and e plot the average number of quakes
$\mathrm{\bar{n}}$ each site per year for $\mathrm{CV}>0$, $\mathrm{CV}<0$ and
$\lvert\mathrm{CV}\rvert$ respectively. The quakes per year are averaged over
a 3-year window time, and $\overline{n}$ is averaged over $N=10512$ sites.
$1/t$ is plotted as a black line. b, d and f plot rescaled average number of
the quakes of a, c and e against the time correspondingly. If
$\mathrm{\bar{n}}$ decays as $1/t$, then the rescaled number should lie in a
horizontal line. The error bars are indicated in black lines. This plot
reveals clearly deviations away from a proposed scaling law from the year 2010
to 2019. The light orange and light green shadings indicate the El Niño and La
Niña periods. We shuffle the dataset for $2000$ times to generate a null time
series. The results are presented in grey with error bars smaller than the
symbol size.

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 $\frac{\mathrm{d}\bar{n}}{\mathrm{d}t}\times
t=\text{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 Niño and La
Niña 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
yang_impacts_2020 , leading to stronger baroclinic disturbances. This can
cause stronger day-to-day time variability machado_influence_2021  and make
extreme weather events harder to predict. In particular, super and extreme El
Niño events may become more frequent under climate change wang_historical_2019
; cai_increasing_2014 . 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 (Coupled Model Intercomparison Project Phase 6)
gmd-9-1937-2016  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.

## Discussion

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 theory 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 gyorgyi2010renormalization ; calvo2012extreme ;
corral2015scaling .

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 system’s 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 Kotz2021b ,
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 a strong bias
in the time statistics of quakes. The result by Kotz2021b  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 for 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 sornette2006critical .

## Methods

### Fitting methods

Fitting with Maximum Likelihood Estimator. The Maximum Likelihood Estimation
(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.

### Goodness of fit test

Kolmogorov–Smirnov goodness of fit test. We apply the Kolmogorov–Smirnov test
(KS test) to the magnitudes, waiting times of the extreme events, and the
quakes to assess the goodness of fit. The Kolmogorov–Smirnov statistic is
defined as $D_{n}=\underset{x}{\sup}\lvert F_{n}(x)-F(x)\rvert$, where
$F_{n}(x)$ and $F(x)$ are the empirical distribution function of the sample
and the cumulative distribution function of the reference distribution
respectively. This statistic quantifies the distance between those two
distributions. By the Glivenko–Cantelli theorem, if the sample comes from
distribution $F(x)$, then $D_{n}$ converges to 0 almost surely in the limit
when the number of samples $n$ goes to infinity. When
$\sqrt{n}D_{n}>K_{\alpha}$, where $\mathrm{Pr}(K\leq K_{\alpha})=1-\alpha$,
the test can reject the null hypothesis (the sample comes from distribution
$F(x)$) at the $1-\alpha$ confidence level.

### Generalised Extreme Value Theory

The probability density function of $\textrm{Generalized Extreme
Value},\leavevmode\nobreak\ \mathrm{GEV}(\mu,\beta,\xi)$ distribution family
is defined for the maximum (or minimum) values of random variables,
mathematically:

$\displaystyle p(z)$ $\displaystyle=\frac{1}{\beta}t(z)^{\xi+1}e^{-t(z)},$ (9)
$\displaystyle t(z)$
$\displaystyle=\left\\{\begin{matrix}(1+\xi(\frac{z-\mu}{\beta}))^{-\frac{1}{\xi}}&\xi\neq
0\\\ e^{-\frac{z-\mu}{\beta}}&\xi=0,\end{matrix}\right.$

and the cumulative density function is:

$G(z)=\left\\{\begin{matrix}\exp\left\\{-\left[1+\xi(\frac{z-\mu}{\beta})\right]^{-1/\xi}\right\\}&\xi\neq
0\\\
\exp\left\\{-\exp\left[(\frac{z-\mu}{\beta})\right]\right\\}&\xi=0,\end{matrix}\right.$
(10)

where $\xi$ is the shape parameter and it can be any real number, $\mu$ and
$\beta$ are the local and scale parameters. When $\xi=0$, the distribution is
the Gumbel distribution which is widely used in climate modelling
lucarini2016extremes .

### Peaks Over Threshold

The Generalised Pareto Distribution is a family distribution $H(y)$ defined
for the excesses $y=\mathbf{X}-u$ of a sequence of independent random
variables $X_{1},X_{2},X_{3}\cdots$ with common distribution family.

Theorem:
$M_{n}=\mathrm{max}\left\\{X_{1},\cdots,X_{n}\right\\},\mathrm{Pr}\left\\{M_{n}\leq
z\right\\}\approx G(z)$, where
$G(z)=\exp\left\\{-\left[1+\xi(\frac{z-\mu}{\beta})\right]^{-1/\xi}\right\\}$,
for some $\mu,\beta>0$ and $\xi\neq 0$ while
$G(z)=\exp\left\\{-\exp\left[(\frac{z-\mu}{\beta})\right]\right\\}$ for some
$\mu,\beta>0$ and $\xi=0$. $G(z)$ represents the GEV family. Then for large
enough $u$, conditional on $X>u$, approximately,

$H(y)=\left\\{\begin{matrix}1-\left(1+\frac{\tilde{\xi}y}{\tilde{\sigma}}\right)^{1/\tilde{\xi}}&\tilde{\xi}\neq
0,\\\
1-\exp\left(\frac{y}{\tilde{\sigma}}\right)&\tilde{\xi}=0.\end{matrix}\right.$
(11)

If convergence between the GEV and GP distributions is realised,
$\xi=\tilde{\xi}$ and $\tilde{\sigma}=\beta+\xi(u-\mu)$. 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 $\xi$ 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 $\sigma=\tilde{\sigma}_{u}-\xi u$.

### Record Statistics and Dynamics

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 $1$, $2$, $\cdots$, $t$. When this time series is
sufficiently long, the number of records found between the first and the $t$th
draw has a Poisson distribution, with an average $\ln t$. That is:

$P_{n}(t)=\frac{(\ln t)^{n-1}}{(n-1)!}\frac{1}{t},$ (12)

where $P_{n}(t)$ is the probability of having $n$ records in $t$ successive
trials, where $1\leq n\leq 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_{k-1})$ are independent random
numbers with an exponential distributions, where $t_{k}$ is the time when the
$k$th record appears. Let $\overline{n(t)}$ be the average number of vents in
time t, we have $\overline{n(t)}=\ln t$. Furthermore, the average number of
events per unit of time decays as:
$\frac{\mathrm{d}\bar{n}}{\mathrm{d}t}\propto\frac{1}{t}$.

## Data availability

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 hersbach2020era5 . ECMWF is the European
Centre for Medium-Range Weather Forecasts. The CMIP6 (Coupled Model
Intercomparison Project Phase 6) model data is provided by ‘The Program for
Climate Model Diagnosis and Intercomparison’ gmd-9-1937-2016 .The data spans
from 2015 to 2100. More details of the models are described in SI.

## Acknowledgments

The authors would like to thank Prof. Jürgen 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 National Natural Science Fundation of
China (Grant No. 12135003) and Liyun Postdoctoral Program of Beijing Normal
University. V.L. acknowledges the support received from the Horizon 2020
project TiPES (grant no. 820970) and from the EPSRC Grant EP/T018178/1.

## 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.

## Additional information

Supplementary Information is available in the online version of the paper.

## Competing interests

The authors declare no competing interests.

## References

  * (1) Ódor, G.  Universality classes in nonequilibrium lattice systems.  _Reviews of modern physics_ 76, 663 (2004).
  * (2) Fan, J. _et al._ Universal gap scaling in percolation.  _Nature Physics_ 16, 455–461 (2020).
  * (3) Munoz, M. A.  Colloquium: Criticality and dynamical scaling in living systems.  _Reviews of Modern Physics_ 90, 031001 (2018).
  * (4) Di Santo, S., Villegas, P., Burioni, R. & Muñoz, M. A.  Landau–ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization.  _Proceedings of the National Academy of Sciences_ 115, E1356–E1365 (2018).
  * (5) Solé, R.  Scaling laws in the drier.  _Nature_ 449, 151–153 (2007).
  * (6) Ross, S. R.-J. _et al._ Universal scaling of robustness of ecosystem services to species loss.  _Nature communications_ 12, 1–7 (2021).
  * (7) Song, C., Koren, T., Wang, P. & Barabási, A.-L.  Modelling the scaling properties of human mobility.  _Nature physics_ 6, 818–823 (2010).
  * (8) Bouchaud, J.-P. & Georges, A.  Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications.  _Physics reports_ 195, 127–293 (1990).
  * (9) Mantegna, R. N. & Stanley, H. E.  Scaling behaviour in the dynamics of an economic index.  _Nature_ 376, 46–49 (1995).
  * (10) Stanley, H. E. _et al._ Scale invariance and universality: organizing principles in complex systems.  _Physica A: Statistical Mechanics and its Applications_ 281, 60–68 (2000).
  * (11) Sornette, D.  _Critical phenomena in natural sciences: chaos, fractals, selforganization and disorder: concepts and tools_ (Springer Science & Business Media, 2006).
  * (12) West, G.  _Scale: The Universal Laws of Growth, Innovation, Sustainability, and the Pace of Life in Organisms, Cities, Economies, and Companies_ (Penguin Press, 2017).
  * (13) Ashkenazy, Y., Baker, D. R., Gildor, H. & Havlin, S.  Nonlinearity and multifractality of climate change in the past 420,000 years.  _Geophysical research letters_ 30 (2003).
  * (14) Fan, J. _et al._ Statistical physics approaches to the complex earth system.  _Physics reports_ 896, 1–84 (2021).
  * (15) Bak, P., Christensen, K., Danon, L. & Scanlon, T.  Unified scaling law for earthquakes.  _Physical Review Letters_ 88, 178501 (2002).
  * (16) Saichev, A. & Sornette, D.  “universal” distribution of interearthquake times explained.  _Physical Review Letters_ 97, 078501 (2006).
  * (17) Corral, Á.  Long-term clustering, scaling, and universality in the temporal occurrence of earthquakes.  _Physical Review Letters_ 92, 108501 (2004).
  * (18) Livina, V. N., Havlin, S. & Bunde, A.  Memory in the Occurrence of Earthquakes.  _Physical Review Letters_ 95, 208501 (2005).
  * (19) Ghil, M. & Lucarini, V.  The physics of climate variability and climate change.  _Reviews of Modern Physics_ 92, 035002 (2020).
  * (20) Franzke, C. L. _et al._ The structure of climate variability across scales.  _Reviews of Geophysics_ 58, e2019RG000657 (2020).
  * (21) Katz, R. W. & Brown, B. G.  Extreme events in a changing climate: variability is more important than averages.  _Climatic change_ 21, 289–302 (1992).
  * (22) Karl, T. R., Knight, R. W. & Plummer, N.  Trends in high-frequency climate variability in the twentieth century.  _Nature_ 377, 217–220 (1995).
  * (23) Easterling, D. R. _et al._ Climate extremes: Observations, modeling, and impacts.  _Science_ 289, 2068–2074 (2000).
  * (24) IPCC.  _Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation_ (Cambridge Univerity Presse, 2012).
  * (25) Allen, M.  Liability for climate change.  _Nature_ 421, 891–892 (2003).
  * (26) Trenberth, K. E., Fasullo, J. T. & Shepherd, T. G.  Attribution of climate extreme events.  _Nature Climate Change_ 5, 725–730 (2015).
  * (27) Otto, F. E.  Attribution of weather and climate events.  _Annual Review of Environment and Resources_ 42, 627–646 (2017).
  * (28) Organization, W. M.  Weather-related disasters increase over past 50 years, causing more damage but fewer deaths.  https://public.wmo.int/en/media/press-release/weather-related-disasters-increase-over-past-50-years-causing-more-damage-fewer (2021).
  * (29) Poumadére, M., Mays, C., Le Mer, S. & Blong, R.  The 2003 heat wave in france: Dangerous climate change here and now.  _Risk Analysis_ 25, 1483–1494 (2005).
  * (30) Shi, L., Kloog, I., Zanobetti, A., Liu, P. & Schwartz, J. D.  Impacts of temperature and its variability on mortality in new england.  _Nature climate change_ 5, 988–991 (2015).
  * (31) Wheeler, T. R., Craufurd, P. Q., Ellis, R. H., Porter, J. R. & Prasad, P. V.  Temperature variability and the yield of annual crops.  _Agriculture, Ecosystems & Environment_ 82, 159–167 (2000).
  * (32) Kotz, M., Wenz, L., Stechemesser, A., Kalkuhl, M. & Levermann, A.  Day-to-day temperature variability reduces economic growth.  _Nature Climate Change_ 11, 319–325 (2021).
  * (33) Albeverio, S., Jentsch, V. & Kantz, H. (eds.) _Extreme Events in Nature and Society_.  The Frontiers Collection (Springer Berlin Heidelberg, 2006).
  * (34) Lucarini, V. _et al._ _Extremes and recurrence in dynamical systems_ (John Wiley & Sons, 2016).
  * (35) Gálfi, V. M., Bódai, T. & Lucarini, V.  Convergence of extreme value statistics in a two-layer quasi-geostrophic atmospheric model.  _Complexity_ 2017, 5340858 (2017).
  * (36) Faranda, D., Messori, G. & Yiou, P.  Dynamical proxies of north atlantic predictability and extremes.  _Scientific Reports_ 7, 41278 (2017).
  * (37) Hochman, A., Alpert, P., Harpaz, T., Saaroni, H. & Messori, G.  A new dynamical systems perspective on atmospheric predictability: Eastern mediterranean weather regimes as a case study.  _Science Advances_ 5, eaau0936 (2019).
  * (38) Gálfi, V. M., Lucarini, V., Ragone, F. & Wouters, J.  Applications of large deviation theory in geophysical fluid dynamics and climate science.  _La Rivista del Nuovo Cimento_ 44, 291–363 (2021).
  * (39) Helbing, D.  Globally networked risks and how to respond.  _Nature_ 497, 51–59 (2013).
  * (40) Leadbetter, M. R., Lindgren, G. & Rootzen, H.  _Extremes and Related Properties of Random Sequences and Processes_.  Springer Series in Statistics (Springer Verlag, 1983).
  * (41) Katz, R. W., Parlange, M. B. & Naveau, P.  Statistics of extremes in hydrology.  _Advances in Water Resources_ 25, 1287–1304 (2002).
  * (42) Ghil, M. _et al._ Extreme events: dynamics, statistics and prediction.  _Nonlinear Processes in Geophysics_ 18, 295–350 (2011).
  * (43) Coles, S., Bawa, J., Trenner, L. & Dorazio, P.  _An introduction to statistical modeling of extreme values_ , vol. 208 (Springer, 2001).
  * (44) Holton, J. R.  _An introduction to dynamic meteorology_.  International Geophysics Series (Elsevier Academic Press,, MA, 2004), 4 edn.
  * (45) Kotz, M., Wenz, L. & Levermann, A.  Footprint of greenhouse forcing in daily temperature variability.  _Proceedings of the National Academy of Sciences_ 118, e2103294118 (2021).
  * (46) Hersbach, H. _et al._ The era5 global reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ 146, 1999–2049 (2020).
  * (47) Gálfi, V. M., Bódai, T. & Lucarini, V.  Convergence of extreme value statistics in a two-layer quasi-geostrophic atmospheric model.  _Complexity_ 2017 (2017).
  * (48) Zahid, M., Blender, R., Lucarini, V. & Bramati, M. C.  Return levels of temperature extremes in southern pakistan.  _Earth System Dynamics_ 8, 1263–1278 (2017).
  * (49) Boers, N. _et al._ Complex networks reveal global pattern of extreme-rainfall teleconnections.  _Nature_ 566, 373–377 (2019).
  * (50) Pruessner, G.  _Self-organised criticality: theory, models and characterisation_ (Cambridge University Press, 2012).
  * (51) Note that in continental and arid areas, day-to-night excursions of more than 30 ${}^{\circ}C$ are far from atypical.
  * (52) Meng, J., Fan, J., Ashkenazy, Y. & Havlin, S.  Percolation framework to describe El Niño conditions.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 27, 035807 (2017).
  * (53) Sibani, P. & Littlewood, P. B.  Slow dynamics from noise adaptation.  _Physical review letters_ 71, 1482 (1993).
  * (54) Sibani, P. & Jensen, H. J.  _Stochastic dynamics of complex systems: From glasses to evolution_ , vol. 2 (World Scientific Publishing Company, 2013).
  * (55) Anderson, P. E., Jensen, H. J., Oliveira, L. & Sibani, P.  Evolution in complex systems.  _Complexity_ 10, 49–56 (2004).
  * (56) Sibani, P. & Jensen, H. J.  Record statistics and dynamics.  _Encyclopedia of Complexity and System Science_ (2009).
  * (57) Yang, M., Li, C., Tan, Y., Li, X. & Chen, X.  Impacts of two types of El-Niño on the winter North Pacific storm track.  _Environmental Research Letters_ 15, 094062 (2020).
  * (58) Machado, J. P., Justino, F. & Souza, C. D.  Influence of El Niño-Southern Oscillation on baroclinic instability and storm tracks in the Southern Hemisphere.  _International Journal of Climatology_ 41, E93–E109 (2021).
  * (59) Wang, B. _et al._ Historical change of El Niño properties sheds light on future changes of extreme El Niño.  _Proceedings of the National Academy of Sciences_ 116, 22512–22517 (2019).
  * (60) Cai, W. _et al._ Increasing frequency of extreme El Niño events due to greenhouse warming.  _Nature Climate Change_ 4, 111–116 (2014).
  * (61) Eyring, V. _et al._ Overview of the coupled model intercomparison project phase 6 (cmip6) experimental design and organization.  _Geoscientific Model Development_ 9, 1937–1958 (2016).
  * (62) Györgyi, G., Moloney, N., Ozogány, K., Rácz, Z. & Droz, M.  Renormalization-group theory for finite-size scaling in extreme statistics.  _Physical Review E_ 81, 041135 (2010).
  * (63) Calvo, I., Cuchí, J. C., Esteve, J. G. & Falceto, F.  Extreme-value distributions and renormalization group.  _Physical Review E_ 86, 041109 (2012).
  * (64) Corral, A.  Scaling in the timing of extreme events.  _Chaos, Solitons & Fractals_ 74, 99–112 (2015).