Atmospheric convection, dynamics and topography shape the scaling pattern of hourly rainfall extremes with temperature globally

Precipitation extremes are expected to intensify under climate change as ground temperature rises with a rate similar to the air’s water holding capacity ~7%/K (Clausius-Clapeyron). Recent studies have been inconclusive on the robustness and global consistency of this behavior. Here, we use hourly weather stations, 40 years of climate reanalysis and two convection permitting models to unravel the global pattern of rainfall extremes scaling with temperature at the hourly scale and identify hotspots of divergence from thermodynamical expectations. We show that in high- and mid-latitudes precipitation extremes closely follow a Clausius-Clapeyron scaling, while divergence occurs over the tropics and subtropics. Local features of atmospheric convection, larger-scale dynamics and orography, affect the dependence of extreme rainfall on surface temperature. In regions with deep convection, persistent large-scale dynamics and complex orography, hourly rainfall extremes diverge from expectations from the atmosphere’s water holding capacity, suggests a global analysis of station data, reanalyses and convection-permitting models.

D ata and climate models show a consistent increasing trend of precipitation extremes (PEx) [1][2][3][4] in a warming climate, accompanied by increases in their frequency 5 . An increase in PEx will impact flooding and ecosystem functioning globally [6][7][8] . Knowledge of precipitation change at hydrologically relevant scales (e.g., hourly and sub-hourly) is thus crucial 7,9 . As atmospheric thermodynamics influence precipitation formation, increases in PEx must be related to the increase in temperature. Since records from weather stations provide data at the ground level, several studies in the last decade have attempted to link surface air temperature (SAT) to PEx [10][11][12][13][14][15][16][17] . The common assumption has a thermodynamic foundation stating that PEx should scale with SAT at a rate similar to air's water holding capacity. According to Clausius-Clapeyron (CC), this rate is 7%/K 10 . Within this context, SAT could be used as a predictor for future changes in PEx. However, whether such dependence is valid globally remains an open question, due to the limited knowledge of the importance of the dynamic component of precipitation formation process 18 , and the limited ground data 2 .
Studies agree toward an emerging pattern of positive~CC scaling rates in the high-and mid-latitudes 13,17 . Diverging patterns have been observed mostly over the tropics and subtropics, including super-adiabatic (super-CC) [10][11][12][13][14][15][16][17]19 and negative rates 13,17,20 , contradicting thermodynamic expectations. Super-CC rates have been attributed to a 'positive feedback' mechanism, where latent heat release during moisture condensation enhances convection and cloud updrafts, intensifying precipitation rates 21 . Super-CC rates have been also interpreted as a statistical artefact due to the mixing of different storm types 11,22 . In some cases, a positive scaling relationship with SAT is exhibited up to a certain temperature threshold (breakpoint), followed by negative scaling rates at higher temperature ranges, thus demonstrating the so-called "peak-like" scaling behavior 23 . This "peak-like" structure and negative rates have been attributed to limited moisture availability at high temperatures diminishing precipitation formation 20,[22][23][24] . Using dewpoint temperature (DPT) as a covariate has been suggested 25 , which eliminates the "peak-like" behavior and strengthens rates globally 13,17 . In some studies, a transition from CC to super-CC rates has also been reported 10,19,21 . Moisture convergence at large scales, variations in temperature lapse rate, changes in vertical velocities, associated with changes in circulations 18,21,26 , and orography effects 11 can decouple ground-level measurements from precipitating systems. Interpretations and comparisons among different studies are not a trivial task, since the methods used for scaling estimation, the spatiotemporal scales examined, and even the definition of extremes are not consistent. Observed divergence from CC and partial disagreement among climate model projections and estimations based on scaling approaches, have raised questions over their applicability in quantitative climate change predictions 15,20,27,28 .
Global scale quantification of PEx changes at fine temporal scales is very challenging. Meteorological data at the hourly or sub-hourly scales with adequate length are scarce and mostly available in the developed world 2 , providing limited information for crucial and vulnerable regions such as the tropics. Atmospheric circulation models have been used to fill this gap. Traditionally their coarse spatial resolution, their inherent parametrization of convection 29 and their limitation in representing cloud microphysics 30 has been known to greatly affect their skill in reproducing PEx 29 . Such models have been used to produce past weather proxies via reanalysis. However, reanalysis products have now increased their spatiotemporal resolution to scales capable of estimating intensification of PEx at hydrologically relevant scales for the first time.
In this study we quantify the scaling of PEx with SAT and DPT globally at the hourly scale for the first time, to our knowledge.
Using a combination of the recent ERA5 global reanalysis product, the convection permitting configuration of the UK Met Office Unified Model (UM), the Weather Research and Forecasting (WRF) convection permitting model (CPM) output over the US and 125,086,356 years of observational data, we validate the capability of atmospheric models to recover the observed dependence of PEx on SAT and DPT, and reveal the global pattern of hourly PEx scaling with SAT and DPT. A consistent CC scaling with SAT is evident thourhout the high-and midlatitudes, while deviations over the tropics and subtropics are exhibited. Divergence from thermodynamical expectations can be linked to the local features of atmospheric convection, larger-scale dynamics experienced in a region and topography.

