# Detecting Extreme Temperature Events Using Gaussian Mixture Models

###### Abstract

Extreme temperature events have traditionally been detected assuming a
unimodal distribution of temperature data. We found that surface temperature
data can be described more accurately with a multimodal rather than a unimodal
distribution. Here, we applied Gaussian Mixture Models (GMM) to daily near-
surface maximum air temperature data from the historical and future Coupled
Model Intercomparison Project Phase 6 (CMIP6) simulations for 46 land regions
defined by the Intergovernmental Panel on Climate Change (IPCC). Using the
multimodal distribution, we found that temperature extremes, defined based on
daily data in the warmest mode of the GMM distributions, are getting more
frequent in all regions. Globally, a 10-year extreme temperature event
relative to 1985-2014 conditions will occur 13.6 times more frequently in the
future under 3.0∘C of Global Warming Levels (GWL). The frequency increase can
be even higher in tropical regions, such that 10-year extreme temperature
events will occur almost twice a week. Additionally, we analysed the change in
future temperature distributions under different GWL and found that the hot
temperatures are increasing faster than cold temperatures in low latitudes,
while the cold temperatures are increasing faster than the hot temperatures in
high latitudes. The smallest changes in temperature distribution can be found
in tropical regions, where the annual temperature range is small. Our method
captures the differences in geographical regions and shows that the frequency
of extreme events will be even higher than reported in previous studies.

JGR: Atmospheres

Deutsches Zentrum für Luft- und Raumfahrt (DLR), Institut für Physik der
Atmosphäre, Oberpfaffenhofen, Germany University of Bremen, Institute of
Environmental Physics (IUP), Bremen, Germany Computational Research Division,
Lawrence Berkeley National Laboratory, Berkeley, CA, USA Center for Climate
Change and Policy Studies, Boğaziçi University, Istanbul, Turkey Aytaç
<EMAIL_ADDRESS>
Extreme temperature events are detected with Gaussian Mixture Models to follow
a multimodal rather than a unimodal distribution.

10-year temperature extremes will occur 13.6 times more frequently under 3.0∘C
future warming.

Colder days are getting warmer faster than hotter days in high latitudes,
whereas it is the opposite for many regions in low latitudes.

## Plain Language Summary

Extreme temperature events are unusual weather conditions with exceptionally
low or high temperatures. Traditionally, the temperature range was determined
by assuming a single distribution, which describes the frequency of
temperatures at a given climate using their mean and variability. This single
distribution was then used to detect extreme weather events. In this study, we
found that temperature data from reanalyses and climate models can be more
accurately described using a mixture of multiple Gaussian distributions. We
used the information from this mixture of Gaussians to determine the cold and
hot extremes of the distributions. We analysed their change in a future
climate and found that hot temperature extremes are getting more frequent in
all analyzed regions at a rate that is even higher than found in previous
studies. For example, a global 10-year event will occur 13.6 times more
frequently under 3.0∘C of global warming. Furthermore, our results show that
the temperatures of hot days will increase faster than the temperature of cold
days in equatorial regions, while the opposite will occur in polar regions.
Extreme hot temperatures will be the new normal in highly populated regions
such as the Mediterranean basin.

##  1 Introduction

Increasing levels of atmospheric carbon dioxide (CO2) concentration
unequivocally transformed the earth’s climate [IPCC (2021)]. This surplus of
CO2 in the atmosphere contributes to the greenhouse effect, and by increasing
the mean and the variability of global temperatures, it amplifies the risk of
high-impact temperature extremes [Baker . (2018)]. The effects of
anthropogenic global warming led to the emergence of heat extremes that would
not have occurred previously [Robinson . (2021)]. This means that
unprecedented heat extremes like the 2010 Russian heatwave or the 2021 Western
North America heatwave would have likely not happened without the warming
effect [Rahmstorf Coumou (2011), Christidis . (2015), Thompson . (2022)]. The
latter was found to be a remarkable four standard deviations away from the
mean [Thompson . (2022)]. The Intergovernmental Panel on Climate Change (IPCC)
Sixth Assessment Report (AR6) concluded that human influence on the climate
system is unequivocal [Eyring . (2021)] and virtually certain to be the main
driver of the changes in hot and cold extremes [Seneviratne . (2021)]. It
introduced more frequent and intense hot extremes since the 1950s on land
areas while a decrease in cold extremes is observed [IPCC (2021)]. Several
studies found that the duration, frequency, and intensity of extreme events
will increase, and extreme events will be introduced at new locations
[Seneviratne . (2012), Rahmstorf Coumou (2011), Kharin . (2013), Sillmann,
Kharin, Zhang . (2013), Sillmann, Kharin, Zwiers . (2013), Pfleiderer .
(2019), Perkins-Kirkpatrick Lewis (2020), Vogel, Hauser Seneviratne (2020),
Raymond . (2020), Seneviratne . (2021), Mallick . (2022)]. As the number of
occurrences of heat extremes like the 2003 European heatwave and their
duration increase, the socio-economic burden of climate change poses a threat
to societies [Meehl Tebaldi (2004), Robine . (2008), García-León . (2021),
Demiroglu . (2020), Perera . (2020), Seneviratne . (2021)].

The warming of the climate causes different changes in different regions.
Tropics, polar regions and the Middle East and North Africa (MENA) region, are
hot spots of notable climate trend shifts [Hao . (2018), Y. Zhang . (2022)].
Iyakaremye2022 have shown that an abrupt shift in the daily maximum
temperatures occurred in Africa in the last two decades compared to the
previous 20 years, which introduced more frequent and intense hot days.
Moreover, regions in Africa will face a higher increase in temperatures
compared to the rest of the globe. Iyakaremye2021 found that the annual
maximum of daily maximum temperatures over Africa is expected to increase by
1.6/2.2∘C in the future, while global temperatures are projected to rise by
1.5/2.0∘C during the same period. In the MENA region, the frequency and
intensity of heatwaves will highly increase by the end of the century under a
business-as-usual pathway scenario, which will affect about half of the MENA
population [Lelieveld . (2016), Zittis . (2021), Ozturk . (2021)]. The number
of occurrences of exceptionally hot summers, which have 2-4∘C hotter
temperatures than the long-term average, has also increased from a single
event between 1951 and 1980 to five events between 2001 and 2010 in Central
and Eastern Europe, where the 2010 heatwave was the hottest and longest event
with the largest geographical extent that ever occurred over Europe [Twardosz
Kossowska-Cezak (2013), Guerreiro . (2018)]. Similarly, other studies also
found that the temperature extremes in Europe will increase 20-fold at the end
of the century, compared to 1961-1990 [Nikulin . (2011), Schär . (2004),
Barriopedro . (2011)]. Over the Americas, the dry and hot extremes showed an
increase both in frequency and spatial scope over the past 122 years [Alizadeh
. (2020), Cai . (2014)].

Correctly characterizing the temperature distributions to analyze extreme
events is a still-continuing issue as extremes are by definition rare events,
and several studies showed that the assumption of distributions or a
stationary climate often underestimates the observed heat records [Benestad
(2004), Schär . (2004), Anderson Kostinski (2010), Fischer Schär (2010),
Barriopedro . (2011), C. Li . (2019), Loikith Neelin (2019)]. Thompson2022
characterized extreme events by calculating a daily extreme index which is the
difference between the daily maximum temperature and mean daily maximum
temperature divided by the standard deviation. With the assumption of a normal
distribution, they found that the 2021 North American heatwave was one of the
most extreme events with 4 standard deviations from the mean. Moreover, the
authors projected that 20% of the weather risk attribution forecast regions
[Stone (2019)] will experience extreme events that are four standard
deviations from the means in the future. Other studies found that hot summers
will be the norm, i.e. mean temperatures exceed the temperature of the
historically hottest summer, within the next 1-2 decades [Mueller . (2016),
Lewis . (2017), Vogel, Hauser Seneviratne (2020), Vogel, Zscheischler .
(2020)].

Common indices to monitor and analyze climate extremes that are used in the
climate community at the moment, such as ETCCDI (the Expert Team on Climate
Change Detection and Indices), are mostly based on daily mean near-surface air
temperature or daily maximum near-surface air temperature [X. Zhang . (2011),
Alexander . (2006)]. Two standard approaches to detect extreme events are the
percentile-over-threshold (POT) and the block maxima method. The block maxima
method groups data into an equal length of blocks, e.g. month, season, or
year, and use the maximum temperature value of each block to fit the data. The
POT method defines a threshold, e.g. percentiles, and uses all temperature
values above this threshold in the analysis. Choosing the percentiles for
defining extremes is not trivial as the temperature extremes have a strong
seasonality and temporal dependence [Huang . (2016)]. The block maxima method
is more commonly used in climate studies because of its simplicity with
monthly, seasonal or annual block periods for fitting generalized extreme
value (GEV) distribution to temperature and precipitation extremes [Kharin .
(2013), Wang . (2016), Paciorek . (2018), Wehner . (2018), Ben Alaya . (2020),
C. Li . (2021), IPCC (2021)]. The block maxima method, however, does not use
all available data, as calculating a single maximum value from a block period
throws out the rest of the data. To be approximated by the GEV distribution,
the blocks are assumed to be long enough and “max-stable”, which means that if
you take the maximum of a group of values selected from a specific GEV
distribution, the result will be GEV distributed with the same shape parameter
[Huang . (2016), Ben Alaya . (2020)]. However, these assumptions might not be
valid for all possible use cases or all possible variables. For example, GEV
is not the best fit for shorter block lengths as the fit improves with
increasing block size [Ben Alaya . (2020), Wang . (2016)]. Alaya2020 argued
that the identically distributed random variables assumption of extreme value
theory might be problematic for extreme precipitation events. They considered
a mixture of GEV distributions to fit precipitation data to demonstrate that
the mixture distribution could be a potential explanation for the instability
of annual maxima. Kollu2012 tested wind speed characteristics using mixture
probability distribution functions (PDF). They found that conventional PDFs
are inadequate to describe wind speed distributions compared to the mixture
distributions that they used in the study. A mixture of Gaussians was used by
Shin2022 to describe the distribution of the daily thermal comfort index in
South Korea, an index that has a strong seasonality. Ice surface temperature
data follows a clear multimodal distribution, according to Clarkson2022. They
also found that a unimodal distribution fit is particularly poor at modelling
the tail probabilities. Probability distributions with one and two components
are called unimodal and bimodal, respectively, whereas distributions with
multiple (two or more) components are called multimodal distributions.

The temperature distributions are expected to move towards warmer temperatures
and to change their shape with changing means and standard deviations [IPCC
(2021)]. Also, the assumption of distribution might not be correct for all
geographical regions as daily weather variables show a distinct non-
Gaussianity [EM. Volodin Yurova (2012), Perron Sura (2013), Kodra Ganguly
(2014), Sardeshmukh . (2015), Linz . (2018), Tamarin-Brodsky . (2019)].
Furthermore, several studies found that daily mean, daily maximum and real
forecast data of 2m temperatures show bimodal features [Grace (1995), Wilks
(2002), Donat Alexander (2012), Cho Jeong (2016), Bertossa . (2021)]. These
changes, shifts and bimodalities in the temperature distributions affect the
probabilities in the tails. As extreme events are rare events that lie in the
tails of a distribution, correctly describing the tails is very important for
extreme event detection. Even though the block maxima method is widely used in
studies which used block sizes large enough to converge asymptotically to GEV
distributions, a GEV distribution is not well suited to describe extreme value
data when the bimodality is apparent or block sizes are short [Sardeshmukh .
(2015), Wang . (2016), Knoben . (2019), Ben Alaya . (2020)]. Therefore, the
properties of the entire probability distribution, i.e. mean, standard
deviation and shape, are needed to get the tail properties right [Sardeshmukh
. (2015)]. A distribution can be described by not only the mean and the
standard deviation but also skewness and kurtosis. Donat2012 found that daily
minimum and maximum temperatures have significantly shifted towards higher
values and skewed towards the hotter part of the distribution. They
highlighted that the changes in extremes are related not only to the means but
also to other parameters of the daily temperature distribution.
Sardeshmukh2009 found a parabolic relationship between kurtosis and skewness
that cause the non-Gaussianity of the observed daily weather anomalies.
Similarly, Tamarin2022 used a mixture model with three Gaussians to describe
the PDF of near-surface atmospheric temperature to analyze the relationship
between kurtosis and skewness, as they are important to explain how the tails
of the distribution change. They found that two- and three-Gaussian models are
useful to explain the relationship between kurtosis and skewness.

In the study presented here, our approach is to utilize the entire temperature
distribution to detect extreme events. We implemented Gaussian Mixture Models
(GMM), which describe the probability distribution function of data points as
a mixture of Gaussian distributions. We determined the number of Gaussian
components in the temperature distribution of each grid cell of 46 land
regions defined by the Intergovernmental Panel on Climate Change (IPCC) using
daily near-surface maximum air temperature data from the historical and future
Coupled Model Intercomparison Project Phase 6 (CMIP6) simulations. This choice
was supported by previous studies which found distinct bimodality in daily
weather variables [Grace (1995), Wilks (2002), Donat Alexander (2012), Cho
Jeong (2016), Bertossa . (2021)] and was verified by applying the same
analysis to the European Centre for Medium-Range Weather Forecasts Reanalysis
5th Generation (ECMWF-ERA5) data for the same historical time period
(1985-2014). The parameters from the determined distribution components,
namely mean, standard deviation and weight, were used to calculate the change
in the return period of extreme temperature events between the historical and
future periods determined by using global warming levels (GWL). In a
stationary climate, the return period of an event describes the average time
between the occurrences of a certain event of a defined size. In this study,
we analysed 1-year, 5-year, 10-year and 20-year events, where an n-year event
has an occurrence probability of 1/n as the climate is not stationary. Hence,
these event magnitudes change as time progresses, where an n-year event means
that the event in question would be expected to occur once in every n years.
We only calculated return periods equal to or less than the available future
data period to prevent overestimating the return periods of extreme events,
since GMM distributions are not bounded. Section 2 presents the climate data
and warming levels used in this study, as well as the analyzed regions, and
explains the methodology of detecting extreme event return periods by using
GMM. Section 3 shows our results obtained using the GMM method for all
analyzed IPCC land regions, and section 4 finalizes the paper with a summary
and discussion.

