Enhanced risk of hot extremes revealed by observation-constrained model projections

Recent decades have witnessed a surge of extreme hot events, some of which are unparalleled in the historical records or even deemed impossible without present-day climate change^{1,2,3,4,5,6,7}. While assessing the human influence on individual events remains challenging, enhanced levels of greenhouse gas (GHG) forcing undoubtedly raise the chance for the extremes to occur^8,9,10,11,12. However, as both observational and model evidence suggest, the probability of hot extremes grows nonlinearly with global warming and highly unevenly in space, reflecting complex interactions between intrinsic regional vulnerabilities and proximate physical drivers^13,14,15. The amplification mechanisms are manifold, varying by region and season, and may be tied to large-scale atmospheric changes (e.g., in circulation and thermal gradients)^{16,17,18,19,20,21} and feedback processes between land and the atmosphere likely triggered by long-term drying^{22,23,24,25,26,27}.

Climate models are essential tools for investigating the Earth system’s response to external forcing and provide insights into extreme event dynamics. Yet, models exhibit broad uncertainties stemming from the range in future emissions, inter-model structural differences, and internal climate variability, and despite progress, several unknowns remain^6,28. Furthermore, they often struggle to accurately capture observed regional climate trends, particularly in extreme events^29,30,31,32, undermining confidence in future projections. Alongside traditional bias correction approaches, the emergent constraint framework offers an effective route to reduce model uncertainties and fix known biases by exploiting insightful relationships between past and future physical observables^33,34,35. In this vein, recent studies have evidenced the potential of present-day thermal variability to constrain projections of hot event probabilities, suggesting their possible underestimation in several regions^36,37. The significance of high-frequency variability for extreme event trajectories is widely recognized, as it can sway the strength and timing of climate change signals by amplifying or buffering the effects of background warming^14,38,39. Additionally, growing evidence suggests that asymmetries in the tails of thermal distributions (skewness) can be remarkable in many places and considerably affect the occurrence frequency of the extremes at the regional scale^40,41,42.

Here we examine in depth the relationships between future changes in the probability of hot extremes and the historical properties of thermal distributions, unravelling the critical roles played by the higher-order moments (i.e., beyond the mean) and the implications for climate change projections. By theoretically modelling the fine structure of distributions and their detailed evolution, we explain the scaling behaviour of hot event probabilities with GHG warming and their heterogeneous changes across the global land area, identifying regions prone to greater risk in the coming decades. Hence, we leverage the tight connections between historical distribution moments and future hot-event probabilities to correct potential model flaws and improve the accuracy of scenario projections, blending insights from observations with global model results. Unlike common constraining strategies, this formalism exploits the relevant model information without relying on model selection or weighting criteria, while ensuring consistency with the observed climate conditions. Despite uncertainties in both simulation and observational data, our results indicate that current models may be overly optimistic, as the growth rate of hot extremes is likely to outpace predictions over much of the land surface, particularly in vulnerable regions. Crucially, due to the nonlinear nature of probability changes, surpassing the critical 2 ^∘C global warming threshold^6,7 could plunge many parts of the world into uncharted climate conditions, where extremes become the norm.

Scaling law of hot event probabilities

Under enhanced GHG forcing, hot extremes are expected to intensify rapidly in many areas. Figure 1a illustrates land changes in hot event probabilities (HEP) at 2 ^∘C of global warming (GW) relative to the early industrial level, as they result from the mean across a set of model projections contributing to the Climate Model Intercomparison Project Phase 6 (CMIP6⁴³, see Methods and Supplementary Table 1 for details). Specifically, for every model and grid point, the fractional changes in probabilities (FCHEP) are given by

where brackets denote the time averages at a fixed GW level (2 ^∘C in Fig. 1a) and over the early industrial period (1851–1900, H) respectively, and P is the total probability for daily TX anomalies (x) to exceed a large historical threshold (the 99th percentile, x_p), empirically estimated from the data (Methods). The analysis is focused on summer (June-August and December-February in the Northern and Southern Hemispheres, respectively), when the more intense TX warming and generally smaller temperature variability compared to other seasons increase the likelihood of hot extremes. Furthermore, since probability changes with the level of GW are relatively insensitive to the GHG emission trajectory¹⁴, the high-end Shared Socioeconomic Pathway (SSP5-8.5) is routinely used as the baseline of future climate projections.

a,b Land-based fractional changes in probabilities at 2 ^∘C of GW in the mean of model projections (a) and their theoretical representations in the full SN scheme (b). c Fraction (R²) of spatial variation of probability changes, in the model mean (red bars) and spread (black lines), theoretically explained with fully varying parameters (full SN), with constant variability and skewness (rigid shift, RS), and with constant variability and zero skewness (Gaussian shift, GS). d–f Multimodel mean of summer TX warming at 2 ^∘C of GW (d), historical TX variability (e) and skewness (f, unitless).

