Extreme rainfall and temperature sampling

We analyzed summer rainfall (that is, rainfall recorded between June and September) for 299 climate stations located along the Alps (Table S3). To examine the extreme rainfall, we applied the unified framework proposed by ref. 40 to identify and sample independent rainfall peaks from the data, following a two-step process. First, independent “rainfall events” are defined as rainfall periods separated by at least 24-h dry intervals, ensuring that the events are temporally isolated and statistically independent. Second, rainfall peaks of duration d are defined as the duration maxima of each rainfall event using a running window approach. The window size d corresponds to the duration of interest, and the time step is set to match the temporal resolution of the rainfall data. This method guarantees that the selected rainfall peaks share the statistical properties of the annual maxima for the same duration, as demonstrated in previous studies (e.g.,17,40,41). This approach provides a robust framework for consistent event identification across varying durations, making it suitable for the analysis of extreme rainfall events across durations—in our case, d was set to hourly and 10-min scales.

Rainfall records have time intervals of 5, 6, and 10 min (Table S3) and were homogenized to 10- or 12-min (for the 5- and 6-min stations, respectively) through aggregation. For each station, years with more than 10% missing records were excluded from the analysis. Additionally, as described by ref. 17, we sampled and averaged the temperature for the 24 h preceding each peak to identify the associated proxy temperature.

Computing empirical scaling rates

The scaling rates between extreme rainfall and temperature are computed using a quantile regression model42 that is fitted to the logarithm of the q percentile of precipitation Pq as a function of near-surface air temperature T:

$$\log ({P}_{q})=\alpha +\beta T,$$

(1)

from which the scaling rate is derived:

$$\frac{\partial {P}_{q}}{\partial T}=({e}^{\beta }-1)\cdot 100.$$

(2)

Estimating rainfall return levels

We estimate rainfall return levels using the TEmperature-dependent Non-Asymptotic statistical model for eXtreme return levels (TENAX)17, which enables assessing changes in short-duration precipitation extremes as a function of shifts in near-surface air temperature during wet days. TENAX separates the dependence of extreme precipitation on temperature from the occurrence of precipitation events and combines this information to estimate precipitation return levels. TENAX is composed of three components: (i) a magnitude component that models the exceedance probability of rainfall peaks as a function of temperature; (ii) a temperature component that models the distribution of temperature associated with the rainfall peaks; (iii) a return level estimation component that estimates return levels based on the first two components. The magnitude x of extreme rainfall peaks for a given duration is modeled using a Weibull distribution, with parameters dependent on near-surface air temperature T:

$$W(x;T)=1-\exp \left[-{\left(\frac{x}{{\lambda }_{0}\cdot {e}^{aT}}\right)}^{{\kappa }_{0}+bT}\right],$$

(3)

where λ0 and a describe the exponential dependence of the scale parameter on temperature (thus conceptualizing based on scaling rates), and κ0 and b describe the temperature-dependent shape parameter. A likelihood ratio test run over all the climate stations showed that in most cases b is not significantly different from zero.

The distribution of temperatures during precipitation events is here modeled using a normal distribution with parameters μ and σ, which describe the mean and standard deviation of near-surface air temperature of the 24 h preceding the peak of the events17. Once the magnitude W(x; T) and the temperature g(T) components are established, a stochastic approach is used to generate a large collection of temperatures Ti with i = 1, …, N sampled from g(T), to obtain a Monte Carlo approximation of the marginal cumulative distribution function \(F(x)={\int}_{{\mathcal{T}}}W(x;T)g(T)\,dT\), where \({\mathcal{T}}\) is the domain of g(T). An estimate of the distribution of annual maxima G(x) can then be obtained as30:

$${G}_{{\rm{TENAX}}}(x)\simeq {\left(\frac{1}{N}\mathop{\sum }\limits_{i = 1}^{N}W(x;{T}_{i})\right)}^{n},$$

(4)

where N is the number of stochastically generated events, and n is the average number of independent events per year. Precipitation return levels are then computed by inverting the above equation17 and41 provide additional information on the TENAX model, its parameterization, calibration, validation procedure, uncertainties, and examples of its application.

Climate projections

The main assumption behind the estimation of future rainfall return levels using TENAX is that the magnitude component described above, which represents the statistics of convective rainfall in the location of interest at a given temperature, is invariant. This means that we assume the influence of other covariates (e.g., aerosol concentration) on precipitation intensity at a given temperature to be negligible. We tested this assumption in hindcast by splitting the record of each climate station with records longer than 30 years (74 stations) into two periods of equal length and estimated W(x; T) separately for the two periods, as presented in Fig. 3a for the Samedan station and in Fig. S1 for nine other stations along the Alps. We then tested whether the null hypothesis of having a unique model estimated for the entire period, as opposed to the alternative hypothesis of having two models, could be rejected using a likelihood ratio test at a 5% significance level. We found that the null hypothesis could be rejected in only 6 out of 299 stations.