##  2 Data and Methodology

###  2.1 Climate Data

For this study, we used multi-year daily near-surface maximum temperatures
from the Coupled Model Intercomparison Project Phase 6 (CMIP6), and for which
both the historical simulations and the simulations for Shared Socioeconomic
Pathways (SSPs) 1-2.6, 2-4.5, 3-7.0 and 5-8.5 scenarios were available
[O’Neill . (2014), Eyring . (2016), O’Neill . (2016)]. Additionally, the
ECMWF-ERA5 dataset was included for the 30-year time period (1985-2014)
[Hersbach . (2018)]. Table 1 shows the list of models and their resolutions.
The 30-year time period from 1985 to 2014 from historical simulations is used
as the base to calculate the return values of extreme temperature events, i.e.
1-year, 5-year, 10-year and 20-year events. The GWL, as introduced in the IPCC
AR6 report, are used to assess the changes in future climate in line with the
warming levels defined in the Paris Agreement which are compared to the pre-
industrial period [IPCC (2021)]. The future period for each model is defined
as a 20-year period between 2015 and 2100 when the central year of the running
window of the global daily near-surface temperature mean of that model first
exceeds 1.5∘C, 2∘C, 3∘C, and 4∘C relative to 1850-1900 global daily near-
surface mean temperatures. We used the same GWL periods defined for and used
in IPCC [IPCC (2021), Hauser . (2022)], similarly to hajat2022 and
Ribeiro2022. Therefore, we obtained the start and end years of 20-year GWL
periods for each CMIP6 simulation from gwl_periods. Here, we used a longer
historical base period (30 years) compared to future GWL periods (20 years)
for the analysis to obtain more robust results. This decision was made based
on the fact that GMM distributions have no bounds. Therefore, we focused our
analysis solely on return periods shorter than our base period. By limiting
our analysis to shorter return periods, we can mitigate the biases and
outliers that may occur beyond the limits of the datasets. As some datasets
did not exceed certain warming levels, they were excluded from the analysis
(e.g NOR-ESM2-MM was not used in calculations for 4∘C warming under SPP5-8.5,
as it did not exceed this level). Figure 1 shows the historical and future GWL
periods for each CMIP6 model used in this study.

Table 1: Reanalysis data and CMIP6 models used in this study to detect extreme temperature events. Climate models with spatial resolutions ranging from 50 to 500 km were used in the analyses. The first available ensemble members were chosen. The Reanalysis dataset that has a resolution of 25km was regridded to 100 km and used for evaluating modality. Model | Variant | Resolution | Reference |
---|---|---|---|---
ECMWF-ERA5 | Reanalysis | 25 km | [Hersbach . (2018)] |
ACCESS-CM2 | r1i1p1f1 | 250 km | [Dix . (2019)] |
ACCESS-ESM1-5 | r1i1p1f1 | 250 km | [Ziehn . (2019)] |
AWI-CM-1-1-MR | r1i1p1f1 | 100 km | [Semmler . (2018)] |
BCC-CSM2-MR | r1i1p1f1 | 100 km | [Wu . (2018)] |
CanESM5 | r1i1p1f1 | 500 km | [Swart . (2019)] |
CNRM-CM6-1 | r1i1p1f2 | 250 km | [Voldoire (2018)] |
CNRM-CM6-1-HR | r1i1p1f2 | 50 km | [Voldoire (2019)] |
CNRM-ESM2-1 | r1i1p1f2 | 250 km | [Seferian (2018)] |
EC-Earth3 | r1i1p1f1 | 100 km | [EC-Earth Consortium (EC-Earth) (20191)] |
EC-Earth3-CC | r1i1p1f1 | 100 km | [EC-Earth Consortium (EC-Earth) (2021)] |
EC-Earth3-Veg | r1i1p1f1 | 100 km | [EC-Earth Consortium (EC-Earth) (20192)] |
EC-Earth3-Veg-LR | r1i1p1f1 | 250 km | [EC-Earth Consortium (EC-Earth) (2020)] |
FGOALS-g3 | r1i1p1f1 | 250 km | [L. Li (2019)] |
GFDL-ESM4 | r1i1p1f1 | 100 km | [Krasting . (2018)] |
HadGEM3-GC31-LL | r1i1p1f3 | 250 km | [Ridley . (20191)] |
HadGEM3-GC31-MM | r1i1p1f3 | 100 km | [Ridley . (20192)] |
INM-CM4-8 | r1i1p1f1 | 100 km | [von . (2019)] |
INM-CM5-0 | r1i1p1f1 | 100 km | [E. Volodin . (2019)] |
IPSL-CM6A-LR | r1i1p1f1 | 250 km | [Boucher . (2018)] |
KACE-1-0-G | r1i1p1f1 | 250 km | [Byun . (2019)] |
MIROC6 | r1i1p1f1 | 250 km | [Tatebe Watanabe (2018)] |
MIROC-ES2L | r1i1p1f2 | 500 km | [Hajima . (2019)] |
MPI-ESM1-2-HR | r1i1p1f1 | 100 km | [Jungclaus . (2019)] |
MPI-ESM1-2-LR | r1i1p1f1 | 250 km | [Wieners . (2019)] |
MRI-ESM2-0 | r1i1p1f1 | 100 km | [Yukimoto . (2019)] |
NESM3 | r1i1p1f1 | 250 km | [Cao Wang (2019)] |
NorESM2-LM | r1i1p1f1 | 250 km | [Seland . (2019)] |
NorESM2-MM | r1i1p1f1 | 100 km | [Bentsen . (2019)] |
UKESM1-0-LL | r1i1p1f2 | 250 km | [Tang . (2019)] |

We extracted daily maximum near-surface air temperature for 30-year historical
and 20-year future periods under GWL for each SSP individually for 46 IPCC
land regions that are shown in Figure 1 [Iturbide . (2020)]. All data
extraction and preprocessing in this study were performed by using the Earth
System Model Evaluation Tool (ESMValTool) version 2.5.0, which is an open-
source software package for analysing and evaluating model simulations [Eyring
. (2020), Lauer . (2020), Righi . (2020), Weigel . (2021)]. We extracted the
daily maximum near-surface air temperature from each model for each region
using shapefiles provided by IPCC [Iturbide . (2020)], converted units from
Kelvin to Celsius, and created a single spatiotemporal Network Common Data
Form (NetCDF) file for each region. The data were then ready to be used in the
diagnostic script written in Python where the extreme events and their return
periods were analyzed.

Figure 1: (Top) We used 46 land regions defined in Iturbide2020. See
Supplementary Material Table S2 for region definitions. (Bottom) Future
periods of the CMIP6 models when the central year of the 20-year running
window exceeds global warming levels relative to 1850-1900 base for the
SSP5-8.5 scenario are extracted using the data from gwl_periods. The colours
in the graph go from light to dark, each colour representing a different level
of warming 1.5∘C, 2∘C, 3∘C, and 4∘C. These levels are expected to be exceeded
around 2026, 2040, 2060, and 2070 respectively. The 30-year historical base
period is indicated in grey. Note that different models have different time
periods when they exceed the GWL. Future periods for other SSP scenarios are
presented in the Supplementary Material Table S3.

###  2.2 Return Period Analyses

For the return period calculation of extreme temperature events, i.e. 1-year,
5-year, 10-year and 20-year events, we defined a temperature threshold for an
event by calculating the standard deviation distance of the event temperature
from the mean temperature in the past, i.e. how many standard deviations away
the event temperature was from the mean. We then applied this temperature
threshold value to the future period but calculated its standard deviation
distance from the mean using the parameters from the future distribution, i.e.
how many standard deviations away the event temperature will be from the mean.
To test the underlying distribution shape of the daily near-surface maximum
temperature distribution, we first analyzed data from individual grid cells of
each climate model. We found that daily maximum near-surface air temperature
data in climate grid cells usually do not follow a unimodal distribution, but
rather follow a bimodal distribution, a probability distribution composed of
two components.

To calculate the return periods of extreme events, we modelled the probability
distribution of multi-year (30 years for the historical base period and 20
years for future GWL periods) daily near-surface maximum temperature data from
a grid cell as mixtures of multiple Gaussian distributions, rather than a
single Gaussian distribution. GMM is a probabilistic model that describes the
data points in a population as a mixture of Gaussian distributions with
unknown parameters which are the mean, standard deviation and weight of each
Gaussian component, five parameters in total for a bimodal distribution. With
this approach, we were able to analyse the change in the distribution, and
accurately model the tails of the data compared to a unimodal distribution.
When a unimodal distribution is fit to multi-year data with bimodality, it is
likely that the resulting distribution will have a larger standard deviation
to encompass both modes. This large standard deviation between a unimodal
distribution fit to bimodal data can have significant implications for
analyses, as the larger standard deviations of unimodal distributions tend to
push the extreme events further away from what would be observed if the
bimodality were properly accounted for. In other words, when one tries to
calculate the threshold of an event as n sigma distance from the mean, this
threshold might be well beyond the maximum value of the distribution.
Furthermore, a unimodal distribution fit will affect the measures of central
tendency when bimodality exists in data [von Hippel (2005)]. An example
goodness-of-fit test for normal distribution, GEV distribution with different
shape parameters and GMM distributions on the daily maximum temperature data
from a random grid cell is presented in Supplementary Material Section 1. We
used an unsupervised machine-learning package, the “GaussianMixture” package
from open-sourced machine-learning library Scikit-learn, to compute the
unknown parameters of the Gaussian components in a mixture that generates all
observed data points [Pedregosa . (2011)]. We applied this package to the
daily maximum near-surface air temperature data in each grid cell of the CMIP6
models. The “GaussianMixture” package first randomly assigns values to the
component parameters and then uses the expectation-maximisation algorithm (EM)
to converge their values. EM algorithm fits GMM to data by alternating between
two steps, Expectation (E) and Maximisation (M). In the E step, it randomly
assumes components and calculates the probability of each point to be
generated by that component. In the M step, the parameters are tweaked to
maximise the likelihood found in the first step. It also uses the Bayesian
Information Criteria (BIC) score, which is used to estimate the goodness-of-
fit of a distribution and which accounts for both the likelihood function and
the number of parameters. Then, the probability distribution function of the
mixture model that was fit to multi-year daily near-surface temperature can be
written as a linear summation of multiple Gaussian components:

$\displaystyle p(x)$
$\displaystyle=\sum_{k=1}^{K}\omega_{k}\mathcal{N}(x\mid\mu_{k},\sigma_{k})$
(1) $\displaystyle\mathcal{N}(x\mid\mu_{k},\sigma_{k})$
$\displaystyle=\frac{1}{\sigma_{k}{\sqrt{2\pi}}}\exp{\left(-\frac{(x-\mu_{k})^{2}}{2\sigma_{k}^{2}}\right)}$
(2) $\displaystyle\sum_{k=1}^{K}\displaystyle\omega_{k}$ $\displaystyle=1$ (3)

where $K$ is the number of Gaussian components in the mixture. $\mu_{k}$,
$\sigma_{k}$ and $\omega_{k}$ are the mean, the standard deviation and the
weight of the $k^{th}$ component, respectively. Implementing Gaussian mixture
models to evaluate multi-year raw daily maximum temperatures allows us to
investigate the long-term characteristics of the individual components. This
method does not consider the temporal changes within one period, as they can
be assumed to be negligible compared to the changes between different time
periods. As shown in other studies, mean temperatures are increasing all over
the globe [IPCC (2021), Robinson . (2021), Eyring . (2020)]. Using the raw
temperatures, we can analyse how the convergence or divergence of the peaks of
the different Gaussian components affect the extremes compared to the used
historical periods. In our analysis, we have disregarded three or more
Gaussian components. This choice was supported by the value of the BIC score
and the fact that increasing the number of components tends to cause
overfitting, even though BIC scores penalise adding more parameters. In some
cases, the BIC scores for the components showed close results for more than
three components (see Figure S2 in the Supplementary Material). For instance,
the lowest BIC score was reached for a mixture with seven Gaussian components
for the distribution of temperatures in a grid cell. However, the highest
change in BIC scores occurred when switching from one component to two
components. Consequently, we used the gradient of BIC scores rather than using
the lowest score. We selected the number of Gaussian components where the
highest gradient change occurs in the BIC scores as the best fit. To further
prevent overfitting, we also applied the following unimodality test after
estimating the BIC scores: If the BIC score returned a bimodal distribution,
then the parameters of the mixture distribution components were used for the
unimodality test. As shown in Equation 4, if the difference between the means
of Gaussian components was less than or equal to twice the minimum of standard
deviations, then unimodal distribution was assumed, otherwise, the bimodal
distribution fit for the data was kept. It is worth noting that this procedure
had a tendency to favour fitting a unimodal distribution. However, after all
these tests and checks, the majority of grid cells showed a clear bimodal
distribution. For a bimodal distribution, hereafter we referred to the right
(left) Gaussian component as “hot (cold) Gaussian” as shown in Figure 2).

$|\mu_{1}-\mu_{2}|\leq 2\min(\sigma_{1},\sigma_{2})$ (4)

Figure 2: Exemplary bimodal distribution of daily maximum temperatures from a
grid cell for the historical 30-year period of 1985-2014 (blue) and future
20-year GWL period (red). Blue and red lines show the corresponding GMM fit
for the historical and future periods, respectively. The shape of the
distribution is determined by the parameters of each Gaussian component, which
are the means, standard deviations and weights. Here, the means of cold and
hot Gaussian peaks are shown with blue(red) dots and squares for the
historical (future) period, respectively. The hot Gaussian component used in
the analysis is shown with a dashed blue (red) line for the historical
(future) period. The bottom two plots show what convergence and divergence of
the peaks mean based on the $\Delta T$ value.