Full size image

As seen in Fig. 1a, future changes in probabilities vary significantly across regions, with the sharpest increases projected over low-latitude land, consistent with previous findings^39,44. Most notably, under 2 ^∘C of GW, parts of the Amazon and subtropical desert areas may experience a more than 30-fold increase in the number of hot extremes compared to historical levels. The roots of differential regional sensitivities lie in the complex interactions among nonuniform TX warming, the historical properties of daily anomalies, and their future changes. To unveil the key amplification factors, we analytically relate the flow of hot event probabilities across a fixed threshold Eq. (1b) to the structure of the underlying thermal distributions, namely

where

and z = (x − ξ)/ω. Here, the density function ({f}_{{{{rm{SN}}}}}) provides a skewed extension of the Gaussian ϕ, known as the skew normal (SN⁴⁵), that effectively captures the detailed properties of daily thermal anomalies^41,46. The factor Φ accounts for data asymmetries through the shape parameter α, while the scale and location (ω and ξ, respectively) result from a combination of the leading distribution moments (i.e., mean, variability, and skewness, as defined by Eq. (5), Methods). Parameters are allowed to vary over time following the evolution of TX anomalies, while the hot event threshold x_p is fixed (the 99th percentile) and determined by data distributions in the reference past.

Using Eqs. (2), (3), future changes in hot event probabilities can be theoretically predicted from those in the leading moments of TX distributions (see Methods). Results at 2 ^∘C of GW are shown in Fig. 1b. Here, the SN parameters are derived from model-simulated early industrial moments and their projected changes, namely the shifts in the TX mean (Fig. 1d), historical variability and skewness for 1851–1900 (Fig. 1e, f, respectively), and their changes (Supplementary Fig. 1). As seen in comparison with Fig. 1a, theoretical predictions closely match empirical probability changes and explain near completely their total spatial variation (R² ~ 96%, see the full SN bar in Fig. 1c and Eqs. (7a–c) in Methods). Similar results hold at higher levels of warming (Supplementary Fig. 2). The accuracy of theoretical predictions clearly proves the effectiveness of the SN-based approach, highlighting the importance of the fine structure of TX distributions. Historical moments play a central role as evidenced by the still large R² achieved ( ~ 92%, RS bar in Fig. 1c) if changes in variability and skewness are disregarded (α and ω in Eqs. (2), (3) held fixed to historical values), namely if TX anomalies are rigidly shifted upward with GW. Differences with theoretical predictions based on fully varying SN parameters (Fig. 1b) are moderate, and primarily stem from neglected increases in TX variability over tropical regions (Supplementary Fig. 1a). If, additionally, skewness is not accounted for (α = 0 throughout and ω constant), the prediction accuracy of Eqs. (2), (3) is degraded, with the explained spatial variation substantially reduced and a wider intermodel spread (R² ~ 65%, GS bar in Fig. 1c). However, due to the obvious inverse relation between past variability and future probability (Fig. 1e, a), simple rigid shifts of Gaussian anomalies can still broadly represent basic large-scale features of probability changes^14,39. For example, future changes are suppressed over the vast land masses of northern mid-to-high latitudes, showing the largest thermal fluctuations. Vice versa, the steepest increases in probability, as previously noted, are projected in the tropics, where daily variability is at its lowest.

On the other hand, the decreased accuracy of Gaussian-based predictions emphasises the relevance of skewness in explaining extreme event trajectories. As shown in Fig. 1f, model-simulated summer TX anomalies exhibit clear-cut deviations from normal behaviour, with left skewness prevailing over tropical and extratropical land areas and a switch to right skewness at higher latitudes. Previous studies have shown that non-Gaussianity is pervasive in atmospheric variables and likely an intrinsic characteristic of geophysical flows^40,47,48,49. While its physical origin remains uncertain, non-zero skewness in thermal anomalies may arise dynamically through horizontal advection mechanisms or from land surface processes, land-ocean contrasts, and complex topography, depending on regions and seasons^21,40,50,51. In practice, neglecting left skewness typically leads to an underestimation of hot event probabilities, as shorter-than-normal tails at the warm end of thermal distributions are associated with higher rates of change, as discussed below.

Role of historical climate conditions

