Future changes in the frequency of temperature extremes may be underestimated in tropical and subtropical regions

In a warming world, temperature extremes are expected to show a distinguishable change over much of the globe even at 1.5 °C warming, and in many regions this change has already been detected in observations. Although many studies predict an increase in heat extreme events, the magnitude of the change varies greatly among different models even for the same mean warming. This uncertainty has been linked to differences in land–atmosphere feedback across models. Here we show that a significant constraint for future projections can be based on the ability of climate models to accurately simulate the present day variability of daily surface maximum temperature. An emergent constraint on Coupled Model Intercomparison Project Phase 5 (CMIP5) and 6 (CMIP6) models, applied to ERA5 reanalysis, indicates that the best estimate in hot extreme changes by the end of the century could be worse than previously estimated, mostly for tropical and subtropical regions as well as South and East Asia. Present-day variability of daily maximum temperature at the surface can serve as an emergent constraint on the frequency of heat extremes in a warmer climate, according to tests with a large multi-model ensemble of climate models.

T emperature extremes are expected to show a distinguishable change over much of the globe 1 , even at 1.5°C warming 2 , and in many regions this change has already been detected in observations [3][4][5][6] . Change in extremes events remain inconsistent between models, even at a same level of mean global warming 7 in CMIP5 8 and CMIP6 9 models. These extremes have a strong impact on society and can have negative consequences on health 10 , agriculture 11 , or water resources 12 . Daily maximum temperature (TX) is often used to measure heat wave intensity. It is governed by many processes, including solar radiation, heat transport, and sensible, and latent heat flux exchange with the surface 13 . In particular, energy used to evaporate surface moisture can mitigate atmospheric warming and thus TX 14 . At any given location, TX variability tends to be larger under drier surface conditions than wetter conditions. Another way to formulate this idea is that soil moisture deficit (and deficit in other surface humidity variables) can lead to amplified TX (and with it, potentially amplified heat waves) and can explain part of TX variability. There is evidence that many current climate models are too dry under the present conditions 15 and we hypothesize that this amplifies TX variability and with it, heat wave frequency 16 , whereas more accurate models may see this amplification in upcoming decades. We postulate that this could lead to large differences between models in terms of heat wave changes under climate warming. Previous work has indicated a reduction in temperature change uncertainties when using prescribed land heat fluxes 17 and has found a systematic increase in interannual variability of summer temperatures over Europe with models that realistically represent variability 18 . Although our main focus here is on surface heat fluxes, many other processes can also impact TX variability locally (as, e.g., dynamics or aerosols) and model parametrization. Thus, we do not expect to explain systematically the temperature signal from heat fluxes alone.