To parameterize TENAX for future rainfall return periods, it is necessary to collect information on the changes in the mean Δμ and standard deviation δσ of temperature during wet days and the change in the total number of annual rainfall events δn17. This information was obtained from the “historical” and “RCP8.5” simulations of 17 regional climate models (summarized in Table S4) covering the Alps. While the use of CMIP5-EURO-CORDEX (GCM, downscaled by a RCM) to assess changes in temperature and precipitation has been well established43, the models exhibit known uncertainties stemming from their physical formulations, the omission of certain atmospheric processes, and the inherited uncertainties from global circulation models43,44. Notably, these models tend to simulate too cold and too wet climates due to factors such as excessive snow accumulation and albedo effects, overestimated cloud cover, underestimated evapotranspiration, and deficiencies in convective parameterization45. For example, summer temperatures can be underestimated by ~0.5 °C due to misrepresented aerosol effects46. In the Alps, the models simulate colder summer temperatures (by 0.8 °C) and a higher frequency of rainfall events (by 14.8%) compared to observations47. We account for these climate model uncertainties in our predictions of changes in extreme short-duration rainfall, as discussed in “Data and model uncertainty”.

We depart from the conventional use of emission scenarios and assess how a fixed increase in near-surface air temperature over the EURO-CORDEX domain affects Δμ, δσ, and δn. This approach isolates temperature change as the primary driver, independent of specific emission pathways (i.e., solely as a function of general warming levels, see ref. 48 for further discussion). Our first step was to calculate μ, σ, and n for each climate model for the historical period 1976–2005 for all grid cells intersecting with one of the climate stations in the Alpine domain. In addition, we computed the historical mean daily near-surface air temperature Ta for each grid. Then, using a moving average window for the projected period 2006–2100, we calculated decadal values of future μF, σF, nF, and \({T}_{a}^{F}\). A power equation was then used to link the change in temperature during wet days Δμ and the shift in near-surface air temperature ΔTa:

$$\Delta \mu ={c}_{1}\cdot \Delta {T}_{a}^{{c}_{2}}+{c}_{3},$$

(5)

where Δμ = μF − μ, \(\Delta {T}_{a}={T}_{a}^{F}-{T}_{a}\), and c1…3 are the regression coefficients. This enabled finding the required value of Δμ for the specific cases of increase in 1 °C, 2 °C, and 3 °C of near-surface air temperature.

The regressions of δσ and δn with ΔTa were also fitted using power relations:

$$\delta \sigma ={c}_{1}\cdot \Delta {T}_{a}^{{c}_{2}}+1,$$

(6)

and

$$\delta n={c}_{1}\cdot \Delta {T}_{a}^{{c}_{2}}+1,$$

(7)

where \(\delta \sigma =\frac{{\sigma }^{F}}{\sigma }\) and \(\delta n=\frac{{n}^{F}}{n}\), and c1 and c2 being the coefficients. We computed the R2 goodness-of-fit for each grid cell for these three regressions and preserved only those fits exceeding 0.75 (52% of all fits). The median changes in Δμ, Δσ, and δn for the 1–3 °C increase in regional temperature are presented in Fig. S2–S4.

Data and model uncertainty

The climate station data we used varies in both temporal scales and recording periods, contributing to uneven uncertainty levels in rainfall return level estimates across stations. While some stations provide extensive records spanning up to 42 years (from Switzerland; over the central Alps), others cover only 15 years (Austria); this is still within a reasonable length of data to apply the TENAX model to estimate return levels, but wider uncertainties are to be expected in some regions17,30. Moreover, the stations recorded data over different time intervals and we treat the longest period covering 1981–2023 as the “present” period. While the TENAX model does not need stationarity assumptions for its training, we assume the climate to be stationary over this period in computing the changes in wet-day air temperature and rainfall occurrence derived from climate models to generate the future rainfall return periods.

The application of the TENAX model for estimating future rainfall return periods involves several layers of uncertainty. First, uncertainties may accumulate throughout the stepwise process of translating climate model outputs into temperature and rainfall frequency change factors, considering both the transition from daily to sub-daily scales (see ref. 17 for a detailed discussion on the topic) and the relatively coarse resolution of the climate models (11 km), which may not adequately represent local climate conditions (particularly in complex terrains). In addition, the model’s reliance on time-invariant temperature-precipitation scaling, while widely supported10 and tested at the station level in hindcast, may not hold exactly in the future considering potential changes in climate dynamics49 (though this seems unlikely26) and that future wet-day temperatures may reach values unexplored in the observations.

Given the potential large uncertainties in climate projections, we have quantified them alongside presenting the change in the multi-model median. First, we calculated the change in rainfall return levels for the individual climate projections, as illustrated for the 17 climate models and three warming scenarios at the Samedan station (Fig. S5a). Next, we calculated the 10–90th percentile range for each return period, defining it as the plausible uncertainty range derived from the ensemble of climate trajectories. We applied this procedure to all stations (see selected 12 in Fig. 5) and found that uncertainties increase with both warming and return period (see Results for details). Further, we summarized the 10–90th percentiles of the ratio between the standard deviation and mean signal change from the ensemble of climate projection per station to illustrate the uncertainties arising from climate model uncertainty in future rainfall return levels as a function of warming scenarios, return periods, and elevation (Fig. S5b).