4.1. Graph conceptualization

Each point within the low-resolution and high-resolution grids corresponds to a specific geographical location, suggesting that they can be naturally modelled as nodes. For convenience, we model both grids as a unified heterogeneous graph featuring two node- and two edge-types:

• • Low nodes: first set of nodes, generated from the points on the low-resolution grid with a spatial resolution of approximately $ 25 $ km.

• • High nodes: second set of nodes, created from the points on the high-resolution grid with a spatial resolution of $ 3 $ km.

• • Low-to-High edges: unidirectional edges, which connect Low to High nodes, ensuring each High node is linked to a fixed number $ k $ of Low nodes. These edges model the downscaling of atmospheric variable information from the Low nodes to the High nodes (Figure 2a).

  

  Figure 2. Graph conceptualisation: Low nodes (blue dots) and High nodes (red dots) and close-up of (a) Low-to-High unidirectional edges (orange), connecting Low nodes to High nodes (b) High-within-High bidirectional edges (red), linking High nodes.

• • High-within-High edges: bidirectional edges that capture relationships between High nodes based on an 8-neighbours approach, ensuring each node is connected to its eight nearest neighbours (Figure 2b).

4.2. Model design

Precipitation data contain a significant amount of zeros, as rain events only occur during a limited number of hours. In the case of the GRIPHO dataset, almost $ 90\% $ of the values fall below the meteorological threshold assumed as $ 0.1 $ mm, effectively rendering them as zeros. In view of this, we explored two different designs for the DL model. In both cases, the model outputs an estimate for the time step $ t $ based on a time series of predictors spanning $ \left[t-L,\dots, t\right] $ , where $ L $ is a hyper-parameter.

In the first approach, we addressed the challenge posed by zero precipitation values by adopting a Hurdle modelling scheme (Cragg, Reference Cragg1971). The method relies on the construction of two distinct models, which are subsequently combined: a Regressor and a Classifier. The Classifier is trained on the entire dataset and discerns between two classes: $ 0 $ , i.e., precipitation below the threshold, and $ 1 $ , i.e., precipitation above the threshold. Conversely, the Regressor is exclusively trained on targets where precipitation values exceed the threshold, and provides a quantitative estimation of hourly precipitation. During the inference phase, predictions from the Regressor and Classifier models are computed independently and then multiplied to yield a singular estimate of the precipitation value. We refer to this model design as RC (Figure 3a). In addition, we considered a second approach where a single Regressor is used. This model is trained using as target the full GRIPHO dataset, i.e., including instances with zero precipitation, and outputs a quantitative estimation of hourly precipitation. We refer to this model design as R-all (Figure 3b). To enhance training efficiency and focus on meaningful targets, time steps containing only values below the threshold were excluded. This resulted in a reduction of approximately 50% in the training set size for the RC Regressor, whereas the impact was negligible for the RC Classifier and the R-all Regressor, with $ 99.9\% $ of time steps retained in both cases. In both the RC and R-all approaches, the models are trained on the complete graph, while the loss is computed exclusively on nodes with valid target values by applying a masking strategy.

Figure 3. Schematic views of (a) RC, designed as a combination of Regressor and Classifier components, (b) R-all consisting of a single Regressor, (c) architecture, composed of four modules: a RNN-based pre-processor, a GNN-based downscaler, a GNN-based processor, and a FCNN-based predictor.

The two alternative designs offer valuable insights due to their differing model complexities and problem formulations. The RC case adopts a two-model structure, separating the classification of precipitation occurrence from the regression of its intensity. This allows for targeted learning but increases the number of parameters and model complexity. Moreover, the classification task remains particularly imbalanced and complex to solve, and even the data obtained by considering only the values larger than the threshold remain very skewed. In contrast, the R-all case adopts a single model, thus resulting in a simpler architecture with fewer parameters but requiring the model to simultaneously handle both occurrence and intensity prediction within a single task, thus increasing the problem complexity. In the experiments, we will analyse both model designs to see whether the added value justifies the cost of using two models instead of one.

4.3. Architecture