First, we grouped grid cells of a region depending on their modality, either
unimodal or bimodal, for each CMIP6 model, and calculated the percentages of
grid modalities among all grid cells of a region for each CMIP6 model. We then
determined the multi-model mean percentages of grid cell modalities of a
region as shown in Figure 3. Additionally, we calculated the global multi-
model mean percentage of grid cell modalities using all regions and CMIP6
models. We found that globally 88.78% of all grid cells follow a bimodal
distribution in the historical period as shown in the white box in the upper
centre part of Figure 3. Furthermore, we analysed the ECMWF-ERA5 dataset for
the same historical time period (1985-2014) to confirm whether bimodality is
also found in data other than model simulations. We regridded the ECMWF-ERA5
data from a 25-km grid to a coarser 1-degree 100-km grid using the nearest
neighbour method to have a similar resolution as many CMIP6 datasets. The
ECMWF-ERA5 reanalysis dataset shows similar results to the CMIP6 models:
Globally 86.95% of all grid cells in the ECMWF-ERA5 reanalysis dataset follow
a bimodal distribution as shown in the white box in the upper centre part of
Figure 4, while only 13.05% of them follow a unimodal distribution. The
temperature distributions in the ERA5 and CMIP6 datasets show predominantly
similar patterns across various regions, although certain exceptions are
observed, particularly in South America. These differences can most likely be
attributed to several factors. First, resampling of ERA5 data from a 25 km
grid to a coarser 1-degree grid introduces a smoothing effect on the data,
which would increase the unimodal grid cells. Additionally, biases in surface
temperature in CMIP6 datasets also contribute to the observed variations from
ERA5 [Bock . (2020)]. Nevertheless, as we aim to evaluate the shape of
temperature distributions, we did not apply a bias correction and used raw
multi-year daily temperature data from CMIP6 models for our analysis.

Figure 3: Multi-model mean percentages of grid modalities for the historical
period in study regions grouped by continents. Dark and light blue bars show
the percentage of grid cells with unimodal or bimodal distribution,
respectively, for the historical period of 29 CMIP6 simulations. Figure 4:
Same as Figure 3 but for ECMWF-ERA5 reanalysis dataset.

Then, the parameters of the hot Gaussian component, $\mu_{hot}^{historical}$,
$\sigma_{hot}^{historical}$ and $\omega_{hot}^{historical}$, were used to
calculate the change in return periods. We only analysed 1-year, 5-year,
10-year and 20-year events, as GMM are unbounded. One should be careful while
calculating the return periods using GMM, as the unbounded tails of the
Gaussian component could overestimate the probabilities of longer return
periods. Therefore, return periods equal to or less than the analysis period
were calculated using GMM. The change in return periods is calculated first in
each grid cell of a region and then averaged together to produce regional
results for each CMIP6 simulation.

For normally distributed data, the expected percentage of the population
inside the $\mu\pm d\sigma$ range is defined as

$E(\mu\pm d\sigma)=\operatorname{erf}\bigg{(}\frac{d}{\sqrt{2}}\bigg{)}$ (5)

where $\operatorname{erf}$ is the error function and $d$ is the standard
deviation distance. The approximate expected frequency, $f$, outside this
range is then defined as the return period of an extreme.

$1\ in\
\frac{1}{1-\displaystyle\operatorname{erf}\bigg{(}\frac{d}{\sqrt{2}}\bigg{)}}\
days$ (6)

The return period of an event describes the average time between the
occurrences of a certain event of a defined size, where an $n-year$ event has
an occurrence probability of 1/n as the climate is not stationary, where
“year” is defined as the number of days covered by the hot Gaussian component.
The reason for this definition is that the entire probability distribution is
composed of both the cold and warm periods of a year, however, our dataset
consists of daily maximum temperature data spanning 30 (or 20) years,
totalling 10,950 (or 7,300) days. Since our analysis specifically aims to
identify extreme values using parameters from the hot Gaussian component, we
need to consider the number of data points generated by this component as the
definition of a ”year”. We determine the length of a ”year” by dividing the
data points falling under the hot Gaussian component by the length of the
analysis period as shown in Equation 7. For example, we can assume that a
symmetrical bimodal distribution results in $\sim$180 days of cold weather and
$\sim$180 days of hot weather in a normal 365-day calendar year. For such a
symmetric case, a 10-year event would then be a temperature event in 1800 days
(10 years$\textstyle\times 180\textstyle\frac{days}{year}$). Since we cannot
assume a symmetric distribution for grid cells of each model, we calculated
the number of days covered by the hot Gaussian component using the component
weights and dataset size.

Let $\mathcal{D}$ denote the number of days in $L$ years. Then, a “year” in
the historical period,
$\mid\mathcal{N}(\mu_{hot}^{historical},\sigma_{hot}^{historical})\mid$ is
defined as

$\mid\mathcal{N}(\mu_{hot}^{historical},\sigma_{hot}^{historical})\mid=\frac{\omega_{hot}^{historical}\mathcal{D}}{L}$
(7)

where $\mu_{hot}^{historical}$ is the mean, $\sigma_{hot}^{historical}$ is the
standard deviation and $\omega_{hot}^{historical}$ is the weight of hot
Gaussian component. The expected frequency of n-year events in the historical
period, $f^{historical}_{n}$, is then calculated by using the length of a
year,

$f^{historical}_{n}=n\times\mid\mathcal{N}(\mu_{hot}^{historical},\sigma_{hot}^{historical})\mid\quad
n=1,5,10,20$ (8)

The standard deviation distance of range, $d^{historical}_{n}$, for an extreme
event in the historical period can be calculated by using Equation 6,

$d^{historical}_{n}=\displaystyle\operatorname{erf}^{-1}\bigg{(}1-\frac{1}{f^{historical}_{n}}\bigg{)}\sqrt{2}$
(9)

where $\operatorname{erf}^{-1}$ is inverse error function. Now, we can
calculate a temperature threshold, $\tau^{historical}_{n}$, for an n-year
event in the historical period.

$\tau^{historical}_{n}=\mu_{hot}^{historical}+d^{historical}_{n}\sigma_{hot}^{historical}$
(10)

Using this temperature threshold from the historical period, we calculate the
standard deviation distance of the temperature threshold of n-year event in
the future, $d^{future}_{n}$, by using the mean $\mu^{future}_{hot}$, and
standard deviation $\sigma^{future}_{hot}$ from the hot Gaussian component of
the future distribution.

$\displaystyle d^{future}_{n}$
$\displaystyle=\frac{\tau^{historical}_{n}-\mu^{future}_{hot}}{\sigma^{future}_{hot}}$
(11) $\displaystyle f_{\dot{n}}^{future}$
$\displaystyle=\frac{1}{1-\operatorname{erf}\bigg{(}\displaystyle\frac{d^{future}_{n}}{\sqrt{2}}\bigg{)}}$
(12)

Finally, the new value of the return period in the future $\dot{n}$, i.e.
$\dot{n}$-year event, is calculated by using Equation 8

$\dot{n}=\frac{f_{\dot{n}}^{future}}{|\mathcal{N}(\mu_{hot}^{future},\sigma_{hot}^{future})|}$
(14)

where $\mathcal{N}(\mu_{hot}^{future},\sigma_{hot}^{future})|$ is length of a
“year” in the future period.

With this method, we can also analyse if and how much the Gaussian components
will shift in the future relative to the historical period. We defined $\Delta
T$, as the change in the difference between the means of cold and hot Gaussian
components as shown in Equation 15:

$\displaystyle\Delta T$ $\displaystyle=\delta T_{cold}-\delta T_{hot}$ (15)
$\displaystyle\delta T_{cold}$
$\displaystyle=\mu_{cold}^{future}-\mu_{cold}^{historical}$ (16)
$\displaystyle\delta T_{hot}$
$\displaystyle=\mu_{hot}^{future}-\mu_{hot}^{historical}$ (17)

In Figure 2, this change in hot and cold Gaussian means is schematically
illustrated. Assuming the future means of Gaussian components are higher than
the historical periods, $\delta T_{cold}$ and $\delta T_{hot}$ will always be
positive. Therefore, a negative $\Delta T$ means that the peaks are diverging
in the future: the hot Gaussian moves toward warmer temperatures faster than
the cold Gaussian, which increases the frequency of hot extremes and induces
an overall warmer climate. A positive $\Delta T$ means that the peaks are
converging: the cold Gaussian moves closer to the hot Gaussian, which
increases the number of days with warmer temperatures in the colder mode.

##  3 Results

First, we checked the change in the percentage of modalities from the present
to the future time periods. For this, we analyzed the modality of the
temperature data from each individual grid cell of an IPCC land region by
counting the number of grid cells with each modality. We found that the
percentages of grid cells with bimodal distributions stay almost the same
under different warming levels. As some of the CMIP6 datasets do not exceed
certain warming levels, the number of datasets are not identical for the
historical and future period and therefore affect the change in percentages.
We analysed modalities of grid cells under different GWL for all SSP scenarios
but we only present SSP5-8.5 results here, as the SSP5-8.5 scenario had data
from 29 CMIP6 models and the GWL is scenario independent. Globally, almost 90%
of all grid cells follow a bimodal distribution as shown in Figure 3 for the
historical period, Figure 4 for the reanalysis data and Figure 5 for GWL 3.0∘C
for different regions grouped by continents (See Supplementary Material Table
S3 for other warming levels). Global averages and the number of datasets are
shown in the white box in the upper centre part of each figure. In the
historical period, the grid cells in tropical and sub-tropical regions have
slightly higher percentages of unimodal distributions compared to higher
latitude regions. However, regions still mostly follow a bimodal distribution
as shown in Figure 3. The multi-model mean percentage of unimodal
distributions does not exceed 50% of grid cells in any of the regions, except
in N.W.South-America (NWS) and South-American-Monsoon (SAM) regions where
51.94% and 50.33% of the grid cells follow a unimodal distribution,
respectively, in the historical period. The higher percentage of unimodal
distributions in lower latitudes is consistent with tropical climate features,
where hot temperatures are observed all year round and the annual temperature
range is small [Richter (2016), Beck . (2018)]. This climate type is therefore
expected to likely experience a temperature distribution close to a single
Gaussian. All grid cells (99.9%) in CMIP6 models follow a bimodal distribution
in the Mediterranean (MED) region in the historical period and under all
future periods. In polar regions, more than 90% of the grid cells follow a
bimodal distribution in the historical period. The percentage of grid cells
with unimodal distributions in polar regions slightly increases under future
global warming levels.

Figure 5: Same as Figure 3 but for future SSP5-8.5 scenario under GWL 3.0∘C.

As previously mentioned in Section 2.2, large values of $\Delta T$ (see
Equation 15) will cause the temperature distribution to change its modality
for future GWL periods with respect to the historical base period of
1985-2014. We analysed all regional grids for all CMIP6 models for the
modality changes under GWL 1.5∘C, 2∘C, 3∘C, and 4∘C. Figure 6 shows the
percentage of changes in grid cell distribution modalities under GWL3.0∘C.
Globally, the percentage of grids changing from a unimodal (bimodal)
distribution in the historical period to a bimodal (unimodal) distribution in
the future periods is between 2.79% (2.26%) and 6.02% (3.88%) for different
scenarios and GWL as shown in Table 2. The change from unimodal to bimodal
distribution in the future period is most prevalent in regions where the
highest percentage of unimodality was observed in the historical period, as
shown in Figure 3. This suggests that regions that were previously
characterized by more consistent temperatures (as indicated by a unimodal
temperature distribution) may experience more variability in temperature in
the future. As our analysis uses the mean and standard deviation of the same
component from the historical and future daily maximum temperature
distributions, we only used the grid cells which have the same modality in the
historical and future periods. We disregarded the grid cells with changing
modalities, i.e. unimodal to bimodal or vice versa, as this will affect the
mean and standard deviation, and hence the return period analysis.

Figure 6: Percentage of changes in grid cell modalities relative to 1985-2014
distribution shape for SSP5-8.5 under GWL3.0∘C. Each cell represents a region
of a CMIP6 model and is divided into 4 quadrants. Each quadrant of squares,
$q_{ij}$, uses index notation, where $i$ represents the modality in the
historical period and $j$ represents the modality in the future period, 1 for
a unimodal distribution and 2 for a bimodal distribution. The top-left
quadrant, $q_{11}$, shows the percentage of grid cells with unimodal
distribution both in the historical and the future periods, i.e. unimodal to
unimodal (UU). The top-right quadrant, $q_{12}$, shows the percentage of grid
cells that change from unimodal distribution in the historical period to
bimodal distribution in the future (UB). The bottom-left quadrant, $q_{21}$,
shows bimodal to unimodal (BU). The bottom-right quadrant, $q_{22}$, shows
bimodal to bimodal distribution (BB). The colour of the quadrants shows the
percentage of grid cells. For several models, East Antarctica (EAN) region is
not included in the analysis because it is composed of many grid cells near
the pole, causing numerical problems. Table 2: Global average percentage of
grid cells with varying distribution modality between the historical and
future periods.

