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. Increasing levels of atmospheric carbon dioxide (CO 2 ) concentration unequivocally transformed the earth's climate (IPCC, 2021). This surplus of CO 2 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 et al., 2018). The effects of anthropogenic global warming led to the emergence of heat extremes that would not have occurred previously (Robinson et al., 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 et al., 2015;Thompson et al., 2022). The latter was found to be a remarkable four standard deviations away from the mean (Thompson et al., 2022). The Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report (AR6) concluded that human influence on the climate system is unequivocal (Eyring et al., 2021) and virtually certain to be the main driver of the changes in hot and cold extremes (Seneviratne et al., 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 et al., 2012;Rahmstorf & Coumou, 2011;Kharin et al., 2013;Sillmann, Kharin, Zhang, et al., 2013;Sillmann, Kharin, Zwiers, et al., 2013;Pfleiderer et al., 2019;Perkins-Kirkpatrick & Lewis, 2020;Vogel, Hauser, & Seneviratne, 2020;Raymond et al., 2020;Seneviratne et al., 2021;Mallick et al., 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 et al., 2008;Garcia-Leon et al., 2021;Demiroglu et al., 2020;Perera et al., 2020;Seneviratne et al., 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 et al., 2018;Y. Zhang et al., 2022). Iyakaremye et al. (2022) 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. Iyakaremye et al. (2021) 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 et al., 2016;Zittis et al., 2021;Ozturk et al., 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 et al., 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 et al., 2011;;Schar et al., 2004;Barriopedro et al., 2011). Over the Americas, the dry and hot extremes showed an increase both in frequency and spatial scope over the past 122 years (Alizadeh et al., 2020;Cai et al., 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;Schar et al., 2004;Anderson & Kostinski, 2010;Fischer & Schar, 2010;Barriopedro et al., 2011;C. Li et al., 2019;Loikith & Neelin, 2019). Thompson et al. (2022) 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 et al., 2016;Lewis et al., 2017;Vogel, Hauser, & Seneviratne, 2020;Vogel, Zscheischler, et al., 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 et al., 2011;Alexander et al., 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 et al., 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 et al., 2013;Wang et al., 2016;Paciorek et al., 2018;Wehner et al., 2018;Ben Alaya et al., 2020;C. Li et al., 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 et al., 2016;Ben Alaya et al., 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 et al., 2020;Wang et al., 2016). Ben Alaya et al. (2020) 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. Kollu et al. (2012) 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 Shin et al. (2022) 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 Clarkson et al. (2022). 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 (E. M. Volodin & Yurova, 2012;Perron & Sura, 2013;Kodra & Ganguly, 2014;Sardeshmukh et al., 2015;Linz et al., 2018;Tamarin-Brodsky et al., 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 et al., 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 et al., 2015;Wang et al., 2016;Knoben et al., 2019;Ben Alaya et al., 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 et al., 2015). A distribution can be described by not only the mean and the standard deviation but also skewness and kurtosis. Donat and Alexander (2012) 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. Sardeshmukh and Sura (2009) found a parabolic relationship between kurtosis and skewness that cause the non-Gaussianity of the observed daily weather anomalies. Similarly, Tamarin-Brodsky et al. ( 2022) 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 6 (CMIP6) simulations. This choice was Wilks, 2002;Donat & Alexander, 2012;Cho & Jeong, 2016;Bertossa et al., 2021) and was verified by applying the same analysis to the European Centre for Medium-5th Generation (ECMWF-ERA5) data for the same historical time period . 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. For this study, we used multi-year daily near-surface maximum temperatures from the 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 et al., 2014;Eyring et al., 2016;O'Neill et al., 2016). Additionally, the ECMWF-ERA5 dataset was included for the 30-year time period (Hersbach et al., 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 et al., 2022), similarly to Hajat et al. (2022) and Ribeiro et al. (2022). Therefore, we obtained the start and end years of 20-year GWL periods for each CMIP6 simulation from Hauser et al. (2022). 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. (Tang et al., 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 et al., 2020). All data extraction and prepro-cessing in this study were performed by using the (ESMValTool) version 2.5.0, which is an open-source software package for analysing and evaluating model simulations (Eyring et al., 2020;Lauer et al., 2020;Righi et al., 2020;Weigel et al., 2021). We extracted the daily maximum near-surface air temperature from each model for each region using shapefiles provided by IPCC (Iturbide et al., 2020), converted units from Kelvin to Celsius, and created a single spatiotemporal (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. 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 et al., 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 S3. 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: where K is the number of Gaussian components in the mixture. k ,  k and  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 et al., 2021;Eyring et al., 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 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). 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 multimodel mean percentages of grid cell modalities of a region as shown in Then, the parameters of the hot Gaussian component,  historical hot ,  historical hot and  historical hot , 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  +/- d range is defined as where 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. The return period of an event describes the average time between the occurrences of a certain event of a defined size, where an nyear 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 180 days of cold weather and 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 yearsx180 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 D denote the number of days in L years. Then, a "year " in the historical period, | N ( historical hot ,  historical hot ) | is defined as where  historical hot is the mean,  historical hot is the standard deviation and  historical hot 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, The standard deviation distance of range, d historical n , for an extreme event in the historical period can be calculated by using Equation 6, where erf 1 is inverse error function. Now, we can calculate a temperature threshold,  historical n , for an n-year event in the historical period. Finally, the new value of the return period in the future n, i.e. n-year event, is calculated by using Equation 8n = where )| 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 T , as the change in the difference between the means of cold and hot Gaussian components as shown in Equation 15: 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, T cold and T hot will always be positive. Therefore, a negative 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 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. 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 (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 We also analysed the movements of the Gaussian components relative to each other using the T definition from Equation 15in grid cells with a bimodal distribution. Figure 7 shows the 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 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, 5year, 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 yearsx188 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 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 observedcompared 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 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 T values as shown in 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 (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, 2019;Walt & Fitchett, 2021;Prodhomme et al., 2022). For example, Qian and Zhang (2015) 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 et al., 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 et al., 2011;Zhao et al., 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 et al., 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 et al., 2016;Wang et al., 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 et al., 2021). Seasonality in temperature extremes adds complexity to the process of selecting percentiles to define extreme temperatures (Huang et al., 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 et al., 2015;Wang et al., 2016;Knoben et al., 2019;Ben Alaya et al., 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 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 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 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: pacal23jgr GaussianMixtureModels Extremes. The regional output files amount to hundreds of GB. The latest release of ESMValTool is publicly at ESMValTool (Andela et al., 2022). 855187). This work used resources of the Deutsches Klimarechenzentrum (DKRZ) We acknowledge the 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 (ESGF) for archiving the data and providing access, and the multiple Hersbach et al. (2018) was downloaded from the (C3S) (2023). The results contain modified 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. manuscript submitted to JGR: Atmospheres 855187). This work used resources of the Deutsches Klimarechenzentrum (DKRZ) We acknowledge the 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 (ESGF) for archiving the data and providing access, and the multiple Hersbach et al. (2018) was downloaded from the (C3S) (2023). The results contain modified 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. manuscript submitted to JGR: Atmospheres