Regressor and Classifier components share the same structure, which consists of four primary modules (see Figure 3c). First, a pre-processor module, which acts at the Low nodes level and handles the predictors’ temporal component through the adoption of a recurrent neural network (RNN), specifically a gated recurrent unit (GRU). The RNN encoder captures the temporal dependencies across time steps, outputting a sequence that is flattened and passed through a fully connected layer to produce a fixed-dimensional latent representation. This encoding serves as input to the graph-based components of the model. The next stage involves the downscaler module, which uses a graph convolution (GC) layer to map the preprocessed atmospheric variables, represented as Low node features, to learned attributes on the High nodes. This transformation is crucial for bridging the different spatial scales within the input and output data. The downscaler incorporates additional high-resolution spatial attributes (elevation and land use data), ensuring that the model is well-informed about local geographical features. Following the downscaling step, the core of the model is a multi-layer processor network comprising several graph attention layers (GAT, Veličković et al., Reference Veličković, Cucurull, Casanova, Romero, Lio and Bengio2018). These layers are designed to dynamically attend to neighbouring nodes, thereby capturing complex spatial relationships. Each GAT layer is followed by batch normalization and ReLU activations to stabilise training and introduce non-linearity, respectively. The use of multiple GAT layers allows the model to progressively refine its understanding of spatial dependencies, essential for accurately predicting local precipitation. The final component of the architecture is the predictor, a fully connected neural network (FCNN) that takes the processed graph features and returns the desired output on each High node. The model is entirely implemented in PyTorch (Paszke et al., Reference Paszke, Gross, Massa, Lerer, Bradbury, Chanan, Killeen, Lin, Gimelshein, Antiga, Desmaison, Kopf, Yang, DeVito, Raison, Tejani, Chilamkurthy, Steiner, Fang, Bai and Chintala2019) and PyTorch Geometric (Fey and Lenssen, Reference Fey and Lenssen2019), using GraphConv (Morris et al., Reference Morris, Ritzert, Fey, Hamilton, Lenssen, Rattan and Grohe2021) and GATv2Conv (Brody et al., Reference Brody, Alon and Yahav2022) as GC and GAT layers, respectively.

The downscaling task requires the specification of the spatial and temporal length scales at which predictors can influence outputs. While this represents an assumption, it is grounded in the underlying physics of the problem and adjusted by empirical evidence. From a spatial point of view, we assume that the relevant region extends to approximately $ 50 $ km around the target location. From a time perspective, we assume that the precipitation prediction at time $ t $ is influenced by the atmospheric variables up to $ 24 $ hours before. These ranges are explicitly considered in designing the graph and the emulator architecture. Respectively, the spatial assumption is realised in choosing a value of $ k=9 $ , as the number of nearest Low nodes to which each High node is connected (see subsection 4.1), effectively representing the downscaling relationship in terms of graph edges. This assumption dictates the graph structure, but has no direct influence on the GNN parameters, as GNNs accept graphs with an arbitrary number of edges. The temporal assumption is instead reflected in choosing $ L=24 $ (see Subsection 4.2), thus designing the predictors’ time series as $ \left[t-24,\dots, t\right] $ , encompassing a total of $ 25 $ time steps. This assumption sets the length of the time series of predictors that are passed as input to the RNN module, thus defining the RNN sequence length parameter.

4.4. Regressor loss function

Mean square error (MSE) is the standard loss function for regression problems, yet it becomes less effective when the target data is highly imbalanced or skewed, as in the case of precipitation data. MSE leads to estimates that are biased towards frequent values, which is detrimental for modelling rare events, thus we chose to use a modified MSE loss function which can explicitly address imbalance and skewness. Multiple studies in the literature use weighted MSE loss for training in such unfavourable conditions (Wang et al., Reference Wang, Wang, Wang, Xue and Wang2022; Scheepens et al., Reference Scheepens, Schicker, Hlaváčková-Schindler and Plant2023). The most sensitive part is in the definition of an optimal weighting strategy, which should be consistent with the target data distribution and training objectives, e.g., giving more weight to the tail of the distribution, in order to estimate the extreme events. Recently, Szwarcman et al. (Reference Szwarcman, Guevara, Macedo, Zadrozny, Watson, Rosa and Oliveira2024) proposed a formulation to quantise the reconstruction loss and improve the synthesis of extreme weather data with variational autoencoders (VAEs). This modification of the MSE loss was designed to address the skewed distribution typical of weather data by giving more weight to rare, extreme values. The idea is to penalise the loss according to the observed values frequency, by quantising the target data and averaging losses for each bin. In this study, we adopt the same loss to train the Regressor, and we call it quantised MSE (QMSE) loss (Equation 4.1). The only parameters of QMSE are the bins into which the target data are quantised. During training, examples are dynamically assigned to the corresponding bins based on their true target values, ensuring that the weights in the QMSE loss reflect the actual value distribution within each batch. This approach is beneficial because batches are formed by randomly selecting time instants, and each batch includes all points in the graph for the chosen time instants. Consequently, the target distribution can vary slightly between batches.