The historical structure of thermal distributions is key to explaining the highly heterogeneous changes in hot-event probabilities across regions. To disentangle the detailed effects of the higher moments, the theoretical growth rates of hot extremes are visually displayed in Fig. 2 as a function of historical TX variability and skewness, at rising levels of GW. Specifically, following Eq. (6) (Methods), higher order changes are disregarded, while the TX mean is assumed to uniformly shift upward at the same pace as GW throughout. As shown in Fig. 2a–c, at 1 ^∘C of GW (i.e., close to current climate conditions), the rates of change in probability are higher for smaller variability and/or more negative skewness, then they steadily increase with warming until a maximum is eventually reached. Once this maximum is exceeded, hot extremes become as common as the mean and their probability continues to rise asymptotically, but at a slower rate. Because of nonlinearity, contour lines get denser for smaller historical moments, indicating a rapid acceleration of hot extremes with both increasing warming and decreasing parameters. This basically explains why, in small-variability, left-skewed distributions, even low levels of GW can make extremes commonplace. Thus, minor differences in model simulations of historical moments may lead to large deviations in probability changes, amplifying the spread of future projections.

a–c Theoretical rates of change in the probability of hot extremes as a function of historical variability and skewness, under the rigid shift approximation of Eq. (6) (Methods) and fixed levels of GW. The TX shifting velocity is held constant and equal to the GW rate for every point of the moment space to disentangle the effects of the higher moments. Symbols denote the pairs of model-simulated historical variability and skewness in the regional mean over pantropical land (30^∘N–30^∘S, grey squares) and northern midlatitude land (40–65 ^∘N, black squares).

Full size image

Figure 2 helps reveal the origins of enhanced regional vulnerabilities, identifying land areas inherently at higher risk under warming conditions. A prominent example is the tropical land region (30^∘N–30^∘S, represented by grey squares in Fig. 2), where the growth rate of hot extremes may even double with a 1 ^∘C increase in GW (Fig. 2a, b). In contrast, over northern midlatitude land (40–65^∘N, black squares) the growth rate is considerably slowed down by the larger historical variability and the near Gaussian behaviour of TX anomalies.

However, reality may be slightly more complex than depicted in Fig. 2, as these results strictly hold in the limit of rigid shifts of TX distributions at the same constant rate. In fact, the shifting velocity varies considerably across regions and, in many places, exceeds the rate of GW (Fig. 1d), amplifying hot event probabilities. Additionally, the higher moments are expected to change in future climate (Supplementary Fig. 1), with substantial increases in summer TX variability over tropical land and parts of midlatitudes, potentially further accelerating the growth of hot extremes. Evidence for such changes over the recent past has also been found in observational data, although their significance is a matter of some debate^52,53,54,55. From a physical standpoint, summer increases in TX variability may be triggered by soil-drying tendencies, through changes in the heat fluxes between land and the atmosphere, and associated feedback mechanisms^{22,25,26,27,56}. Moreover, advection-based arguments suggest that regional strengthening of summer mean temperature gradients, due to uneven warming of land and oceans, may lead to stronger TX anomalies and thus to increasing variability in some areas^20,21. Hence, higher-order changes, while of lesser importance compared to historical moments (Fig. 1c), may still provide regionally relevant contributions to future hot-event probabilities.

Constraining model projections

As argued above, the higher moments of historical TX distributions play a crucial role in shaping the increase in hot extremes, either amplifying or suppressing the effects of regional background warming (Fig. 2). How climate models represent real-world distributions may thus significantly impact future projections. Unfortunately, observational data coverage across the Earth’s surface is unequal^15,57, with limited and discontinuous data prior to 1950 in tropical regions, such as parts of the Amazon, continental Africa, and South Asia, complicating the evaluation of model performance. Here, we use available land-based observations to diagnose model biases in the early industrial past and improve future projections of hot event probabilities, with a focus on highly vulnerable regions. We primarily rely on the Berkeley Earth Surface Temperature (BEST) dataset⁵⁸, which provides regularly spaced TX observations at the daily scale, dating back to 1880. However, all data up to 1950 are included in the historical climate analysis to balance the gaps affecting the earliest record period and gain insights into tropical regions (see Methods). Thus, to ensure consistent comparison with observational results, model TX anomalies and all derived quantities – including historical moments, hot event thresholds, and probabilities for every model and grid point – are recalculated relative to the same 1880–1950 baseline used for observations. Given the near-stationarity prior to 1950, this forward adjustment to the model reference period has minimal impact. In particular, model-simulated moments for 1880–1950 (Supplementary Fig. 3a,b) and projected changes in probabilities relative to this period (Supplementary Fig. 4) are very similar, in the multimodel mean, to their equivalents based on 1851–1900 (Fig. 1e, f, a, and Supplementary Fig. 2a, respectively).