Results
Scaling rates at stations. A valid global quantification of the scaling rate of PEx with either SAT or DPT using reanalysis requires that (a) a statistically significant scaling exists; (b) it is constant in time; (c) and reanalysis truthfully represents the ground truth (i.e., the scaling estimated from weather stations).
To test whether (a) and (b) hold true we first used station records exclusively ( Supplementary Fig. 1). A robust scaling relation between PEx and both SAT or DPT was found, with uncertainties less than ±2%/K for 75% and 78% of all stations respectively (Fig. 1a). The scaling rates with SAT for two subperiods in the record (before and after 1980) were highly correlated (R 2 = 0.57) despite the increase in average temperature under climate change (Fig. 1b). The bias was not strong (slope = 0.71) and the distributions of scaling rates between the two subperiods were not statistically significantly different (p-value = 0.25 for a two sample Smirnov-Kolmogorov test).
Scaling rates with both SAT and DPT, based on station data and reanalysis for different climates were estimated (Fig. 2a). ERA5 closely captures the climatic signal observed by stations in most climates. Station data in semi-arid, oceanic, humid continental, and humid-subtropical climates exhibit a positivẽ CC scaling with SAT, while Mediterranean climates exhibit lower, yet still positive rates~4%/K. ERA5 exhibits similar results, except in semi-arid and oceanic climates, where the estimates are consistently lower. Differences could be partially attributed to sample size, as reanalysis captures the entire geographic area of each climate, while station data, a small proportion of it.
In all the climates examined except humid-subtropical, scaling with SAT is lower than DPT for station data. In ERA5 rates with SAT are higher than with DPT in Mediterranean, almost equal in oceanic and humid-continental and lower for semi-arid regions. Super-CC occur mostly at semi-arid, oceanic, and humid-subtropical areas primarily for station data, while negative rates rarely exist in using either datasets.
To validate (c) we compared the scaling rates with SAT and DPT using station data and ERA5 for the same period at 1460 locations. Estimates using ERA5 are 1.88 and 3.2%/K lower than stations for SAT and DPT respectively. Despite the bias, the rates between the two datasets were correlated (R 2 = 0.12, p-value < 0.01 for SAT and R 2 = 0.295, p-value < 0.01 for DPT) (Fig. 2b, c) and ERA5 captures the general spatial pattern over the US ( Supplementary Fig. 2), the only area in our study where station density was adequate for the quantification of large-scale patterns. Large scatter between estimates can be attributed partially to the uncertainty of the estimate itself, and the mismatch between the spatial scale of the data sets.
Global pattern. Using ERA5, distinct spatial features of the scaling rate between hourly PEx and both SAT and DPT were revealed (Fig. 3). A~CC rate with SAT is evident for most parts of the world 13,17 , except for the tropics and subtropics exhibiting rates spanning from negative up to >2CC (Fig. 3a). Negative rates were eliminated using DPT for most parts of the world and stronger super-CC rates appeared in the tropics (Fig. 3d) 13,17 . Scaling relationship with SAT and DPT was robust globally (Fig. 3e, f), with uncertainty less than ±2%/K over 48% and 55% of gridpoints respectively. Scaling with SAT exhibited larger uncertainties than DPT in the Sahel, north Australia and over the west Indian peninsula (Fig. 3e, f). A seasonal investigation provided similar, yet strongly more uncertain results .
Scaling with SAT exhibited a~CC rate over most of the highand mid-latitudes (Fig. 3a). In the tropics, super-CC rates appeared over the Andes, northwestern South America, parts of the Sahel, the Congo river basin and the northern part of the Indian and Indochinese peninsulas. Negative rates were pronounced over the drier parts of eastern Brazil, coastal northwest Africa, coastal central east Africa, parts of the eastern Indian and Indochinese peninsulas, western central America, and the Coastal mountains in Canada. "Peak-like" scaling was exhibited in the Sahel, tropical Australia, tropical and subtropical South Africa, the Mediterranean coast of Africa, the Middle East, central Asia, subtropical and west India, the south and southeastern US and over parts of the North American Cordillera (Fig. 3b, c).
A similar~CC pattern over the high-and mid-latitudes was revealed with DPT, with 52% of gridpoints exhibiting greater rates than SAT over the temperate zone (Fig. 3d). Over the tropics a consistent super-CC rate was apparent. Negative rates appeared over the southeastern US, Bangladesh, central South America, the northern coasts of Gulf of Guinea, the eastern coast of the Red Sea, the eastern coast of the Andes and the Coast mountains in Canada. In the subtropics and tropics, 34% and 70% of gridpoints respectively exhibited greater rates than with SAT.
The analysis has also been performed for the 0.95 th quantile and at the daily scale as well . Using the ERA5 dataset, the general global scaling pattern remains roughly similar throughout both quantiles and temporal scales. Over the high-and mid-latitudes, aggregating at the daily scale yields lower scaling rates for both quantiles (Supplementary Figs. 11 and 12), while changes due to decreasing quantile levels alone are less prominent (Supplementary Figs. 9 and 10). However, this behavior does not hold over the tropics and subtropics, where at the daily scale greater super-CC rates are evident with DPT, while a strongly negative band with SAT is exhibited over tropical Africa.
Convection-permitting models. Most of the regions with rates strongly diverging from thermodynamical expectations (e.g., ≫ CC or <0) experience intense convection, orographic enhancement of precipitation or persistent larger-scale weather patterns. Given the well-documented biases of coarse resolution models in simulating deep convection and orographic effects 29 , we examined whether this divergence could be attributed to ERA5's~31 km resolution. Comparing ERA5 and station data, events of strong convective nature are more strongly underestimated than non-convective ones, both in terms of total rainfall amounts ( Fig. 4a) and in terms of extremes, quantified here as the heaviness of the tail of their distribution ( Fig. 4b) (see "Methods"). In order to investigate the impact of this underestimation in representing the scaling behavior, we compared ERA5 estimates and station data with outputs from the WRF (~4 km) CPM and the CPM employed in the UKCP18 project (~2.2 km) (see "Methods"). CPMs outperform ERA5 in simulating rainfall (Fig. 4c, d), and perform better than ERA5 in estimating scaling rates, even though a positive bias exists (Fig. 4e, f).
Over the US, both ERA5 and WRF can adequately capture the main features of the observed scaling pattern ( Fig. 5 and Supplementary Fig. 13), despite existing model biases in representing rainfall extremes. As expected, WRF is more reilable (Fig. 5a-f and Supplementary Fig. 14), providing stronger scaling estimates than ERA5 over most regions, except for the North American Cordillera, the Pacific Coast and parts of the south ( Supplementary Fig. 14). WRF estimates are consistently larger than ERA5 when DPT is used over the Great Plains, Midwestern States and the Southeast ( Supplementary Fig. 13). In terms of model performance, WRF outperforms ERA5 over the Southeast, the Midwestern States and the Northeast with SAT (Fig. 5g). Performance over those regions improves when DPT is used, while WRF performs better over the Great Plains as well (Fig. 5h). WRF's deviation from reality in terms of PEx scaling with SAT is more evident over the North American Cordillera, its' Western Foothills, the Pacific Northwest and Southwest, while with DPT is biased over the Cordillera (Supplementary Fig. 15). ERA5's b Density plot of the observational versus reanalysis scaling rate with SAT for each station (R 2 = 0.12, p-value < 0.01). c Same as b, but with DPT (R 2 = 0.295, p-value < 0.01). The black line indicates the 1-1 line. In cases when a peak-like behavior is exhibited (see methods), only the ascending branch of the bilinear model is presented.
deviation is more widely spread over the domain, especially with DPT ( Supplementary Fig. 15). Over the UK a west/north-east/ south dipole is exhibited, whereby~CC rates over the first and super-CC rates over the latter are evident in both models in agreement with observations ( Supplementary Fig. 16). The CPM performs better than ERA5 when SAT is used, however both models exhibit larger errors with DPT ( Supplementary Fig. 16).