Experiment | GWL | Unimodal$\rightarrow$Unimodal | Unimodal$\rightarrow$Bimodal | Bimodal$\rightarrow$Unimodal | Bimodal$\rightarrow$Bimodal
---|---|---|---|---|---
SSP1-2.6 | 1.0 ∘C | 11.01% | 2.79% | 2.26% | 83.94%
SSP1-2.6 | 2.0∘C | 10.31% | 3.53% | 2.45% | 83.71%
SSP2-4.5 | 1.5∘C | 11.02% | 2.78% | 2.26% | 83.95%
SSP2-4.5 | 2.0∘C | 10.24% | 3.56% | 2.70% | 83.50%
SSP2-4.5 | 3.0∘C | 8.79% | 4.71% | 3.08% | 83.42%
SSP2-4.5 | 4.0∘C | 7.21% | 6.02% | 3.42% | 83.35%
SSP3-7.0 | 1.5∘C | 10.82% | 2.95% | 2.40% | 83.83%
SSP3-7.0 | 2.0∘C | 10.15% | 3.62% | 2.89% | 83.34%
SSP3-7.0 | 3.0∘C | 8.92% | 4.68% | 3.44% | 82.96%
SSP3-7.0 | 4.0∘C | 7.80% | 5.50% | 3.81% | 82.89%
SSP5-8.5 | 1.5∘C | 11.05% | 2.85% | 2.31% | 83.78%
SSP5-8.5 | 2.0∘C | 10.32% | 3.58% | 3.04% | 83.06%
SSP5-8.5 | 3.0∘C | 9.14% | 4.76% | 3.78% | 82.32%
SSP5-8.5 | 4.0∘C | 8.21% | 5.47% | 3.88% | 82.45%

We also analysed the movements of the Gaussian components relative to each
other using the $\Delta T$ definition from Equation 15 in grid cells with a
bimodal distribution. Figure 7 shows the $\Delta T$ results for all analysed
regions for SSP5-8.5 under 3.0∘C warming (see Supplementary Material Figure S7
to S12 for other warming levels). Changes in distribution peaks are smaller
for the lower warming levels. This is consistent with the fact that the time
periods for exceeding warming levels are very close to the historical period
as shown in Figure 1. For the future 3.0∘C warming scenario, we observed that
the mean temperatures are increasing in all regions. Temperature distributions
for the European regions have negative $\Delta T$ values, -0.42 degrees on
average. This will cause already bimodal peaks in the historical period to
separate further from each other in the future, while the whole distribution
moves towards higher temperatures. Divergence of peaks will result in more
extreme hot temperatures in Europe, as the hot Gaussian moves faster. This
result is in agreement with findings from the IPCC AR6 report, in which
temperatures in Europe are reported to increase faster than the rest of the
globe [IPCC (2021)]. Polar regions, Northern America and parts of Northern
Asia have positive $\Delta T$ values, i.e. converging peaks in grid cells with
bimodal distributions. The distribution shape shifts to warmer temperatures
and approaches a unimodal distribution as the cold Gaussian part of the
distribution moves toward the warmer temperatures faster than the hot Gaussian
part. This convergence is also consistent with the slight increase in the
percentage of unimodal distribution in polar regions as shown in Figure 5.
This will cause polar regions to have more days with warmer temperatures also
in the colder mode while also having an overall warmer climate. The
convergence of peaks in three polar regions (EAN, WAN, GIC) and three northern
regions (RAR, NEN and NWN) becomes clear when the regions are sorted by the
mean temperature of cold Gaussian component as shown in Figure 8. High $\Delta
T$ values in polar regions are also supported by previous studies reporting
that Arctic regions are warming faster than the global average [Taylor .
(2022)]. The lowest $\Delta T$ values are in MED and SAM regions, -0.90 and
-1.21 degrees respectively, which will cause both bimodal peaks to diverge
from each other while both are moving towards warmer temperatures. Regions in
Oceania, Central- and parts of South-America have $\Delta T$ values close to
zero, i.e. the cold and hot Gaussian peaks shift toward the warmer
temperatures at the same rate. This will cause these regions to have warmer
cold and hot periods under future global warming levels compared to the
historical period. When all regions are considered, we observe that the
extreme temperature events will increase everywhere, as the mean temperatures
increase in all regions compared to the historical distributions. The fact
that the peaks are converging only in cold climate regions while diverging in
other regions shows that shifts in the Gaussian components with respect to
each other are essential for extreme temperature event analyses as these
changes affect the overall distribution shape and extent. Also, these results
are consistent with the change in skewness in temperature distribution as
shown in previous studies [Tamarin-Brodsky . (2020), Skelton . (2020)].
Skelton2020 found an abrupt change in skewness in Europe. Tamarin2020 found
that changes in skewness in winter and summer months will cause cold anomalies
in Southern Europe, while warm anomalies intensify in Northeastern Europe.
They emphasize the importance of analysing the shape of temperature
distributions.

Figure 7: Multi-model peak mean change of region temperature distributions
from bimodal grid cells for SSP5-8.5 under GWL3.0∘C. Blue (red) dots and
squares are the means for cold (hot) peaks of the historical (future) period,
respectively. They are plotted on the left y-axis. Green bars describe $\Delta
T$, the change in the difference between the means of cold and hot Gaussian
components, and are plotted on the right y-axis. The upward shift in markers
represents the overall warming (see Supplementary Material Figure S7 to S9 for
other warming levels). Figure 8: Multi-model peak mean change of region
temperature distributions sorted by cold Gaussian mean temperatures (blue
dots) for SSP5-8.5 under GWL 3.0∘C. Blue (red) dots and squares are the means
for cold and hot peaks of the historical (future) period, respectively. They
are plotted on the left y-axis. Green bars describe $\Delta T$, the change in
the difference between the means of cold and hot Gaussian components, and are
plotted on the right y-axis. The colder regions have positive $\Delta T$
values and their absolute values are higher than the other regions. The upward
shift in blue dots shows that the temperature of cold days is getting warmer
and this increase is faster in polar regions compared to the rest of the world
(see Supplementary Material Figure S10 to S12 for other warming levels).

After analysing the distribution shapes and peak movements, we calculated the
return periods -the average time between the occurrences of a certain event-
of 1-year, 5-year, 10-year and 20-year events using only the grid cells with
constant modalities, i.e. unimodal or bimodal both for the historical and
future periods, as described in Equation 14. Instead of analyzing extreme
temperatures within specific time blocks, our analysis focused on the extremes
in the region’s probability distribution of 30(20)-years of daily maximum
temperatures. Since we used the hottest component in the mixture of Gaussian
components to define n-year events, we considered the number of data points
falling under the Gaussian component to define year-length according to
Equation 7. For example, globally a 10-year event was a temperature event once
in every 1880 days (10 years$\textstyle\times 188\textstyle\frac{days}{year}$)
(for bimodal distributions) in the historical period, but it will occur once
in every 643, 355, 138, and 63 days under GWL 1.5∘CC, 2.0∘C, 3.0∘C and 4.0∘C
scenarios, as shown in the plot showing global results in Figure 9 (also in
Figure 10), respectively. In other words, historical 10-year events will be
3.42-year, 1.89-year, 0.73-year and 0.34-year events under the future GWL
1.5∘C, 2.0∘C, 3.0∘C and 4.0∘C scenarios, respectively. After calculating the
frequency of extreme events using the temperature distributions in each grid
cell individually for an IPCC land region, we averaged the results for the
whole region for a single model. The global map with box plots in Figure 9
shows multi-model 10-year event frequencies of each region for SSP5-8.5
scenario under different GWL, where the boxes from light to dark shades of red
represent 1.5∘C, 2.0∘C, 3.0∘C and 4.0∘C. Results for 1-year, 5-year, and
20-year events are left out for simplicity and presented in the Supplementary
Material Figure S13 to S27. The length of a “year” in each region that is used
for return period calculations, i.e. the number of days in 10 years, is shown
on the top right corner of each sub-plot in Figure 9.

As shown in Figure 9, return periods of extreme temperature events are getting
shorter for all regions under all GWL scenarios as the median of each box is
smaller than the historical period. The frequency of extreme events is higher
in lower latitudes compared to higher latitudes. For example, the return
periods are getting prominently shorter in regions around the equator -where a
higher percentage of unimodal grid cells was observed- compared to the other
regions. Furthermore, CMIP6 models show narrower boxes and shorter whiskers in
lower latitudes compared to wider boxes and longer whiskers in higher
latitudes for all analyzed GWL. Among all analysed regions, the Caribbean
(CAR) region has the highest increase in the frequency of a 10-year event,
from once in 1910 days for the historical period to once in every 137.3,
35.32, 5.5 and 2.0 days under GWL 1.5, 2, 3, and 4∘C, respectively. Regions
around the equator (namely CAR, NSA, NWS, NES, SEA, SCA, SAM, MDG, WAF, and
SEAF regions) are the top 10 regions with the highest increase in the
frequency of extreme events under all GWL. The frequency of a temperature
event equivalent to a 10-year event (historically once in every 1610 days) in
the Mediterranean (MED) region increases to once in 405.6, 215.7, 72.4, and
30.6 days in the future under GWL 1.5, 2, 3, and 4∘C, respectively. Within the
European continent, the West&Central Europe (WCE) region has a higher increase
in the frequency of extreme events compared to the Eastern Europe (EEU) and
the North Eastern Europe (NEU) regions, where the latter two regions are among
the regions with the least increase in extreme temperature event frequency.
The smallest increase in the frequency of hot extremes is observed in the
Western Antarctica (WAN) region, where the return periods of 10-year events
will decrease from once in 1790 days to once in 1070.1, 827.6, 542.7 and 338.7
days under GWL 1.5, 2, 3, and 4∘C, respectively. High latitude regions, such
as WAN, NEU, EAN, NWN, ESB, GIC, RAR, SSA, TIB, and NEN regions are the 10
regions with the smallest decrease in return periods of extreme hot
temperature events. Some of these regions are polar regions with positive
$\Delta T$ values as shown in Figure 8. This will cause more days with warmer
temperatures in the colder mode of these regions while having an increase in
hot extremes.

Figure 9: Multi-model median of event frequencies for 10-year hot temperature
events compared to the 1985-2014 period under GWL 1.5, 2, 3 and 4∘C relative
to 1850-1900 baseline for SSP5-8.5 scenario. The boxes from light to dark
shades red represent 1.5∘C, 2.0∘C, 3.0∘C and 4.0∘C, respectively. The orange
lines inside the boxes show the CMIP6 multi-model median, and the boxes extend
between the first quartile (Q1) to the third quartile (Q3) of the data, i.e.
inter-quartile range (IQR). The vertical lines, i.e. whiskers, stretch out 1.5
IQR from the box. The circles represent the models outside of the
interquartile range, i.e. outliers. The length of the hot period used for
return period calculations, i.e. number of days in 10 years, is shown in the
top right corner of each plot. The global return periods are shown on the
left. The more outlier points in the global box plot are because of the
differences between regional return periods (See Supplementary Material Figure
S13 to S27 for other return periods).

Figure 10: Global multi-model median of event frequencies for 10-year
temperature events under 1.5, 2, 3 and 4∘C warming levels for a) SSP1-2.6, b)
SSP2-4.5, c) SSP3-7.0 and d) SSP5-8.5 scenarios. The orange lines inside the
boxes show the CMIP6 multi-model median, and the boxes extend between the
first quartile (Q1) to the third quartile (Q3) of the data, i.e. inter-
quartile range (IQR). The vertical lines, i.e. whiskers, stretch out 1.5 IQR
from the box. The circles represent the models outside of the interquartile
range, i.e. outliers. The length of the hot period used for return period
calculations, i.e. number of days in 10 years, is shown in the top right
corner of each plot. The number of datasets is given in parenthesises. All
plots show similar results for different SSP scenarios as the GWL are
scenario-independent.

##  4 Summary and Discussion

Detection of extreme events is important to mitigate their impact on natural
and anthropogenic systems. Future projections suggest that the mean and
standard deviations of maximum surface temperature will increase. This change
in the shape of maximum surface temperature distributions increases the
intensity and frequency of extreme events in the future. However, not only the
shift to warmer temperatures but also the modality of temperature distribution
affects the parameters of the entire distribution which is important to
calculate the return periods as shown in this study.

GMM are a promising method for calculating the return periods of extreme
events, and additionally determining the shape of the entire distribution for
daily maximum temperature data. GMM can provide information on different
climate features in different regions such as cold and hot periods, and their
changes. We showed that bimodality is a prominent characteristic observed in
multi-year daily near-surface maximum temperature data. To understand the
underlying factors of this bimodal pattern, we analyzed temperature
distributions from grid cells with distinct bimodality across different
months, seasons and 6-month running windows. We observed that the winter and
summer seasons emerged as the primary contributors to the peaks observed in
the bimodal distribution. In grid cells of different regions with distinct
bimodal distributions, the transition from winter to summer occurs swiftly,
leading to a more distinct separation of the temperature modes. Consequently,
the distributions during transitional seasons, such as spring and autumn,
appeared to be wider (covering a broader value range) compared to the more
distinct distributions observed during winter and summer (covering a very
small value range). Furthermore, analyses of 6-month running windows also
showed an agglomeration of similar temperatures around winter(summer) months
from November(May) to April(October) that creates the peaks in the bimodal
distribution (See Supplementary Material Figure S2 for the distributions of
seasons and months.). Here, the advantage of GMM becomes evident. For analyses
to uncover the origins of bimodality, we had to select certain seasons or
months. Seasonal periods are commonly used in previous studies to analyse
extreme events [Qian Zhang (2015), Qian Zhang (2019), Walt Fitchett (2021),
Prodhomme . (2022)]. For example, Qian2015 found that the seasonality is
weakening in the northern high-latitude regions and East Asia while
strengthening in the Mediterranean. This can also be seen in Figure 8, as the
northern regions and east/central Asia regions have converging peaks which
means that these regions will have a distribution closer to a unimodal
distribution. Meanwhile, diverging peaks in MED will introduce more distinct
cold and warm periods. However, onsets and length of seasons are predicted to
change with climate change [Wang . (2021)]. Therefore, the definition of
current seasonal periods or months will not necessarily be valid for future
climates. One can utilize GMM to determine the hot Gaussian component of a
region to define the length of the analysis period instead of using fixed
seasonal definitions. Moreover, the bimodality analysis also shows how peaks
are changing in the future, effectively changing the expected climate of the
area.