Results
Metrics and techniques. Many indices of TX can be used to describe hot events (some as defined by the Expert Team on Climate Change Detection and Indices). We chose a simple derived index that can be applied easily at the global scale, namely the number of days above the 98th percentile (TX98p, see "Methods" for detailed computation). We focus only on the warmest season (June to August for North Hemisphere, December to February for South Hemisphere, and all year for the 15°S-15°N tropical area). TX98p indicates for each location the number of days that are considered as extremely hot (relative to the 1995-2005 daily climatology of TX at this location) and we evaluate its change in climate projections (i.e., the change in the number of days above the climatological threshold, see "Methods" for details). We also define a metric to quantify the historical variability of TX at each location, ΔTX. This metric indicates, at each grid point and for each calendar day, the distance between mean TX and the 95th percentile of TX in°C. ΔTX gives an indication of the temperature difference between a hot day compared to the climatology. It is used to evaluate models against a global reference dataset, the ERA5 reanalysis 19 . This difference between hot and average days has been found to be too high in some climate models 20 . Computation of ΔTX implies that we ignore any bias in the mean TX of a model (compared to ERA5) and focus only on TX variability. It is noteworthy that our results are not sensitive to the exact definition of the heat metric. If using another threshold, e.g., the 95th percentile instead of 98th percentile (see Supplementary Fig. 9), results are very similar.
The role of surface humidity. Previous studies have shown that soil moisture deficits enhance surface temperature extremes [21][22][23] as demonstrated, e.g., for European extremes 24 . Here we focus on daily temperatures, and due to limited availability of humidity model outputs at daily temporal resolution (especially integrated soil moisture is not available at daily timescales for CMIP5 outputs), we use the surface heat fluxes (latent and sensible) as an indicator of land-atmosphere interaction. Latent heat flux can be a good indicator of evaporation during hot days. We verify that latent heat flux anomalies during days above the 98th percentile exhibit a negative correlation with ΔTX over 80% of the land ( Supplementary Fig. 1a), i.e., models with high ΔTX have less evaporation. Sensible heat flux ( Supplementary Fig. 1) has positive correlations almost everywhere (indicating that hot days occur under warmer land surface conditions), except over the driest deserts. In these regions, the air may be so warm that the sensible heat flux is actually reduced. It is also verified that models with the largest change in TX98p have stronger decreases in latent heat flux ( Supplementary Fig. 1c) over tropical/subtropical areas and South and East Asia, i.e., they are drying more compared to other models between the present time and the end of the century. This relationship is not observed (or even reverses) for midlatitudes, corresponding to regions where ΔTX-TX98p is not valid as an emergent constraint (EC) and masked in the EC analysis. Other constraints, e.g., related to precipitation change, may be valid there 24 . These results confirm the relationship between surface humidity and TX variability during the baseline period, especially over tropical areas. A noticeable exception are the desert regions of North Africa and Middle East, where we also observe a strong ΔTX-TX98p change relationship but it is apparently not related to land-drying processes. In these regions, surface shortwave radiation (related to aerosol, cloud cover, or dynamics) could play a major role to explain ΔTX (Supplementary Fig. 2).
An EC to improve projection in extreme temperature events. We verify where ΔTX is correlated with the TX98p change for different warming targets (end of century or +1.5 and +2°C warming periods). Over most of land, the relationship between ΔTX and TX98p change is negative (Fig. 1) and significant ( Supplementary Fig. 3), indicating that in regions with presentday overestimated variance for hot days, the future change in TX98p is smaller on average. Thus, this simple metric is justified, both physically and statistically, to constrain model projections. In the following we apply the EC only where the ΔTX-TX98p correlation is significant. This is the case for tropical and subtropical areas and South and East Asia mainly. Mid-and highlatitude show less significant correlations.
The EC methodology requires understanding and accounting for observational and model variability and uncertainties to check consistency 25 . We use the internal variability of a large multimember historical ensemble (HAPPI 26 ), which is forced with observed sea surface temperatures to estimate the uncertainty in the ΔTX variability associated with weather variability at each location (model internal variability/error). As an estimate of observational error, we use the difference between ERA5 and the blended satellite surface CHIRTS 27 dataset, both of which provide high-resolution daily output for TX. We then consider these two pieces of information as an uncertainty range for ΔTX based on ERA5 and to evaluate when models fit within this range (with multi-member models having narrower uncertainty, see "Methods") and select only these models to simulate future change in heat wave metrics.
Regional and global constraint. We find that changes in TX98p are larger than estimated by an unconstrained ensemble over a large part of tropical and subtropical regions when using ΔTX to constrain climate projections by selecting at the regional scale (see "Methods") the models within the observed constraint (Fig. 2). Africa, South and Central Asia, and South America have a particularly strong signal, locally above 50% increase in the number of exceedances of the 98th percentile, i.e., twice as many hot days as in unconstrained predictions. This means models that represent ΔTX more accurately during the baseline period (and thus, hypothetically, humidity feedbacks) tend to warm hot extremes faster compared to the other models. Similar relationships are found for all climate warming targets (Supplementary Figs. 10 and 11), although the area with significant correlation is reduced for a 1.5°C target. Thus, the influence of our EC persists through different warming targets and is confirmed robustly by several sensitivity tests. It is also able to improve the prediction of any climate model's future change in the heat wave metric based on the ensemble, passing a key test for ECs (perfect model test; see ref. 25 and "Methods"). As an uncertainty, we note that EC results are the strongest where only a few models are selected as realistic ( Supplementary Fig. 4a). Thus, the magnitude of the amplification of heat wave projections may be less robust over these area, as it relies on a more limited number of models. Applying an EC based on global mean ΔTX (i.e., selecting single models based on a global mean relationship only) leads to slightly weaker, but still valid, amplification with seven models selected as more realistic (Fig. 3). Using a regional constraint to select the best models at regional scale seems more appropriate, as no model is considered realistic in hot day variability everywhere ( Supplementary Fig. 4).
The constrained TX98p signal (either by the local or global method) suggests that the level of increase in frequency of hot days that was previously estimated to occur by the end of the century could be reached by 2060 instead, i.e., 40 years earlier. All these results are verified to be independent of model selection by performing sensitivity tests where one model is removed randomly from the ensemble (Fig. 3). The regional constraint remains highly skillful in this sensitivity test. The global constraint is still consistent but slightly more sensitive to model  (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) plotted against a variability metric for daily maximum temperature (ΔTX) during the historical period (x axis, in°C). ΔTX measures the difference between daily TX95p and mean TX in a model compared to that observed. The solid black line is the linear regression between ΔTX and TX98p, and dashed black lines show the 95% confidence interval for the slope line. Light gray shading represents the ΔTX difference between ERA5 and CHIRTS (observation uncertainty), and dark gray shading shows uncertainties estimated from the HAPPI ensemble (model internal atmospheric variability) combined. Acronyms refer to AR5 region definitions and numbers refer to models in Supplementary Table 1. selection (due to the small size of n models that fall near the global variability uncertainty range). Similar spatial patterns of increased frequency in extremes to those shown for the regional constraint are also found if applying the global constraint ( Fig. 4), except over desert areas where the signal is not present. This again suggests distinct processes between desert regions and other parts of the globe.
Physical constraints. To further evaluate the EC, we evaluate the physical mechanism linking change in TX98p and land drying. The relationship between the future change in TX98p and latent heat flux is overall negative where the EC is significant and has the strongest influence on results ( Supplementary Fig. 1), indicating larger temperature variability for weaker fluxes (drier soils), supporting our hypothesis. Larger sensible heat flux (heat transfer from land to atmosphere) during hot extremes (Supplementary Fig. 1c) also highlights the importance of land-atmosphere interaction during hot days for our EC. We find a decrease in latent and an increase in sensible heating that correlates with an increase in the frequency of hot extremes. We note that over some regions constrained models do not indicate an increase in TX98p, especially over northern parts of America and generally in mid-latitudes. These correspond to areas with weak correlation between ΔTX and TX98p at present. Other processes may be more dominant in these regions (including change in rainfall and atmospheric dynamics) and drying of soil may be not a factor in high latitudes. In addition, permafrost land-atmosphere exchanges and humidity processes are different there.