(4.1) $$ \mathrm{QMSE}=\sum \limits_j^B\frac{1}{\mid {\Omega}_j\mid}\sum \limits_{i\in {\Omega}_j}{\left({y}_i-{\hat{y}}_i\right)}^2 $$

In Equation 4.1, $ j $ represents the bin index, defined based on a histogram of the training data, $ {\Omega}_j $ is the set of target indices whose values fall within bin $ j $ , thus $ \mid \Omega \mid $ is the observed frequency and $ 1/\mid \Omega \mid $ weighs the loss inversely to the frequency. The quantities $ {y}_i $ and $ {\hat{y}}_i $ are the true and estimated target values, respectively. Finally, in the training, we used a combined MSE-QMSE loss (Equation 4.2) with a coefficient $ \overline{\alpha} $ which accounts for the different scale and balances the contribution of the two terms.

(4.2) $$ L=\mathrm{MSE}+\overline{\alpha}\cdot \mathrm{QMSE}\hskip0.1em $$

4.5. Classifier loss function

The Classifier is trained using focal loss (FL) (Lin et al., Reference Lin, Goyal, Girshick, He and Dollár2020), specifically designed to address class imbalance during training. In this setting, standard cross entropy (CE) loss tends to focus on minimising errors for the majority class, often leading to poor performance on the minority class (Johnson and Khoshgoftaar, Reference Johnson and Khoshgoftaar2019). FL introduces a modulating term $ \gamma $ in the CE formulation to dynamically scale the CE loss. Thanks to this scaling factor, the contribution of easy examples is down-weighted and the model is guided to focus on hard examples. Additionally, a hyper-parameter $ \alpha $ helps to handle the class imbalance. The formulation of the FL loss is in Equation 4.3, where $ p $ is function of the Classifier output, i.e. depends on the input data, and $ y\in \left\{0,1\right\} $ is the ground-truth class.