ETCCDI indices are commonly used in extreme event analysis as they offer a
simple and concise way to define extremes [X. Zhang . (2011), Zhao . (2021),
Vogel, Hauser Seneviratne (2020)]. ETCCDI indices use block maxima methods
such as TXx (Monthly maximum value of daily max temperature), TNx (Monthly
maximum value of daily min temperature) or percentile-over-threshold (POT)
methods such as TX90p (Percentage of time when daily max temperature
$>90^{th}$percentile), TN90p (Percentage of time when daily min temperature
$<90^{th}$percentile) [X. Zhang . (2011)]. These exceedances can be modelled
with GEV distributions or generalized Pareto distribution (GPD). However, GEV
distributions are a better fit for longer block sizes than for shorter blocks
like daily data. If the available dataset is short, the longer block sizes
will produce fewer data which can increase the variability in parameter
estimation [Huang . (2016), Wang . (2016)]. For example, if there is more than
one extremely hot day in the block (month, season or year), e.g. several
consecutive days, block maxima methods consider only the hottest, and hence
only one day in a block, while GMM considers all days hotter than the
threshold. Assuming that a heat wave lasts usually days to a few weeks, a
substantial number of hot days might not be seen by block maxima methods as
long as they fall into the same block. Percentile-over-threshold methods
together with count-day indices such as WSDI (Warm spell duration indicator)
are useful for analysing the durations of events. However, the derivation of
percentiles is strongly affected by the choice of the base period, a right
shift in the distribution will result in a higher threshold and erroneously
reduce the frequency of extreme events [Yosef . (2021)]. Seasonality in
temperature extremes adds complexity to the process of selecting percentiles
to define extreme temperatures [Huang . (2016)]. The advantage of GMM is that
the model analyses the distribution of temperatures without any previous
assumption and learns the hot periods from the data. Also, GMM uses all
available data in contrast to block maxima methods, which makes it useful if
the available data is short or bimodality exists [Sardeshmukh . (2015), Wang .
(2016), Knoben . (2019), Ben Alaya . (2020)].

However, since the Gaussian components of GMM are not bounded, it is important
to only calculate the return periods of extreme events equal to or less than
the study period when applying GMM. Additionally, we only used grid cells
which have the same number of Gaussian components in their temperature
distribution, i.e. unimodal or bimodal distribution, both for the historical
and future periods. Grid cells with changing distribution shapes, e.g.
transforming from a bimodal distribution in the historical period to a
unimodal distribution in the future or vice versa, were found in less than 10%
of the grid cell for each GWL as shown in Table 2, and were disregarded in the
analysis as calculating the temperature thresholds becomes problematic with
the abrupt change in means and standard deviations.

For the first time, the IPCC AR6 Report includes a new dedicated chapter on
weather and climate extreme events [IPCC (2021)]. This emphasizes the
importance of robust methods of extreme event detection to be able to mitigate
the impact of such events. IPCC AR6 reports that the return periods of 10-year
events will increase around the world, with the highest changes projected to
happen in some mid-latitude and semi-arid regions. Our findings are in
agreement with these results. Furthermore, IPCC AR6 projects the warming rate
in mid-latitudes to be higher than the average global warming rate. GMM might
explain why these regions are projected to have higher warming, as we observed
that grid cells in these regions predominantly follow a bimodal distribution
in the historical (future) period as shown in Figure 3 (5). Furthermore, these
regions have diverging peaks as shown in Figure 8, i.e. mode for warm
temperatures moving towards warmer temperatures faster than the mode for
colder temperatures. These diverging bimodal peaks will create distinct
Gaussian components in the entire multi-year daily maximum temperature, which
in turn results in a higher increase in extremes in these regions. For
example, almost all grid cells in the Mediterranean region follow a bimodal
distribution, and the peaks of bimodal distribution will diverge in the future
as shown in Figure 7. Mediterranean region is identified as one of the most
responsive regions to climate change and a hot spot of climate extremes [IPCC
(2021), Feng . (2022)]. Similarly, Arctic regions are projected to have the
highest increase in temperature of the coldest days [IPCC (2021), C. Li .
(2021)]. Our results are also consistent with these increases as shown in
Figure 7, where diverging bimodal peaks in mid-latitude regions will shift the
mode for warm temperatures, i.e. hot Gaussian, to the higher temperature
ranges. This shift in the Gaussian components of temperature distribution will
cause those land regions to have warmer temperature extremes and can explain
the higher average warming rate than the global average. Likewise, converging
peaks in polar regions as shown in Figure 7 will move the cold Gaussian part
toward warmer temperatures, thereby introducing higher warming on the coldest
days.

According to our analyses, 10-year events will increase almost 3-fold under
GWL 1.5∘C compared to the historical period for all SSP scenarios as shown in
Figure 10 when looking at the whole globe. This means a temperature event that
occurs once in every 10 years (1880 days) will be expected to occur 2.9 times
in every 10 years under GWL 1.5∘C. 10-year extreme temperature events will
become even more frequent globally under GWL 2∘C, 3∘C and 4∘C; 5.3, 13.6, and
29.5 times every 10 years, respectively. In other words, current 10-year
events will be 3.42-year, 1.89-year, 0.73-year and 0.34-year events in the
future under GWL 1.5∘C, 2∘C, 3∘C and 4∘C, respectively. Our results show a
higher increase compared to the IPCC AR6 report, where the frequency of
10-year events is projected to increase approximately 3, 4, 5.5 and 9-fold
under GWL 1.5∘C, 2∘C, 3∘C and 4∘C, respectively [IPCC (2021)], using a block
maxima method for determining the extreme events. The higher increase in our
method compared to IPCC AR6 can most likely be explained by the fact that we
used GMM to model the distribution of temperatures and GMM considers all days
hotter than the threshold, while the block maxima method only uses the maximum
of a block. Another important point deduced from the analyses of different
regions for several CMIP6 models is that the ensemble of analyzed CMIP6 models
shows coherent results for regions as shown in the regional box plots in
Figure 9. Most of the individual model results fall within the first and third
quartile, and only a few models fall outside this range. The higher number of
outlier points in the global box plot in Figure 9, and also shown for
different SSP scenarios in Figure 10, are caused by the differences between
regional return periods. All SSP scenarios show similar results with each
other as the return periods are calculated for GWL which have the same forcing
on climate.

Return periods of extreme events become shorter in every region, which means
that the frequency of extreme temperature events increases. This will become
larger with increasing global warming levels. Some climate models have already
exceeded GWL 1.5∘C with respect to the 1850-1900 period as shown in Figure 1.
This fact further emphasises the importance of robust methods to detect
extreme events. Even though there is a delay in taking the necessary
precautions to reduce the speed of the warming of the climate, as time goes
by, tomorrow’s projections become today’s reality.

## Code and data availability

The recipes to extract regional data from CMIP6 models using ESMValTool,
python scripts to analyse extreme events and to produce all figures of this
manuscript are accessible in the following GitHub repository:
https://github.com/EyringMLClimateGroup/pacal23jgr_GaussianMixtureModels_Extremes.
The regional output files amount to hundreds of GB.

The latest release of ESMValTool is publicly at
https://github.com/ESMValGroup/ESMValTool [Andela . (2022)].

###### Acknowledgements.

Funding for this study was provided by the European Research Council (ERC)
Synergy Grant “Understanding and Modelling the Earth System with Machine
Learning (USMILE)” under the Horizon 2020 research and innovation programme
(Grant Agreement No. 855187). This work used resources of the Deutsches
Klimarechenzentrum (DKRZ) granted by its Scientific Steering Committee (WLA)
under project ID bd1083. We acknowledge the World Climate Research Programme,
which, through its Working Group on Coupled Modelling, coordinated and
promoted CMIP6. We thank the climate modelling groups for producing and making
available their model outputs, the Earth System Grid Federation (ESGF) for
archiving the data and providing access, and the multiple funding agencies
that support CMIP6 and ESGF. ERA5_data was downloaded from the C3S. The
results contain modified Copernicus Climate Change Service information 2020.
Neither the European Commission nor ECMWF is responsible for any use that may
be made of the Copernicus information or data it contains. We would like to
thank Dr. Pauline Bonnet for her valuable comments and suggestions to improve
the manuscript. We would like to extend our sincere gratitude to the anonymous
reviewers whose invaluable feedback and constructive comments significantly
contributed to the improvement and quality of this work.