Discussion
Our results indicate that climatological bias in the difference between hot and average days in climate models leads to an  Regionally constrained model results are shown in red. A 9-year running mean of TX98p is first applied at each grid point then globally averaged to obtain yearly global mean values. The solid back line shows the mean of a sub-ensemble (seven models) where EC is based on globally averaged ΔTX (instead of applying EC regionally). Gray shading highlights the baseline period used to compute the TX98p threshold and the constraint. Red (and black) dashed lines show sensitivity tests where one model is removed from the ensemble before computation of regional (and global) emergent constraints (test repeated for each model of each ensemble). For each model and each year, TX98p is first normalized by the change in mean Tas then linearly scaled to the multi-model ensemble mean increase in Tas (see "Methods"; note EC results are not sensitive to this, Supplementary Fig. 6). underestimate of the frequency of unusually hot days in the future over many low latitude regions. If the EC is based on the latent heat flux (i.e., the ability of models to reproduce correctly the land-atmosphere humidity exchanges) rather than temperature variability, results support findings over regions where the EC is statistically strongest (tropical regions and South-East Asia). This suggests that changes in hot extremes are related to land surface humidity processes. Results over other regions are less clear, suggesting that other physical processes may dominate the changes in temperatures variability there. When using the regional constraint, the limited number of models selected over some areas may question the robustness of the actual EC signal. However, results using a global constraint are consistent with regional constraints and support the key findings of the latter methodology.
We further note that we focus on daily dry-bulb temperature as metric for hot extremes. In the tropical and subtropical regions, where our constraint is most valid, dangerous heat stress for humans specifically will also be influenced by humidity, which is not considered here. High humidity can lead to deadly conditions when combined with high temperatures and could be a major threat in the future 28,29 . As our EC is physically linked to heat and humidity exchange between land and atmosphere, the wetbulb temperature signal could be investigated with a similar methodology.