Summer TX variability and skewness from BEST data over 1880–1950 are shown in Supplementary Fig. 5a and Supplementary Fig. 6a, respectively, and compared in the zonal mean with their model counterparts in Fig. 3a, b. Models are found to considerably overestimate observed variability and, in some places, to underestimate skewness, overall suggesting that future increases in hot extremes might be faster than predicted over much of the global land. Hence, by relying on the SN analytical scheme (Eqs. (2), (3) with fully varying parameters), we combine the historical insights from observations with model projected changes. Specifically, for every model and grid point within the observed area, the higher-order moments are debiased by correcting historical mean values while retaining the future changes, as schematically shown in Fig. 3c. For instance, variability becomes

where brackets denote the historical values derived from BEST data and models (GCM) respectively, and δTX_SD are the moment’s model changes over time. Likewise for skewness, whereas for the TX mean the historical values are vanishing by construction in both models and observations. For simplicity, projected changes in the leading moments are preserved as they are, even though addressing potential model biases in simulated TX warming could in principle reduce future uncertainties and further improve results⁵⁹. The SN parameters are then computed from debiased moments (Eq. (4)) for each year of simulation, and the associated hot event probabilities are obtained, as with bare model moments (Fig. 1b), via Eqs. (2), (3), where the exceedance thresholds (x_p) are now drawn from the observed distributions consistently with historical moments (see Methods).

a, b Zonal-mean summer TX variability (a) and skewness (b) from BEST data (orange lines) and model simulations for the common historical period (1880–1950). Model results are shown as the mean (MM, thick black lines), total spread (thin black), and interquartile range (IQR, grey-shaded bands). c Conceptual drawing of the observation-constrained flow of thermal distributions with GW, compared to the unconstrained case. d, e BEST-constrained fractional changes in hot event probabilities at 2 ^∘C of GW over the observed area (d) and compared in the zonal mean (e, orange line) to bare model projections (shown as the MM, total spread, and IQR as above). Hot event thresholds and probability changes are referenced to 1880–1950. Red boxes in panel d highlight regions with enhanced increases.

Full size image

The changes in observation-constrained probabilities at 2 ^∘C of GW are shown as the model mean for each grid point in Fig. 3d, where blank areas denote insufficient real data coverage. As expected, since observed variability is systematically smaller than that simulated (Fig. 3a), probability increases are globally enhanced relative to bare model projections (Supplementary Fig. 4a), as evidenced by the zonal comparison in Fig. 3e and the grid-point differences in Supplementary Fig. 7.

Land areas showing a considerably heightened risk are highlighted in Fig. 3d (red boxes). For these regions, the full probability evolution with the increase in GSAT is displayed in Fig. 4a–f, where model-mean BEST constrained projections are compared with bare model results and their unconstrained SN representations. Constrained probabilities exhibit a steeper rise than projected in all regions, however the most impressive enhancement is foreseen over tropical America (TAM, Fig. 4d), particularly in the Amazon basin (Fig. 3d). This is primarily due to historical TX variability, whose model simulated value is so small that any further reduction can lead to a giant amplification of probability changes (Fig. 2). Based on historical observations, this region may thus experience the most dramatic increase in the number of hot extremes, with a nearly continuous heatwave state at less than 2 ^∘C GW (Fig. 4d), despite the slower rate of summer warming compared to northern midlatitudes. These findings have important implications for the Amazon forest system, since sustained hot daytime conditions can have detrimental effects on carbon stocks and fluxes, increasing the risk of large biome loss and forest dieback, potentially further exacerbating regional climate change^60,61.

a–f Regional aggregation of BEST constrained probabilities over the selected areas (Fig. 3d), as a function of changes in GSAT. For each region, the model mean of BEST constrained projections (orange line) is compared to the mean and spread of bare model results (thick black and thin black lines, respectively) and their unconstrained SN representation (in the model mean, red line). The latter shows an almost perfect match with bare model projections in all regions. Stochastic fluctuations in probabilities are filtered out using 40y running averages. Dashed lines at 0.5 probability mark the overshoot of critical climate thresholds, when hot extremes are projected to become the norm.

Full size image

Over the dry areas of the Sahara and Arabian desert (SAD, Fig. 4e) BEST corrections are milder and constrained probabilities, in the model mean, fall within the spread of bare model projections. This is because the imposed reduction in historical variability is partly offset by the larger, either near-vanishing or positive skewness, over much of this region (Supplementary Figs. 5a, 6a, 7). Furthermore, unlike bare model results, constrained probabilities grow faster in the TAM than the SAD region, consistent with observed changes in the recent past⁵⁷. Both models and observations, however, indicate the SAD region as one of the most jeopardized under GW⁶, where summer TX warming is increasing at the highest rates (e.g., Fig. 1d), substantially contributing to the acceleration of hot extremes.