Discussion
Although climate change impacts on rainfall extremes have already been documented 2-5 , the scaling relationship remains roughly constant over the historical record, indicating that it could remain roughly unchanged in the future climate as well. This validates the findings by Prein et al. 24 , who found that scaling rates over future projections with WRF remain roughly similar. Model projections confirm the relevance of thermodynamics in constraining future PEx, even though confidence on straightforward extrapolations to the future using simplified scaling approaches remains low 23,24,31,32 . However, despite the potential of using a thermodynamically consistent~CC scaling as a future PEx intensification predictor, its behavior is highly divergent in many parts of the world.
A "peak-like" behavior, which has been attributed to moisture limitations at the higher temperature ranges 20,22-24 has been found in many regions ( Fig. 3c and Supplementary Fig. 17). According to a recent hypothesis, an additional transition regime from CC to super-CC rates exists 10,19,21,33 , however such a behavior was not prominent in both the station and model output data that we analyzed. Studies have been inconclusive regarding a full physical explanation of such a transition 11,21,34 and, hence a more complex model that would also account for such a behavior has not been examined, in order to reduce uncertainties in our analysis.
Our results demonstrated a consistent~CC scaling in most parts of the high-and mid-latitudes for both SAT and DPT. In the tropics and subtropics, scaling with SAT spans from negative to >2CC rates. However, the scaling pattern with DPT over those regions is vastly different, since mostly super-CC rates are exhibited, in agreement with previous findings 12,13,15,17 . This global pattern is roughly consistent over all quantile levels and temporal scales examined. In the observational dataset, results comply with previous findings 10, 11,15,16,25,35 , exhibiting lower rates as the temporal scale becomes coarser and the quantile level decreases (Supplementary Fig. 8). However, it should be noted that when aggregating to the daily scale strong short-duration convective events can be masked, larger-scale events can become dominant 34,35 , while the feedback of rainfall itself on temperature is inevitably lumped in the analysis as well. A similar pattern has been obtained when daytime and nighttime events are separately examined ( Supplementary Fig. 18), however deviation between the two can be strong in regions with a pronounced diurnal signal in the occurrence of hourly rainfall extremes, such as the tropics.
Super-CC rates have also been exhibited and could be attributed to the contribution of the dynamic component of precipitating systems strengthening PEx 18,26 . Negative rates with SAT imply that the most extreme events do not occur when SAT is highest. In the tropics, given the low variation in solar insolation throughout the year, ground air temperature is linked with land surface feedbacks. High SAT occurs during the dry season when cloudiness is lower and thus net radiation on the ground is higher, which can increase substantially SAT during the day, especially if soil water content is low 20 . This is a clear indication of the decoupling of integrated air temperature that impacts the air's water holding capacity and the temperature at the surface. This can be further supported by the fact that negative rates are fewer when DPT is used. In the dry season when SAT is high, due to increased sensible heat at the ground level, relative humidity is typically low, and thus the corresponding DPTs are low.
However, there are a few regions exhibiting negative rates even with DPT. Such a negative patch is exhibited over the coast of the Guinea Highlands, where the hotspot of rainfall activity during the West African Monsoon (WAM) is located 36 , due to the orographic enhancement of precipitating systems reaching the Highlands. However, the complexity of WAM's atmospheric dynamics and the topography allow only for limited confidence over ERA5's estimations, since even Regional Climate Models with finer resolutions have been found to fail in simulating extreme rainfall over the region 36 . Over Bangladesh, a patch of strongly negative rates with DPT and a mixture of strongly negative and super-CC rates with SAT emerges ( Supplementary  Figs. 19 and 20). Most of the extreme rainfall during the monsoon period is associated with the break of the Indian monsoon's active phase and the shifts of the monsoon trough 37,38 . When the trough reaches the foot of the Himalayas, increased rainfall occurs. During this break period, westerly/southwesterly winds are dominant, transporting moisture from the Ganga Basin 37,38 . Upon reaching the Shillong Plateau in northeastern India and the Chittagong Hill Tacts in west Myanmar, the incoming moisture gets stagnated, large-scale moisture convergence occurs and orographically enhanced convective systems are formed 37 , even though the local ground temperatures can be low.
Southeastern US, along the Gulf of Mexico also exhibits negative scaling with DPT ( Fig. 3d and Supplementary Fig. 20), while station data contradict this finding (Fig. 5 and Supplementary Fig. 2). Northward atmospheric rivers 39 are often detected and intense PEx associated with tropical cyclones 40 commonly occur in late Summer and Fall, when SAT and DPT are usually lower. Deep convection associated with Mesoscale Convective Systems (MCS) is also experienced throughout the year 41 .
Over parts of the extratropical east coast of the Atlantic and the Pacific Ocean, a weakened or negative scaling pattern with both SAT and DPT is evident ( Supplementary Figs. 19 and 20). Ground data confirm this behavior for the west coast of the US (Supplementary Fig. 2; i.e., the only region we could validate it robustly). In those regions, atmospheric rivers transporting moisture from the tropics eastwards, generate heavy rainfall, mostly during the winter 39 and, hence, lower temperature ranges are occupied with more extreme events.
In South America, along the central and northern parts of the La Plata river basin and northern Bolivia, a well-defined band of negative and weakened scaling rates is evident ( Supplementary  Figs. 19 and 20). Its trajectory mostly coincides with that of the South American Low Level Jet (SALLJ) that gets diverted southeastward by the Andes 39,42,43 . Severe MCS form over those regions during the warm season, which can be associated with deep convective activity in the South American Convergence Zone (SACZ) over parts of southeastern Brazil 43 . MCS are fed with moisture from the SALLJ, carrying warm and moist air masses from the tropics 42,43 .
The true nature and magnitude of larger-scale weather patterns with their associated mechanisms, like atmospheric rivers (occurring higher in the troposphere 39 ) and LLJs (occurring lower in the troposphere, roughly at the 850 hPa level 42 ) cannot be fully captured by the variability of ground measurements. Hence, these moisture flows feeding clouds decouple ground  39,44 . The effects of topographic heterogeneity, altitude and the presence of orographic barriers, which are known to affect precipitating systems locally 44 also increase the uncertainty of the scaling relationship (e.g., uncertainty is larger over large mountain ranges such as the Andes, the Rocky Mountains, and the Himalayas than their surrounding areas (Fig. 3e, f and Supplementary Fig. 7)). The predictive skill of any scaling relationship over such regions for future rainfall intensification is therefore limited.
Overall, the regions diverging from CC are either tropical/ subtropical experiencing deep convection, or with pronounced influence of topography, or affected persistently by larger-scale weather patterns and their associated mechanisms. Here we demonstrated that in ERA5, events of strong convective nature are underestimated more strongly than non-convective. The differential underestimation between storm types occupying different temperature ranges can lead to underestimating the scaling rates, thus limiting the applicability of coarse reanalyses as good data proxies for this task. It should also be noted that generally model estimates, especially non convection permitting ones, are less reliable over the tropics 18,45 .
Over the UK, during the warm months convection occurs throughout the island, but mostly over the southeastern parts 46 . During the winter, the westerly low pressure systems and atmospheric rivers landfall upon reaching the Isles and get orographically enhanced, producing more rainfall over the western/ northern parts of the country 46 . This explains the resulting dipole showing weaker rates in the western/northern regions where enhanced winter PEx occur.
Over the US CPMs perform better than ERA5 and can represent the observed scaling pattern more faithfully, even though ERA5 can still exhibit a roughly similar, yet negatively biased pattern (Fig. 5 and Supplementary Figs. 13-15). Performance is better especially over areas experiencing deep convection, such as the Southeast and the Midwestern Sates 41 (Fig. 5g, h and Supplementary Fig. 15), where the coarse spatial scales of current reanalysis products can be very limiting. Performance is improved over the Northeast as well, which experiences enhanced convective activity throughout the summer and a few MCS propagating eastwards as well 41 . However, WRF does not outperform ERA5 over the Central and Northern Great Plains markedly, even though MCS are also known to occur there 41 . This could be attributed to the uncertainty of the estimation itself or the limiting length of the WRF simulation. The documented WRF strong underestimation of strong MCS over the region 47 , could also be biasing the results, even though many stations in the area are located towards the eastern parts of the Plains that do not experience such frequent MCS occurrence 41 . What is most striking is WRF's inability to outperform ERA5 over the Cordillera ( Fig. 5 and Supplementary Fig. 15) where the topography is complex and one would reasonably expect otherwise. Simulation length and the uncertainty of the estimation alone probably cannot explain this behavior and there are possibly other underlying reasons. Further research is therefore needed to reveal WRF's inconsistencies over the region. However, it should be noted that the station network over the Cordillera is relatively sparse and located predominantly in valleys and hence model validation is not a trivial task.
Even though WRF's deviation from reality scatters roughly around 0 with a small positive bias for both SAT and DPT (Fig. 4e, f), it is more distinctly organized spatially when DPT is used (Supplementary Fig. 15). In ERA5, deviation from the ground truth is notably higher with DPT (Fig. 4e, f). This can be related to the peak-like behavior of PEx scaling with SAT. If a model tends to underestimate convective events, which naturally occur when temperatures are higher, a peak-like behavior that can be a numerical artifact arises. In these cases, the ascending branch of the scaling mostly reflects rainfall events of nonconvective nature. Contrarily, when DPT is used, model discrepancies are more likely to be revealed, especially when the mismatch in the extremes is large. This can be supported by the fact that the areas where a peak-like behavior was found, mostly coincides with hotspots of model disagreement with the ground truth ( Supplementary Figs. 15 and 17). As a result, even though the proposed bilinear model is theoretically more consistent, it is not fully informative when direct comparisons between model and reality are to be made and scaling with DPT is considered to be more reliable in this context. However, accounting for moisture limitations is important when studying observational datasets, in order to reveal the true nature of PEx sensitivity to SAT. Another possible explanation, especially over regions like the UK, where a bilinear behavior is not exhibited, is that such a deviation could be associated with model biases regarding humidity estimations.
The overall better performance of CPMs for both the US and the UK, highlights the need for global scale application of such models focusing on PEx 24,32 , particularly over the tropics and orographically complex regions. Such global scale use of CPMs has been prohibitive so far, primarily because of computational constraints, however now for the first time the resources to achieve that task are becoming available. Since data availability is constrained mostly within the boundaries of the developed world, model validation and improvement could only be achieved with a simultaneous large-scale effort to establish high-resolution measuring networks globally, especially in the tropics and sub-tropics. Recent development of reliable and affordable sensors and citizen science projects open new possibilities towards a substantial increase of data volume 48 .
However, even if all model biases get remedied and scaling rates can be confidently estimated globally, future predictions based on thermodynamic scaling alone may be problematic over the tropics and subtropics. Future PEx changes over these regions will be mostly determined by changes in the InterTropical Convergence Zone (ITCZ), which is expected to get narrower and exhibit a weaker mean ascent 49 , affecting the formation of MCS. Weakened vertical velocities might lead to negative trends in PEx 26 , while changes in PEx associated with shifts in monsoonal activity will most likely be jointly determined by aerosol forcing, land-use changes and shifts in circulation anomalies 50,51 .
In conclusion, our findings support the relevance of thermodynamics in investigating future PEx and the applicability of scaling approaches to quantitative climate change predictions for most regions in the higher and mid-latitudes. However, the results are much more complex when strong convection is present, local topography influences precipitation formation or larger-scale weather systems consistently affect an area. Any climate change assessment of PEx and their influence on ecosystem functioning 7,52 and natural hazard management must take into account the detailed atmospheric dynamics, at the proper spatial and temporal scales. To achieve that for hydrological applications, the use of CPMs should be expanded, and given the inherent uncertainties in PEx 53 , ensemble simulations must be considered. This is a computationally challenging task, which in our opinion is necessary for proper climate change impact studies.

