Projected Changes in the Probability Distributions, Seasonality, and Spatiotemporal Scaling of Daily and Subdaily Extreme Precipitation Simulated by a 50‐Member Ensemble Over Northeastern North America

Global warming is expected to produce modifications in the intensity, as well as in the seasonality and spatiotemporal structure of extreme precipitation. In the present study, the temporal evolution of simulated daily and subdaily precipitation extremes was analyzed to assess how they respond to climate warming over different time horizons. Pooling series from the recent 50‐member Canadian Regional Climate Model v5 Large Ensemble, the probability distributions, date and time of occurrences, and spatiotemporal structure of simulated Annual Maxima (AM) precipitation were analyzed at various spatial scales and for durations between 1 hr and 3 days. In agreement with previous studies, the results underline the large increases in AM precipitation quantiles, especially for the shortest durations and for the more extreme events (i.e., longest return periods), and modifications in their spatiotemporal scaling properties and annual and diurnal cycles. For instance, subdaily AM extremes are expected to occur later in the evening, while, no matter the duration, the extremes are expected to occur over a wider period of the year in future climate. Finally, the analysis of projected AM probability distributions showed that heavy‐tail Generalized Extreme Value (GEV) distributions will most likely be observed in the future climate, with some model grid boxes experiencing a significant increase of GEV shape parameters. These results may have major consequences in terms of the occurrence and possible impacts of the most extreme precipitation events.