In the other regions (Fig. 4a–c, f), amplification of probabilities over model projections is primarily driven by reductions in historical variability due to BEST constraints, with the largest corrections observed in the Mediterranean zone (MED) and Southeast Asia (SEA). In these regions, constrained results fall close to the upper edge of the bare model spread, indicating a nearly 20-fold increase in the number of hot extremes already under 2 ^∘C of GW, which is about twice as large as model-based expectations. However, as noted above, the limited observational coverage over tropical areas of Asia and South America makes constrained projections more uncertain. In Western North America (WNA) and Australia (AUS), the enhancement of constrained probabilities is smaller, yet still remarkable.

Except in tropical regions, historical skewness is broadly captured by models or, at worst, slightly underestimated (e.g., in the MED zone), and its corrections therefore generally produce minor changes.

Over northern high-latitude land, the effects of BEST constraints are mostly buffered by the large intrinsic variability and prevailing right skewness (Supplementary Fig. 3a,b and Fig. 2). Additionally, at high latitudes (above  ~ 50 ^∘N), observational uncertainty, as discussed further below, is too wide to draw any firm conclusions.

Assessing uncertainties

Model uncertainty in future hot event probabilities depends to varying degrees on the range in TX warming and the spread in the historical simulations of the higher moments, which may combine in complex ways. Observational constraints, by removing the latter spread, lead to uncertainty reductions in tropical areas where, as argued above (Fig. 2), adjustments to historical climate conditions have major effects on future probabilities and thus on their overall uncertainty. This is illustrated for the TAM and the SAD region in Fig. 5a, b, at relevant levels of GW. Conversely, at higher latitudes, owing to larger variability and a null-to-right skewness, removing the model spread in the historical moments has a lesser impact, and the contribution of the TX warming range is predominant, as seen in the MED zone (Fig. 5c). This range is particularly large over Australia and several areas of northern midlatitudes³⁷, inevitably limiting the accuracy of hot event projections, as shown in Supplementary Fig. 8 for other regions (WNA, SEA and AUS) and the global land.

a–c Uncertainties in BEST-constrained projections of probability changes (FCHEP, green boxes) compared to those in bare model results and their unconstrained SN representations (white and grey boxes, respectively), for selected regions and future levels of GW. Boxplots represent the 10th to 90th percentile range of model projections, the inside bars indicate the median, and whiskers the total spread.

Full size image

In turn, constrained model results closely depend on the reliability of observed distribution moments, which becomes increasingly challenging to assess as we move further back into the early industrial era. Nevertheless, some qualitative insights into the level of observational uncertainty can be gained from comparisons of TX moments over the 1880–1950 period (referred to as early H henceforth) with those over the later 1951–2000 period (late H), when daily data availability is greater. Specifically, Fig. 6 compares zonal-mean summer variability and skewness from BEST data with those from the NOAA/CIRES/DOE 20th Century Reanalysis (20CRv3⁶²) over both the early and late H periods, as well as with late H moments from the HadGHCND gridded observations⁶³ (Methods). Additional results for the late H period from the JRA-3Q⁶⁴ and ERA5⁶⁵ reanalyses are displayed in Supplementary Fig. 9a,b.

Zonal-mean summer TX variability (a) and skewness (b) from the BEST, HadGHCND, and 20CRv3 datasets (orange, red, and violet lines, respectively), compared over the late (1951–2000, solid lines) and early (1880–1950, dashed lines) historical periods. HadGHCND data exhibit gaps in coverage at tropical latitudes (e.g., Supplementary Fig. 5d). Model-simulated moments for the early historical period are also shown as the multimodel mean (MM, thin black lines) and interquartile range (IQR, grey-shaded bands). Observational and model results were remapped to a common resolution (N32) before averaging over latitude circles.

Full size image

As shown in Fig. 6a and Supplementary Fig. 9a (and Supplementary Fig. 5c–g for grid point comparisons), in the late H period, TX variability patterns are generally consistent across datasets over midlatitudes and the northern tropics, where observational coverage is denser, less so at higher latitudes. In particular, BEST and HadGHCND data exhibit close agreement across most regions, while reanalyses often show larger discrepancies, with either positive or negative biases relative to the observations (Supplementary Fig. 9a), in line with previous findings⁶⁶.