Methods
Data. We applied the following quality control rules to derive our station records. Only stations with at least 30 years of data and at least 25 years of non-missing values were considered. All data flagged as spurious were treated as missing, while a visual inspection for possible outliers was also performed. The final dataset consists of 1460 stations, located in the US, Netherlands, UK, Switzerland, Japan, and New Zealand (Supplementary Fig. 1). For every station hourly precipitation and daily average SAT data are extracted, while daily mean DPT measurements were available for the 619 stations. The ECMWF ERA5 reanalysis dataset 54 , has a spatial resolution of 0.25°× 0.25°(~31 km) spanning from 1979 to 2019 and operates at the hourly temporal scale. The UK Climate Projections UKCP18 dataset 55 has a spatial resolution of 2.2 km × 2.2 km spanning from 1980 to 2000. UKCP18 is a historical run of the convection permitting UM. Eleven ensemble members have been used. The high-resolution WRF simulation of the current climate of North America 56 , has a spatial resolution of 4 km × 4 km, operates at hourly temporal scale and consists of a single retrospective simulation spanning from October 2000 to September 2013, driven by ERA-Interim reanalysis.
Model. Scaling estimation is performed using quantile regression which has been proven to be more robust than the widely used binning technique 57 (Supplementary Fig. 21). For each statistically independent precipitation event the peak intensity was extracted. SAT and DPT are aggregated at the daily scale since we consider it more representative of the conditions leading to an event. In order to avoid possible feedbacks on temperature from precipitation itself (e.g., evaporative cooling) daily aggregation is not made based on calendar days, but using a 24-h time window prior to the corresponding PEx. The scaling rate corresponds to the 0.99 th quantile regression between peak intensities and the corresponding SAT or DPT. Statistically independent events are separated by a 2-h dry interval. In order to account for possible moisture limitations in the upper range of SAT, the scaling rate is estimated by using a single linear model or a combination of two linear models (i.e., a change in the scaling regime). Each linear model is estimated with quantile linear regression. The breakpoint where the scaling changes in the bilinear model is estimated as the inflection point of a quadratic equation fitted on the 0.99 th quantiles of overlapping 1°C. Choice between a single linear or a bilinear model was done based on the Akaike Information Criterion. Contrary to other similar approaches 33 , continuity is not enforced at the breakpoint and, since the bilinear model is meant to account for moisture limitations at high temperatures, it can only have an ascending branch at the lower SAT range (prior to the breakpoint) and a descending branch at the higher SAT range. Therefore, a bilinear model cannot result to two positive or two negative scaling rates and cannot have a descending branch over the lower SAT range. In these cases, the linear model is chosen instead. When DPT is used, only the single linear model was chosen. The model's performance in identifying the true underlying behavior of a sample has been evaluated with a Monte Carlo simulation scheme as in Van de Vyver et al. 33 ( Supplementary Fig. 22).
Analysis. The uncertainty range of the scaling estimate was defined by bootstrapping the PEx and their corresponding daily SAT and DPT. We defined as "uncertainty" the 0.025-0.975 interquartile range of all estimates. Bootstrapping is performed 200 times, yielding errors~±0.5(%/K) in the uncertainty estimations at the 0.99 th quantile level. Since usually a small number of values lies within the higher temperature range, thus increasing uncertainty, the bilinear model is considered only when at least 5% of the events lie within that range. However, uncertainty might still remain high and, for that reason, in cases where moisture limitations are identified, we refrain from making interpretations based on the magnitude of negative scaling rates. In direct comparisons between observations and models, we only use the stations that have enough values during the model reference period and make estimations based on this part of the full record. Analysis has also been conducted by aggregating rainfall at the daily scale and for the 0.95 th quantile as well. The observed scaling pattern (Fig. 5) is obtained by using an inverse distance weight function to interpolate values.
We also estimated the scaling behavior of observational records for the two subperiods up to and from 1980. For this purpose we only used stations with at least 25 years of non-missing values for each subperiod and also excluded stations where inhomogeneities and changes in the record's depth resolution between the two subperiods were found.
Events possibly having a strong convective nature were defined based the Convective Available Potential Energy data (CAPE ≥ 500 J/kg) obtained from ERA5. ERA5 is a notable improvement of its predecessor ERA-Interim 54 , whose estimations of convective parameters have been evaluated and extensively used throughout the literature, in order to infer the real-world state of convection and identify environments favoring thunderstorm generation 58 . The Weibull distribution was fitted on the tails (upper 5% of non-zero values) of rainfall distributions by minimizing the probability ratio mean square error as in Papalexiou et al. 59 and the Γ distribution using maximum likelihood on all nonzero values.
All model output data are freely available. ECMWF ERA5 reanalysis dataset (https:// www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5) is obtained by the Copernicus Climate Change Service (C3S) Climate Data Store (https://cds.climate. copernicus.eu). UK Climate Projections UKCP18 is obtained by the Centre for Environmental Data Analysis (CEDA) Archives (https://catalogue.ceda.ac.uk/uuid/ ad2ac0ddd3f34210b0d6e19bfc335539). The High Resolution Weather Research and Forecasting (WRF) Simulations of the Current and Future Climate of North America is obtained by the Research Data Archive (RDA) of the National Center for Atmospheric Research (NCAR) (https://rda.ucar.edu/datasets/ds612.0/).