We use vector generalized additive/linear models to super-resolve hourly precipitation GEV parameters from 13.2, 26.4, and 52.8km to 2.2km. We define a climate change robustness gap to explain how our generative superresolution model generalizes via its splines. We propose an upper bound for super-resolution factors using spatial auto-and cross-correlation functions of precipitation and elevation. Climate adaptation requires projecting high-impact weather events at local scales, notably extreme precipitation due to its impacts on ecosystems and infrastructure (Fowler et al., 2007;Gimeno et al., 2022). General circulation models (GCMs) are too coarse (with horizontal grid spacing around 100 km) to explicitly simulate such extremes (Maraun, 2016;Benestad, 2004). This limitation warrants the use of expensive regional climate models that can only be run selectively, especially at convection-permitting scales (Schar et al., 2020), potentially under-sampling projection uncertainty (Hawkins & Sutton, 2011). Empirical downscaling, including statistical (Maraun & Widmann, 2018) and machine learning (ML; Rampal, Hobeichi, et al. (2024)) methods, offers a promising complement or alternative by directly predicting relevant variables at local scales when suitable training data are available. Recent advances in ML-based super-resolution (Wang et al., 2022) and generative modeling (Yang et al., 2023) have further fueled the rapid development of empirical downscaling for precipitation (Rampal et al., 2025;Srivastava et al., 2024;Vandal et al., 2019). However, their application to climate change remains limited by challenges in understanding their robustness-i.e., how well they extrapolate to warmer climates (Hernanz et al., 2022). This is especially true for extremes, where stationarity is difficult to assess from historical performance alone (Dixon et al., 2016). This motivates pseudo-reality experiments (also called "model as truth" or "perfect model"), in which outputs from a dynamical regional climate model are treated as pseudo-observations for empirical downscaling, enabling direct benchmarking of generalization capabilities (Maraun et al., 2015). While such experiments have helped identify best practices (e.g., optimal predictor sets) for improving the robustness of simple statistical downscaling algorithms (Charles et al., 1999;Vrac et al., 2007;Dayon et al., 2015), they remain under-used for more sophisticated methods. Notable exceptions include Legasa et al. (2023), who showed that a posteriori random forests targeting parameters of the precipitation gamma distribution (Legasa et al., 2022) generalize better to a warmer climate than generalized linear models and convolutional neural networks; Bano-Medina et al. (2024), who found that deep learning emulators generalize more reliably when GCM predictors are bias-adjusted to the upscaled regional model; and Rampal, Gibson, et al. (2024), who showed that generative adversarial networks outperform deterministic convolutional neural networks in projecting warming-driven precipitation extremes. An emerging challenge is that these results remain scattered-typically evaluated per case or quantile-and we still lack a simple framework to explain what drives robustness across varying degrees of extremeness. Here, we aim to address this gap for the task of super-resolving precipitation extremes, i.e., improving their spatial resolution. We focus on this task because it is tractable: (1) although fine-resolution and coarse-resolution extremes do not occur simultaneously, they typically fall within the same temporal block, making them weakly paired; and (2) super-resolution is generally not treated on a per-GCM basis. This contrasts with the full empirical downscaling pipeline required for local climate projections, where global and sometimes even regional climate model outputs must be bias-corrected (Francois et al., 2020;Cannon, 2018;Vrac & Friederichs, 2015) before being brought to the desired spatial scale. We emphasize that super-resolution is not trivial, especially for precipitation extremes, where analog-based and deterministic methods often fail (Sachindra et al., 2018;Reddy et al., 2023). We therefore introduce the concept of super-resolving extremal distributions, where we learn to increase the spatial resolution of a distribution's parameters. Friederichs (2010) previously downscaled precipitation block maxima by targeting GEV parameters from large-scale circulation features. However, we focus here on cases where convection-permitting extremes (1-5 km) are at least partially resolved at coarser resolutions (12-50 km), consistent with the HighResMIP2 range targeted by nextgeneration of Earth system models (Roberts et al., 2025), allowing a super-resolution setup. As illustrated in Figure 1, our generative super-resolution model combines regressionbased modeling with extreme value theory (Vrac & Naveau, 2007). It learns the parameters of the generalized extreme-value (GEV) distribution governing fin-resolution precipitation maxima from those governing their coarse-resolution counterparts. To ensure interpretability, we use (VGAM) and (VGLM), which allow us to specify a distribution family and examine how features influence predicted distributions through spline terms. This approach offers three key advantages: low training costs, reliable probability estimates, and transparent model behavior. To quantify our models' robustness to climate change, we adopt a pseudo-reality framework over Switzerland, where complex topography-ranging from alpine areas above 4,000 m to low-lying regions below 200 m-drives strong spatial variability in precipitation. Observations show a positive trend in annual maximum daily precipitation at 91% of Swiss stations, with 10-minute summer precipitation intensities increasing at 5.7% per decade (Bauer & Scherrer, 2024;Scherrer et al., 2016). Convection-permitting model en- sembles reproduce these changes more faithfully than coarser models and project a 6-7% increase in heavy summer precipitation intensity per degree of warming, despite an overall summer mean decrease (Estermann et al., 2025;Ban et al., 2021). As summer precipitation is dominated by convective events, it remains particularly challenging to downscale without underestimating variability (Zubler et al., 2014;Schmidli et al., 2007), though this can be partly mitigated by explicitly incorporating temperature dependence in precipitation scaling (Moraga et al., 2024). Motivated by the challenge of super-resolving extreme summer precipitation, we base our pseudo-reality experiment on the present-future convection-permitting simulation pair of Hentgen et al. (2019) (Section 2). We then formalize the task of super-resolving distributions and define the "robustness gap" to quantitatively assess model generalization across climates (Section 3). Applied to the super-resolution of the GEV using VGAMs and VGLMs (Section 4), our diagnostic framework explains generalization errors across quantiles, interprets model behavior via splines, and identifies performance limits when the super-resolution factor becomes too large (Section 5). We conclude in Section 6. We super-resolve hourly precipitation model output over Switzerland for 11 European summer (JJA) seasons in both historical (1999-2009) and projected (2079-2089) climates.
Swiss summers are characterized by weak synoptic forcing, frequent afternoon thunderstorms and convective rainfall. Data are derived from simulations using the regional weather and climate model COSMO (Consortium for Small-scale Modeling), run over a European domain of approximately 3000x3000 km at 2.2 km grid spacing. Full simulation details are available in Hentgen et al. (2019), with a brief overview provided in Section 2.1. Section 2.2 describes the construction of the coarse-resolution field from the fine-resolution data. The elevation statistics employed are detailed in Section 2.3, and the training, validation, and test set configurations are outlined in Section 2.4. We calibrate and test our method using regional climate simulations at 2.2 km horizontal resolution. Kilometer-scale climate modeling has become increasingly common in recent years (Ban et al., 2021;Stevens et al., 2019), offering several advantages over coarser-resolution simulations. The finer representation of topography and land surface enables more realistic precipitation patterns in complex terrain such as the Alps. In addition, vertical air motion is explicitly resolved by the governing equations, bringing the model formulation closer to physical first principles. Compared to coarser resolutions that require convective parameterization, explicit convection improves the realism of the hydrological cycle, particularly for extreme precipitation and its associated mechanisms and feedbacks (Prein et al., 2015;Schar et al., 2020;Lenderink et al., 2025). To assess generalization, we use model output from a simulation based on the pseudoglobal warming (PGW) method (Rasmussen et al., 2011;Adachi & Tomita, 2020;Schar et al., 1996). PGW aims to simulate a warmer climate by preserving the spatiotemporal structure of historical weather patterns. A regional simulation is first performed using reanalysis-based boundary conditions, referred to as the control simulation (CTRL). This serves as the baseline for comparison with the PGW scenario. A physically consistent climate change signal ("climate delta") is then applied to these boundaries, and the simulation is repeated. The deltas are derived from multi-year climatological means, representing differences between the MPI-ESM-LR global climate model's future (RCP 8.5) and historical scenarios (Kroner et al., 2017), and vary with latitude, longitude, elevation, and month. The PGW simulation hence retains the sequence of weather events from CTRL, adjusted by the imposed climate signal. As such, the PGW simulation primarily captures thermodynamic changes (e.g., increases in temperature and moisture) while preserving large-scale circulation patterns from the historical record (Hall et al., 2024). This reduces confounding effects from circulation biases, providing a dataset of intermediate complexity with greater control for evaluating the generalization of statistical downscaling methods to warmer climates. To ensure comparability across climates, we adopt an idealized super-resolution framework in which the coarse-resolution inputs are coarsened versions of the fine-resolution targets. Hourly precipitation data are aggregated to coarser grid spacings via mean pooling, i.e., using a square spatial filter of prescribed length on the native COSMO grid. To assess robustness across not only climates but also super-resolution factors, we apply pooling lengths of 13.2 km, 26.4 km, and 52.8 km. We then compute monthly maxima on the coarsened low-resolution data (features) and on the fine-resolution data (targets). As shown in Figure 2, decreasing spatial resolution leads to substantial information loss. Precipitation values become more uniform, with reduced means and variances of extremes consistent with the geostatistical concept of "change of support" (Park, 2013;Onibon et al., 2004). The highest extremes are particularly affected, with distribution tails shortening as resolution decreases. This is due to the localized nature of extremes: coarse-resolution cells blend multiple events within a grid block, smoothing out their intensity and variability. Future climates exhibit broader, flatter distributions with heavier tails and higher standard deviations, indicating increased dispersion. These patterns remain consistent across all resolutions in both climates. To incorporate elevation as a covariate, we calculate spatial statistics from the digital elevation model used by COSMO. Specifically, we compute two statistics from the 2.2 km-resolution elevation field h by defining a circular neighborhood of radius R centered on each grid point: the mean elevation h m and the standard deviation h s , both calculated over the values of h within the circle. The quantities h, h m , and h s are all considered during feature selection, and the radius R is treated as a model hyperparameter. To prevent overfitting and to ensure objective model evaluation, we define spatially separated splits (Valavi et al., 2019;Brenning, 2012), shown in To better understand the super-resolution of extremes and its sensitivity to spatial resolution and climate change, this section presents novel tools to characterize the distribution of extreme events at high spatial resolution, using information from coarser-resolution distributions and auxiliary variables. We begin by introducing and formalizing the concept of super-resolving parametric distributions, which enables the inference of fine-scale extremes from aggregated inputs. We then examine the spatial scales at which coarse-resolution predictors lose their capacity to inform fine-scale extremes. Finally, we propose strategies for assessing the robustness of data-driven models in the context of non-stationary climate conditions. Although super-resolution of physical fields has been widely studied and applied, the concept of super-resolution of distributions is, to the best of our knowledge, novel. Access to fine-resolution distributional information is essential for risk assessment and impact studies at fine spatial scales, as it enables the computation of various distributional summaries-such as return periods and exceedance probabilities-that are critical for decision-making under uncertainty. Let X denote a physical field-such as precipitation intensity-defined over a spatial domain S. In statistical terms, X is a random field, and the value of the field at any location s  S is a random variable X(s), whose cumulative distribution function (CDF) is given by For simplicity, we assume that this distribution is parametric and characterized by a parameter vector  s  R p . Consider two spatial grids, G 1 and G 2 , such that the resolution of G 1 is coarser than that of G 2 . We assume that the distribution of X is known at each location s 1  G 1 , but unknown at the finer-resolution locations s 2  G 2 .
The goal of super-resolution of distributions is to infer the distribution of X(s) at each fine-resolution location s  G 2 , using the known distributions at coarse-resolution locations in G 1 , possibly along with auxiliary covariates such as topography. A concrete example for G 1 and G 2 consists of regular two-dimensional grids with mesh sizes  1 > 0 and  2 > 0 and side length m, respectively: (2) where  1 >  2 , so that G 2 m has higher spatial resolution than G 1 m . Let  1 denote the concatenation of all parameter vectors  s for s  G 1 . The task of super-resolving distributions then consists in estimating, for each s  G 2 , where  s represents a set of features specific to location s (e.g., topographic or land-use characteristics), and f is a function-possibly learned from data-that maps coarse-resolution distributional information and local covariates to fine-resolution distribution parameters. As spatial resolution decreases, fine-resolution information about the variable of interest is progressively lost due to the homogenization of values across larger grid cells. This degradation may cause super-resolution techniques to fail. The objective of this section is to introduce a general and easy-to-implement methodology for identifying the resolution threshold at which super-resolution becomes ineffective-specifically, when the distribution of the variable of interest at low resolution provides little to no insight into its distribution at higher resolution. For clarity of exposition, we assume the existence of two regular grids, G 1 m and G 2 m , as introduced in Section 3.1, with  1 >  2 > 0, and where  1 is an integer multiple of  2 . While this assumption simplifies the presentation, our approach naturally extends to more general settings involving arbitrary grids. A natural approach for addressing the question posed in this section is to examine the spatial correlation function of the field of interest (e.g., precipitation intensity), denoted by X. This function captures how similar the values of the field are at different locations, depending on their spatial separation-providing insight into the scale and structure of spatial patterns, such as whether events are localized or spread out over larger regions. In particular, it allows us to assess whether the value of the field at a distant site still carries meaningful information about the value at a reference point, which is crucial for understanding the spatial coherence of the field and identifying the resolution beyond which super-resolution techniques may no longer yield significant benefits. However, in most cases, the field X is not stationary, and we therefore standardize it so that it has zero mean and unit variance at each grid point, resulting in a field denoted by X. In practice, we compute, for each s  G 2 m , the temporal mean and temporal standard deviation and define X as This standardization has the advantage that the resulting field X can reasonably be treated as second-order stationary, meaning it has a constant mean across space and a covariance structure that depends only on the relative displacement between locations, not their absolute positions. We denote by  X (r), for r  G 2 m , the spatial correlation function of the transformed field. Under the assumption of isotropy-which we adopt here-this function depends only on the distance between locations, and is thus written  X (r). Let N denote the number of time points (equal to 33 in our application) and M the number of fine-resolution grid points. We define as the set of all pairs of sites whose pairwise distances fall within a ball of radius r centered at r, and denote by M r the number of such pairs. Then the spatial correlation function of X is defined by where Xi (s j ), std are the spatial mean and the spatial variance of the field Xi observed at time i. In (7), the term within the sum over i corresponds to the empirical estimator of the spatial correlation function at distance r for the i-th temporal replicate of X. These terms are then averaged over time. Let B 1 denote a generic block of side length  1 , and let B 2 i represent the i-th subblock of side length  2 contained within B 1 . One way to define a resolution threshold below which super-resolution becomes ineffective is to identify the smallest  1 , denoted  l , such that the spatial correlation between the most distant sub-blocks within B 1 falls below a predefined threshold  (e.g., 0.1), which may depend on the variable of interest. Formally, we define where d(, ) denotes the distance between the barycenters of the respective blocks. The term involving the maximum captures the largest separation between any two sub-blocks of size  2 within a block of size  1 . If this separation is large enough for the correlation between the blocks to drop below , it indicates that the distribution at the coarser resolution  1 no longer provides meaningful information about the finer-scale distribution at resolution  2 . Since the spatial correlation function  X is typically continuous,  l can often be characterized as the solution to the equation: Up to this point, we have defined a resolution threshold for the breakdown of superresolution based on the spatial correlation function of the field of interest. A complementary perspective involves identifying the resolution at which auxiliary variables-such as topography-become more informative than the coarse-resolution representation of the target variable itself. In the case of rainfall, for example, relevant auxiliary variables may include topography, land use, or the dot product between wind vectors and topographic slope (to account for the orientation of terrain relative to prevailing circulation). These variables often retain fine-scale spatial structure by capturing terrain-driven features that strongly influence the underlying physical processes. As a result, they may explain a substantial portion of the variability in extremes of the target variable and, in some cases, offer greater predictive power than the coarse-resolution distribution of the variable of interest.
We denote by Y the field corresponding to the alternative variable (e.g., topography). We use X and Y to denote the potentially normalized versions of X and Y , respectively. Depending on the context, normalization may or may not be applied to these fields. When the alternative variable Y is time-invariant-as is the case for topography at the considered timescales-we do not normalize it, and thus set Y = Y . Similarly, when X represents precipitation and Y is topography, as in our application, we avoid temporal normalization of X in order to preserve systematic relationships with elevation, and therefore set X = X. To assess the spatial relationship between X and Y , we compute their spatial cross-correlation function, defined, similarly as in (7), by This function quantifies how strongly the values of the two fields are related across space, depending on their separation. To determine the resolution at which the alternative variable Y becomes more informative than the coarse-resolution representation of X, we define This expression seeks the smallest block size  1 such that the spatial correlation between distant sub-blocks of X becomes lower than the cross-correlation between X and Y at the same spatial scale. In other words, it identifies the resolution beyond which the alternative variable Y provides more useful information about the fine-resolution structure of X than X itself at low resolution. As before, since  X, Y is continuous,  l is typically the solution of The overall goal of this section is to (i) formalize the notions of generalization and robustness abilities of a model in a climate change context; (ii) propose a concrete solution tailored to estimating the quantiles of fine-resolution distributions. To do this, we consider three strategies: (i) a model trained on present-day data and applied to present conditions; (ii) a model trained on present-day data and applied to future conditions; and (iii) a model both trained and evaluated on future data. Let y denote an arbitrary characteristic of the distribution under consideration (e.g., the mean, a quantile, etc.) at a given grid point. Throughout this section, to maintain notational simplicity, we omit explicit reference to grid-point specificity and use lowercase letters even for random variables. Let y P P be the predicted quantity for the present climate using a model trained on present-day data, y P F be the predicted quantity for the future climate using a model trained on present-day data, and y F F be the predicted quantity for the future climate using a model trained on future data. Moreover, let y P and y F denote the observed (true) quantities in the present and future climates, respectively. Let l be a generic point-wise loss function that measures the discrepancy between predicted and true values, and let E denote the expectation operator with respect to the spatial dimension. The use of this expectation enables the integration of pointwise losses over the entire grid. A concrete example of loss function is the squared distance, defined by l(y, y) = (y  y) 2 . We can now define the generalization gap (that can also be called extrapolation gap in our context), as Generalization Gap def = E l(y F , y P F )  E l(y P , y P P ) , which quantifies the performance gap between a model trained and evaluated on presentday data versus the same model evaluated on future data. It captures both the inadequacy of parameters learned under current climate conditions when applied to future climates, and the impact of shifts in the model's covariates over time. Thus, this generalization gap can be challenging to interpret, as a high value-indicating poor performance under future climate conditions-may primarily result from shifts in covariates, reflecting fundamentally different environmental conditions, rather than a failure of the model itself. To overcome this, we introduce the notion of robustness gap, which is defined by This metric compares two models evaluated on the same future data, effectively removing the influence of covariate shifts. It quantifies how well a model trained on presentday data performs under future conditions, using a model trained on future-data as reference, thereby serving as an indicator of its transferability across climates. A low value indicates that the model parameters are robust to climate shifts, whereas a high one suggests that a model trained on present-day data may not be directly applicable under future climate conditions. We now focus on the specific case where the previously defined quantity y corresponds to a quantile at a generic level   (0, 1). Quantiles are particularly relevant in the context of risk assessment, as they are directly linked to the concept of return periods commonly used in environmental sciences, especially hydrology. The T -year return level z T , defined as the level expected to be exceeded on average once every T temporal units, corresponds to the -quantile q  of the fitted distribution, with  = 11/T . Return levels provide an interpretable and widely used risk metric; for example, the 100year return level represents the rainfall intensity expected to be exceeded, on average, once per century. Let q P F, denote the predicted quantile for the future climate obtained using a model trained on present-day data, and q F F, the predicted quantile obtained using a model trained on future data. The observed (true) quantiles in the present and future climate are denoted by q P  and q F  , respectively. In this context, the robustness gap introduced above becomes Robustness Gap = E l(q F  , q P F, )  l(q F  , q F F, ) . ( From now on, to explicitly link the robustness gap with quantile generalization errors, we use the pinball loss (also known as quantile loss) to evaluate quantile predictions at level   (0, 1). Given a predicted quantile q and an observed (true) value q  , the loss is defined as where I {} denotes the indicator function that equals 1 if the condition in the subscript is satisfied and 0 otherwise.
Although widely used-particularly in machine learningthis loss function may appear less intuitive than alternatives such as the squared loss, which could also have been employed. Nevertheless, this asymmetric loss function is particularly well-suited for evaluating quantile predictions, as it imposes different penalties for over-and under-predictions depending on the quantile level . Specifically, underestimations are penalized more heavily for high quantiles, while overestimations incur greater penalties for low quantiles. In addition, this loss function facilitates a straightforward decomposition of the robustness gap, enabling a more interpretable analysis of model performance. We now introduce two interpretable quantities that will naturally appear in the decomposition of the pointwise robustness gap (see below). Let  F, = q F F,  q F  , which is the fit bias of the model trained on the future climate for the quantile at level , and   = q P F, q F F, , as it quantifies the sensitivity of the predicted quantile to the choice of training data (whether from the present-day or future climate). We expect  F, to be small if the model is adequate and has been properly trained. With these notations, it is straightforward to obtain that and the term within the expectation in ( 13) is this difference simplifies to the sum of two terms: one involving the product of   with a number in (0, 1) and the other involving the product of  F with an indicator function difference that takes values in {1, 0, 1}. Provided the model is suited and well calibrated,  F is small and so is the second term. Two interesting limiting cases emerge. If   = 0 (i.e., q F P, = q F F, ), then the pointwise robustness gap is zero for any fit bias  F . If  F = 0 (i.e., the future-trained model is perfectly calibrated), then showing that the degradation in robustness is directly proportional to the quantile shift and exhibits asymmetry with respect to . In the case where the object of interest is a parametric distribution, it is useful to investigate which parameters are primarily responsible for the model's lack of robustness under a climate shift. This can be achieved by expressing   as below. Let us assume that we have a parametric distribution with vector  = ( 1 , . . . ,  p )  , where  denotes transposition. Then, the quantile robustness gap   can be decomposed in terms of the individual gaps coming from each parameter: where q  denotes the predicted -quantile of the distribution, and q  / k reflects the sensitivity of the quantile to variations in that parameter. When incorporated into ( 16), this decomposition establishes a general framework for analyzing the contribution of individual parameter shifts to the overall quantile robustness gap. It maintains the generality of the expression, ensuring applicability across a wide range of parametric distributions. Overall, our framework provides a comprehensive understanding of the robustness gap-clarifying not only when the super-resolution model generalizes effectively, when its performance begins to degrade, and when it ultimately fails, but also uncovering the underlying factors driving these behaviors. This section introduces the methodological tools necessary to apply the theoretical framework developed in Section 3 to rainfall extremes. We implement the super-resolution framework described in Section 3.1 within the context of the parametric Generalized Extreme-Value (GEV) distribution, which is well-suited for the statistical modeling of maxima. Section 4.1 formally presents the GEV distribution. We then define, in Section 4.2, the super-resolution function f appearing in (4), which maps the characteristics of the fine-resolution distributions from the parameters of the coarse-resolution distributions and a set of auxiliary features. Finally, Section 4.3 outlines the procedures used for feature selection and hyperparameter tuning to ensure optimal performance of the superresolution model. To characterize extreme hourly precipitation, we adopt the Generalized Extreme-Value framework from Extreme-Value Theory (EVT)-a branch of statistics focused on modeling the behavior of distribution tails. Unlike classical methods that describe central tendencies such as means and variances, EVT provides asymptotically justified models for block maxima or threshold exceedances, making it particularly well suited for assessing the risk of rare, high-impact events like extreme hourly rainfall (Coles, 2001). Let M n = max{Z 1 , Z 2 , . . . , Z n } denote the sample maximum of a sequence of random variables Z 1 , Z 2 , . . . , Z n . The subscript may represent, for example, the time index in hours, with each variable corresponding to the rainfall amount measured during the preceding hour. In this context, M n represents the maximum hourly rainfall observed over a period of n hours. Under fairly mild conditions, it is known that, for sufficiently large n, the sample maximum M n approximately follows the GEV distribution, whose CDF is given by and defined on the set {z : 1 + (z  )/ > 0} with   R,  > 0 and   R. The location parameter  shifts the distribution along the real line, the scale parameter  controls the dispersion, and the shape parameter  governs the heaviness of the tail. High hourly precipitation events lie in the far right tail of the rainfall distribution, where data are sparse and extremes dominate impacts. This key result underpins the block maxima method, which partitions a time series into non-overlapping blocks, extracts the maximum value from each block, and fits a GEV distribution to the resulting maxima-typically using maximum likelihood estimation; see Appendix A for the expression of the log-likelihood. In this study, the GEV distribution is independently fitted at each grid point for both fine-resolution and coarsened (coarse-resolution) datasets. The fine-resolution GEV distribution serves as the reference, while the coarse-resolution distributions provide the baseline for evaluating our approach. A key limitation is the relatively short data record, with only 11 years available for both present-day and future climate scenarios.
To increase the number of maxima considered in the block-maxima approach, we treat the monthly maxima from June, July, and August as independent realizations of extremes representative of a typical summer month, thereby increasing the number of extreme values per grid point from 11 to 33. To ensure temporal independence among these monthly maxima, a minimum separation of five days between events is enforced. While we adopt a block-maxima framework, alternative approaches-such as the r -largest order statisticscould be explored in future work. We now specify the form of the function f introduced in (4), adapted to our specific setting in which the underlying parametric distribution is the GEV distribution. To flexibly capture how the features (parameters of the GEV distributions at coarse-resolution and auxiliary features) influence each parameter of the GEV distribution, we use the (VGAM) framework (Yee, 2015). VGAM extends the familiar Generalized Additive Model (GAM) by allowing for multiple response variables. In our case the vector of response variables is composed of the GEV parameters. Each GEV parameter is linked to features x i = (x i1 , . . . , x ip )  , i = 1, . . . , m via an additive predictor in the VGAM framework (Yee, 2015): where log denotes the natural logarithm,   ,   and   are the intercepts for the parameters ,  and ; f ,k (), f ,k () and f ,k () are possibly smooth functions of the k-th feature, typically represented using basis expansions (e.g. splines). The log-link on  i ensures positivity of the scale parameter. In the case where all f .,k are linear, the class of VGAMs reduces to the subclass of so-called (VGLMs). We typically model  and log  as smooth functions of the features, but-for stability and identifiability-take  to be constant (i.e. f ,k  0 for all k). Indeed  estimation is notoriously imprecise even with large samples of block maxima and allowing  to vary with covariates often yields unstable fits (Coles, 2001;Davison & Smith, 1990). Moreover, extreme events are inherently rare, so there is typically insufficient data within each covariate "slice" to support a reliable smooth trend in , and attempts to do so can lead to over-parameterization and degraded predictive performance (Eastoe & Tawn, 2009). Finally, empirical comparisons indicate that introducing time-or covariate-dependence in the location (and sometimes scale) parameters captures the bulk of observed non-stationarity, while varying  provides only marginal improvements at the cost of substantially increased uncertainty (Katz et al., 2002;Coles, 2001). We estimate all unknown quantities by maximizing a penalized log-likelihood (see (A2) in Appendix A2), which balances model fit and smoothness. The penalization discourages overfitting by controlling the complexity of the smooth functions, ensuring that the estimated relationships remain interpretable and generalize well to unseen data. Covariate selection is performed using a forward selection procedure, guided by model performance on the hold-out validation set described in Section 2.4. Since our objective is to accurately model the fine-resolution distribution, model fit is evaluated using a statistical distance between the modeled and empirical distributions at each grid point. We adopt the Cramer-von Mises (CVM) distance, a robust metric that quantifies the discrepancy between two CDFs. For a given grid point, let F n denote the empirical CDF derived from the observed sample x 1 , . . . , x n , and let F (x, ) represent the theoretical CDF parameterized by . The CVM statistic is defined as It integrates the squared difference between empirical and modeled CDFs across the entire distribution, providing a sensitive and stable assessment of overall distributional fit. In practice, we employ a computationally efficient approximation: which assumes continuity of the theoretical distribution F (x, ). Here, the term (2i 1)/(2n) corresponds to the expected CDF values under a uniform distribution, serving as a reference for comparison. The squared differences quantify the deviation of the theoretical CDF from this uniform benchmark at the observed data points. Since we work within the GEV family, this formula can be adapted by explicitly using the known analytical form of the cumulative distribution function (Equation 19). The total CVM score is computed by summing the individual CVM distances across all grid points in the validation set and, finally, covariates yielding the lowest aggregate CVM score are selected. A preliminary step is conducted to assess the most informative topographic covariates. For this purpose, elevation-based features are calculated from the digital elevation model using circular neighborhoods of varying radii R  {10, 20, . . . 100} km. For each radius, spatial statistics such as the mean and standard deviation of elevation are extracted and forward selection is applied exclusively to this subset. The results indicate that covariates derived using a 50 km radius yield the best model fit as measured by the CVM criterion. Therefore, these covariates are retained for inclusion in the full feature selection procedure. We first investigate on present-day data the performance of our super-resolution models from 13.2 km to 2.2 km trained on present-day data. We then study how the performance of our super-resolution models varies when decreasing the resolution of the input data to which they are trained (26.4 and 52.8 km, instead of 13.2), and explain the observed behavior using the tools introduced in Section 3.2. Finally, we investigate, using the methodology developed in Section 3.3, the robustness of our super-resolution models in a warmer climate. Our best-performing model for super-resolving precipitation distributions from 13.2 km to 2.2 km is a VGAM using eight features. These include two topographic features -elevation (h) and its local average (h m )-and six features derived from the coarse-resolution GEV distributions: the location ( 1 ), scale ( 1 ), and shape ( 1 ) parameters of the nearest block, along with the location ( 2 ) and scale ( 2 ) of the second-nearest block. The VGAM yields interpretable expressions (see ( 20)) for the parameters of the fine-resolution GEV distribution: where x represents the covariates at a specific grid point.
Therefore, at each grid point, the parameters of the target GEV distribution (, , and ) are obtained by evaluating each spline at the corresponding covariate value at that point, summing all contributions, and adding an intercept. Coarse-resolution features are dominant in this setup: they account for 81% of the model's explanatory power, as determined by the AIC drop when each variable is removed . Among the coarse precipitation features,  1 and  1 are most informative, while  2 and  2 contribute less, suggesting that most of the information comes from the nearest coarse grid cell. The second-nearest block serves primarily as a correction, which is an expected and reassuring result. This configuration corresponds to a canonical superresolution setting, where most of the predictive skill stems from the coarse field of interest rather than from external covariates. The advantage of our trained models is 2 km-resolution precipitation). Both VGLM and VGAM outperform the coarse-resolution baseline, but the VGAM achieves lower errors across the training, validation, and test sets in the reference climate. We attribute this gain to the model's greater 0 1000 2000 3000 4000 5000 6000 showing each feature's explanatory power for the target GEV's location (), scale (), and shape () parameters when super-resolving from 13.2 km (a) and 52.8 km (b). Coarse-resolution features (gray) become less informative than elevation statistics (blue) as the super-resolution factor increases from 6 (left) to 24 (right). ability to capture nonlinear relationships, underscored by VGAM's consistent outperformance of VGLM. One of the key advantages of VGAMs lies in their interpretability through smooth additive functions, which describe how each feature contributes to the predicted parameters of the target distribution. Figure 6 shows the learned splines for each feature in (24), revealing meaningful relationships in the present climate (brown lines). In the top two rows, which govern the location parameter , the spline f ,1 (top left) is monotonically increasing, indicating that higher coarse-scale location values predict higher fine-scale extremes. The spline f ,2 (top middle) also increases, but with a lower slope, suggesting that the second-nearest block provides only a secondary correction. The spline f ,h (center left) rises with elevation before flattening at higher altitudes, consistent with orographic enhancement up to a threshold. The average elevation spline f ,hm (center) refines this relationship by incorporating topographic context over a broader scale. In the bottom row, which models the logarithm of the scale parameter log(), similar patterns emerge. The spline f ,1 (bottom left) increases nearly linearly, reinforcing the importance of coarse-scale scale values. The spline f ,2 (bottom center) also increases but shows a threshold effect, indicating that secondary features affect the spread only beyond a certain magnitude. The function f ,1 (bottom right) shows a sharp, monotonic increase, revealing that regions with higher location values also exhibit greater variability. To assess model generalization across super-resolution factors, we train VGAMs on 13.2 km, 26.4 km, and 52.8 km input data, keeping the feature set fixed to that selected at 13.2 km. At 26.4 km, the model coefficients and splines remain similar, but performance degrades , indicating information loss from spatial coarsening . The 52.8 km model diverges further: it selects fewer features and relies more on topography than on coarse GEV parameters. This shift marks a departure from a canonical superresolution regime, where predictive skill stems from the coarse-resolution target field, toward a topographic downscaling setup. This is reflected in the 52.8 km model equations: where we remind the reader that h m and h s denote the mean and standard deviation of local elevation over the circular neighborhood of radius R. Over the Swiss Alps, the 52.8 km model using GEV features underperforms compared to an elevation-only model, reflecting the loss of predictive content in coarse precipitation fields. At such resolutions, alternative downscaling methods-such as perfectprognosis approaches that incorporate temperature or humidity-or deep learning architectures exploiting broader spatial contexts may improve performance. To better understand the observed performance drop, we apply the methdology described in Section 3.2 to the the 2.2 km precipitation and elevation fields. More precisely, we investigate the behavior of the spatial correlation and autocorrelation functions defined in ( 7) and (10). As illustrated in Figure 7a, the spatial correlation of precipitation decreases approximately exponentially with distance, whereas the spatial cross-correlation between precipitation and elevation remains nearly constant. These behaviors can be modeled using the following functional forms, as shown in Figure 7b: 7) and ( 10), respectively). Gray labels on the x-axis indicate block sizes matching the maximum radial distances for 13.2 km, 26.4 km, and 52.8 km grids. (b) Fitted exponential decay (spatial correlation) and constant fit (spatial cross-correlation). Their intersection at a block size of 30.8 km marks the limit beyond which spatial detail is insufficient for super-resolution. where r is radial distance, and  P P and  P H are the spatial correlation and cross-correlation functions, respectively. Interpreting r as the maximum radial distance between two blocks of side length 2.2 km within a block of side length x (maximum distance appearing in (9)), the curves intersect at r  38 km, or x  30.8 km. Beyond this, elevation becomes more informative than precipitation, marking the end of the super-resolution regime. This provides a heuristic threshold: once block size exceeds 30.8 km, coarse precipitation loses predictive advantage, and topographic features dominate. In flatter terrain, the heuristic could be adapted using alternative covariates or by identifying a correlation cutoff (e.g., 0.1) below which super-resolution is no longer effective; see details in Section 3.2. Finally, we use our distribution super-resolution framework to assess challenges in generalization/robustness across climates. Climate change is treated here as a domain shift in model features, affecting the location, scale, and shape parameters of the coarseresolution GEV distributions. As model errors increase in the future climate, where precipitation extremes intensify, we evaluate robustness by comparing the performance on future data of reference models (trained on present-day data) to that of models trained directly on future climate data .
This comparison enables an assessment of both spatial generalization and generalization across climates. When recalibrated on future data, the VGAM retains qualitatively similar splines as shown in Figure 6, but the intercepts for both the scale and shape parameters increase, with the shape parameter rising from 0.19 to 0.25. To diagnose the source of generalization errors, we analyze the robustness gap introduced in Section 3.3, defined via the pinball loss, which quantifies the discrepancy between predicted quantiles for future climates when using models trained in present and future climates. In the case of the GEV distribution, we decompose the quantile gap   (difference between the quantiles for future climate, computed from a model trained on present-day data and a model trained on future data) at a specific level   (0, 1) us- ing a first-order Taylor expansion: where R  is the residual, and the partial derivatives of the GEV quantile function (see (B1) in Appendix B) are given by for the detailed derivation of these terms, see Appendix B. By incorporating this decomposition into (16), we obtain the decomposed form of the pointwise robustness gap G in the GEV case: where Interestingly, for lower quantile levels (below the 30 th percentile), the VGLM model shows better generalization, as evidenced by consistently smaller robustness gaps. Looking at panel (b), we see that the term of the location parameter, T  , plays a dominant role. For the VGAM, its contribution increases earlier and more sharply than in the VGLM model, as shown by the difference between the solid and dashed blue lines. In such cases, the location parameter accounts for nearly all of the robustness behavior, while the contributions from scale and shape remain close to zero. In this work, we introduced a framework for super-resolving distributions within interpretable statistical models, focusing on risk-relevant quantities such as high quantiles (return levels). This approach is particularly valuable for impact modeling, as it enables the estimation of risk measures at resolutions suited to downstream applications. We developed a method to identify the resolution scale at which super-resolution begins to fail and proposed the novel concept of a robustness gap, which we analyzed in the context of quantile estimation using the pinball loss. Applied in a pseudo-reality experiment over Switzerland, our framework showed that the quantile-wise robustness gap is an effective diagnostic for evaluating how well models trained on present-day data generalize to warmer climates. Leveraging a GEV distribution parameterized via VGAMs, we could pinpoint which model components-parameters and splines-contribute to robustness failures, offering insights into model limitations and guiding future improvements. Overall, the proposed methodology enables a tractable decomposition of spatial generalization and robustness across climates. Several directions could further improve this framework. First, our use of mean pooling to approximate coarse-resolution fields is a first-order simplification; future work could explore alternative filters (e.g., Gaussian) or incorporate temporal and physical biases through explicit bias correction. Second, we assumed a constant shape parameter for the GEV distribution, which may not reflect how extremes respond to climate shifts. Using the r-largest values approach could stabilize tail estimates and support state-dependent shape parameter modeling, helping address an open question: which covariates most strongly influence the shape parameter ? Third, since the original pseudo-global warming simulation used here (Hentgen et al., 2019), the pseudo-global warming method has been refined (Brogli et al., 2023;Heim et al., 2023), and future work should update the analysis using the latest available simulations. Finally, tailored versions of the GEV, such as the blended GEV (Vandeskog et al., 2022), may improve the stability and realism of tail estimates. Overall, our results suggest that combining probabilistic super-resolution with quantifiable generalization/robustness diagnostics provides a principled framework for modeling extremes under distributional shifts, with potential relevance to other fields that study rare events in non-stationary systems. A static version of the code, along with the data necessary to reproduce the results presented in this manuscript, is publicly available at the following DOI: .org/10.5281/zenodo.15863970. This study was conducted by researchers based in Switzerland (LL, EK, DL, GM, VC, TB) with contributions as follows: Conceptualization, Methodology, Resources, Validation, and Writing -review & editing: all authors; Data curation: LL, DL; Formal analysis, Investigation, Software, Visualization, Writing -original draft: LL, EK, TB; The work relied exclusively on publicly available global and regional datasets, with no fieldwork, permits, or engagement with local or Indigenous communities required. The authors are committed to equity, transparency, and responsible research partnerships in line with The TRUST Code of Conduct for Equitable Research Partnerships. We thank Marine Berthier and Leo Micollet for preliminary work that laid useful foundations for this study, Linda Mhalla and Shivanshi Asthana for feedback that improved our manuscript, and and Christoph Schar, for performing the kilomter-scale PGW simulations. These simulations were conducted on the Piz Daint supercomputer of the (CSCS) under project ID pr144. EK would like to thank the Expertise Center for Climate Extremes (ECCE) at UNIL for financial support. TB acknowledges partial support from the (SNSF) under 10001754 ("RobustSR" project). -24-manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) We thank Marine Berthier and Leo Micollet for preliminary work that laid useful foundations for this study, Linda Mhalla and Shivanshi Asthana for feedback that improved our manuscript, and and Christoph Schar, for performing the kilomter-scale PGW simulations. These simulations were conducted on the Piz Daint supercomputer of the (CSCS) under project ID pr144. EK would like to thank the Expertise Center for Climate Extremes (ECCE) at UNIL for financial support. TB acknowledges partial support from the (SNSF) under 10001754 ("RobustSR" project). -24-manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Appendix A Log-likelihoods for the GEV parameters and VGAM/VGLM A1 Estimation of GEV parameters Assume that we have extracted m block-maxima denoted by M 1 , . . . , M m .
Under the independence and identical distribution assumption, the log-likelihood of the GEV parameters (, , ) is and the maximum-likelihood estimates , ,  are obtained by numerically maximizing l. Let us recall that in the VGAM framework, the GEV parameters are expressed as where   ,   and   are the intercepts for the parameters ,  and ; f ,k (), f ,k () and f ,k () are possibly smooth functions of the k-th covariate. We estimate all unknowns by maximizing a penalized log-likelihood over the parameter vector where b ,k (respectively b ,k and b ,k ) collects the basis coefficients for f ,k (respectively f ,k and f ,k ). Writing l() for the GEV log-likelihood (A1), the penalized objective is where S is a m x m penalty matrix (e.g. second-derivative penalty for splines), m being the number of basis functions used for that smoothing and  ,k  0 (respectively  ,k  0 and  ,k  0) is the smoothing parameter controlling the wiggliness of f ,k (respectively f ,k and f ,k ). We solve for  by a penalized Fisher scoring (or quasi-Newton) algorithm, iterating until convergence. Smoothing parameters  ,k ,  ,k and  ,k are chosen automatically, e.g. by minimizing a generalized cross-validation (GCV) criterion or by approximate restricted maximum likelihood (REML). Standard errors for  j (x) (j = 1, 2, 3) and for any function of the GEV parameters are obtained from the inverse penalized Fisher information. Appendix B Partial Derivatives of the GEV Quantile Function The quantile function at level   (0, 1) of the GEV distribution with location, scale and shape parameters ,  and , is In this appendix, we focus on the case   = 0. While  may approach zero in some regions, the GEV quantile function for   = 0 continuously converges to the  = 0 (Gumbel) case as   0. Therefore, using the   = 0 expression provides a unified and consistent formulation across all locations without requiring a separate derivation. In practice, this approximation is accurate and numerically stable even for small values of || (Coles, 2001). For the partial derivatives with respect to  and  we immediately obtain B3 Partial derivative with respect to  Let A = ( log )  for brevity. Then we have Applying the product and chain rules, we get where Putting everything together yields which completes the derivation of the three partial derivatives of the GEV quantile function used in the robustness gap analysis. Appendix A Log-likelihoods for the GEV parameters and VGAM/VGLM A1 Estimation of GEV parameters Assume that we have extracted m block-maxima denoted by M 1 , . . . , M m . Under the independence and identical distribution assumption, the log-likelihood of the GEV parameters (, , ) is and the maximum-likelihood estimates , ,  are obtained by numerically maximizing l. Let us recall that in the VGAM framework, the GEV parameters are expressed as where   ,   and   are the intercepts for the parameters ,  and ; f ,k (), f ,k () and f ,k () are possibly smooth functions of the k-th covariate. We estimate all unknowns by maximizing a penalized log-likelihood over the parameter vector where b ,k (respectively b ,k and b ,k ) collects the basis coefficients for f ,k (respectively f ,k and f ,k ). Writing l() for the GEV log-likelihood (A1), the penalized objective is where S is a m x m penalty matrix (e.g. second-derivative penalty for splines), m being the number of basis functions used for that smoothing and  ,k  0 (respectively  ,k  0 and  ,k  0) is the smoothing parameter controlling the wiggliness of f ,k (respectively f ,k and f ,k ). We solve for  by a penalized Fisher scoring (or quasi-Newton) algorithm, iterating until convergence. Smoothing parameters  ,k ,  ,k and  ,k are chosen automatically, e.g. by minimizing a generalized cross-validation (GCV) criterion or by approximate restricted maximum likelihood (REML). Standard errors for  j (x) (j = 1, 2, 3) and for any function of the GEV parameters are obtained from the inverse penalized Fisher information. Appendix B Partial Derivatives of the GEV Quantile Function The quantile function at level   (0, 1) of the GEV distribution with location, scale and shape parameters ,  and , is In this appendix, we focus on the case   = 0. While  may approach zero in some regions, the GEV quantile function for   = 0 continuously converges to the  = 0 (Gumbel) case as   0. Therefore, using the   = 0 expression provides a unified and consistent formulation across all locations without requiring a separate derivation. In practice, this approximation is accurate and numerically stable even for small values of || (Coles, 2001). For the partial derivatives with respect to  and  we immediately obtain B3 Partial derivative with respect to  Let A = ( log )  for brevity. Then we have Applying the product and chain rules, we get where Putting everything together yields which completes the derivation of the three partial derivatives of the GEV quantile function used in the robustness gap analysis.