Results from the longest datasets, BEST and 20CRv3, indicate an overall increase in TX variability during the late H period compared to the early H, albeit with differences in the magnitude of changes. Whether this increase reflects a true climate signal or improved data coverage over time remains to be assessed. Notably, despite this increase, TX variability in the late H period, from both observations and reanalyses, still falls mostly below the range of model simulations for the early H period. This picture again suggests a systematic overestimation of variability by climate models. For instance, over midlatitudes and the northern tropics, BEST late-H variability, though larger than its early-H counterpart, is about 0.5 ^∘C smaller, in the zonal mean, than the multimodel early-H result (Fig. 6a). In contrast, model TX variability remains constant throughout the historical simulations (Supplementary Figs. 3a, c, 9c).

The observed variability estimates, taken together, provide a comprehensive uncertainty range, with late H values representing plausible upper bounds for early H variability. This range can be used to assess the impact of observational uncertainty on constrained probability changes. As shown in Fig. 7 (and Supplementary Fig. 10), since larger variability implies smaller growth rates of hot extremes (Fig. 2), using observed moments from the late rather than the early H period to constrain projections generally results in slower probability changes. Yet, these changes remain faster than model-based expectations in most areas. For example, in the TAM region (Fig. 7d), hot extremes may become commonplace at GW levels well below model predictions, according to both BEST and 20CRv3 data. Similarly, in several regions, such as MED and SEA, constrained results based on multiple observational datasets consistently suggest enhanced changes in probabilities. BEST- and HadGHCND-constrained projections show a particularly close match in areas with full data coverage, such as WNA and AUS. Future changes, however, are affected by substantial uncertainty, which further widens when including results from additional reanalysis data, as shown in Supplementary Fig. 10. Nonetheless, constrained probability changes generally exceed bare multimodel projections, indicating that models likely underestimate the increasing rates of hot extremes.

a–f Scaling of constrained hot-event probabilities with GW, based on multiple observational estimates of TX moments, for the same regions as in Fig. 4. Orange, red, and violet lines represent the model mean of probabilities constrained by BEST, HadGHCND, and 20CRv3 data, respectively, while solid and dashed lines indicate the observational period (1951–2000 and 1880–1950, respectively) for TX moments, as in Fig. 6. Grey shading denotes the uncertainty range outlined by constrained probabilities based on BEST late and early historical moments. HadGHCND-based results are displayed only for regions with full (WNA and AUS) and near-full data coverage (MED). Bare model projections are shown as the mean (MM, black lines) and total spread (thin grey lines).

Full size image

The origins of model biases in historical variability, and thereby in hot event probabilities, are currently unclear. Misrepresentation of land surface processes could play a role in tropical regions³⁶, although several factors may contribute at the global scale, including horizontal grid resolution, physical parametrizations, and biases in simulated atmospheric circulation patterns, among others, needing further investigation^15,31. At high latitudes (above  ~ 50 ^∘N and over Antarctica), as previously noted, observational discrepancies in variability become more pronounced, reducing the reliability of historical data for constraining future projections.

Additionally, historical TX skewness exhibits a complex regional structure (Fig. 6b and Supplementary Figs. 3b, d, 6, 9b, d). Both models and observations similarly capture certain distinctive features, such as the reversal of the skewness sign at northern latitudes and the overall prevalence of left skewness down to the tropics, although considerable local differences may affect probability changes. For instance, 20CRv3 data show particularly low skewness in regions like the Sahara-Arabian area and Australia, throughout the historical period (Supplementary Fig. 6). In these regions, the associated constrained probabilities increase at the highest rates (Fig. 7), further highlighting the relevance of skewness (Fig. 2). In contrast, the agreement between HadGHCND and BEST results for the late H period is notably high (Fig. 6b), except in the northern tropics due to the lack of HadGHCND data in the Sahara region. Here, BEST data suggest a substantial decrease in skewness from the early to the late H period, which is not corroborated by 20CRv3 results and may reflect undetected flaws in observational data.

Caution is thus required when interpreting constrained projections, as they may inherit potential deficiencies from the reference historical dataset. Nevertheless, extensive assessment of observations and reanalysis data clearly reveals model misrepresentations of historical climate conditions, strongly suggesting that hot event probabilities could rise more rapidly than expected in the coming decades and surpass critical thresholds sooner in highly vulnerable regions.

On the theoretical side, errors arising from SN related calculations (namely, due to approximations in parameter fitting and numerical integrations as well as to inherent SN limitations, see Methods) are virtually negligible compared to both model and observational uncertainties. Indeed, as shown in Fig. 1a–c, Supplementary Fig. 2, and Fig. 4, the SN-based approach offers compelling theoretical representations of hot event probabilities under GHG warming, both regionally and globally, outperforming simple Gaussian-based approximations. This framework thus provides an effective paradigm for addressing historical model biases and constraining future projections with observations.