## References

  * Alexander . (2006) Alexander2006Alexander, LV., Zhang, X., Peterson, TC., Caesar, J., Gleason, B., Klein Tank, AMG.Vazquez-Aguirre, JL.  2006\.  Global observed changes in daily climate extremes of temperature and precipitation Global observed changes in daily climate extremes of temperature and precipitation.  Journal of Geophysical Research: Atmospheres111D5.  10.1029/2005JD006290
  * Alizadeh . (2020) Alizadeh2020Alizadeh, MR., Adamowski, J., Nikoo, MR., AghaKouchak, A., Dennison, P. Sadegh, M.  2020\.  A century of observations reveals increasing likelihood of continental-scale compound dry-hot extremes A century of observations reveals increasing likelihood of continental-scale compound dry-hot extremes.  Science Advances639eaaz4571.  10.1126/sciadv.aaz4571
  * Andela . (2022) Andela_ESMValTool_2022Andela, B., Broetz, B., de Mora, L., Drost, N., Eyring, V., Koldunov, N.Winterstein, F.  2022\.  ESMValTool ESMValTool [software].  https://github.com/ESMValGroup/ESMValTool/ 10.5281/zenodo.3401363
  * Anderson Kostinski (2010) Anderson2010Anderson, A. Kostinski, A.  2010\.  Reversible Record Breaking and Variability: Temperature Distributions across the Globe Reversible record breaking and variability: Temperature distributions across the globe.  Journal of Applied Meteorology and Climatology4981681–1691.  10.1175/2010JAMC2407.1
  * Baker . (2018) Baker2018Baker, HS., Millar, RJ., Karoly, DJ., Beyerle, U., Guillod, BP., Mitchell, D.Allen, MR.  2018\.  Higher CO2 concentrations increase extreme event risk in a 1.5${}^{\circ}C$ world Higher CO2 concentrations increase extreme event risk in a 1.5${}^{\circ}C$ world.  Nature Climate Change87604–608.  10.1038/s41558-018-0190-1
  * Barriopedro . (2011) Barriopedro2011Barriopedro, D., Fischer, EM., Luterbacher, J., Trigo, RM. García-Herrera, R.  2011\.  The Hot Summer of 2010: Redrawing the Temperature Record Map of Europe The Hot Summer of 2010: Redrawing the Temperature Record Map of Europe.  Science3326026220-224.  10.1126/science.1201224
  * Beck . (2018) Beck2018Beck, HE., Zimmermann, NE., McVicar, TR., Vergopolan, N., Berg, A. Wood, EF.  2018\.  Present and future Köppen-Geiger climate classification maps at 1-km resolution Present and future Köppen-Geiger climate classification maps at 1-km resolution.  Scientific Data51180214.  10.1038/sdata.2018.214
  * Ben Alaya . (2020) Alaya2020Ben Alaya, MA., Zwiers, F. Zhang, X.  2020\.  An Evaluation of Block-Maximum-Based Estimation of Very Long Return Period Precipitation Extremes with a Large Ensemble Climate Simulation An evaluation of block-maximum-based estimation of very long return period precipitation extremes with a large ensemble climate simulation.  Journal of Climate33166957–6970.  10.1175/JCLI-D-19-0011.1
  * Benestad (2004) Benestad2004Benestad, R.  2004\.  Record-values, nonstationarity tests and extreme value distributions Record-values, nonstationarity tests and extreme value distributions.  Global and Planetary Change441–411–26.  10.1016/j.gloplacha.2004.06.002
  * Bentsen . (2019) NorESM2-MMBentsen, M., Oliviè, DJL., Seland, y., Toniazzo, T., Gjermundsen, A., Graff, LS.Schulz, M.  2019\.  NCC NorESM2-MM model output prepared for CMIP6 CMIP historical NCC NorESM2-MM model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.8040
  * Bertossa . (2021) Bertossa2021Bertossa, C., Hitchcock, P., DeGaetano, A. Plougonven, R.  2021\.  Bimodality in ensemble forecasts of 2 m temperature: identification Bimodality in ensemble forecasts of 2 m temperature: identification.  Weather and Climate Dynamics241209–1224.  https://wcd.copernicus.org/articles/2/1209/2021/ 10.5194/wcd-2-1209-2021
  * Bock . (2020) Bock2020Bock, L., Lauer, A., Schlund, M., Barreiro, M., Bellouin, N., Jones, C.Eyring, V.  2020\.  Quantifying Progress Across Different CMIP Phases With the ESMValTool Quantifying Progress Across Different CMIP Phases With the ESMValTool.  Journal of Geophysical Research: Atmospheres12521e2019JD032321.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JD032321 https://doi.org/10.1029/2019JD032321
  * Boucher . (2018) IPSL-CM6A-LRBoucher, O., Denvil, S., Levavasseur, G., Cozic, A., Caubel, A., Foujols, MA.Cheruy, F.  2018\.  IPSL IPSL-CM6A-LR model output prepared for CMIP6 CMIP historical IPSL IPSL-CM6A-LR model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.5195
  * Byun . (2019) KACE-1-0-GByun, YH., Lim, YJ., Sung, HM., Kim, J., Sun, M. Kim, BH.  2019\.  NIMS-KMA KACE1.0-G model output prepared for CMIP6 CMIP historical NIMS-KMA KACE1.0-G model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.8378
  * Cai . (2014) Cai2014Cai, W., Borlace, S., Lengaigne, M., van Rensch, P., Collins, M., Vecchi, G.Jin, FF.  2014\.  Increasing frequency of extreme El Niño events due to greenhouse warming Increasing frequency of extreme El Niño events due to greenhouse warming.  Nature Climate Change422111–116.  10.1038/nclimate2100
  * Cao Wang (2019) NESM3Cao, J. Wang, B.  2019\.  NUIST NESMv3 model output prepared for CMIP6 CMIP historical NUIST NESMv3 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.8769
  * Cho Jeong (2016) Hongyeon2016Cho, H. Jeong, S.  2016\.  Estimation and Comparative Analysis on the Distribution Functions of Air and Water Temperatures in Korean Coastal Seas Estimation and Comparative Analysis on the Distribution Functions of Air and Water Temperatures in Korean Coastal Seas.  Journal of Korean Society of Coastal and Ocean Engineers28171-176.  10.9765/KSCOE.2016.28.3.171
  * Christidis . (2015) Christidis2015Christidis, N., Jones, GS. Stott, PA.  2015\.  Dramatically increasing chance of extremely hot summers since the 2003 European heatwave Dramatically increasing chance of extremely hot summers since the 2003 European heatwave.  Nature Climate Change51146–50.  10.1038/nclimate2468
  * Clarkson . (2022) Clarkson2022Clarkson, D., Eastoe, E. Leeson, A.  2022\.  Melt probabilities and surface temperature trends on the Greenland ice sheet using a Gaussian mixture model Melt probabilities and surface temperature trends on the Greenland ice sheet using a Gaussian mixture model.  The Cryosphere1651597–1607.  https://tc.copernicus.org/articles/16/1597/2022/ 10.5194/tc-16-1597-2022
  * Copernicus Climate Change Service (C3S) (2023) C3SCopernicus Climate Change Service (C3S).  2023\.  ERA5 hourly data on single levels from 1940 to present. ERA5 hourly data on single levels from 1940 to present.  Copernicus Climate Change Service (C3S) Climate Data Store (CDS).  10.24381/cds.adbb2d47
  * Demiroglu . (2020) Demiroglu2020Demiroglu, OC., Saygili-Araci, FS., Pacal, A., Hall, CM. Kurnaz, ML.  2020\.  Future Holiday Climate Index (HCI) Performance of Urban and Beach Destinations in the Mediterranean Future Holiday Climate Index (HCI) Performance of Urban and Beach Destinations in the Mediterranean.  Atmosphere119911.  10.3390/atmos11090911
  * Dix . (2019) ACCESS-CM2Dix, M., Bi, D., Dobrohotoff, P., Fiedler, R., Harman, I., Law, R.Yang, R.  2019\.  CSIRO-ARCCSS ACCESS-CM2 model output prepared for CMIP6 CMIP historical CSIRO-ARCCSS ACCESS-CM2 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4271
  * Donat Alexander (2012) Donat2012Donat, MG. Alexander, LV.  2012\.  The shifting probability distribution of global daytime and night-time temperatures The shifting probability distribution of global daytime and night-time temperatures.  Geophysical Research Letters3914.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2012GL052459 https://doi.org/10.1029/2012GL052459
  * EC-Earth Consortium (EC-Earth) (20191) EC-Earth3EC-Earth Consortium (EC-Earth).  20191\.  EC-Earth-Consortium EC-Earth3 model output prepared for CMIP6 CMIP historical EC-Earth-Consortium EC-Earth3 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4700
  * EC-Earth Consortium (EC-Earth) (20192) EC-Earth3-VegEC-Earth Consortium (EC-Earth).  20192\.  EC-Earth-Consortium EC-Earth3-Veg model output prepared for CMIP6 CMIP historical EC-Earth-Consortium EC-Earth3-Veg model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4706
  * EC-Earth Consortium (EC-Earth) (2020) EC-Earth3-Veg-LREC-Earth Consortium (EC-Earth).  2020\.  EC-Earth-Consortium EC-Earth3-Veg-LR model output prepared for CMIP6 CMIP historical EC-Earth-Consortium EC-Earth3-Veg-LR model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4707
  * EC-Earth Consortium (EC-Earth) (2021) EC-Earth3-CCEC-Earth Consortium (EC-Earth).  2021\.  EC-Earth-Consortium EC-Earth-3-CC model output prepared for CMIP6 CMIP historical EC-Earth-Consortium EC-Earth-3-CC model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4702
  * Eyring . (2020) Eyring2020Eyring, V., Bock, L., Lauer, A., Righi, M., Schlund, M., Andela, B.Zimmermann, K.  2020\.  Earth System Model Evaluation Tool (ESMValTool) v2.0 – an extended set of large-scale diagnostics for quasi-operational and comprehensive evaluation of Earth system models in CMIP Earth System Model Evaluation Tool (ESMValTool) v2.0 – an extended set of large-scale diagnostics for quasi-operational and comprehensive evaluation of Earth system models in CMIP.  Geoscientific Model Development1373383–3438.  10.5194/gmd-13-3383-2020
  * Eyring . (2016) Eyring2016Eyring, V., Bony, S., Meehl, GA., Senior, CA., Stevens, B., Stouffer, RJ. Taylor, KE.  2016\.  Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization.  Geoscientific Model Development951937–1958.  10.5194/gmd-9-1937-2016
  * Eyring . (2021) IPCC_2021_WGI_Ch_3Eyring, V., Gillett, NP., Achuta Rao, KM., Barimalala, R., Barreiro Parrillo, M., Bellouin, N.Sun, Y.  2021\.  Human Influence on the Climate System Human influence on the climate system [Book Section].  V. Masson-Delmotte . (), Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change Climate change 2021: The physical science basis. contribution of working group i to the sixth assessment report of the intergovernmental panel on climate change ( 3).  Cambridge, United Kingdom and New York, NY, USACambridge University Press.
  * Feng . (2022) Feng2022Feng, X., Qian, C. Materia, S.  2022\.  Amplification of the Temperature Seasonality in the Mediterranean Region Under Anthropogenic Climate Change Amplification of the Temperature Seasonality in the Mediterranean Region Under Anthropogenic Climate Change.  Geophysical Research Letters4920e2022GL099658.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022GL099658 https://doi.org/10.1029/2022GL099658
  * Fischer Schär (2010) Fischer2010Fischer, EM. Schär, C.  2010\.  Consistent geographical patterns of changes in high-impact European heatwaves Consistent geographical patterns of changes in high-impact European heatwaves.  Nature Geoscience36398-403.  10.1038/ngeo866
  * García-León . (2021) Garcia2021García-León, D., Casanueva, A., Standardi, G., Burgstall, A., Flouris, AD. Nybo, L.  2021\.  Current and projected regional economic impacts of heatwaves in Europe Current and projected regional economic impacts of heatwaves in Europe.  Nature Communications1215807.  10.1038/s41467-021-26050-z
  * Grace (1995) Grace1995Grace, W.  1995\.  A binormal model of frequency distributions of daily maximum temperature A binormal model of frequency distributions of daily maximum temperature.  Australian meteorological magazine42151-161.
  * Guerreiro . (2018) Guerreiro2018Guerreiro, SB., Dawson, RJ., Kilsby, C., Lewis, E. Ford, A.  2018\.  Future heat-waves, droughts and floods in 571 European cities Future heat-waves, droughts and floods in 571 European cities.  Environmental Research Letters133034009.  10.1088/1748-9326/aaaad3
  * Hajat . (2022) hajat2022Hajat, S., Gampe, D., Sarsour, A. Abuzerr, S.  2022\.  Climate Change and Diarrhoeal Disease Burdens in the Gaza Strip, Palestine: Health Impacts of 1.5C and 2C Global Warming Scenarios Climate change and diarrhoeal disease burdens in the gaza strip, palestine: Health impacts of 1.5c and 2c global warming scenarios.  International Journal of Environmental Research and Public Health198.  https://www.mdpi.com/1660-4601/19/8/4898 10.3390/ijerph19084898
  * Hajima . (2019) MIROC-ES2LHajima, T., Abe, M., Arakawa, O., Suzuki, T., Komuro, Y., Ogura, T.Tachiiri, K.  2019\.  MIROC MIROC-ES2L model output prepared for CMIP6 CMIP historical MIROC MIROC-ES2L model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.5602
  * Hao . (2018) Hao2018Hao, Z., Hao, F., Singh, VP. Zhang, X.  2018\.  Changes in the severity of compound drought and hot extremes over global land areas Changes in the severity of compound drought and hot extremes over global land areas.  Environmental Research Letters1312124022.  10.1088/1748-9326/aaee96
  * Hauser . (2022) gwl_periodsHauser, M., Engelbrecht, F. Fischer, EM.  2022\.  Transient global warming levels for CMIP5 and CMIP6 Transient global warming levels for CMIP5 and CMIP6 [dataset].  Zenodo.  https://doi.org/10.5281/zenodo.7390473 10.5281/zenodo.7390473
  * Hersbach . (2018) ERA5_dataHersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J.Thépaut, JN.  2018\.  ERA5 hourly data on single levels from 1940 to present. ERA5 hourly data on single levels from 1940 to present.  Copernicus Climate Change Service (C3S) Climate Data Store (CDS).  10.24381/cds.adbb2d47
  * Huang . (2016) Huang2016Huang, WK., Stein, ML., McInerney, DJ., Sun, S. Moyer, EJ.  2016\.  Estimating changes in temperature extremes from millennial-scale climate simulations using generalized extreme value (GEV) distributions Estimating changes in temperature extremes from millennial-scale climate simulations using generalized extreme value (GEV) distributions.  Advances in Statistical Climatology, Meteorology and Oceanography2179–103.  10.5194/ascmo-2-79-2016
  * IPCC (2021) IPCC2021IPCC.  2021\.  Climate Change 2021: The Physical Science Basis Climate Change 2021: The Physical Science Basis (V. Masson-Delmotte ., ).  Cambridge University Press (In Press).
  * Iturbide . (2020) Iturbide2020Iturbide, M., Gutiérrez, JM., Alves, LM., Bedia, J., Cerezo-Mota, R., Cimadevilla, E.Vera, CS.  2020\.  An update of IPCC climate reference regions for subcontinental analysis of climate model data: definition and aggregated datasets An update of IPCC climate reference regions for subcontinental analysis of climate model data: definition and aggregated datasets.  Earth System Science Data1242959–2970.  10.5194/essd-12-2959-2020
  * Iyakaremye . (2022) Iyakaremye2022Iyakaremye, V., Zeng, G., Ullah, I., Gahigi, A., Mumo, R. Ayugi, B.  2022\.  Recent Observed Changes in Extreme High‐Temperature Events and Associated Meteorological Conditions over Africa Recent Observed Changes in Extreme High‐Temperature Events and Associated Meteorological Conditions over Africa.  International Journal of Climatologyjoc.7485.  10.1002/joc.7485
  * Iyakaremye . (2021) Iyakaremye2021Iyakaremye, V., Zeng, G. Zhang, G.  2021\.  Changes in extreme temperature events over Africa under 1.5 and 2.0C global warming scenarios Changes in extreme temperature events over africa under 1.5 and 2.0c global warming scenarios.  International Journal of Climatology4121506–1524.  10.1002/joc.6868
  * Jungclaus . (2019) MPI-ESM1-2-HRJungclaus, J., Bittner, M., Wieners, KH., Wachsmann, F., Schupfner, M., Legutke, S.Roeckner, E.  2019\.  MPI-M MPI-ESM1.2-HR model output prepared for CMIP6 CMIP historical MPI-M MPI-ESM1.2-HR model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.6594
  * Kharin . (2013) Kharin2013Kharin, VV., Zwiers, FW., Zhang, X. Wehner, M.  2013\.  Changes in temperature and precipitation extremes in the CMIP5 ensemble Changes in temperature and precipitation extremes in the CMIP5 ensemble.  Climatic Change1192345–357.  10.1007/s10584-013-0705-8
  * Knoben . (2019) Knoben2019Knoben, WJM., Woods, RA. Freer, JE.  2019\.  Global bimodal precipitation seasonality: A systematic overview Global bimodal precipitation seasonality: A systematic overview.  International Journal of Climatology391558-567.  https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/joc.5786 https://doi.org/10.1002/joc.5786
  * Kodra Ganguly (2014) kodra_asymmetry_2014Kodra, E. Ganguly, AR.  2014\.  Asymmetry of projected increases in extreme temperature distributions Asymmetry of projected increases in extreme temperature distributions.  Scientific Reports415884.  [2023-01-08]https://www.nature.com/articles/srep05884 10.1038/srep05884
  * Kollu . (2012) Kollu2012Kollu, R., Rayapudi, SR., Narasimham, S. Pakkurthi, KM.  2012\.  Mixture probability distribution functions to model wind speed distributions Mixture probability distribution functions to model wind speed distributions.  International Journal of Energy and Environmental Engineering3127.  10.1186/2251-6832-3-27
  * Krasting . (2018) GFDL-ESM4Krasting, JP., John, JG., Blanton, C., McHugh, C., Nikonov, S., Radhakrishnan, A.Zhao, M.  2018\.  NOAA-GFDL GFDL-ESM4 model output prepared for CMIP6 CMIP historical NOAA-GFDL GFDL-ESM4 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.8597
  * Lauer . (2020) Lauer2020Lauer, A., Eyring, V., Bellprat, O., Bock, L., Gier, BK., Hunter, A.Zechlau, S.  2020\.  Earth System Model Evaluation Tool (ESMValTool) v2.0 – diagnostics for emergent constraints and future projections from Earth system models in CMIP Earth System Model Evaluation Tool (ESMValTool) v2.0 – diagnostics for emergent constraints and future projections from Earth system models in CMIP.  Geoscientific Model Development1394205–4228.  10.5194/gmd-13-4205-2020
  * Lelieveld . (2016) Lelieveld2016Lelieveld, J., Proestos, Y., Hadjinicolaou, P., Tanarhte, M., Tyrlis, E. Zittis, G.  2016\.  Strongly increasing heat extremes in the Middle East and North Africa (MENA) in the 21st century Strongly increasing heat extremes in the Middle East and North Africa (MENA) in the 21st century.  Climatic Change1371–2245–260.  10.1007/s10584-016-1665-6
  * Lewis . (2017) Lewis2017Lewis, SC., King, AD. Perkins-Kirkpatrick, SE.  2017\.  Defining a New Normal for Extremes in a Warming World Defining a new normal for extremes in a warming world.  Bulletin of the American Meteorological Society9861139–1151.  10.1175/BAMS-D-16-0183.1
  * C. Li . (2019) Li2019Li, C., Zwiers, F., Zhang, X. Li, G.  2019\.  How Much Information Is Required to Well Constrain Local Estimates of Future Precipitation Extremes? How much information is required to well constrain local estimates of future precipitation extremes?  Earth’s Future7111–24.  10.1029/2018EF001001
  * C. Li . (2021) Li2021Li, C., Zwiers, F., Zhang, X., Li, G., Sun, Y. Wehner, M.  2021\.  Changes in Annual Extremes of Daily Temperature and Precipitation in CMIP6 Models Changes in Annual Extremes of Daily Temperature and Precipitation in CMIP6 Models.  Journal of Climate3493441 - 3460.  10.1175/JCLI-D-19-1013.1
  * L. Li (2019) FGOALS-g3Li, L.  2019\.  CAS FGOALS-g3 model output prepared for CMIP6 CMIP historical CAS FGOALS-g3 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.3356
  * Linz . (2018) Linz2018Linz, M., Chen, G. Hu, Z.  2018\.  Large‐Scale Atmospheric Control on Non‐Gaussian Tails of Midlatitude Temperature Distributions Large‐scale atmospheric control on non‐gaussian tails of midlatitude temperature distributions.  Geophysical Research Letters45179141–9149.  10.1029/2018GL079324
  * Loikith Neelin (2019) Loikith2019Loikith, PC. Neelin, JD.  2019\.  Non-Gaussian Cold-Side Temperature Distribution Tails and Associated Synoptic Meteorology Non-gaussian cold-side temperature distribution tails and associated synoptic meteorology.  Journal of Climate32238399–8414.  10.1175/JCLI-D-19-0344.1
  * Mallick . (2022) Mallick2022Mallick, J., Islam, ARMT., Ghose, B., Islam, HMT., Rana, Y., Hu, Z.Ismail, ZB.  2022\.  Spatiotemporal trends of temperature extremes in Bangladesh under changing climate using multi-statistical techniques Spatiotemporal trends of temperature extremes in Bangladesh under changing climate using multi-statistical techniques.  Theoretical and Applied Climatology1471–2307–324.  10.1007/s00704-021-03828-1
  * Meehl Tebaldi (2004) Meehl2004Meehl, GA. Tebaldi, C.  2004\.  More Intense, More Frequent, and Longer Lasting Heat Waves in the 21st Century More Intense, More Frequent, and Longer Lasting Heat Waves in the 21st Century.  Science3055686994–997.  10.1126/science.1098704
  * Mueller . (2016) Mueller2016Mueller, B., Zhang, X. Zwiers, FW.  2016\.  Historically hottest summers projected to be the norm for more than half of the world’s population within 20 years Historically hottest summers projected to be the norm for more than half of the world’s population within 20 years.  Environmental Research Letters114044011.  10.1088/1748-9326/11/4/044011
  * Nikulin . (2011) Nikulin2011Nikulin, G., Kjellström, E., Hansson, U., Strandberg, G. Ullerstig, A.  2011\.  Evaluation and future projections of temperature, precipitation and wind extremes over Europe in an ensemble of regional climate simulations Evaluation and future projections of temperature, precipitation and wind extremes over Europe in an ensemble of regional climate simulations.  Tellus A63141–55.  10.1111/j.1600-0870.2010.00466.x
  * O’Neill . (2014) ONeill2014O’Neill, BC., Kriegler, E., Riahi, K., Ebi, KL., Hallegatte, S., Carter, TR.van Vuuren, DP.  2014\.  A new scenario framework for climate change research: the concept of shared socioeconomic pathways A new scenario framework for climate change research: the concept of shared socioeconomic pathways.  Climatic Change1223387-400.  10.1007/s10584-013-0905-2
  * O’Neill . (2016) Oneill2016O’Neill, BC., Tebaldi, C., van Vuuren, DP., Eyring, V., Friedlingstein, P., Hurtt, G.Sanderson, BM.  2016\.  The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6 The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6.  Geoscientific Model Development993461–3482.  10.5194/gmd-9-3461-2016
  * Ozturk . (2021) Ozturk2021Ozturk, T., Saygili-Araci, FS. Kurnaz, ML.  2021\.  Projected Changes in Extreme Temperature and Precipitation Indices Over CORDEX-MENA Domain Projected Changes in Extreme Temperature and Precipitation Indices Over CORDEX-MENA Domain.  Atmosphere1255622.  10.3390/atmos12050622
  * Paciorek . (2018) Paciorek2018Paciorek, CJ., Stone, DA. Wehner, MF.  2018\.  Quantifying statistical uncertainty in the attribution of human influence on severe weather Quantifying statistical uncertainty in the attribution of human influence on severe weather.  Weather and Climate Extremes2069-80.  https://www.sciencedirect.com/science/article/pii/S2212094717300841 https://doi.org/10.1016/j.wace.2018.01.002
  * Pedregosa . (2011) scikit-learnPedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O.Duchesnay, E.  2011\.  Scikit-learn: Machine Learning in Python Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research122825–2830.
  * Perera . (2020) Perera2020Perera, ATD., Nik, VM., Chen, D., Scartezzini, JL. Hong, T.  2020\.  Quantifying the impacts of climate change and extreme climate events on energy systems Quantifying the impacts of climate change and extreme climate events on energy systems.  Nature Energy52150-159.  10.1038/s41560-020-0558-0
  * Perkins-Kirkpatrick Lewis (2020) PerkinsKirkpatrick2020Perkins-Kirkpatrick, SE. Lewis, SC.  2020\.  Increasing trends in regional heatwaves Increasing trends in regional heatwaves.  Nature Communications1113357.  10.1038/s41467-020-16970-7
  * Perron Sura (2013) Perron_Sura_2013Perron, M. Sura, P.  2013\.  Climatology of Non-Gaussian Atmospheric Statistics Climatology of Non-Gaussian Atmospheric Statistics.  Journal of Climate2631063–1083.  10.1175/JCLI-D-11-00504.1
  * Pfleiderer . (2019) Pfleiderer2019Pfleiderer, P., Schleussner, CF., Kornhuber, K. Coumou, D.  2019\.  Summer weather becomes more persistent in a 2 °C world Summer weather becomes more persistent in a 2 °C world.  Nature Climate Change99666–671.  https://doi.org/10.1038/s41558-019-0555-0 10.1038/s41558-019-0555-0
  * Prodhomme . (2022) Prodhomme2022Prodhomme, C., Materia, S., Ardilouze, C., White, RH., Batté, L., Guemas, V.García-Serrano, J.  2022\.  Seasonal prediction of European summer heatwaves Seasonal prediction of European summer heatwaves.  Climate Dynamics587-82149–2166.  [2022-05-02]https://link.springer.com/10.1007/s00382-021-05828-3 10.1007/s00382-021-05828-3
  * Qian Zhang (2015) Qian2015Qian, C. Zhang, X.  2015\.  Human Influences on Changes in the Temperature Seasonality in Mid- to High-Latitude Land Areas Human Influences on Changes in the Temperature Seasonality in Mid- to High-Latitude Land Areas.  Journal of Climate28155908 - 5921.  https://journals.ametsoc.org/view/journals/clim/28/15/jcli-d-14-00821.1.xml https://doi.org/10.1175/JCLI-D-14-00821.1
  * Qian Zhang (2019) Qian2019Qian, C. Zhang, X.  2019\.  Changes in Temperature Seasonality in China: Human Influences and Internal Variability Changes in Temperature Seasonality in China: Human Influences and Internal Variability.  Journal of Climate32196237 - 6249.  https://journals.ametsoc.org/view/journals/clim/32/19/jcli-d-19-0081.1.xml https://doi.org/10.1175/JCLI-D-19-0081.1
  * Rahmstorf Coumou (2011) Rahmstrof2011Rahmstorf, S. Coumou, D.  2011\.  Increase of extreme events in a warming world Increase of extreme events in a warming world.  Proceedings of the National Academy of Sciences1084417905-17909.  10.1073/pnas.1101766108
  * Raymond . (2020) Raymond2020Raymond, C., Matthews, T. Horton, RM.  2020\.  The emergence of heat and humidity too severe for human tolerance The emergence of heat and humidity too severe for human tolerance.  Science Advances619.  10.1126/sciadv.aaw1838
  * Ribeiro . (2022) Ribeiro2022Ribeiro, AFS., Brando, PM., Santos, L., Rattis, L., Hirschi, M., Hauser, M.Zscheischler, J.  2022jun.  A compound event-oriented framework to tropical fire risk assessment in a changing climate A compound event-oriented framework to tropical fire risk assessment in a changing climate.  Environmental Research Letters176065015.  https://dx.doi.org/10.1088/1748-9326/ac7342 10.1088/1748-9326/ac7342
  * Richter (2016) Richter2016Richter, M.  2016\.  Temperatures in the Tropics Temperatures in the tropics.  L. Pancel M. Köhl (), Tropical Forestry Handbook Tropical forestry handbook ( 343–361).  Berlin, HeidelbergSpringer Berlin Heidelberg.  10.1007/978-3-642-54601-3_38
  * Ridley . (20191) HadGEM3-GC31-LLRidley, J., Menary, M., Kuhlbrodt, T., Andrews, M. Andrews, T.  20191\.  MOHC HadGEM3-GC31-LL model output prepared for CMIP6 CMIP historical MOHC HadGEM3-GC31-LL model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.6109
  * Ridley . (20192) HadGEM3-GC31-MMRidley, J., Menary, M., Kuhlbrodt, T., Andrews, M. Andrews, T.  20192\.  MOHC HadGEM3-GC31-MM model output prepared for CMIP6 CMIP historical MOHC HadGEM3-GC31-MM model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.6112
  * Righi . (2020) Righi2020Righi, M., Andela, B., Eyring, V., Lauer, A., Predoi, V., Schlund, M.Zimmermann, K.  2020\.  Earth System Model Evaluation Tool (ESMValTool) v2.0 – technical overview Earth System Model Evaluation Tool (ESMValTool) v2.0 – technical overview.  Geoscientific Model Development1331179–1199.  10.5194/gmd-13-1179-2020
  * Robine . (2008) Robine2008Robine, JM., Cheung, SLK., Le Roy, S., Van Oyen, H., Griffiths, C., Michel, JP. Herrmann, FR.  2008\.  Death toll exceeded 70,000 in Europe during the summer of 2003 Death toll exceeded 70,000 in Europe during the summer of 2003.  Comptes Rendus Biologies3312171-178.  10.1016/j.crvi.2007.12.001
  * Robinson . (2021) Robinson2021Robinson, A., Lehmann, J., Barriopedro, D., Rahmstorf, S. Coumou, D.  2021\.  Increasing heat and rainfall extremes now far outside the historical climate Increasing heat and rainfall extremes now far outside the historical climate.  npj Climate and Atmospheric Science4111–4.  10.1038/s41612-021-00202-w
  * Sardeshmukh . (2015) Sardeshmukh2015Sardeshmukh, PD., Compo, GP. Penland, C.  2015\.  Need for Caution in Interpreting Extreme Weather Statistics Need for caution in interpreting extreme weather statistics.  Journal of Climate28239166–9187.  10.1175/JCLI-D-15-0020.1
  * Sardeshmukh Sura (2009) Sardeshmukh2009Sardeshmukh, PD. Sura, P.  2009\.  Reconciling Non-Gaussian Climate Statistics with Linear Dynamics Reconciling non-gaussian climate statistics with linear dynamics.  Journal of Climate2251193 - 1207.  https://journals.ametsoc.org/view/journals/clim/22/5/2008jcli2358.1.xml 10.1175/2008JCLI2358.1
  * Schär . (2004) Schar2004Schär, C., Vidale, PL., Lüthi, D., Frei, C., Häberli, C., Liniger, MA. Appenzeller, C.  2004\.  The role of increasing temperature variability in European summer heatwaves The role of increasing temperature variability in European summer heatwaves.  Nature4276972332-336.  10.1038/nature02300
  * Seferian (2018) CNRM-ESM2-1Seferian, R.  2018\.  CNRM-CERFACS CNRM-ESM2-1 model output prepared for CMIP6 CMIP historical CNRM-CERFACS CNRM-ESM2-1 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4068
  * Seland . (2019) NorESM2-LMSeland, Ø., Bentsen, M., Oliviè, DJL., Toniazzo, T., Gjermundsen, A., Graff, LS.Schulz, M.  2019\.  NCC NorESM2-LM model output prepared for CMIP6 CMIP historical NCC NorESM2-LM model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.8036
  * Semmler . (2018) AWI-CM-1-1-MRSemmler, T., Danilov, S., Rackow, T., Sidorenko, D., Barbi, D., Hegewald, J.Jung, T.  2018\.  AWI AWI-CM1.1MR model output prepared for CMIP6 CMIP historical AWI AWI-CM1.1MR model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.2686
  * Seneviratne . (2012) SREXch3Seneviratne, S., Nicholls, N., Easterling, D., Goodess, C., Kanae, S., Kossin, J.Zhang, X.  2012\.  Changes in climate extremes and their impacts on the natural physical environment Changes in climate extremes and their impacts on the natural physical environment.  C. Field . (), Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation ( 109-230).  Cambridge, UK, and New York, NY, USACambridge University Press.
  * Seneviratne . (2021) IPCC_2021_WGI_Ch_11Seneviratne, S., Zhang, X., Adnan, M., Badi, W., Dereczynski, C., Di Luca, A.Zhou, B.  2021\.  Weather and Climate Extreme Events in a Changing Climate Weather and Climate Extreme Events in a Changing Climate [Book Section].  V. Masson-Delmotte . (), Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change ( 11).  Cambridge, United Kingdom and New York, NY, USACambridge University Press.  10.1017/9781009157896.013
  * Shin . (2022) Shin2022Shin, JY., Kang, M. Kim, KR.  2022\.  Outdoor thermal stress changes in South Korea: Increasing inter-annual variability induced by different trends of heat and cold stresses Outdoor thermal stress changes in South Korea: Increasing inter-annual variability induced by different trends of heat and cold stresses.  Science of The Total Environment805150132.  https://www.sciencedirect.com/science/article/pii/S0048969721052074 https://doi.org/10.1016/j.scitotenv.2021.150132
  * Sillmann, Kharin, Zhang . (2013) Sillmann2013aSillmann, J., Kharin, VV., Zhang, X., Zwiers, FW. Bronaugh, D.  2013\.  Climate extremes indices in the CMIP5 multimodel ensemble: Part 1. Model evaluation in the present climate Climate extremes indices in the CMIP5 multimodel ensemble: Part 1. Model evaluation in the present climate.  Journal of Geophysical Research: Atmospheres11841716–1733.  10.1002/jgrd.50203
  * Sillmann, Kharin, Zwiers . (2013) Sillmann2013bSillmann, J., Kharin, VV., Zwiers, FW., Zhang, X. Bronaugh, D.  2013\.  Climate extremes indices in the CMIP5 multimodel ensemble: Part 2. Future climate projections Climate extremes indices in the CMIP5 multimodel ensemble: Part 2. Future climate projections.  Journal of Geophysical Research: Atmospheres11862473–2493.  10.1002/jgrd.50188
  * Skelton . (2020) Skelton2020Skelton, A., Kirchner, N. Kockum, I.  2020\.  Skewness of Temperature Data Implies an Abrupt Change in the Climate System Between 1985 and 1991 Skewness of Temperature Data Implies an Abrupt Change in the Climate System Between 1985 and 1991.  Geophysical Research Letters4723e2020GL089794.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020GL089794 https://doi.org/10.1029/2020GL089794
  * Stone (2019) WRAF2019Stone, DA.  2019\.  A hierarchical collection of political/economic regions for analysis of climate extremes A hierarchical collection of political/economic regions for analysis of climate extremes.  Climatic Change1554639–656.  10.1007/s10584-019-02479-6
  * Swart . (2019) CanESM5Swart, NC., Cole, JN., Kharin, VV., Lazare, M., Scinocca, JF., Gillett, NP.Sigmond, M.  2019\.  CCCma CanESM5 model output prepared for CMIP6 CMIP historical CCCma CanESM5 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.3610
  * Tamarin-Brodsky . (2019) Tamarin2019Tamarin-Brodsky, T., Hodges, K., Hoskins, BJ. Shepherd, TG.  2019\.  A Dynamical Perspective on Atmospheric Temperature Variability and Its Response to Climate Change A dynamical perspective on atmospheric temperature variability and its response to climate change.  Journal of Climate3261707–1724.  10.1175/JCLI-D-18-0462.1
  * Tamarin-Brodsky . (2020) Tamarin2020Tamarin-Brodsky, T., Hodges, K., Hoskins, BJ. Shepherd, TG.  2020\.  Changes in Northern Hemisphere temperature variability shaped by regional warming patterns Changes in northern hemisphere temperature variability shaped by regional warming patterns.  Nature Geoscience1366414–421.  10.1038/s41561-020-0576-3
  * Tamarin-Brodsky . (2022) Tamarin2022Tamarin-Brodsky, T., Hodges, K., Hoskins, BJ. Shepherd, TG.  2022\.  A Simple Model for Interpreting Temperature Variability and Its Higher-Order Changes A simple model for interpreting temperature variability and its higher-order changes.  Journal of Climate351387–403.  10.1175/JCLI-D-21-0310.1
  * Tang . (2019) UKESM1-0-LLTang, Y., Rumbold, S., Ellis, R., Kelley, D., Mulcahy, J., Sellar, A.Jones, C.  2019\.  MOHC UKESM1.0-LL model output prepared for CMIP6 CMIP historical MOHC UKESM1.0-LL model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.6113
  * Tatebe Watanabe (2018) MIROC6Tatebe, H. Watanabe, M.  2018\.  MIROC MIROC6 model output prepared for CMIP6 CMIP historical MIROC MIROC6 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.5603
  * Taylor . (2022) Taylor2022Taylor, P., Boeke, R., Boisvert, L., Feldl, N., Henry, M., Huang, Y.Tan, I.  2022\.  Process Drivers, Inter-Model Spread, and the Path Forward: A Review of Amplified Arctic Warming Process drivers, inter-model spread, and the path forward: A review of amplified arctic warming.  Frontiers in Earth Science9.  10.3389/feart.2021.758361
  * Thompson . (2022) Thompson2022Thompson, V., Kennedy-Asser, AT., Vosper, E., Lo, YTE., Huntingford, C., Andrews, O.Mitchell, D.  2022\.  The 2021 western North America heat wave among the most extreme events ever recorded globally The 2021 western North America heat wave among the most extreme events ever recorded globally.  Science Advances818eabm6860.  10.1126/sciadv.abm6860
  * Twardosz Kossowska-Cezak (2013) Twardosz2013Twardosz, R. Kossowska-Cezak, U.  2013\.  Exceptionally hot summers in Central and Eastern Europe (1951–2010) Exceptionally hot summers in Central and Eastern Europe (1951–2010).  Theoretical and Applied Climatology1123617–628.  10.1007/s00704-012-0757-0
  * Vogel, Hauser Seneviratne (2020) Vogel2020Vogel, MM., Hauser, M. Seneviratne, SI.  2020\.  Projected changes in hot, dry and wet extreme events’ clusters in CMIP6 multi-model ensemble Projected changes in hot, dry and wet extreme events’ clusters in CMIP6 multi-model ensemble.  Environmental Research Letters159094021.  10.1088/1748-9326/ab90a7
  * Vogel, Zscheischler . (2020) Vogel2020aVogel, MM., Zscheischler, J., Fischer, EM. Seneviratne, SI.  2020\.  Development of Future Heatwaves for Different Hazard Thresholds Development of future heatwaves for different hazard thresholds.  Journal of Geophysical Research: Atmospheres1259e2019JD032070.  10.1029/2019JD032070
  * Voldoire (2018) CNRM-CM6-1Voldoire, A.  2018\.  CMIP6 simulations of the CNRM-CERFACS based on CNRM-CM6-1 model for CMIP experiment historical CMIP6 simulations of the CNRM-CERFACS based on CNRM-CM6-1 model for CMIP experiment historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4066
  * Voldoire (2019) CNRM-CM6-1-HRVoldoire, A.  2019\.  CNRM-CERFACS CNRM-CM6-1-HR model output prepared for CMIP6 CMIP historical CNRM-CERFACS CNRM-CM6-1-HR model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4067
  * E. Volodin . (2019) INM-CM5-0Volodin, E., Mortikov, E., Gritsun, A., Lykossov, V., Galin, V., Diansky, N.Emelina, S.  2019\.  INM INM-CM5-0 model output prepared for CMIP6 CMIP historical INM INM-CM5-0 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.5070
  * EM. Volodin Yurova (2012) Volodin2012Volodin, EM. Yurova, AY.  2012\.  Summer temperature standard deviation, skewness and strong positive temperature anomalies in the present day climate and under global warming conditions Summer temperature standard deviation, skewness and strong positive temperature anomalies in the present day climate and under global warming conditions.  Climate Dynamics405-61387–1398.  https://doi.org/10.1007/s00382-012-1447-4 10.1007/s00382-012-1447-4
  * von . (2019) INM-CM4-8von, E., Mortikov, E., Gritsun, A., Lykossov, V., Galin, V., Diansky, N.Emelina, S.  2019\.  INM INM-CM4-8 model output prepared for CMIP6 CMIP historical INM INM-CM4-8 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.5069
  * von Hippel (2005) gmm_mean_medianvon Hippel, PT.  2005\.  Mean, Median, and Skew: Correcting a Textbook Rule Mean, Median, and Skew: Correcting a Textbook Rule.  Journal of Statistics Education132null.  https://doi.org/10.1080/10691898.2005.11910556 10.1080/10691898.2005.11910556
  * Walt Fitchett (2021) Walt2021Walt, AJ. Fitchett, JM.  2021\.  Trend analysis of cold extremes in South Africa: 1960–2016 Trend analysis of cold extremes in south africa: 1960–2016.  International Journal of Climatology4132060–2081.  10.1002/joc.6947
  * Wang . (2021) Wang2021Wang, J., Guan, Y., Wu, L., Guan, X., Cai, W., Huang, J.Zhang, B.  2021\.  Changing Lengths of the Four Seasons by Global Warming Changing lengths of the four seasons by global warming.  Geophysical Research Letters486e2020GL091753.  10.1029/2020GL091753
  * Wang . (2016) Wang2016Wang, J., Han, Y., Stein, ML., Kotamarthi, VR. Huang, WK.  2016\.  Evaluation of dynamically downscaled extreme temperature using a spatially-aggregated generalized extreme value (GEV) model Evaluation of dynamically downscaled extreme temperature using a spatially-aggregated generalized extreme value (GEV) model.  Climate Dynamics479-102833–2849.  https://doi.org/10.1007/s00382-016-3000-3 10.1007/s00382-016-3000-3
  * Wehner . (2018) Wehner2018Wehner, M., Stone, D., Shiogama, H., Wolski, P., Ciavarella, A., Christidis, N. Krishnan, H.  2018\.  Early 21st century anthropogenic changes in extremely hot days as simulated by the C20C+ detection and attribution multi-model ensemble Early 21st century anthropogenic changes in extremely hot days as simulated by the C20C+ detection and attribution multi-model ensemble.  Weather and Climate Extremes201-8.  https://www.sciencedirect.com/science/article/pii/S2212094717301159 https://doi.org/10.1016/j.wace.2018.03.001
  * Weigel . (2021) Weigel2020Weigel, K., Bock, L., Gier, BK., Lauer, A., Righi, M., Schlund, M.Eyring, V.  2021\.  Earth System Model Evaluation Tool (ESMValTool) v2.0 – diagnostics for extreme events, regional and impact evaluation, and analysis of Earth system models in CMIP Earth System Model Evaluation Tool (ESMValTool) v2.0 – diagnostics for extreme events, regional and impact evaluation, and analysis of Earth system models in CMIP.  Geoscientific Model Development1463159–3184.  10.5194/gmd-14-3159-2021
  * Wieners . (2019) MPI-ESM1-2-LRWieners, KH., Giorgetta, M., Jungclaus, J., Reick, C., Esch, M., Bittner, M.Roeckner, E.  2019\.  MPI-M MPI-ESM1.2-LR model output prepared for CMIP6 CMIP historical MPI-M MPI-ESM1.2-LR model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.6595
  * Wilks (2002) Wilks2002Wilks, DS.  2002\.  Smoothing forecast ensembles with fitted probability distributions Smoothing forecast ensembles with fitted probability distributions.  Quarterly Journal of the Royal Meteorological Society1285862821-2836.  https://rmets.onlinelibrary.wiley.com/doi/abs/10.1256/qj.01.215 https://doi.org/10.1256/qj.01.215
  * Wu . (2018) BCC-CSM2-MRWu, T., Chu, M., Dong, M., Fang, Y., Jie, W., Li, J.Zhang, Y.  2018\.  BCC BCC-CSM2MR model output prepared for CMIP6 CMIP historical BCC BCC-CSM2MR model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.2948
  * Yosef . (2021) Yosef2021Yosef, Y., Aguilar, E. Alpert, P.  2021\.  Is it possible to fit extreme climate change indices together seamlessly in the era of accelerated warming? Is it possible to fit extreme climate change indices together seamlessly in the era of accelerated warming?  International Journal of Climatology41S1E952-E963.  https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/joc.6740 https://doi.org/10.1002/joc.6740
  * Yukimoto . (2019) MRI-ESM2-0Yukimoto, S., Koshiro, T., Kawai, H., Oshima, N., Yoshida, K., Urakawa, S.Adachi, Y.  2019\.  MRI MRI-ESM2.0 model output prepared for CMIP6 CMIP historical MRI MRI-ESM2.0 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.6842
  * X. Zhang . (2011) Zhang2011Zhang, X., Alexander, L., Hegerl, GC., Jones, P., Tank, AK., Peterson, TC.Zwiers, FW.  2011\.  Indices for monitoring changes in extremes based on daily temperature and precipitation data Indices for monitoring changes in extremes based on daily temperature and precipitation data.  WIREs Climate Change26851-870.  10.1002/wcc.147
  * Y. Zhang . (2022) Zhang2022Zhang, Y., Li, Q., Ge, Y., Du, X. Wang, H.  2022\.  Growing prevalence of heat over cold extremes with overall milder extremes and multiple successive events Growing prevalence of heat over cold extremes with overall milder extremes and multiple successive events.  Communications Earth and Environment3173.  10.1038/s43247-022-00404-x
  * Zhao . (2021) Zhao2021Zhao, Y., Qian, C., Zhang, W., He, D. Qi, Y.  2021\.  Extreme temperature indices in Eurasia in a CMIP6 multi-model ensemble: Evaluation and projection Extreme temperature indices in Eurasia in a CMIP6 multi-model ensemble: Evaluation and projection.  International Journal of Climatology41115368-5385.  10.1002/joc.7134
  * Ziehn . (2019) ACCESS-ESM1-5Ziehn, T., Chamberlain, M., Lenton, A., Law, R., Bodman, R., Dix, M.Druken, K.  2019\.  CSIRO ACCESS-ESM1.5 model output prepared for CMIP6 CMIP historical CSIRO ACCESS-ESM1.5 model output prepared for CMIP6 CMIP historical [dataset].  Earth System Grid Federation.  10.22033/ESGF/CMIP6.4272
  * Zittis . (2021) Zittis2021Zittis, G., Hadjinicolaou, P., Almazroui, M., Bucchignani, E., Driouech, F., El Rhaz, K.Lelieveld, J.  2021\.  Business-as-usual will lead to super and ultra-extreme heatwaves in the Middle East and North Africa Business-as-usual will lead to super and ultra-extreme heatwaves in the Middle East and North Africa.  npj Climate and Atmospheric Science4120.  10.1038/s41612-021-00178-7