(4.3) $$ {\displaystyle \begin{array}{c} FL\left({p}_t\right)=-{\alpha}_t{\left(1-{p}_t\right)}^{\gamma}\cdot \log \left({p}_t\right)\\ {}{p}_t=\left\{\begin{array}{ll}p& \mathrm{if}\;y=1\\ {}1-p& \mathrm{otherwise}\end{array}\right.\hskip2em {\alpha}_t=\left\{\begin{array}{ll}\alpha & \mathrm{if}\;y=1\\ {}1-\alpha & \mathrm{otherwise}\end{array}\right.\end{array}} $$

4.6. Training and evaluation

As mentioned in Section 2, the GRIPHO observational data are available for a period of $ 16 $ years, with spatial coverage of the entire Italian territory. To evaluate the emulator’s capability of generalisation in spatial domains not used during training, we chose to train only on northern Italy and use the whole peninsula for the inference phase (see Figure 4a). In this setting, the geographical area considered for training is approximately $ \mathrm{120,000} $ km $ {}^2 $ with around $ 400 $ Low nodes and $ \mathrm{14,000} $ High nodes. Instead, the evaluation area is approximately $ \mathrm{300,000} $ km $ {}^2 $ , with around $ 1000 $ Low nodes and $ \mathrm{33,000} $ High nodes. The time range spanned by the GRIPHO dataset is rather limited, which is usually the case when dealing with high-resolution observational data, which are quite difficult to find. We decided to use the first $ 15 $ years for training ( $ 2001 $ – $ 2006 $ and $ 2008 $ – $ 2015 $ ) and validation ( $ 2007 $ ) and leave the last available year ( $ 2016 $ ) to test the GNN4CD model in the reanalysis to observation downscaling task. Considering the number of time instants and High nodes, we obtain a ratio of approximately $ 75-12.5-12.5 $ for the train-validation-test datasets. In the RCM emulation setting, we considered the area of the Italian peninsula covered by RegCM and three different time slices of the RegCM simulations: historical ( $ 1996 $ – $ 2005 $ ), mid-century ( $ 2041 $ – $ 2049 $ ) and end-of-century ( $ 2090 $ – $ 2099 $ ). All projections were performed under the RCP8.5 scenario (Riahi et al., Reference Riahi, Rao, Krey, Cho, Chirkov, Fischer, Kindermann, Nakicenovic and Rafaj2011; IPCC, 2014), which represents a high-emissions pathway associated with the most pronounced climate change signal. This choice increases the difficulty of the emulation task, particularly for capturing changes in the distribution of extreme precipitation events.

Figure 4. (a) training (northern Italy) and inference (entire Italy) areas, (b) locations of original stations used to create the GRIPHO dataset, and (c) percentage of valid time steps for each station.

The RC Regressor and Classifier components and the R-all Regressor component were all trained separately for $ 50 $ epochs, i.e. approximately $ 24 $ hours each on $ 4\times $ NVIDIA Ampere GPUs on Leonardo, the new pre-exascale Tier-0 EuroHPC supercomputer, currently hosted by CINECA in Bologna, Italy (Turisini et al., Reference Turisini, Amati and Cestari2023). The trained emulator needs only few minutes to compute the hourly precipitation estimates for an entire year over Italy, much less than the dynamical downscaling of a CP-RCM, which needs approximately a couple of days on an equivalent high-resolution grid.

We used the validation year to empirically tune the hyper-parameters that appear in the model architecture and loss functions. Considering the computational cost of training the model, we opted for manual hyper-parameters tuning to limit the resource usage to what was strictly necessary, while still achieving quick convergence to good values. Systematic hyper-parameter tuning is expected to further improve performance. Thus, we view this not as a limitation, but as a promising direction for future enhancement, contingent on the availability of additional computational resources and time. For the Regressor components, we empirically tested the number of bins in QMSE and we concluded that it does not have a strong impact on the training, provided that the coefficient $ \overline{\alpha} $ is properly tuned. We chose to use logarithmically equispaced bins with a bin size of $ \log \left(0.5+1\right) $ . Empirical findings showed a trade-off when training the emulator with the combined MSE-QMSE loss, with lower values of $ \overline{\alpha} $ that favour average results, while higher values of $ \overline{\alpha} $ lead to improved accuracy in the tail of the distribution. We chose $ \overline{\alpha} $ $ =0.025 $ for the RC Regressor component and $ \overline{\alpha} $ $ =0.005 $ for the R-all model, both of which lean toward the latter behaviour. This choice reflects our interest in accurately capturing the full distribution of precipitation values, but with particular emphasis on the extremes. For the FL, we started from the parameter values suggested in Lin et al. (Reference Lin, Goyal, Girshick, He and Dollár2020) ( $ \alpha =0.25 $ and $ \gamma =2 $ ) and we operated a manual grid search to adjust these values to our specific case. We found that the FL loss used to train the RC Classifier component works best with parameters $ \alpha =0.75 $ and $ \gamma =2 $ .

6.1. Reanalysis to observation downscaling

The aim of this initial evaluation is to provide a comprehensive assessment of the GNN4CD behaviour in a setting that closely resembles the training environment.

We begin by assessing the PDFs and diurnal cycles (Figures 5 and A1, respectively, for GNN4CD RC and R-all). The PDFs computed over the entire Italian domain show a good agreement between the estimates and the observational reference (panels a). GNN4CD tends to slightly overestimate precipitation amounts (in the order of a few millimetres per hour) while higher values, from the p $ 99 $ onward and in the tail of the distribution, are more accurately captured. The p $ 99 $ and p $ 99.9 $ computed on the same aggregated data are also close to the GRIPHO reference, especially for the RC model (Table 2). Panels c-d-e display the diurnal cycles of hourly precipitation average, frequency (percentage of rainy hours), and intensity, respectively. Each panel is further divided by seasons and presents the daily evolution at one-hour intervals. Overall, GNN4CD provides a good match with GRIPHO observations, particularly in terms of average precipitation and frequency in the RC model configuration. Both RC and R-all exhibit a larger bias in average JJA precipitation compared to other seasons. In the RC case, this arises predominantly due to too high precipitation intensities. In the R-all case, too frequent precipitation is the main contributor. Nonetheless, in both cases, GNN4CD is well able to capture the evolution of the precipitation with a very good timing of the precipitation peak in the late afternoon (around 17:00-18:00 hours).

Figure 5. Results in the reanalysis to observation downscaling setting and comparison with GRIPHO observations for the testing year $ 2016 $ for the PDF of hourly precipitation [mm/h] with bin size of $ 0.5 $ mm for (a) Italy (I) (b) northern Italy (N) and central-south Italy (C-S); the insets provide a magnified view of the tail of the distribution; (c) average [mm/h], (d) frequency [%] and (e) intensity [mm/h] seasonal diurnal cycles for Italy (I).

Table 2. Extreme percentiles computed for GRIPHO and the GNN4CD RC and R-all model designs for Italy (I), northern Italy (N), and central-south Italy (C-S)

Next, we examine the spatial maps and the corresponding maps of percentage bias (Figures 6 and A2 for GNN4CD RC and R-all, respectively). In the case of average precipitation (panels a), GNN4CD generally leads to overestimation, with the largest bias happening in areas of complex topography. Systematic biases across the estimates have been identified in the Apulia region, as well as along the Tyrrhenian and Adriatic coasts of the Tuscany and Marche regions. These biases are more pronounced for the RC model. In the case of p $ 99 $ and p $ 99.9 $ (panels b-c), GNN4CD still tends to overestimate extreme precipitation in regions characterised by complex topography, although this overestimation is less pronounced relative to the bias observed in average precipitation. Conversely, a clear underestimation is evident in plain and hilly areas. The RC model exhibits overestimation in Apulia and along the Tyrrhenian coastlines, also for the extreme percentiles. The overestimation in regions of complex topography may be likely linked to the well-known issue of gauge undercatch, as the stations used to create the reference observational dataset are primarily located in valleys. Regions of complex topography are instead rarely covered, leading to poorer interpolation and thus affecting the DL model learning (Figure 4b). The temporal coverage of station data is also very diverse and may have negatively influenced the quality of the GRIPHO dataset in the less covered locations (Figure 4c).

Figure 6. Results in the reanalysis to observation downscaling setting and comparison with GRIPHO observations for the testing year $ 2016 $ for (a) average precipitation [mm/h] and percentage bias [%], (b) p $ 99 $ [mm/h] and percentage bias [%], (c) p $ 99.9 $ [mm/h] and percentage bias [%].

Similarly, we analyse the seasonal spatial maps for GRIPHO and the corresponding percentage bias for GNN4CD RC and R-all. The evaluation focuses on average precipitation and extreme percentiles, presented in Figures A3, A4, and A5. Figure A6 shows instead the seasonal PDFs. Both RC and R-all models tend to be wetter in JJA and drier in MAM. The seasonal PDFs suggest that the principal cause might be an overestimation of the occurrence of low-to-moderate events for the RC model and low events for the R-all model. The tails are instead better represented, in line with the p $ 99.9 $ maps showing the lowest biases.

Moreover, we evaluate the PCC for the entire Italian peninsula, the north and central-south areas (Table 3). In five cases out of nine, the correlation coefficients computed for the R-all model are higher than the RC values. Specifically, R-all performs better in all the metrics for the central-south area. Nevertheless, the RC model also reports positive correlation coefficients for all the cases investigated, with the highest values observed for the north of Italy, as expected. However, the systematic biases that we observe in some parts of the central-south area (e.g., the Apulia region) may have a detrimental influence on the aggregated spatial correlations values, where the gap with northern Italy is much more pronounced than in other metrics we assessed.

Table 3. Spatial correlation between the reference GRIPHO maps and the GNN4CD RC and R-all estimated maps for precipitation average, p $ 99 $ and p $ 99.9 $ ; results are shown for Italy (I), northern Italy (N), and central-south Italy (C-S)

GNN4CD is further evaluated in representing the total precipitation for $ 10 $ flood episodes occurring within the GRIPHO time span. All flood events exceed the $ 99 $ th percentile of the precipitation distribution in the affected area, except one case, which remains above the $ 90 $ th percentile, and seven events that also surpass the $ 99.9 $ th percentile. For this specific application, both the RC and R-all models have been retrained by excluding the time steps of the floods from the training set, in order to allow a fair evaluation. Figure 7 displays the comparison with the GRIPHO observational reference for all the events, for both the RC and R-all model designs. The results are promising, as all flood events are captured in terms of both spatial extent and severity. The overestimation of precipitation amounts is consistent with the patterns observed in previous evaluations, except for the second event, where the overestimation is more pronounced. This case requires a more in-depth study of the physical event to understand whether the event is particularly out-of-distribution compared to the cases encountered during training, and will be the subject of further studies.

Figure 7. Total precipitation [mm] for $ 10 $ flood events in Italy. Events $ 1 $ , $ 4 $ , $ 5 $ , $ 6 $ , $ 8 $ , $ 10 $ took place in northern Italy (N), events $ 2 $ , $ 7 $ in central Italy (C), events $ 3 $ , $ 9 $ in southern Italy (S).

Overall, GNN4CD demonstrates consistent performance in the reanalysis to observation downscaling task over the grid points of the Italian peninsula, indicating a degree of spatial transferability across the precipitation distribution. This includes the DL model’s capacity to represent both average and extreme spatial patterns, as well as the characteristics of individual flood events across the entire domain. A comparison of the PDFs relative to northern Italy and central-south Italy further illustrates the spatial transferability of GNN4CD, although some regional differences remain evident (Figures 5b and A1b). Indeed, notable degradation in performance is observed in certain specific regions where systematic biases persist and worsen the aggregate results. In future work, we aim to address these issues by investigating the underlying causes of the observed discrepancies and developing targeted solutions to improve the spatial transferability and overall performance of the model in the affected areas.

6.2. RCM emulation

In the RCM emulation setting, GNN4CD is used as a proper emulator to downscale RCM simulations.

We begin by comparing the PDFs (Figure 8). Panels a and d show the precipitation PDFs estimated by GNN4CD RC and R-all for the historical period. We compare them to the original climate model output and to the $ 10 $ -year subset of the GRIPHO observational dataset. The PDFs estimated by the RC and R-all models tend to yield a higher frequency of low precipitation values while underestimating the tail of the distribution, relative to the observational data. In contrast, RegCM tends to overestimate precipitation in the distribution tail when compared to observations. As expected, GNN4CD estimates are closer to the observed PDFs than to the RegCM distribution, especially for the RC case. The R-all model shows instead a more evident underestimation. Nevertheless, results are generally promising and in line with the reanalysis to observation downscaling case. Panels b-c and e-f present a comparison between the PDFs shown in panels a and d and those generated by the GNN4CD and RegCM models for the mid-century and end-of-century projections. The results indicate that the emulator generally captures the climate change signal exhibited by the RegCM model, reflected in a shift of the precipitation distribution towards more frequent and intense events, when moving to the future time slices. The magnitude of the shift differs slightly between the RC and R-all model, with the latter being closer to the RegCM shift and the former slightly under-representing the change between mid-century and end-of-century. These findings are quite remarkable, considering that GNN4CD was not trained on any precipitation scenario data, yet it shows the surprising ability to capture general trends even beyond the climate regimes where it was trained.

Figure 8. PDF of hourly precipitation [mm/h] with bin size of $ 0.5 $ mm for Italy (I); comparison of GRIPHO $ 10 $ -years (grey) and RegCM historical (black) with (a) historical GNN4CD RC (blue), (b) mid-century RegCM (pink) and GNN4CD RC (orange), (c) end-of-century RegCM (magenta) and GNN4CD RC (dark orange), (d) historical GNN4CD R-all (blue), (e) mid-century RegCM (pink) and GNN4CD R-all (orange), (f) end-of-century RegCM (magenta) and GNN4CD R-all (dark orange); the insets provide a magnified view of the tail of the distribution.

Next, we compare the spatial change projected by GNN4CD RC and R-all, and RegCM. The results when moving from historical to end-of-century are shown in Figure 9, while Figure A7 shows the case of moving from historical to mid-century. Panels a-b show the reference maps in terms of average precipitation for the historical period and the corresponding changes, computed consistently for each of the models. The emulator’s projections show an end-of-century dry change signal in the central regions towards the Tyrrhenian coast and in the Sardinia island, in line with the projections of RegCM. The intensification of the signal from mid-century to end-of-century, both wet and dry, is also confirmed by the projections of GNN4CD. The emulator also agrees in representing the wet change signal over the Alpine chain for the mid-century time slice. The same agreement is observed for the precipitation increase in the Apulia, Basilicata, Veneto and Friuli-Venezia-Giulia regions, mainly evident in the end-of-century estimates. However, there are multiple cases in which the projections of GNN4CD and RegCM disagree. For instance, GNN4CD projects a wetter climate signal in the Padania region (central-northern Italy) for the mid-century, more evident for the R-all model. The same behaviour is observed over the Alpine chain in the end-of-century change. Panels c-d and e-f show the same results (historical reference and corresponding change) in terms of p $ 99 $ and p $ 99.9 $ . When looking at the p $ 99 $ end-of-century change, the disagreement between the projections of GNN4CD and RegCM is even more evident. In this case, both RC and R-all project a wetter climate in almost all the peninsula, whereas the RegCM projections show a dry signal along the Tyrrhenian coast. Notably, the projections of GNN4CD and RegCM for the p $ 99.9 $ change show an overall alignment. The emulator generally agrees with the RegCM change signal sign, although it generally continues to project a wetter climate.

Figure 9. Maps for GNN4CD RC, GNN4CD R-all, and RegCM showing (a) historical average hourly precipitation [mm/h] and (b) end-of-century average change [%]; (c)-(d) the same for p $ 99 $ and (e)-(f) the same for p $ 99.9 $ .

Additionally, we derive box-plots from the spatial maps of precipitation average, p $ 99.9 $ , and percentage of rainy hours over Italy (Figure 10). These quantities are displayed for each time slice for GNN4CD RC and R-all and for RegCM. Similarly, box plots are also displayed for the relative percentage bias maps. The boxes span from the first to the third quartile of the data, with a horizontal line indicating the median value. The whiskers span from the edges of the box to the most extreme data point that falls within 1.5 times the interquartile range from the lower and upper quartiles. The average precipitation map box-plots highlight the opposite behaviour of the RC and R-all models. The former leads to a much higher median than RegCM, whereas the latter is much closer, with a slightly lower value. Accordingly, the corresponding relative bias map box-plot shows a median value very close to zero for the R-all case. For the p $ 99.9 $ map statistics, both RC and R-all models lead to lower median values, smaller in the case of R-all. The median percentage of rainy hours is instead higher in both cases, again with the larger difference produced by the R-all model. For the average precipitation map, the spread in the estimates of the two model designs is comparable. Instead, in the p $ 99.9 $ case the RC model exhibits a significantly larger spread. Finally, in the rainy hours case, the R-all model exhibits the larger spread. When the same statistics are derived for northern Italy and central-south Italy (Figure A8), similar conclusions can be drawn. However, the case of Northern Italy exhibits much greater spreads, which could have worsened the aggregated statistics relative to Italy.

Figure 10. Box-plots for RegCM (red) and GNN4CD RC (green) and R-all (blue), derived for Italy (I) from the spatial maps of (a) average precipitation [mm/h], (b) p $ 99.9 $ [mm/h] and (c) percentage of rainy hours [%]; the lower panels show the box plots for the relative bias maps of the same quantities.

The observed performance in projecting the precipitation change signal across both spatial and temporal dimensions suggests that the proposed emulator is potentially capable of projecting precipitation changes associated with global warming, as well as it may also possess a degree of spatial transferability in the RCM emulation setting.