Introduction
The impact of climate change on precipitation is a major issue in hydrological and climate science due to the potential impacts these changes may have on natural ecosystems and socio-economic activities (Fischer & Knutti, 2016). Modifications in the probability distributions of precipitation are expected to occur at regional and local scales as a result of global warming partly because of the increased moisture holding capacity of a warmer atmosphere (Trenberth, 2011;Westra et al., 2014). According to these theoretical considerations, more intense and frequent precipitation extremes are also expected, especially at subdaily and subhourly temporal scales Berg et al., 2013;Wasko et al., 2015). Additionally, changes in the seasonality and spatial distribution of extreme rainfall events may be expected based on both thermodynamic considerations and changes in small and large-scale dynamics (Dhakal et al., 2015;Pfahl et al., 2017;Skeeter et al., 2018;Touma et al., 2018).
Using historical observational records and simulations in historical climate, numerous studies supported these theoretical arguments, showing increases in the frequency and/or the intensification of daily and multidaily precipitation extremes over the past decades at global and regional scales (e.g., Barbero et al., 2017;Fischer & Knutti, 2016;Kendon et al., 2018). Some evidence of changes in the duration and spatial characteristics of extreme rainfall events, such as the spatial extent and autocorrelation structure, have also been reported (e.g., Wasko et al., 2016, for observed precipitation series; Li et al., 2015, andGuinard et al., 2015, for Regional Climate Model (RCM) simulations under future conditions; and Prein et al., 2017, for convection permitting RCM projections). This has important implications for security issues, resource management, and infrastructure design since watershed and ecosystem response to extreme rainfall events highly depends on their spatial and temporal features (Mallakpour & Villarini, 2017). However, while there is a significant and growing literature on extreme precipitation intensity, fewer studies have jointly examined the future evolution of extreme precipitation characteristics such as duration, seasonality, timing, and spatial features (Wasko et al., 2016). As a result, the possible modifications of the spatiotemporal structure of extreme rainfall and the mechanisms driving these changes in a warming climate are not completely understood (Dwyer & O'Gorman, 2017;Mallakpour & Villarini, 2017;Touma et al., 2018).
Although there is a broad consensus that climate change particularly affects extremes at short duration and for long return periods (Kharin et al., 2018;Westra et al., 2014), modifications of subdaily rainfall extremes are generally difficult to assess due to their high temporal and spatial variability (e.g., Barbero et al., 2017;Kendon et al., 2018) and the paucity of high-quality recorded rainfall series (Westra et al., 2014).
Sparse gauge networks, short records, and measurement errors limit the possibility of assessing the spatiotemporal characteristics of extreme precipitation from observed rainfall series (Tapiador et al., 2017). Moreover, uncertainties associated with dynamical models (e.g., structural modeling and tuning uncertainties, including those related to the convective parametrization; Kendon et al., 2017;Tebaldi & Knutti, 2007) radiative forcing (e.g., scenario uncertainties; Hawkins & Sutton, 2009), and the possible biases arising from the relatively coarse resolution of state-of-the-art RCMs (e.g., Kendon et al., 2014;Prein et al., 2015) should be accounted for when inferring rainfall extreme characteristics from climate model simulations. Finally, the natural climate variability associated with the chaotic and nonlinear nature of the climate system may hide the temporal changes in precipitation extreme statistics due to climate warming, especially at fine spatiotemporal scales and for the most extreme events (Hawkins, 2011). However, the impacts of all these uncertainties on the spatiotemporal structure of rainfall extremes are difficult to evaluate, partly because they emerge at different time horizons and spatiotemporal scales (Fatichi et al., 2016, and references therein). For instance, whereas model and scenario uncertainties play an important role for most climate variables for projections at the global scale and/or for centennial lead times, internal variability is expected to be the major source of uncertainty for precipitation, especially for the most extreme events and short time horizons (e.g., one decade or two; Hawkins & Sutton, 2009;Hingray & Said, 2014).
One way to assess the impact of internal variability on rainfall extreme estimation is done by using initial-condition ensembles from a single climate model and a given radiative forcing scenario (e.g., Flato et al., 2013;Sanderson et al., 2018). Ensemble members are thus created by applying minor perturbations to the initial state of the model simulation so that the different climate trajectories are surrogate representations of the natural climate variability (Deser et al., 2012). Initial-condition ensembles have been therefore used for detecting and discriminating climate change signals from natural variability in global and regional climate extremes (e.g., Fyfe et al., 2017;Martel et al., 2018). In this study, a recent 50-member initial-condition large ensemble produced using the Canadian Regional Climate Model Version 5 (CRCM5 Martynov et al., 2013;Separovic et al., 2013) for the period 1950-2100 , is used to investigate the multiscale characteristics of extreme precipitation changes in a warming climate. The analysis focuses on the assessment of temporal changes in the extreme rainfall statistical properties, such as depth quantiles and date or time of extreme occurrence, at various spatial scales and for durations ranging from 1 hr to 3 days. In this context, the salient question is to examine the effects of sampling errors on the spatiotemporal structure of simulated extremes (Li et al., 2019;Thompson et al., 2017). By pooling series from various CRCM5 members over short periods of different length, the objective is to investigate how the large ensemble intermember variability can be used to improve the estimation of Annual Maxima (AM) precipitation statistics.
The paper is organized as follows. Section 2 introduces the CRCM5 large ensemble (CRCM5-LE), and section 3 describes the extraction AM precipitation series at various spatiotemporal scales. Sections 3.1 and 3.2 outline the methodology used for pooling AM series from the various CRCM5-LE members and to assess the impacts of climate change on rainfall extremes. The statistical characterization of AM probability distributions, spatiotemporal structure, and the annual and diurnal cycles of AM occurrences are described in section 3.3 to 3.5. Results are presented and discussed in section 4, while section 5 provides a summary and presents the perspective for future works. Table 1 provides the definitions of recurrent acronyms used in the text.

Data and Study Area
The study is based on a 50-member ensemble simulated at the 0.11°resolution (≈12 km) by the CRCM5 for the period 1950-2100 over northeastern North America (Figure 1; Leduc et al., 2019).
The initial-condition CRCM5-LE was generated using the CanESM2-LE (Sigmond & Fyfe, 2016;Fyfe et al., 2017). The 50 independent CanESM2-LE simulations were simulated using the Canadian Fourth generation Atmospheric Global Climate Model (CanESM2; von; Arora et al., 2011;von Salzen et al., 2013) by applying random perturbations to the initial state of cloud-overlap parameters while all other simulation settings (e.g., forcing scenario and model parameters) remained unchanged (Fyfe et al., 2017). For each member, the radiative forcing was prescribed as the observed concentrations of historical greenhouse gases until 2005 and through the Representative Concentration Pathway 8.5 (RCP8.5) scenario thereafter (Meinshausen et al., 2011).
A spin up period of four years was discarded for each CRCM5-LE member series, resulting in 146 years of hourly precipitation available at 280×280 grid points over the 1954-2099 period. The reader can refer to Leduc et al. (2019) for additional details on the experimental setup and validation of CRCM5-LE monthly temperature and precipitation statistics. Innocenti et al. (2019) evaluated the Canadian RCM performances in representing the extreme rainfall characteristics (including diurnal and annual cycles) at several spatiotemporal scales by comparing CRCM5 simulations driven by the ERA-Interim reanalysis (Dee et al., 2011) to other simulated and observational gridded dataset. Their results showed an overall good agreement between ERA-Interim CRCM5 simulations and station quantiles at hourly time scale, while the model generally overestimated rainfall extreme quantiles at daily and longer durations. Interestingly ERA-Interim CRCM5 simulations also showed comparable performance to a convection-permitting RCM (a WRF v3.4.1 simulation, Liu et al., 2017) and two observation-based datasets (the CMORPH, Xie et al., 2017, andthe multisource MSWEP dataset, Beck et al., 2017) for representing observed precipitation extreme annual and diurnal cycles. Finally, Innocenti et al. (2019) highlighted that minor differences were found between the ERA-Interim CRCM5 simulations and CRCM5-LE diurnal cycles, while the large ensemble showed lower probabilities of AM occurrences during fall. CRCM5-LE subdaily AM quantiles were also smaller than the corresponding ERA-Interim CRCM5 estimates in the south of the domain and larger in the northeast regions.

Methods
For each CRCM5-LE member, a moving window was applied to aggregate grid box hourly precipitation series at durations d= 1, 2, 3, 4, 6, 12, 18, 24, 36, 48, 60, and 72 hr. AM precipitation values were then extracted for each model grid box and spatiotemporal scales (r,d). As in Innocenti et al. (2019), grid box hourly precipitation series were also aggregated at six spatial scales using a fixed window in space (r=12,24,…,72 km). Grids at scales r>12 km were defined by considering nonoverlapping grid boxes starting at the southwest corner of the native CRCM5 grid and moving toward the northeast corner of the domain. Grid boxes associated with the ocean and water bodies were removed, and series at coarser spatial scales were discarded if they included less than 75% native land grid boxes. Finally, AM precipitation series aggregated at each spatiotemporal scale (r,d) were extracted.
The date and time (starting hour at local time, UTC-5) of occurrence of the AM were also extracted for each spatiotemporal scale, year, and grid box.
For the analysis of AM statistics, each native CRCM5 grid box was associated with the overlapping grid box at coarser spatial scales. Also, to reduce the computational time, the CRCM5 grid was undersampled with a ratio of 1/4, resulting in the selection of ≈11,800 grid boxes randomly distributed over the CRCM5 North American domain. This strategy also reduced the number of neighboring grid points and therefore the spatial autocorrelation of computed statistics.

Pooling CRCM5-LE Series
The probability of detecting changes in the time series of a given climatic variable X depends on the signal-tonoise ratio (Hawkins & Sutton, 2012), which is a function of the magnitude of the signal, the series length, and the variability of X. Since time series of short-duration AM rainfall usually display large variability, the signal-to-noise is generally low at local spatial scales but can be enhanced by pooling AM series from several spatial and/or temporal scales (e.g., Innocenti et al., 2017;Li et al., 2019), from various spatial locations (e.g., Shephard et al., 2014), and/or from several simulated series (e.g., Martel et al., 2018).
CRCM5-LE members represent 50 equiprobable climate realizations over the 1954-2099 period providing 50 independent AM, and their corresponding date and time of occurrence, for each year, spatiotemporal scale, and grid box. In order to reduce sampling uncertainty on estimated AM characteristics, the 50 CRCM5-LE series were pooled over short sub-periods of 1, 3, and 7 years, hereinafter referred to as 1SP, 3SP, and 7SP. Figures 2a-2b show the steps of this procedure in a schematic diagram for the 3SP case.

Trend Analysis for AM Statistic Series
Consider a time series Y t , where Y is a generic statistic computed from AM series (e.g., Y may represent the empirical quantile of AM probability distribution; see sections 3.3-3.5 for a description of the statistics defined in this study). The subscript t indicates the median year of a given 1SP, 3SP, or 7SP. For instance, t=1955 for the first 3SP, as shown by the example in Figure 2c.
The existence of significant monotonic trends for Y t was assessed for each grid box using the rank-based Mann-Kendall (MK) test with the variance correction proposed by Yue and Wang (2004) to account for the serial correlations in Y t . To take into consideration the spatial dependence between CRCM5 grid box series, the field significance of local (grid box) MK tests was then evaluated through a False Discovery Rate (FDR) approach (Wilks, 2006). The field significance at the global (field) significance level α glo = 0.1 was computed using the Benjamini and Hochberg (1995) procedure.
The magnitude of the significant temporal Y t trends was then estimated using the following model: where the subscript Y for α Y and β Y refers to the considered AM statistic.
The log-linear model in equation (1) adequately fits the observed trends (see, e.g., Figures 4a and 4b) and was chosen as it allows a straightforward interpretation of β Y slopes. For small values of β Y , we have e β Y ≈1 þ β Y . Hence, the slope of equation (1) represents the relative change in the Y t statistics expected for a unit increase in t, namely, the mean annual relative change of Y t . Equivalently, b y =10 3 β Y approximates the mean decadal percent increase of Y t and will be used as reference measure for the trend magnitude.
For statistics that can assume negative values (e.g., the shape GEV parameter; see section 3.4), classical linear regression models were used to estimate the absolute magnitude of expected Y t trends. In these cases, the regression slope β Y represents the absolute variation of Y t expected for a unit increase in t, that is, the expected annual variation of the statistic.

Statistical Characterization of AM Rainfall
Let x t ¼ x 1;t ; x 2;t ; …; x i;t ; …; x n;t À Á be the sample of n AM pooled from the 50 CRCM5-LE members for the subperiod t, a given grid box s, and the spatiotemporal scale (r,d). To simplify notation, grid box and spatiotemporal scale subscripts are omitted, unless otherwise stated. Empirical AM quantiles, x t,q , were estimated for each x t and various return periods q (years) based on the Hazen plotting position (Cunnane, 1978).
The spatial and temporal structure of the sampled extremes was investigated considering the following loglog linear model : lnfx t;q ðr; dÞg ¼ α r;t þ ðh 0;t þ h 1;t rÞlnðdÞ; (2) where x t,q (r,d) represents q-year AM quantiles extracted at the spatiotemporal resolutions (r,d) for the subperiod t and the intercepts α r,t represent estimates of lnfx t;q ðr; 1hÞg for each r=12,24,…,72 km. Note that, for simplicity, the index q for quantile orders have been omitted for α r,t , h 0,t , and h 1,t parameters. Equation (2) is a concise expression that includes both the temporal scaling of AM quantiles at the point scale (Casas-Castillo et al., 2018;Menabde et al., 1999; through the h 0,t coefficient, hereinafter the temporal scaling parameter) and the sensitivity of AM temporal scaling to changes in the spatial scale (through h 1,t , hereinafter the spatiotemporal scaling parameter).
Previous investigations on the CRCM5-LE and other gridded datasets justify the choice of equation (2)   and showed how it extends the classical simple scaling formulations (e.g., Burlando & Rosso, 1996;Gupta & Waymire, 1990) of Depth-Duration-Frequency (DDF) and Intensity-Duration-Frequency (IDF) curves, widely used in hydrological and engineering applications (Koutsoyiannis et al., 1998). Innocenti et al. (2019) also showed that h 1,t corresponds to the change of AM quantile Areal-Reduction-Factors (ARFs; Svensson & Jones, 2010, and references therein) with (r,d) for each subperiod t, that is, where ARF t,q (r,d)=x t,q (r,d)/x t,q (0,d) and ARF t,q (r,1h)=x t,q (r,1h)/x t,q (0,1h) respectively represent the q-year quantile ARFs for durations d and 1h at the spatial scale r and for the subperiod t.
As a result, the future evolution of the extreme precipitation spatiotemporal structure simulated by the CRCM5 can be assessed through the projected h 0,t and h 1,t values. For instance, previous studies stressed that temporal scaling parameters might change with time, questioning the adequacy of stationary frameworks classically used to estimate DDF and IDF slopes (e.g., Casas-Castillo et al., 2018 for observed rainfall extremes; Cannon & Innocenti, 2019 for daily and subdaily simulated AM). Finally, note that some authors underlined the presence of two temporal scaling regimes for short (e.g., hourly and subhourly) and long (e.g., daily and longer) durations (e.g., Ceresetti, 2011;Eggert et al., 2015). This is explained by the fact that different weather regimes and processes drive precipitation extremes at different spatiotemporal scales (Eggert et al., 2015;Innocenti et al., 2017). Hence, the scaling properties of CRCM5-LE AM were separately investigated for Short Duration (SD, 1 hr ≤d< 6 hr) and Long Durations (LD, 6 hr ≤d≤ 72 hr; see Innocenti et al., 2019, for details).

SS-GEV Model
AM probability distributions were also assessed through the GEV distribution (Coles, 2001). Considering X t to be the random variable for AM in the x t sample, the GEV cumulative distribution function for the subperiod t and a generic spatiotemporal scale can be written as where −∞<x<+∞. The parameters μ t ∈R, σ t >0, and ξ t , respectively, represent the location, scale, and shape parameters of the distribution for the subperiod t. Since the shape parameter characterizes the distribution tail, an accurate estimation of ξ t is critical to adequately assess extreme precipitation quantiles.
Uncertainties on ξ t are large for short records (Koutsoyiannis, 2004a(Koutsoyiannis, , 2004bPapalexiou et al., 2013), and various strategies to increase the sample size have been proposed to improve GEV parameter estimation (e.g., regional analysis, Hosking & Wallis, 1997; or approaches combining AM series at different temporal scales, Blanchet et al., 2016). In particular, temporal scaling relationships have been shown to increase the robustness of GEV parameter estimation for station series (Innocenti et al., 2017). Extending typical simple scaling GEV formulations (e.g. Blanchet et al., 2016;Panthou et al., 2014), the following expressions were thus considered for the GEV parameters at each spatiotemporal scale (r,d): where H t,r (r=12,24,…,72 km) are generally referred to as simple scaling exponents and μ Ã t;r , σ Ã t;r , and ξ Ã t;r are the GEV parameters for the spatial scale r, the subperiod t, and the d * =1h being used as reference duration. Distribution in equation (6) is hereinafter referred to as SS-GEV model.
Note that, according to equation (2), the SS-GEV simple scaling exponents should verify H t,r = h 0,t +h 1,t r. However, while equation (2) may result in spatiotemporal scaling parameters that depend on the quantile return period q, the proposed SS-GEV model is based on unique h 0,t and h 1,t values estimated for each grid box and subperiod t (i.e., SS-GEV assumes scaling parameters, and consequently IDF and DDF slopes, to be independent of the return period).
Several methods can be considered for the estimation of the parameters ðH t;r ; μ Ã t;r ; σ Ã t;r ; ξ Ã t;r Þ of the SS-GEV model over each subperiod t. The Nelder-Mead numerical approximation of the Maximum Likelihood (ML) estimate was used, considering a Generalized Linear Model formulation of the SS-GEV. To this end, the duration d was introduced as a model covariate in SS-GEV probability distribution function formulations (Coles, 2001), as in temporal SS expressions used in previous studies (e.g., Mélèse et al., 2018, and references therein). Likelihood-Ratio (LR) tests were also used to test the null hypothesis H 0 : ξ t = 0 for each subperiod t (Gumbel vs. GEV distributions).
SS-GEV distribution parameters were estimated independently at each grid box and spatial scale r, assuming that AM extracted at different d are independent. This last assumption is not fully respected since AM occurring during the same year at different spatiotemporal scales may be associated with the same precipitation event or system (Innocenti et al., 2017). However, previous studies noticed that the GEV estimation is robust to model misspecification due to the dependence of AM at various durations, with little improvement when this dependence is explicitly modeled (e.g., Blanchet et al., 2016;Mélèse et al., 2018). Also, to reduce computation time and improve convergence for Generalized Linear Model ML optimizations, only durations d=1, 2, 3, 6, 24, and 72 hr were considered for GEV and SS-GEV estimation.

Characterization of AM Annual and Diurnal Cycles
Descriptive statistics for annual and diurnal cycles of AM occurrences were assessed at the native CRCM5-LE spatial resolution (r=12 km) within the circular statistics framework (Pewsey, 2004), specifically developed for the analysis of periodic random variables. Within this framework, the date and time of AM occurrence at each spatiotemporal scale (r,d) are represented by angular variables (in radians) as (Berens, 2009) where z i,t and w i,t respectively represent the ordinal day and the hour at which the ith AM of the x t sample begins. Subscript for duration d has been omitted to simplify notation. A visual representation and an example for these variables are shown in Figure 3. Equation (7) can then be used to measure the sample mean, variance, and shape of the AM annual and diurnal cycles. For instance, the mean date of AM occurrence can be determined for the subperiod t by averaging the cosine and sine components of the angular vectors (δ 1,t ,…,δ i,t ,…,δ n,t ): The mean date of AM occurrence can then be computed as the mean angle,¯δ t , using the four quadrant inverse tangent function (Berens, 2009) and inverting equation (7) to estimate the corresponding¯z t (days).
The mean time of occurrence of AM can be estimated in a similar way considering mean vectors¯v t and corresponding angles¯η t , as shown in the examples of Figure 3b.
For nonperfectly symmetrical (e.g., bimodal) variables, the length of the mean resultant vectors can also be interpreted as a measure of the circular dispersion of δ i,t and η i,t values around their mean values,¯δ t and¯η t , for each subperiod t. The closer ‖¯u t ‖ (‖¯v t ‖) value is to one, the less dispersed are the AM occurrences around the mean date (time) of occurrence, while ‖¯u t ‖ ¼ 0 and ‖¯v t ‖ ¼ 0 for δ i,t or η i,t uniformly distributed on the unit circle (i.e., for maximum variance). Finally, a similar approach can be used to define measures of circular symmetry, skewness, and higher-order moments of circular probability distributions (Pewsey, 2004). The reader can refer to Berens (2009) for further details.

Results
The projected temporal evolution of extreme precipitation properties is presented by showing the expected decadal percent variation of AM probability distributions (empirical quantiles and GEV parameters) and annual and diurnal cycles (mean date and time of occurrence and corresponding variability). The expected trends in extreme precipitation spatiotemporal structure (scaling and SS-GEV parameters) are thenconsidered.
Unless otherwise stated, only results for 3SP are presented since similar conclusions were reached for the 1SP and 7SP (the supporting information provides some examples). Figure 4 summarizes the main results for projected changes in empirical extreme precipitation distributions. Figures 4a and 4b show that hourly and daily AM quantiles increase for virtually all grid boxes. MK significant trends were found for more than 98% of grid boxes at the field significance level α glo =0.1 (not shown). Also, while significant increases appear during the 1990s for both hourly and daily 100-year quantiles (Figures 4a and 4b), trends were relatively stronger for subdaily AM than for daily and longer durations ( Figure 4c).

Projected Changes in AM Probability Distributions
The decadal percent variations b xq (equation (1)) of empirical AM quantiles were particularly large for d≤3 hr, namely, larger than 3.3% per decade for more than half of CRCM5 grid boxes at the hourly scale ( Figure 4c). Estimated trends also displayed the largest spatial variability for d=1 hr, as shown by the 10th to 90th and 25th to 75th quantile ranges in Figure 4c. Finally, similar results were observed for shorter return periods but weaker trends were observed for q=2 years and q=25 years compared to 100-year AM (Figure 4d), meaning that more extreme AM are expected to experience larger relative increases in future climate. Figure 5 shows that the largest decadal percent increases were generally obtained in the northeastern portion of the domain, while a decreasing gradient was observed along the NE to SW axes, especially at the hourly time scale. These findings are consistent with the prominent changes of mean and extreme temperatures

10.1029/2019JD031210
Journal of Geophysical Research: Atmospheres and precipitation expected in Canada at high latitudes which are partly enhanced by regional climate feedbacks (e.g., the snow/ice albedo and cloud feedbacks which are particularly relevant at the highest latitudes; Zhang et al., 2019). For all durations and return periods, estimated trends over the Great Lakes region and the southern Atlantic coast were generally small compared to south inland regions. Differences in trends between coastal and inland grid boxes were locally observed in some areas over the northeastern coasts for daily and longer durations (e.g., along the St. Lawrence Gulf coast; Figure 5, second and third columns). Due to higher sampling uncertainties, larger spatial variability was observed for b xt;q for 100-year quantiles compared to shorter return periods, resulting in noisier maps and less obvious spatial structures, especially for daily and longer-duration estimates.
GEV distribution parameters estimated through equations (4) and (5) are presented in Figure 6 for each 3SP and 7SP and for d = 1, 3, 24, and 72 hr (to improve the readability of Figure 6a, the 72-hr curve is not shown). The comparison between 3SP and 7SP aims at evaluating the uncertainty on estimated GEV parameters and the impact of CRCM5-LE member pooling over subperiods of various lengths. The LR-tests assessing the significance of the ξ t parameters showed that the Gumbel distribution (ξ t =0) adequately represents the 3SP distributions of AM for most of grid boxes (Figure 6a), especially for d≥24h. Also, the fraction f ξ t <0 of grid boxes with significant and negative shape parameters (bounded tail distributions) was negligible (f ξ t <0 <6ñ10 −3 ) for all considered durations and pooling subperiods (not shown). Moreover, the fraction of model grid boxes rejecting the LR hypothesis due to significantly nonzero ξ t clearly decreases with increasing t for all d>1h (Figure 6a). This result was confirmed for 7SP, with f ξ t ¼0 values being lower in the second half of the 21 st century than during the 20th century (e.g., f ξ t ¼0 was greater than 0.24 before 2000 and f ξ t ¼0 ≈0:2 after 2050 for d=3 hr AM). This suggests that the large proportions of grid boxes with Gumbel distribution found for 3SP were mainly due to the large statistical uncertainty of GEV shape estimation. This uncertainty can in principle be reduced using larger samples (e.g., using 7SP). Medians over the CRCM5 grid boxes of the estimated GEV parameter values are shown in Figures 6b-6d for 3SP and 7SP and several durations. For shape parameters, Figure 6b, only considers grid boxes with ξ t >0. The 3SP systematically displayed higher shape values than 7SP, with positive median ξ t values distributed around ξ t ≈ 0.17. This result is consistent with typical GEV shape parameter values estimated from long records (e.g., more than 100 years) and reported in the literature (e.g., Ragulina & Reitan, 2017, and references therein). For instance Koutsoyiannis (2004a) and Koutsoyiannis (2004b) reported positive ξ values lower than ≈0.23 for various locations around the world. Finally, note that for both 3SP and 7SP, median ξ t values remain almost constant over time for all durations. Owing to the high uncertainty on estimated shape parameters, the fraction of grid boxes with statistically significant trends for ξ t was generally negligible (e.g., lower than 1% for most durations, not shown). Conversely, most of the grid boxes displayed statistically significant trends when considering statistics related to the shape characteristics of the AM distribution that are less impacted by the sample size. In particular, for all durations, more than 75% of grid boxes displayed significant positive trends for the normalized dispersion coefficient σ t /μ t ( Figure S1 in the supporting information), suggesting that the increase in more extreme quantiles could be relatively larger than for less extreme quantiles.
A substantial shift of AM distributions toward higher values is observed for all durations due to increase in median location parameters with t, while the increase of grid box median scale parameters corresponds to higher variability of AM precipitation in future climate (Figures 6c and 6d). According to the MK tests, μ t and σ t displayed significant increases for almost all grid boxes and durations ( Figure S1 in the supporting information). For both parameters, decadal percent variations showed low local variability, with spatial patterns similar to those of the 2-year AM quantiles (e.g., Figure S2 in the supporting information).
Finally, note that 3SP and 7SP estimates of μ t and σ t did not present any appreciable difference (Figures 6c  and 6d), indicating that sample size has no impact on estimated values, contrary to shape parameters.

. Projected Changes in Annual and Diurnal Cycles
The temporal evolution of the basic statistics for dates of occurrence of hourly and daily AM is shown in Figure 7. The distribution over the CRCM5 grid boxes of the mean date of occurrence is represented for each 3SP considering the spatial angular medians of¯δ t values, namely, the diameter of the unit circle that divides the¯δ t values into two groups containing the same number of grid boxes. The procedure is then repeated for each of these two groups to estimate pseudo-quartiles of the¯δ t distribution, that is, angular measures that can be considered equivalent to 25th and 75th percentiles of the distribution over CRCM5 grid boxes of the mean date of occurrence of the AM. Figure 7a shows that hourly AM are generally concentrated in July, with the pseudo-quartiles of the mean date of occurrence between the end of June and the beginning of August for all 3SP. At the daily scale, the mean dates of occurrence are more variable over summer months (Figure 7c). This result is expected at midlatitudes, where convection-driven short-duration extremes more likely occur during the warmest season of the year, while daily and longer-duration AM are more evenly distributed over the year, as they are more influenced by large-scale dynamics and less constrained by specific thermodynamical conditions than subdaily AM. As expected, clear spatial distributions also emerged for¯δ t in each subperiod t, with the obvious influence of water bodies and oceans and marked latitudinal gradients (see, e.g., Figure S3 of the supporting information). The median¯δ t remained almost unchanged over time for both d=1 hr and d=24 hr. However, at the end of the simulation period hourly and daily AM more frequently occur during nonsummer months for some grid boxes, as the¯δ t pseudo-quartile intervals become wider and asymmetrical with time, especially

10.1029/2019JD031210
Journal of Geophysical Research: Atmospheres for d=24 hr. In fact, at the end of the 21st century, daily and longer-duration AM more frequently occur in winter and fall in southern coastal areas and close to the Great Lakes and the Atlantic coast, while in southern continental areas the annual peaks of daily AM are projected to occur more frequently during spring in future climate (see, e.g., Figure S3 in the supporting information). These results agree with the projected changes in daily extreme (Rx1day) seasonality estimated by Marelle et al. (2018) and Bronnimann et al. (2018) in the Northern Hemisphere at middle to high latitudes. In particular, Marelle et al. (2018) found that, in northeastern North America, the daily AM peak of occurrence could shift from midsummer to late spring or fall and early winter by the end of the 21st century. Furthermore, according to Marelle et al. (2018), the seasonal changes projected for temperatures drive the extreme seasonality changes at high latitudes (e.g., for the Arctic region), while the contribution of dynamical mechanisms these changes is larger at lower latitudes, as also discussed in Pfahl et al. (2017).
Hourly extremes showed stronger seasonality than daily AM and annual cycles were more peaked around their means. In fact, larger ‖¯u t ‖ values were observed for d=1 hr compared to d=24 hr (Figures 7b-7d), which indicates that AM occurrence dates have less dispersed distributions for any given subperiod t. Also, the local temporal variability of mean dates decreased with time for both hourly and daily AM, with significant ‖¯u t ‖ trends detected at more than 81% of the grid boxes for all durations and pooling strategies (not shown). This means that the annual cycles of AM occurrence are expected to be more dispersed around their mean values in a future warmer climate and for most grid boxes AM are expected to occur over a wider period of the year. Preliminary analyses on the circular skewness and kurtosis of δ i,t values further supported this conclusion, as they showed a global tendency to more symmetric and less peaky grid box annual cycles, especially for sub-daily durations (not shown).
Large decadal percent variations were estimated for ‖¯u t ‖ and were generally negative, showing an increase in the local temporal variability of the mean dates. However, these trends of the annual cycle dispersion displayed different spatial patterns for different durations. For daily and longer AM, significant increases in local temporal variability of the mean dates (i.e., b ‖¯ut‖ <0) were estimated in northern regions west of the Great Lakes and for Quebec and Labrador, while both increasing and decreasing trends (i.e., increasing or decreasing variability of the mean dates) were found for ‖¯u t ‖ in the south and in the eastern regions of the domain (e.g., Figure S4 in the supporting information). For most grid boxes, however, the occurrence times, η i,t , displayed large variability within each subperiod, making it difficult to identify a clear peak in the diurnal cycle of the AM. This resulted in relatively small ‖¯v t ‖ values (smaller than 0.5 for almost all grid boxes) that further decrease with t for the shortest durations (e.g., for d≤3 hr), especially in southeast areas (not shown). Hence, the times of occurrence Figure 10. Spatial distribution of the decadal percent variation over 3SP estimated at the native spatial resolution (r= 12 km) for the temporal (h 0,t first row; decreasing significant trends) spatiotemporal (h 1,t , second row; increasing significant trends) scaling parameters for 2-year quantiles for (a) SD and (b) LD. Smaller points indicate grid boxes with no statistically significant trend. The percentage f H1 of grid boxes with significant trend are indicated on each panel.

10.1029/2019JD031210
Journal of Geophysical Research: Atmospheres of subdaily AM were relatively evenly distributed over the 24-hr circle within each subperiod t and¯η t is likely a poor approximation of the diurnal cycle mode, especially for the future climate.

Projected Changes in the Spatiotemporal Structure of AM Precipitation
Distributions of the estimated SD and LD scaling parameters (equation (2)) over CRCM5 grid boxes are shown in Figure 9 for each 3SP. Confirming previous investigations on the spatiotemporal scaling of gridded datasets , SD displays stronger scaling regimes (i.e., higher h 0 and h 1 values) than LD. This implies that changes in AM quantiles related to changes in the spatiotemporal scale (r,d) are more pronounced for SD because short-duration extremes are generated by weather systems that are more spatially and temporally heterogeneous than LD.
Moreover, the temporal scaling parameters decreased with time for both SD and LD (Figure 9, first row). This indicates that the changes of point-scale AM quantiles across durations are expected to be larger in future climate than those observed in past periods. Similar results were also obtained in previous studies (e.g., and references therein , which reported an increase in time of SS exponents for rainfall intensity distributions or, equivalently, decreasingly negative slopes of IDF curves in future climate. Note, in fact, that point scale IDF slopes can be expressed as h int 0 ¼ h 0 −1, where h int 0 is the temporal scaling parameters for AM precipitation intensity. The projected decreases of h 0 in time thus correspond to a relatively stronger intensification of short-duration extremes compared to long duration extremes. Increases of the spatiotemporal scaling parameter, h 1,t , were observed for the shortest return period, while trends are less obvious for q≥25 years, especially for SD (Figure 9, second row). Hence, the heterogeneity of the spatial characteristics of precipitation events producing AM is generally expected to increase in future decades, as higher h 1 values corresponds to larger variations of the ARFs (see equation (3) for a specific change in the spatiotemporal scale (r,d). Also, h 1,t estimates were more spatially heterogeneous (not show) and displayed larger short-term temporal variability than h 0,t for q≥25 years (see, e.g., median curves in Figure 9), suggesting that large uncertainty affects the spatiotemporal scaling estimation, as previously found by Innocenti et al. (under review, 2019).
The analysis of the spatial distribution of scaling parameter trends revealed that large decadal percent variations were observed for both h 0,t and h 1,t , especially in northeastern areas, as shown in Figure 10. For h 0 , significant negative trends were observed for more than 98% of grid boxes, with the largest relative changes found in northern North Atlantic regions (Figure 10, first row). Significant trends were found for a smaller fraction of grid boxes (between 50% and 60% of grid boxes for 3SP and q=2 years) for h 1 and some differences emerged between northern (positive) and southern (negative) b h1;t estimates ( Figure 10, second row). This result is in agreement with previous studies that project a change in the size of weather systems producing extreme precipitation in future climate (e.g., Prein et al., 2017). For a given spatial scale r, in fact, higher h 1 values correspond to higher variations of ARF values when changing the duration. Accordingly, field significance was found for the h 1,t decadal percent increases over the whole CRCM5 domain, while b h1;t values were field significant only when the False Discovery Rate approach was restricted to southern regions. For instance, for grid boxes located south of the 40°parallel, h 1,t temporal trends were generally negative, especially for the longest return periods, which indicates a decreasing sensitivity of ARFs to changes in (r,d) and a more homogeneous spatiotemporal structure of most extreme short-duration events.
Note that, although the decadal percent variations for longer return periods were generally larger and characterized by similar spatial distributions than those presented in Figure 10, the fraction f H1 of field significant trends for q=100 year was generally low (e.g., lower than 1%) for h 0 SD estimates and for h 1 (for both SD and LD). Moreover, f H1 decreased when increasing the number of years pooled within each subperiod (not shown). This underlines that changes in the spatiotemporal structure of projected AM distributions are more difficult to assess for more extreme quantiles, especially for the shortest durations, while the use of longer pooling periods (i.e., 3SP and 7SP) limit the possibility of detecting trends for h 1 as shorter series are considered for the MK tests.
Estimated SS-GEV parameters are presented in Figure 11 for each 3SP considering the reference duration d * =1 hr for SD and d * =24 hr for LD at the native CRCM5 spatial resolution (r=12 km). This figure confirms the statistically significant decrease of the temporal scaling parameter (Figure 11a) and the increase of the scale and location SS-GEV parameters (Figures 11c and 11d) with t that were suggested by results in section 4.1 and those for the h 0 and h 1 parameters. Results for coarser spatial scale (i.e., r>12 km) were also consistent with these estimates (not shown). Figure 11 highlights three important results. First, the difference between short and long duration scaling regimes is substantial, as SD scaling parameter values are much larger than the corresponding LD values (Figure 11a).
Second, the SS-GEV distributions were mostly heavy tailed, with 3SP shape parameters being significantly positive for most of grid boxes for all t and for both SD and LD. The proportions of 3SP Gumbel distributions (ξ Ã t;r ¼ 0) were considerably lower than the corresponding proportions estimated for nonscaling GEV models (e.g, f ξ Ã t;r ¼0 ≤0:21 for SD and f ξ Ã t;r ¼0 ≤0:31 for LD), and f ξ Ã t;r ¼0 decreased with time for LD, while it remained almost unchanged for SD (see Figure S5 in the supporting information). Moreover, only few grid boxes (e.g., less than 1%) displayed significant trend for SD shape parameters values for all pooling strategies, while 50.5% of grid boxes had significant increasing trends for LD (Figure 11b for 3SP, while 50.7% and 42.0% of significantly increasing trends were detected for 1SP and 7SP, respectively). These results are particularly important as they imply that crucial modifications can be expected for the characteristics of AM probability distribution tails, which have huge impacts on the statistical properties of more extreme AM quantiles.
Third, SD shape parameter values were larger than corresponding LD values for each 3SP for more than 61% of the grid boxes, as also shown by the SD and LD ξ Ã t;r distributions ( Figure 11b; note that only significant ξ Ã t;r >0 are considered in this figure). Accordingly, SD extremes displayed heavier tailed distributions than longer-duration AM, which is consistent with the hypothesis that links the shape parameter values to the characteristic precipitation systems generating AM at different temporal scales (Ragulina & Reitan, 2017).

Journal of Geophysical Research: Atmospheres
Finally, note that H r,t values estimated for r=12,24,…,72km verified, for each subperiod, the linear relationship H t,r =h 0,t +h 1,t r presented in equation (2) for empirical quantile scaling. Equally important, the SS-GEV scale and location parameters were also linear functions of r (e.g., Figure S6 in the supporting information). These linear relationships between the reference duration μ Ã t;r and σ Ã t;r values and the spatial scales were statistically significant for most of grid boxes. However, higher statistical uncertainty characterized the spatiotemporal scaling of the SS-GEV shape estimates. It was therefore not possible to define a unique function describing the dependence of ξ Ã t;r on the spatial scale r for all grid boxes (not shown). Figure 12 presents the spatial distribution of the temporal trends estimated at the native CRCM5 resolution for ξ Ã t;r , σ Ã t;r , and μ Ã t;r . Decadal percent variations are shown for the scale and location parameters for SD (d * =1 hr) and LD (d * =24 hr), while the annual (absolute) variation β ξ Ã t;r (year −1 ) is considered for the SS-GEV shape for LD (β ξ Ã t;r are not shown for SD since only a few grid boxes displayed significant trends). Trends estimated for H r,t are not shown, as they were consistent with results presented for the empirical quantiles in Figures 9 and 10 (Figures S7 and S8 in the supporting information).
As expected, the projected annual variations of SS-GEV shape parameters did not display coherent spatial patterns, as opposed to the location and scale parameters. Decadal percent increases for σ Ã t;r , and μ Ã t;r were large in northeastern regions, while smaller increases (e.g. ≈2%/10 year) were estimated for northwest interior areas and in the southeast region. Also, Figure 12b show some spatial patterns due to the effects hr) at the native spatial resolution (r= 12 km) for the shape, scale, and location parameters (first to third column, respectively). Smaller points represent grid boxes with no statistically significant trend. Only LD estimations are shown for ξ Ã t , as SD virtually no significant trend was detected.
of the Appalachian mountains, and the presence of lower σ Ã t;r , and μ Ã t;r increases along the northern St. Lawrence Gulf coasts for LD.
Finally, it can be noted that grid point trends for σ Ã t;r were larger than corresponding μ Ã t;r trends in some areas (e.g., over the Great Lakes region for SD; Figure 12a). Standardized measures of AM distribution dispersion, such as the normalized dispersion coefficient σ Ã t;r =μ Ã t;r , are thus expected to increase in time over these regions, as suggested by Cannon and Innocenti (2019). In fact, significant trends of the normalized dispersion coefficients were observed for a large majority of grid boxes (e.g., Figure S9 in the supporting information). However, the estimation of σ Ã t;r =μ Ã t;r trends strongly depends on the temporal evolution of the SS-GEV shape parameter, as high cross correlation exists between ξ Ã t;r and σ Ã t;r . Hence, the changes projected for the AM distribution shape in future climate corresponding to different assumptions for the evolution of ξ Ã t;r parameters in time should be further investigated, as the increase in most extreme precipitation quantiles strongly depend on the combined changes in ξ Ã t;r and σ Ã t;r .
For coarser spatial scale (i.e., r>12km), the magnitude of the trends for SS-GEV parameters decreased with increasing r (e.g., Figures S8 and S9 in the supporting information).

Discussion and Conclusion
The projected changes in the statistical properties of subdaily and daily precipitation extremes were investigated over northeastern America using the recently available initial-condition CRCM5-LE. More specifically, the impacts of climate change on AM probability distributions (empirical AM quantiles and GEV distribution parameters), the annual and diurnal cycles (date and time of AM occurrences), and the spatiotemporal structure (AM quantile and SS-GEV scaling) were analyzed.
The CRCM5-LE is the downscaled version (0.11°resolution) of the 50-member CanESM2-LE which was simulated over the 1954-2099 period under the RCP8.5 scenario. Simulated grid box precipitation series were pooled over short subperiods of 1, 3, and 7 years (1SP, 3SP, and 7SP, respectively) to create large samples and estimate AM statistics at spatial scales between 12 and 72 km and durations between 1 and 72 hr. The analysis of the AM spatiotemporal scaling was performed separately for short and long duration intervals: SD (1 hr ≤d<6 hr) and LD (6 hr ≤d≤ 72 hr). The field significance and the magnitude of temporal trends were assessed for the characteristics of the AM probability distributions, the annual and diurnal cycles of AM occurrences, and the AM spatiotemporal statistical structure.
Differences between 1SP, 3SP, and 7SP estimates were detected for AM statistics that are substantially affected by sampling errors (e.g., GEV shape and spatiotemporal scaling parameter estimates), showing the benefit of pooling the information from different members of the CRCM5-LE. However, the major conclusions of this study do not differ between 1SP, 3SP, and 7SP. Results globally confirmed the theoretical arguments pointing at the intensification of extreme precipitation over the study domain, with relatively stronger trends for short-duration AM and more extreme events. This conclusion is consistent with previous analyses over North America (e.g., Kharin et al., 2018;Mailhot et al., 2012) and was highlighted by results on AM empirical quantile, as well as from the analysis of their spatiotemporal scaling.
These findings are in agreement with the theoretical argument that, while extreme precipitation will increase with global warming at the same rate as the atmospheric moisture-holding capacity (i.e., according to the Clausius-Clapeyron relation by ≈ 7 % per degree of temperature rise), relatively large increases are expected for the most severe short-duration precipitation events, as more moisture may be brought into storms (Trenberth, 2011;Westra et al., 2014). For longer-duration events, instead, smaller increases can be expected due to the importance of dynamical processes involved in the generation of rainfall events, which intensity is not only a function of the local atmospheric moisture content, but also on the larger-scale moisture advected from surrounding areas , and references therein).
Increasing trends in AM empirical quantiles were observed for more than 99% of model land grid boxes for all considered return periods. Globally, the decadal percent increases of short-duration AM were larger than corresponding daily and longer AM estimates. Relative changes were also larger for more extreme AM quantiles and in the northeastern part of the domain.

Journal of Geophysical Research: Atmospheres
Adjusting GEV models to AM series allowed to further examine changes in extreme distribution properties. Significant increases in the location and scale parameters were obtained for all durations and grid boxes. The fraction of grid boxes with heavy-tailed AM distributions (i.e., with positive GEV shape parameters) was found to increase with time for d≥3 hr, and the estimation of the GEV shape parameter was greatly affected by the length of AM series. For instance, 3SP systematically displayed a lower fraction of grid boxes with heavy-tail distributions. Estimated 3SP shape parameter values were higher than corresponding 7SP estimates, since with larger samples it is possible to assess the significance of smaller nonzero shape parameters.
The impact of climate change on annual and diurnal cycles was then investigated. No substantial change was observed for the annual cycle in terms of the mean date of occurrence of subdaily precipitation extremes. Conversely, 24 hr and longer-duration AM will more likely occur earlier in the year (before July) in future climate and their seasonal variability is projected to increase for all durations, as the dates of AM occurrences are more dispersed around their summer peaks at the end of the simulation. Hence, more frequent warm meteorological conditions are likely to increase the probability of occurrence of intense precipitation events over an extended period of the year (e.g., Vincent et al., 2018). The projected changes in annual cycle statistics showed a well-defined spatial structure over the study domain. For instance, daily and longer-duration AM more frequently occur during fall and winter in future climate for grid boxes along the North Atlantic and Great Lakes coasts, and in the southern region East of the Appalachians. This suggests that important dynamical changes may also affect the large-scale precipitation systems that are likely to produce extremes in Northeast United States and Atlantic regions (e.g., Touma et al., 2018, and references therein), while thermodynamical considerations might explain the projected changes in the annual cycles for the interior regions. Finally, modification in the temperatures and ice cover of water bodies may be involved in these projected changes in coastal extreme precipitation.
A small shift toward later mean times of occurrence was projected for the diurnal cycles of subdaily AM in the future climate, suggesting that short-duration extremes will more likely occur later in the evening. Basic considerations on projected temperature increases (e.g., higher evening and night temperatures or possible changes in the surface heat budget) could partly explain this result. However, various processes and local to large-scale land-atmosphere feedbacks, may also be involved (e.g., and references therein Lo et al., 2017).
The temporal evolution of the spatiotemporal structure of precipitation extremes was analyzed based on the scaling properties of empirical AM quantiles and SS-GEV distributions. The empirical temporal scaling parameters h 0,t displayed decreasing trends in time, meaning that point-scale IDF curves will become steeper (negative slopes) and the variations of AM intensity quantiles across durations will increase in future climate. The decadal percent variations of h 0,t were in fact statistically significant and negative for most grid boxes for 2-year AM quantiles for both SD and LD intervals, while for longer return periods field significant trends were only found for LD. The stronger intensification of short-duration extremes suggests that major changes may occur in the type of systems producing AM at different temporal scales as well in storm dynamics, as higher h 0,t values imply that the statistical structure of precipitation extremes will be more heterogeneous across the different durations.
The trends estimated for the spatiotemporal scaling parameters h 1,t highlight that changes may be expected in the spatial structure of empirical AM quantiles. In particular, the variation of ARF values across the spatiotemporal scales (r,d) is projected to increase, as h 1,t increased with time, especially in northern regions. Therefore, the size of the weather systems generating AM at some of the considered durations is expected to change, since higher h 1,t values imply a stronger sensitivity of ARFs to the duration for a given spatial scale. However, large uncertainty affects h 1,t estimation and projected trends were not significantly different from 0 for most of the grid boxes when considering return periods q≥25 year and field significance over the entire spatial domain. Interestingly, however, the h 1,t trends showed a clear spatial structure with positive values found in the North and West of the domain and decreasing temporal trends found in the southeast. SD decadal percent variations of h 1,t were larger than for LD but a smaller fraction of grid boxes had significant trends.
The estimation of SS-GEV models confirmed that a significant intensification of the temporal scaling regimes of precipitation extremes is projected at the native CRCM5 spatial resolution. Moreover, increasing fractions of grid boxes with significant and positive SS-GEV shape parameters were found for both SD and LD.
Accordingly, the assumption of Gumbel distributions (ξ Ã t ¼ 0) is often inadequate in future climate for subdaily AM and could lead to the underestimation of long return period quantiles.
Some evidence for statistically significant trends of the values of the SS-GEV shape parameter was also observed, with approximately half of the grid boxes displaying significant increases of ξ Ã t for LD. This result is important as it suggests that the assumption of a stationary GEV distribution shape parameter, often made because of the large uncertainties affecting its estimation (Katz, 2013), may be inadequate for the analysis of daily and subdaily AM. At the same time, however, weak statistical evidence was found for trends in ξ Ã t values for SD. Instead, changes in longer return period quantiles were mainly observed through the significant increases of the scale-location parameter ratios, σ Ã t;r =μ Ã t;r . It is therefore more difficult to provide a robust assessment of extreme precipitation characteristics and their changes for the shortest durations, even when data from various temporal scales and from different model members are pooled.
Considering the ensemble of these results, caution can also be advised when using very long precipitation series for the assessment of climate change impacts on AM characteristics, as also suggested by De Gaetano and Castellano (2018). In fact, while there is a clear advantage in using large samples of precipitation extremes for improving the estimation of most uncertain AM statistics (e.g., GEV shape distribution parameters or spatiotemporal scaling statistics), substantial nonstationarity characterizes the AM probability distributions (e.g., changes in GEV location and scale parameters) over very short periods of time (e.g., 10 years).
Future works should thus consider the possibility of introducing temporal trends in SS-GEV models (e.g., for spatiotemporal scaling and/or probability distribution parameters) and further explore the relationships between the SS-GEV parameters and the spatial scale at which precipitation extremes are extracted. Further investigations of severe precipitation characteristics should also consider alternative definitions of extremes, such as event-based definitions of heavy rainfall (e.g., Dunkerley, 2008).
Moreover, besides that the presented results are globally consistent with previous studies based on simulations from other RCMs and GCMs (e.g., Kharin et al., 2018;Prein et al., 2017), it would be important to extend this study to higher-resolution climate model simulations to assess the impact of convection representation on the reported conclusion, especially those related to the diurnal cycle and short-duration AM (Coppola et al., 2018). The use of a large ensemble from a convection parametrized RCM may in fact represent a strong constraint on the detectable changes for precipitation characteristics.
Similarly, modifications in timing and seasonality of subdaily extremes, as well as their relationships with other variables (e.g., temperature and general circulation patterns) should be more thoroughly analyzed to better understand the weather processes and climate feedbacks underlying the projected changes , while other large ensembles should be considered for evaluating the impacts of structural model uncertainties on these results.
As a final remark, the practical implications of the presented results also need to be further examined, as the consequences of the highlighted changes in extreme precipitation characteristics may be critical in terms of long-and middle-term infrastructure design and engineering practice, as well as for the monitoring of natural resources. In this context, the use of the presented scaling relationship can help the development of robust statistical tools for the projection of the IDF and IDAF curves based on precipitation series simulated by climate models (Agilan & Umamahesh, 2017;Cheng & AghaKouchak, 2014).