Comprehensive analyses of observational data indicate that climate models consistently overestimate high-frequency thermal variability over the recent past across all latitudes, while showing minor regional biases in skewness. The factors behind these flaws remain unclear and warrant further investigation. Alongside skewness, historical variability is critical for realistically predicting the growth rates of hot extremes, and its misrepresentation can lead to a significant underestimation of future risks.

To address this issue, we designed an SN-based theoretical approach whereby changes in hot event probabilities with GW and their uneven spread across regions can be accurately explained through the leading moments of thermal distributions. Historical climate conditions, particularly variability, play a key role in shaping the increasing rates of hot extremes, either amplifying or suppressing the effects of regional background warming. Based on universal analytical relationships between historical distribution moments and future hot-event probabilities, the SN framework provides a robust pathway to correct model errors using real-world data, resulting in more reliable projections. Furthermore, it leverages the wealth and heterogeneity of climate model information by retaining all model outcomes and constraining them with observations.

Observationally constrained projections suggest that, in the coming decades, hot event probabilities may grow faster than models imply, over much of the global land. For instance, in several regions like the Mediterranean and Southeast Asia, future increases could exceed bare model predictions by nearly twofold, although considerable model and observational uncertainties remain. Moreover, results provide valuable insights into highly vulnerable areas, such as the Amazon, where even a 2 ^∘C rise in GW could lead to near-continuous heatwave conditions. In tropical regions, however, limited data coverage adds further uncertainty to constrained projections. Progress in the development of long-term observational datasets with global reach may help reduce uncertainties and draw more robust conclusions.

Given these findings, it is crucial to address the risks associated with underestimating future increases in hot extremes, as their frequency could become unbearably high in large areas even under low levels of GW. Urgent action to curb GHG forcing is essential to prevent critical climate shifts and protect the most vulnerable regions from the harmful impacts of climate change.

Data sources and processing

This study is based on a set of near-surface maximum temperatures (TX) at the daily scale, taken from the latest global climate simulations (CMIP6) under the highest forcing scenario SSP5-8.5. The models are listed in Supplementary Table 1, alongside the approximate timings for global-mean surface air temperature (GSAT) to surpass future levels of global warming (GW) relative to the early industrial era (1851-1900). The GW timings were estimated for each model using 20y running averages.

We focus on the global land area to allow for comparison between model simulations of historical climate and results from regularly spaced daily observational datasets. Among these, the Berkeley Earth Surface Temperature dataset (BEST⁵⁸) serves as our primary reference, offering extensive spatiotemporal coverage, with land-based TX observations since 1880 interpolated on a regular 1-degree grid. Given the experimental stage of the daily BEST dataset, we assess observational uncertainty using additional gridded products, namely the Met Office Hadley Center daily dataset (HadGHCND⁶³), providing land-based TX observations since 1950, and, as proxies for real data, the latest NOAA/CIRES/DOE 20th Century Reanalysis (20CRv3⁶²), with daily TX dating back to 1836, the JRA-3Q⁶⁴, and ERA5 reanalyses⁶⁵.

We initially removed the annual cycle from both model and observational data at each grid point by calculating TX anomalies relative to the local day-of-the-year normals over the historical period. For models, the reference period was set to the early industrial era (1851–1900) to align with common standards^6,7 and facilitate comparisons with related results, such as projected changes in hot event probabilities at fixed GW levels (Fig. 1a and Supplementary Fig. 2a). However, for evaluation of model performance, both modelled and observed TX anomalies as well as all derived quantities are referenced to the 1880–1950 period, where BEST data provide adequate coverage. Given the substantial stationarity prior to 1950, this forward adjustment to the reference period has minimal impact on model results.

From the BEST dataset, we retained all grid points except those with less than 5% of the total days across the baseline period (Supplementary Fig. 5a). In fact, most of the selected area has high or near-full data coverage, except for certain low latitude regions (e.g., the Amazon basin, southern Arabian Peninsula, and tropical Asia), where data are scarcer, and results are therefore more uncertain. A similar approach was used to handle HadGHCND data for uncertainty assessment, but with reference to a later period (1951–2000).

The behaviour of daily TX anomalies in both models and observations is characterized by the leading statistical moments, namely

where E[  ] denotes the expectation operator. Besides the mean μ₁, distribution variability (standard deviation) and skewness are derived from Eq. (5), by letting (,sigma =sqrt{{mu }_{2}},) and (,gamma ={mu }_{3}/{mu }_{2}^{3/2}), respectively. Common empirical estimators were then used to compute moments from the data⁶⁷.