Methods
Definition and computation of indices. Our analysis focuses on TX extremes (TX98p). We define TX98p as the number of days above the daily climatological 98th percentile during the warm season (see below). The latter is computed for each location and each calendar day by pooling together all days within a ±15 days window of this calendar day during the 1995-2005 period and selecting the 2% highest values.
We also define a metric, ΔTX, to evaluate the variability of TX during the baseline period over the warm season. It is defined by first calculating the mean and 95th percentile of the temperature distribution for each calendar day at each location (by pooling 15 days around each calendar day as for 98th percentile described above). The distance between the 95th percentile and the mean gives an indication of TX variability for each day and each location. It is computed for each model (ΔTX model ) and for the reference dataset (ΔTX ref ; the ERA5 reanalysis) and the difference between the two defines our metric: ΔTX = ΔTX model − ΔTX ref . Figure 5 illustrates the computation of these indices.
We only focus on the warm season, when hot extremes are likely to happen (June-August for the Northern Hemisphere, December-February for Southern Hemisphere, and all year for the 15°S-15°N tropical areas). Positive values of ΔTX mean a model overestimates the TX variability compared to the reference (i.e., it overestimates high values of TX), negative values indicate an underestimate. For the metric, we choose the 95th percentile to ensure reasonably good sampling of the variability across the base period (as it is used to constrain models), whereas for future changes we focus on the 98th percentile, which correspond to more extreme values. We verified that EC results are not very sensitive to the choice of threshold by doing a sensitivity test using the 95th percentile as threshold instead 98th ( Supplementary Fig. 9).
Each index is computed individually for each model and each ensemble member on their native grid. Results are then interpolated on a common 1°grid before being averaged across all models. As temperature extremes are relatively large scale and grids vary only between 1°and 2.5°latitude/longitude across models, results should not be sensitive to the order of operation.
Datasets CMIP models. An ensemble of 27 individual models from CMIP5 8 and 7 from CMIP6 9 is used. Some only have a single member available while some provide a multi-member ensemble. In the latter case, multi-member results are always computed individually and then averaged to provide one mean result for a single model. We consider a reference period as the historical 1995-2005 decade (being the last decade of CMIP5 historical forcing). Climate projections are investigated using the RCP4.5 30 and SSP245 31 pathways for CMIP5 and CMIP6 models, respectively. Both scenarios are expected to be similar, although each model leads to different mean temperature increases ( Supplementary Fig. 5).
Three climate projection targets are considered as follows: End-of-century, by selecting the 2091-2100 decade for each model. A 1.5°C and +2°C warming above the pre-industrial mean. For these two, we follow a similar approach as in previous study 32 and select for each member of each model the first decade when the average atmospheric surface temperature (Tas) of each year of the decade is above the corresponding threshold ( Supplementary  Fig. 5). As we use 1995-2005 as a baseline, the actual threshold (relative to the baseline) is chosen as +0.7°C and +1.2°C for targets +1.5°C and +2°C above pre-industrial, respectively, as in the HAPPI experiment design 26 . Although the exact definition of these levels can be sensitive to how the baseline period is defined 33 , for this work the main point is that each model or member should reach a similar magnitude of warming. A few members and models do not meet the condition for reaching the +2°C target before the end of the century. For these cases, we select instead the last projection decade 2091-2100. If the mean increase in Tas over this decade is above the threshold (+1.2°C), then we keep the model or member. Otherwise we do not include it in the analysis for this projection target. This leads us to discard four members.
For each climate projection target, results of each member or model are normalized by their respective mean change over the decade (relative to our baseline) in Tas and then averaged to provide ensemble mean results. Thus, for the three specific target projections, all results are shown normalized to an overall warming of +1°C above the baseline.
When showing time series of TX98p (Fig. 3), normalized results for individual models are re-scaled to the multi-model ensemble mean increase in temperature for each year (to include more explicitly the global warming trend). We tested the sensibility of the results by using raw results (without normalization) for each model but both methods lead to very close results in terms of EC amplification ( Supplementary Fig. 6), although raw results have larger uncertainties. Thus, we largely focus on normalized results in the body of the paper.
For most of the models we could get daily TX data for both historical and projection periods while heat flux daily data are more limited. Supplementary Table 1 provides details about outputs used for each variable.
HAPPI ensemble and ΔTX uncertainty. To evaluate the model uncertainties on ΔTX during the baseline period, we use several atmosphere only model simulations driven with observed SST patterns from the historical HAPPI ensemble 26 . Each model provides daily output for the 1995-2005 decade. We select 5 models with a 100 or more members and compute ΔTX for each member (same method as for CMIP models). Then, using the internal variability of each model (multi-member ensemble SD, σ), we estimate ΔTX uncertainties for each location and calendar day ( Supplementary Fig. 7). One model has a mean bias that is much larger than other models (CanAM4) and we thus exclude it (although, as its internal variability is close to other models, we found similar results when including it too). For other models, the ΔTX internal variability is consistent, so we use the mean of four remaining model variabilities (i.e., averaging the four internal standard deviations, STDs) as a measure of ΔTX uncertainties (σ HAPPI ).
The sensitivity of this choice is also tested by using individual model STD instead of ensemble mean (Supplementary Fig. 8). It shows that results stay consistent for each case. We note that the uncertainty so described is that of atmospheric variability only. However, both the HAPPI ensemble and the ERA5 reanalysis are driven by the same SSTs hence this choice is conservative to characterize observational uncertainty.
Internal variability in the climate models used is reduced by ensemble averaging. To take into account the specific number of members for each individual model, the uncertainty between observation, OBS, (σ 2 OBS , see "ERA5" section II The distribution of these temperatures is used to compute the mean (TX mean ) and the 95th (TX 95 ) and 98th (TX 98 ) percentiles. The difference between TX 95 and TX mean is computed (ΔTX data ). III ΔTX and TX98p are then computed from the previous results.
below) and models is expressed as the quadrature sum: (σ 2 OBS + (σ 2 HAPPI /N)) 1/2 with N the number of members of a model. When the absolute value of ΔTX fits within that range then a model (eventually the multi-members ensemble mean) is considered as consistent with OBS.
ERA5. The ERA5 reanalysis 19 is available for the full satellite observation period (1979-present). It provides hourly timescales data at 0.25°resolution on a reduced Gaussian grid, from which we computed daily TX for the 1995-2005 period.
We evaluated the variability of TX in ERA5 against two dense regional observational datasets (Supplementary Fig. 13): a network of 756 homogenized station measurements for China, provided by the Chinese Meteorological Administration 34 , and gridded 0.25°E-OBS v19.0 dataset for Europe 35 . Chinese observations are first gridded on the same regular grid as ERA5 by linear interpolation. Although the TX variability tends to be weaker in ERA5 than in observations, differences are within the range of uncertainties estimated from the HAPPI ensemble variability (Supplementary Fig. 7) for both regions; hence, we consider ERA5 sufficient.
Another observational dataset is used to evaluate ERA5 globally: The Climate Hazards Center InfRared Temperature with Stations CHIRTS 27 ( Supplementary  Fig. 14). Some differences are visible locally but are within the range of HAPPI variability. The difference between ERA5 and CHIRTS is used as an observation error (σ 2 OBS ) and added (in quadrature) to HAPPI variability when applying the EC (see above).
EC method. To decrease model projection uncertainties, model weighting based on individual model performances can be applied, providing an accurate knowledge of the single model skill 36 . This remains a challenge in the presence of regional dynamics 37 . Here we use an EC method based on model selection. This method is illustrated by Fig. 6. ΔTX is our target metric, i.e., we select models that are able to reproduce the width of the TX distribution (the distance between the 95th percentile and the median) within observational and natural variability uncertainty, and select those models for prediction. It has already been demonstrated that TX variability is a useful indicator to weight models and reduce projection uncertainties 38 . To do this here, CMIP models are evaluated against ERA5 during the 1995-2005 period and selected to agree with it within atmospheric internal variability. We use variability from the HAPPI ensemble to characterize variability uncertainty for better sampling, combined with an estimate of observational uncertainty (see above). Models (ensemble mean in case of multi-members model) within the range of error (as described above) are considered reasonably realistic and selected for use in the constrained climate projections. Comparing constrained against unconstrained ensemble projections provides an estimate of the potential current bias in climate forecasts.
Constraints may be applied using global or regional processes 25 . Here we use a regional constraint to take advantage of model spatial information. We first apply a spatial smoothing of 5°on ΔTX over land (to improve sampling and avoid spatial discontinuity; although results are not very sensitive to the scale of the spatial smoothing, see sensitivity test discussed below) then select the models that comply with the constraint within uncertainty at each grid point. Over most of the regions, the number of selected models is between 5 and 10, except in central Africa where it is below 5. This is mainly due to very narrow observational variability over this region ( Supplementary Figs. 4 and 7). Most of the models contribute to the projection over some part of the land. Applying an EC at a global scale instead (Figs. 3 and 4) leads to similar patterns with slightly weaker amplification of future projections.
We also tested the sensitivity of EC results to different choices of uncertainty around the observational distribution and different spatial smoothing ( Supplementary  Fig. 8). Using narrower (wider) range of variability leads to slightly different results with less (more) models selected. This corresponds to a noisier but more intense signal for narrow smoothing (and opposite for larger smoothing). However, global patterns are still consistent with main results. Weaker spatial smoothing (3°) leads to slightly nosier results, while using a larger smoothing area (11°) leads to a large masked area (because we use only land grid points or alternatively to large variation in actual applied smoothing). Thus, 5°smoothing is a good compromise.
Following previous recommendations 25 , we first confirm the strong statistical relationship between ΔTX and TX98p ( Fig. 1 and Supplementary Fig. 3). We then use a resampling method (by removing randomly a model from the ensemble) to test the robustness of the constraint (Fig. 3). Finally, the physical mechanism hypothesis linking land-atmosphere interactions, ΔTX and TX98p is evaluated ( Supplementary Fig. 1) and a perfect model test is conducted (see below).
Validity of the EC method. First, to avoid selection bias 39 and to verify that results are reflective of a physical constraint and roughly independent from the choice of the metric this constraint draws on to determine our EC, we performed a similar analysis with a set of other indicators, all potentially related to Tmax variability, as follows: the interannual variability of Tmax, the diurnal temperature range (DTR) variability, the surface latent heat flux variability, and the surface sensible heat flux variability. Indeed, previous studies have shown the link between heat fluxes variability and temperature variability 40 . All indicators are computed in similar way to ΔTX (except that the fifth percentile is used for latent heat flux as it corresponds to drier conditions). Results of EC on TX98 using these indicators are shown in Supplementary Fig. 12. All indicators lead to very similar results to those using TX variability. It confirms that the EC applies to most of the tropical/ subtropical areas and South and East Asia (i.e., humid regions). They also confirm that mid-and high latitudes do not show similar results (as it is already the case when using ΔTX); thus, different processes are involved in these areas. The main difference occurs over North African and Middle Eastern arid regions when using DTR as a constraint, with a decrease in projection compared to the increase using our EC. This may be due to the fact that DTR variability is also related to Tmin which is expected to have a larger variation over dry regions. Thus, large DTR variability may be an indication of models cooling down too quickly during the night or warming up too quickly during the day, which makes results less reliable here and DTR less suitable for an EC. Results may also be influenced by shortwave radiation and with it clouds (Supplementary Fig. 12).
Second, to verify the validity of our regional constraint, we used a perfect model test as used in previous studies 41 . We select a model as a reference instead of observations, and apply the EC using this new reference, and then compare the projection in TX98p (after excluding the reference model from the ensemble projections) between the constrained models versus all models to the target of prediction (the prediction by the selected model). We first chose the model which showed the largest fraction of grid points consistent with ERA5 (model 34 in Supplementary Table 1, IPSL-CM6A-LR). Results of this test are shown in Supplementary Fig. 15. It clearly indicates that the error in the TX98p projection for this model (difference between the reference model and other models) is reduced for constrained models in the tropical and subtropical areas (i.e., where we also see the strongest signal with EC using ERA5). It is worth noting that even over desert regions, the error is reduced despite the less clear mechanism there. We also tested the same method using a model that compares less favorably to the observed constraint as target (model 1 in Supplementary Table 1, BCC-CSM1) and results are shown in Supplementary Fig. 16. Improvements are also very clear, indicating that "bad models" tend to attract each other and are consistent in their hot extreme projections. Finally, we repeated the process for all individual models and confirmed that errors are always reduced for constrained ensembles versus all models mean ( Supplementary  Fig. 17). This is a powerful confirmation of the validity of our EC 25 .

Data availability
The authors declare that all data that support the findings in the main article are available. All model data are publicly accessible via the Earth System Grid Federation node (https://esgf-node.ipsl.upmc.fr/). ERA5 data can be downloaded from ECMWF website (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5). CHIRTS data can be downloaded from https://www.chc.ucsb.edu/data/chirtsdaily. Fig. 6 Schematic illustration of the emergent constraint (EC) methodology. The projection of a given variable (TX98p in our analysis) is first estimated by an ensemble of individual models (a, red and yellow circles), providing a mean and uncertainty in the projection (right panel, red symbol). A metric (ΔTX in our analysis) is then computed in both models and a reference (ref) during the historical period. It must be verified that this metric is reasonably correlated to the projected variable and that a physical link can be explained. If so, models within the error range (gray shading, corresponding to the observational error and model internal variability) are selected as good models (yellow circles). The projection using only these selected models (b, yellow symbol) is computed and can be compared to the first estimated projection. An efficient EC will usually provide a more accurate projection (reduce uncertainties) and eventually a different estimate of the mean forecast.

Code availability
Scripts used to generate the main results and other data and codes that support the figures in the Supplementary Information are available from the corresponding author on request.
Received: 3 August 2020; Accepted: 9 December 2020; Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.