For models, calculations were performed separately for each grid point on the native grid and for each summer subsample of TX anomalies across all years of historical and future simulations, to track the detailed evolution of thermal distributions with GW. We parallelly defined the probability of hot extremes (Eq. 1b) for each model grid point and year, as the empirical frequency for summer TX anomalies to exceed a fixed climatological threshold, given by the 99th percentile of the historical distributions. Broadly speaking, the selected events are expected to occur nearly once per summer in a stationary climate, with magnitude varying widely across regions.

For BEST data, grid-point historical distributions and the associated hot-event thresholds were obtained from the summer TX anomalies within 1880–1950, consistent with climate normals. For models, hot event probabilities were calculated with reference to distributions and thresholds over both the early-industrial 1851–1900 and the common 1880–1950 periods. Similarly, historical variability and skewness for both models and observations were derived from Eq. (5) using all TX anomalies within the respective reference periods.

Finally, model-projected patterns of change in moments and probabilities at fixed GW levels were computed from long-term averages (20y) centred around the model-specific years associated with the given increase in GSAT (Supplementary Table 1). Native grid-point changes as well as historical moments were later interpolated on a common coarse-grained grid (N32) to allow for averaging over models and latitude circles.

Hot event probabilities in the SN-based approach

Theoretical hot-event probabilities given by Eqs. (2), (3) rely on the continuous SN distribution⁴⁵, which provides an analytical representation of the structure of daily TX anomalies, including deviations from normality. The SN probability density incorporates non-zero skewness through the real parameter α, while retaining the mathematical tractability of the normal density, encompassed as a special case (α = 0). Parameters (ξ, ω and α) were fitted to model data at the native grid points and for each simulation year, based on empirical estimates of the leading moments (μ₁, σ and γ) of summer TX anomalies, obtained as described above. The skills of the SN density and, in particular, the benefit of accounting for non-zero skewness in a theoretical representation of thermal distributions were thoroughly assessed in previous studies^41,46.

The flow of hot event probabilities across the historical fixed threshold (x_p) was then calculated for every model and grid point using Eqs. (2), (3), i.e., by integrating the time-varying SN densities associated with summer TX distributions. The algorithm is the same for observation-constrained model probabilities, except that the sample estimates of distribution moments were debiased according to Eq. (4) before computing the SN parameters, and the hot event thresholds were consistently drawn from the historical observations. The observed moments and thresholds were remapped to each model native grid to allow for bias correction via Eq. (4). As for bare model projections, time-varying theoretical probabilities (both constrained and unconstrained) were averaged over regions and models, and finally re-expressed as a function of changes in GSAT (Fig. 4). Patterns of change at fixed levels of GW were obtained as described above for empirical probabilities, after remapping grid point results at the common resolution (N32).

Due to the SN flexibility, the scaling behaviour of hot event probabilities with GW can be explored across a range of climate conditions and changes, unravelling the relative roles of TX moments. In particular, if changes in the higher moments are disregarded, the growth rate of theoretical probabilities (Eq. (2)) becomes

where z = (x_p − ξ)/ω_H, α_H and ω_H are the historical shape and scale parameters respectively, erf is the error function, and T the global temperature change. Eq. (6) provides a valuable analytical approximation for investigating the nonlinear relationships between historical TX moments and future hot-event probabilities and was used to generate the results presented in Fig. 2.

The accuracy of the SN formalism is assessed through its performance in explaining the global patterns of change in model hot-event probabilities, at fixed levels of warming. Following previous work¹⁴, we used the coefficients of determination defined by

where y_i denotes the model probability change at grid point i, y_i,th the theoretical change (in the SN scheme), and w_i the cell area fraction, with the index i running over the total number of grid points. S_tot is the total spatial variation of model changes around the global mean (bar{y}) and ({S}_{{{{rm{res}}}}}) is the fraction unexplained. A higher R² indicates better accuracy of theoretical predictions (see Fig. 1c).

Finally, there are some caveats to the current approach, stemming from inherent limitations of the SN density. Specifically, the SN kurtosis excess is positive by construction (fixed by absolute skewness), and only moderate levels of skewness, falling within (-1,1), can be represented. However, TX anomalies rarely exhibit skewness outside this range (Fig. 1f and Supplementary Fig. 6). Moreover, any slight discrepancies between empirical and theoretical kurtosis at certain grid points have minimal effects on hot event probabilities and their changes (Fig. 1c). These potential shortcomings along with theoretical uncertainties, arising from parameter fits and numerical integrations, appear negligible in this context, thus making the SN density a suitable analytical representation of TX anomalies.