The evaluation of possible climate change consequence on extreme rainfall has significant implications for the design of engineering structure and socioeconomic resources development. While many studies have assessed the impact of climate change on design rainfall using global and regional climate model (RCM) predictions, to date, there has been no comprehensive comparison or evaluation of intensity-duration-frequency (IDF) statistics at regional scale, considering both stationary versus nonstationary models for the future climate. To understand how extreme precipitation may respond to future IDF curves, we used an ensemble of three RCMs participating in the North-American (NA)-CORDEX domain over eight rainfall stations across Southern Ontario, one of the most densely populated and major economic region in Canada. The IDF relationships are derived from multi-model RCM simulations and compared with the stationbased observations. We modeled precipitation extremes, at different durations using extreme value distributions considering parameters that are either stationary or nonstationary, as a linear function of time. Our results showed that extreme precipitation intensity driven by future climate forcing shows a significant increase in intensity for 10-year events in 2050s (2030-2070) relative to 1970-2010 baseline period across most of the locations. However, for longer return periods, an opposite trend is noted. Surprisingly, in term of design storms, no significant differences were found when comparing stationary and nonstationary IDF estimation methods for the future (2050s) for the larger return periods. The findings, which are specific to regional precipitation extremes, suggest no immediate reason for alarm, but the need for progressive updating of the design standards in light of global warming. Extreme precipitation events have significant societal consequences and important implications for urban and rural development, public infrastructure, agriculture and human health [Rosenzweig et al., 2001;Mailhot and Duchesne, 2009;Rosenberg et al., 2010;Hatfield et al., 2011;Murray and Ebi, 2012]. Historical climate records indicate increasing trend in annual precipitation over Canada since 1950s, often statistically significant [Kunkel et al., 1999;Shephard et al., 2014]. The short-duration (sub-hourly to daily) extreme rainfall is expected to intensify due to global warming [Lehmann et al., 2015;Fischer and Knutti, 2016;Prein et al., 2016;Bao et al., 2017;Pfahl et al., 2017]. Furthermore, observed long-term trends (1951Furthermore, observed long-term trends ( -2005) ) in daily annual maximum (AM) precipitation over Northern Hemisphere land areas are attributed to anthropogenic influences, and the magnitude of increase is consistent with an increase in global mean surface temperature [Zhang et al., 2013]. Overall, the mean precipitation in Canada has increased by about 12% (5 to 35% in Southern Canada) during the period 1950 -1998, although annual snowfall has been significantly decreased over Southern region since 1950s [Ligeti et al., 2006]. Daily and sub-daily precipitation are extensively used in engineering design and infrastructural planning. The design of hydraulic infrastructure, in general, relies on short-duration extreme precipitation statistics obtained from T-year return periods at different durations, which are in turn represented in the form of intensity-duration-frequency (IDF) curves [Bonnin et al., 2006]. Although there is clear evidence that climate is changing, the traditional IDF analyses assume stationarity of the climate, means trend in extremes will not vary significantly over time. The concept of nonstationarity in water resources [Milly et al., 2008[Milly et al., , 2015]], poses a key question to the scientific community whether it is significant enough at the scale relevant to design, operation, and maintenance of hydraulic infrastructures. Further, it is an open question to stakeholders and decision makers how the changing climate information will be translated into the development of future IDFs. In Canada, the nationwide historical IDF information is produced based on short-duration groundbased rainfall records and archived at Environment Canada (EC) Engineering database [EC, 2014]. At regional and local scales, one of the fundamental challenges is the scale mismatch between climate information provided by the climate models and the impact relevant metrics [McNie, 2007;Bonnin, 2013]. The global climate models (GCMs), which are typically designed to simulate largescale atmospheric processes (large scale models; horizontal grid spacing of 150-300 km), do not directly simulate all extreme rainfall processes, such as convection [Zhang et al., 2017]. To capture fine-scale regional information for impact, adaptation and vulnerability (IAV) studies, different downscaling methods, such as statistical [Hayhoe et al., 2004[Hayhoe et al., , 2008;;Dibike and Coulibaly 2006;Stoner et al., 2013] and dynamical [Giorgi and Mearns, 1991;Giorgi, 2006] methods have been developed. The statistical method is less computationally demanding and is based on finding statistical relationships between a set of predictors and predictands [Wilby et al., 1998]. However, its performance largely relies on good quality long-term observational records and assumes that the observed statistical relationship would remain unchanged in the future [Wilby and Wigley, 1997;Wilby et al., 1998]. Furthermore, statistical downscaling approaches are generally applied to aggregate (seasonal or monthly time step) rather than daily time scales. At a daily time scale, the perfect prognosis assumption of statistical downscaling makes them quite susceptible to GCM biases [Ines and Hansen, 2006;Maraun, 2016]. On the other hand, the dynamical downscaling approach uses nested modeling in which GCM boundary conditions drive regional climate models (RCMs) and involve explicit solving of process-based physical dynamics of the system. Although dynamical downscaling is computationally expensive, it allows simulation of topographic complexity and fine-scale atmospheric dynamics and is not limited to either locations or length of records. Hence dynamically downscaled climate variables, in general, respond more physically consistent ways to external forcing [Giorgi and Mearns, 1991;Wang and Kotamarthi, 2015] as compared to their statistically downscaled counterparts. The added value [Feser et al., 2011;Gutmann et al., 2012] of dynamical downscaling method in improving spatial resolution typically by the factor of five -ten times as compared to the host GCMs, allows modelers to widely use the downscaled climate information in impact assessment studies [Gao et al., 2012;Mailhot et al., 2012;Glotter et al., 2014;Ganguli and Ganguly, 2016a;Yang et al., 2016;Yoo and Galewsky, 2016;Li et al., 2017]. However, a common problem associated with the use of RCM outputs directly in hydrologic impact assessment is that the simulated precipitation exhibits significant biases [Teutschbein andSeibert, 2010, 2012].
RCM biases stem from the intrinsic RCM physics or from the errors in the lateral boundary conditions provided by the reanalysis data or by a GCM [Schoetter et al., 2012]. Therefore, bias correction is often required to match RCM simulated rainfall to observations for impact assessment studies. The wide use of dynamically downscaled climate information for impact assessment has accelerated the growth of computational resources. Currently, many large collaborative projects are generating downscaled climate information, which is archived and publicly available. One such effort includes the (NARCCAP) since 2006, which generated high resolution (50 km) climate projections for the United States, Canada and Northern Mexico [Mearns et al., 2012]. Two of such programs are available for Europethe Prediction of Regional Scenarios and Uncertainties for and Effects (PRUDENCE) [Christensen and Christensen, 2007] and ENSEMBLES [Van der Linden and Mitchell, 2009]. Over North America, NARCCAP offers a suite of climate models consist of six different RCMs driven by four different GCMs ( which helps in estimating climate uncertainties in multi-model ensembles on a regional scales. However, the NARCCAP simulations are available only at a limited 30-year period [1970 -2000; historical and 2040 -2077, the end term future runs; Mearns et al., 2012]. Furthermore, to project future climate scenarios, the host GCMs are forced with the SRES (Special Report on Emission Scenarios) 'A2' emission scenario as applied in the III (CMIP3). In this context, the latest generation global climate models that are part of the fifth phase of (CMIP5; Taylor et al., 2012), involve over 20 modeling groups and more than 40 climate models, offers improved horizontal resolutions and physical processes in many models over its earlier version (CMIP3) [Sillmann et al., 2013a[Sillmann et al., , 2013b;;Kumar et al., 2014]. (CORDEX) is a program started by the (WCRP) to develop fine-scale regional climate projections within the CMIP5 modeling framework in the IPCC AR5 [Intergovernmental Panel of 5; Stocker et al., 2013] timeline. The CORDEX North American domain involves several universities and partner institutions to simulate high-resolution regional climate information (horizontal resolution of ~ 50 km), driven by several global climate models ( and archived at Earth System Grid (ESG) data portal for regional impact assessment studies. So far, most efforts to develop future IDF curves are limited to individual city or regional scale assuming stationarity in precipitation extremes [Coulibaly and Shi, 2005;Prodanovic and Simonovic, 2007;Hassanzadeh et al., 2013;Rana et al., 2013;Rodriguez et al., 2014;Srivastav and Simonovic, 2014;Srivastav et al., 2014;Chandra et al., 2015;Elshorbagy et al., 2015;Lima et al., 2016;Li et al., 2017;Vu et al., 2017]. Mirhosseini et al., [2013] analyzed the impact of climate change in the IDF curves at the city of Auburn, Alabama using NARCCAP RCMs. Although their analysis suggests, future precipitation pattern tends towards less intense rainfall for short duration events, for longer durations (> 4-hour) the results are not consistent across the models. Wang et al. [2014] developed projected IDF curves for the city of Toronto using dynamically downscaled precipitation data from the regional climate modeling system (PRECIS) driven by a version 3 forced with SRES A1B emission scenario. Kuo et al. [2015] assessed climate change impact on future IDFs over the city of Edmonton in central Alberta using MM5 (the Pennsylvania State University/National Center for Atmospheric Research numerical model) regional climate models with boundary conditions from four different global climate models (CGCM3, ECHAM5, CCSM3, and MIROC3.2) forced with A2, A1B, and B1 emission scenarios from CMIP3. Their analysis suggests overall; design storm intensities would gradually increase from 2020 to 2080s for short duration AM intensities of less than 1-hour and return period less than 25-year. Simonovic et al. [2016] assessed the impact of climate change on future IDFs across 567 meteorological stations across Canada using a web-based IDF_cc tool [Simonovic et al., 2016]. The results for the end term period (the year 2100) was based on the second generation (CanESM2) regional climate model and the multi-model ensemble median of 24 global climate models. The CanESM2-based simulations indicated a reduction in AM intensities in central regions of Canada and increased in other regions, whereas ensemble median GCMs simulated lower value of design storm intensity than the RCM. At the regional scale, using third generation (CRCM) Mailhot et al., [2007] performed regional frequency analysis for grid boxes over Southern Quebec region. Comparison of regional frequency estimates of control (1961-1990) versus future (2041-2070) climatic condition reveals return periods corresponding to 2-and 6-hour events would be approximately halved in the projected scenario, whereas they will decrease by a third for 12-and 24-hour events. In an another study, Mailhot et al., [2012] compared return periods of sub-daily, daily, 3-, and 5-day AM precipitation intensity for the future (2041-2070) versus historical periods using NARCCAP simulations over fourteen climatic regions of Canada and part of Northern United States. Their analysis suggested for 2-to 25-year return periods, inland regions, especially Southern Ontario, the Prairies, and Southern Quebec would experience the largest increase in AM intensities while the coastal regions, the least. Switzman et al., [2017] characterized variability among an ensemble of future IDF curves generated using a combination of 5-different regional and global climate models, 3-climate change scenarios, and 2-downscaling methods. For 2050s, the authors found statistically significant variability in the direction of change and magnitude in the projected IDFs at the Greater Toronto Area (Oshawa WPCP, Toronto, ), with some ensemble members showing increase and decrease in rainfall intensity within the ensemble members with no specific trend. On the other hand, over Windsor regions (Windsor Airport, Chatham WPCP, and Harrow), a statistically significant increase in storm intensity is noted for the future. However, there is an evidence of variability in the magnitude of change among ensemble members. Recently, [Agilan andUmamahesh, 2016, 2017], compared nonstationary IDFs for the future (2015-2056 and 2057-2098) versus historical (1961 -1990) period using an ensemble of 24 global climate model simulations over the city of Hyderabad, India.
The results of the study indicate that the estimated return periods using nonstationary IDFs with the trend as a covariate can reasonably cope the effect of climate change for at least next fifty years. From existing literature on future IDF development, it may be conjectured that to date more scientific efforts are being towards on characterizing the uncertainty in climate change projections than on developing robust adaptation strategies from a range of plausible climate outcomes [Wilby and Dessai, 2010]. Secondly, till date, no study has been carried out comparing stationary and nonstationary approaches to design storms in a future climatic condition considering climate uncertainty at a regional scale. Third, a robust assessment of post-processing of climate risk information is required for translating climate change scenarios into decision-relevant metrics through appropriate bias correction schemes taking into account locations, durations and ensemble of climate model runs. Effective management of climate-induced risks requires robust characterization of the probability of extremes accounting both historical nonstationarity and the likelihood of future changes [Mote et al., 2016;Diffenbaugh et al., 2017]. The cities across Southern Ontario, the prominent economic hubs [Bourne and Simmons, 2003;Partridge et al., 2007], are especially vulnerable to extreme climatic events given their larger population concentration, property, and aging infrastructure than any other part of Canada [Kling et al., 2003;Ligeti et al., 2006;Hayhoe et al., 2010;Henstra and Thistlethwaite, 2017]. Based on regional climate model experiments, two of the recent literature has indicated an increase in frequency and magnitude of rainfall extremes in Southern Ontario by the 2050s [Mailhot et al., 2012;Deng et al., 2016]. Given these challenges, this study has three preliminary objectives:  To evaluate the credibility of the suite of regional climate models participating in the NA-CORDEX program in attributing extreme rainfall statistics at different durations and return periods with reference to ground-based historical observations. To select appropriate bias correction methods based on various skill measures for translating RCM information to hydrometeorological impact assessment. To evaluate the impact of climate change on IDF curves by comparing design storms of historical versus projected climate assuming both stationary and nonstationary conditions. Thus we performed a three-way comparison between historical and projected IDFs for impact assessment: stationary (future) versus stationary (historical), nonstationary (future) versus stationary (historical) and nonstationary (future) versus nonstationary (historical). To address these research objectives, we analyze the trend and frequency of AM precipitation and temperature data of Southern Ontario region from the station-based observational record. We characterize extreme precipitation frequency at different durations using the concept of "T-year return period" (i.e., likelihood of the particular event to occur is once in every T-year) and develop stationary and nonstationary IDF relations for the historical and projected time periods. The typical planning horizon of 30 ~ 40 years is used by climatologist to filter out natural variability as well as a tradeoff between data length required to perform nonstationary analysis [Palutikof et al., 1999;WMO, 2009]. In this study, estimates of 50-year return period of AM hourly and daily intensity are considered as the representative of rare and extreme events. This return period length and the durations are in agreement with previous studies [Palutikof et al., 1999;Kharin et al., 2007;Mailhot et al., 2012;Wehner, 2013]. Since RCM simulations are available at a daily time scales, the RCM model output is temporally downscaled using a multiplicative cascade-based disaggregation technique [Olsson, 1995[Olsson, , 1998;;Guntner et al., 2001;Rana et al., 2013]. The downscaled climate information is corrected by employing the quantile-based parametric and nonparametric bias correction methodologies [Panofsky and Brier, 1965]. First, the statistical characteristics of ground-based extreme rainfall events for the historical period, 1970 -2010 (centered around 1990), are compared with extreme rainfall statistics from an ensemble of regional climate models. Then, comparative analyses of the future climate, 2030 -2070 (centered around 2050), relative to the past are performed to estimate expected changes in extreme precipitation statistics for the resilient design of water resources, reliability and resources planning [Mills, 2005;Milly et al., 2008]. Since the 2050s represent a reasonable planning horizon for multisector stakeholders' perspective [Hall et al., 2014;Ganguly et al., 2015;Watts et al., 2015], we consider RCM simulations from 2030 -2070 as projected planning period. We use RCM outputs to project future precipitation for the 2050s considering the high emission scenario, Representative Concentration Pathways 8.5 (RCP8.5 hereafter; Meinshausen et al., [2011]) to maximize adaptation and mitigation strategies. We also explore the effect of bias-correction taking into account the duration and location of observation on model performance for historical simulation and projected scenario. A flowchart of this study is shown in Figure S1. The organization of the paper is as follows: Section 2 describes study area, observational data, regional climate models, the bias correction method and evaluation statistics; Section 3 presents historical evaluation at sub-daily and daily scales at different return periods; Section 4 discusses relative changes in precipitation intensities under stationary and nonstationary conditions. The summary and discussion follow in Section 5. The study is conducted across eight rain gauge locations over Southern Ontario, specifically the Southwest-northeast transect, in the southernmost region of Canada. The Southern Ontario is surrounded by three Great Lakes, Erie, Huron, and Ontario (Jien and Gough, 2013). The study sites include eight rain gauge locations of the Windsor-Kingston corridor (in the order from southwest to northeast): Stratford wastewater treatment plant (WWTP), Shand Dam in Fergus on the (WPCP), Trenton Airport (Figure S2). (DEM) of the study area indicates a shallow slope with a maximum altitude of 670 m above mean sea level (MSL). The proximity to Great Lakes and topographic effect, especially in areas to the lee of and the Georgian Bay significantly modifies the climate in the region. Convective storms and thunderstorms primarily modulate the summer rainfall, but fall rainfall is dominated by reduced convective activity and increased lake effect precipitation [Lapen and Hayhoe, 2003;Baldwin et al., 2011].
Since a recent observation-based study has indicated a steady increase in the global warming trend from the 1970s onwards [Rahmstorf et al., 2017], we selected the baseline period for the current analysis from 1970 to 2010. Station-based historical (1970Station-based historical ( -2010) ) AM precipitation time series at durations (d = 1-, 2-, 6-, 12-, 24-hour) are obtained from ( for the eight rain gauge locations over Southern Ontario. The available data is thoroughly quality controlled [Shephard et al., 2014] and have been previously analyzed for the assessment of national extreme rainfall trends [Burn and Taleghani, 2013;Shephard et al., 2014;Simonovic et al., 2016;Switzman et al., 2017]. The extent of missing values in the sub-daily AM rainfall time series ranges between 2 and 20% (average ~ 13%) with the least being in Toronto (only AM rainfall in the year 2005 is missing) and the highest are in Hamilton (1970( , and 2004( -2010( are missing) and Stratford (1973( , 1999( , 2005( -2010) ) respectively. We obtained daily and hourly rainfall records and daily maximum air temperature data from the EC Historical database ( and (TRCA; We infilled missing values and updated the extreme precipitation records till 2010 by successively disaggregating daily rainfall values to hourly and sub-hourly time steps using multiplicative random cascade (MRC)based disaggregation tool [Olsson, 1995[Olsson, , 1998;;Guntner et al., 2001]. The details of the disaggregation procedure of historical precipitation time series are in the supplementary section. We used archived prediction data from three RCMs available at NA-CORDEX domain. The specific models include, fourth generation of the (CanRCM4) driven by the second generation of the (CanESM2); fifth generation of the (CRCM5) driven by CanESM2; and Regional Climate Model version 4 (RegCM4) nested in version 2 (HadGEM2-ES) global climate model. The choice of RCMs are based on their extensive use of the current and previous versions over North America for high-resolution multi-decadal climate change simulations [Ashfaq et al., 2010;Separovic et al., 2013;Singh et al., 2013;Naz et al., 2016;Whan and Zwiers, 2016;Jalbert et al., 2017]. All models' outputs are available in a grid mesh of 0.44 deg horizontal resolution at a daily time step for the historical and projected scenario, except CanRCM4-CanESM2, for which historical simulations are available at an hourly time step. Use of common horizontal resolution removes large sources of variability between RCMs and helps us to evaluate how the differences in the configuration of the RCMs influence simulation of extremes [Whan and Zwiers, 2016]. Both CORDEX generation Canadian RCMs (CanRCM4;[von Salzen et al., 2013;Scinocca et al., 2016] and CRCM5; [Martynov et al., 2013]) share the common dynamical core, however, differ in nesting strategy employed, and land-surface and physics schemes [Zadra et al., 2003;Whan and Zwiers, 2016]. For computational purposes, all climate model outputs are regridded to a common grid point of 0.5 deg latitude/longitude resolution using bilinear interpolation scheme from the climate data operators [CDO, 2017]. Further, we consider multi-ensemble approach (multi-model median and associated bounds, defined by minimum and maximum simulation of AM series) to take into account the internal variability of the climate system, which is particularly suitable for near-term impact assessment [Hawkins and Sutton, 2009;Meehl et al., 2009;Hawkins and Sutton, 2011]. While the ensemble median (multi-model median, hereafter MM-Med) values represent the most likely case, the ensemble minima (multi-model minimum, hereafter MM-Min) and maxima (multimodel maximum, hereafter MM-Max) are considered as the best and worst case scenarios, which indicate the spread of climate model [Ganguli and Ganguly, 2016b;Vousdoukas et al., 2017]. Grid-based RCM simulations are downloaded at the four nearest neighbor values of station-based observation, and a distance weighted average remapping was employed [CDO, 2017]. Since RCMs, on average simulates a small amount of precipitation on regular time steps, as a threshold for discrimination between wet and dry days a value of 0.1 mm is chosen [Wehner, 2013;Frei et al., 2003]. For bias correction of RCM output, we employ quantile mapping (QM) for historical (1970 -2010) RCM simulations. For applying bias correction to climate model simulated rainfall, different options are available. For example, bias correction can be applied to entire time series, which then can be used to extract AM series for the IDF development. On the other hand, since only AM rainfall is required for the IDF development, it is also possible to correct AM precipitation. [Li et al., 2017] found that bias correcting AM rainfall based on empirical distribution followed by frequency analysis yields design storm closest to the observations. Hence, we employ quantile mapping bias correction on RCM simulated historical AM time series. In quantile mapping, a quantile of the present day simulated distribution is replaced by the same quantile of the historical observed distribution [Maraun, 2016]. Given a precipitation time series x, the method is formulated as [Ines and Hansen, 2006;Maraun, 2016;Vu et al., 2017]: Where, x is the bias-corrected precipitation,   condition. Since RCMs simulate too many wet days (the 'drizzle effect'), the QM is automatically able to adjust the number of wet days [Gutowski Jr et al., 2003;Hay and Clark, 2003]. Based on distributional choices, QM can be both parametric and nonparametric. However, for high quantiles, where sampling noise is high, nonparametric QM may produce noisy results and applies random correction [Maraun, 2016] Where h K is the kernel function at a bandwidth 'h'. Following a previous study [McGinnis et al., 2015] We employ Gaussian kernel function and Silverman's rule of thumb for 'h' calculation. However, the basic assumption of QM is that the future distribution properties (such as variance and skew) remain similar to the reference period, and only mean changes. However, owing to nonstationarity of the climate system, this assumption may not hold true in the future [Milly et al., 2008;Li et al., 2010]. Hence, following previous studies [Li et al., 2010;Srivastav et al., 2014;Vu et al., 2017], we applied Equidistant Quantile Mapping (EQM) to bias correct projected AM, which is presented as follows: However, Eq.
( 3) cannot guarantee that the corrected AM series will have positive value (for example, Multimodel minimum 24-hour AM series in London and Multimodel maximum 12-hour AM series in Toronto when KDE-based bias correction is applied). Hence, in these cases, following an earlier study [Wang and Chen, 2014], we employed equiratio QM to the AM series using following mathematical expression: We used non-parametric Mann-Kendall test with correction for ties and autocorrelation [Hamed and Rao, 1998;Reddy and Ganguli, 2013] to detect signature of monotonic trends in the observed AM historical time series at different durations (d = 1-, 2-, 6-, 12-, 24-hours). In addition, we estimate the trend using Theil Sen's slope estimator [Sen, 1968;Gilbert, 1987]. To understand the extent of warming we also present results of trend estimates for AM temperature anomaly records of historical observations. Significance of trends are detected at 10% significance level (i.e., pvalue < 0.10). (GEV) distribution is a combination of three different distributions, viz., Frechet ( 0   ), Weibull ( 0   ) and the Gumbel ( 0   ) depending on the sign of the shape ( ) parameter. The cumulative distribution function (CDF) of stationary GEV model is given as [Coles et al., 2001]   Where, max ,0 yy   with '+' sign indicates positive part of the argument,  is a location parameter,  is a scale parameter and the shape parameter,  determines the symmetry and heaviness of the tail [El Adlouni et al., 2008]. For nonstationary GEV model, we incorporate timevarying covariates into GEV location (GEV I), to describe trends as a linear function of time (in years), i.e.,   10 tt    . Since we limit our planning horizon to 40 years, and modeling temporal changes in shape and scale parameters requires long-term records, following an earlier study [Cheng et al., 2014a], we assume these two parameters as constant. For estimation of GEV parameters, a Bayesian inference is performed combined with (DE-MC) Monte Carlo (MC) simulation as suggested by [Martins et al., 2000;Cheng et al., 2014b;Cheng and AghaKouchak, 2014]. For the AM time series, the parameters are derived by computing 50 th (median), 5 th and 95 th (the lower and the upper bounds) of the DE-MC sampled GEV parameters. For stationary GEV model, the design storm intensity, p q associated with T-year return period is obtained using following expression [Coles and Tawn, 1996]     The computation of nonstationary design storm intensity is similar to the stationary model except the inclusion of time-varying location parameter [Cheng et al., 2014b;Cheng and AghaKouchak, 2014]. We perform the calculations following Cheng and AghaKouchak, [2014] using an MATLAB-based software package, (NEVA), Version 2.0]. The evaluation of observation and models during present-day (1970 -2010) climatic condition is performed using following test statistics: 1. Time-series plot of AM series of observed versus multi-model ensemble climate models before and after bias correction. 2. Taylor diagrams [Taylor, 2001] of AM series to assess pattern error over eight locations. To better assess models' performance, besides Taylor diagrams, we employ a skill score based on the model's spread and correlation reference to the observation [Taylor, 2001]: Where R is the pattern correlation, f  is the normalized standard deviation, 0 R is the maximum attainable correlation calculated from the maximum of inter-ensemble correlation. 8. The percentage of change between future (calculated from MM-Med ensemble run) and historical period can be attributed to changes in precipitation extremes. Further, the (statistically) significant differences of future versus baseline design intensity is computed using standard z-statistics, which is given by [Mikkelsen et al., 2005;Madsen et al., 2009] Where zT future is the T-year intensity obtained from the MM-Med climate model ensemble from a stationary/nonstationary model, zT baseline describes the same but with the present-day condition. The denominator indicates predictive uncertainty; var{zT future } and var{zT baseline } are the estimated variance obtained from the T-year event estimates and associated confidence interval (i.e., 5 th , and 95 th quantiles). The z-statistic can be interpreted as statistically equivalent to quantiles of standard normal distribution, i.e., z = +/- 1.96 and +/- 1.64 correspond to 5% and 10% significance levels. The null hypothesis of the test assumes that T-year event estimate in the future planning period is significantly different from the baseline value. Therefore, we graphically depicted bias corrected models' performance using Taylor diagrams [Taylor, 2001]. A Taylor diagram provides comprehensive statistical summary of how well models agree with observation in the form of a radial plot where distance from the origin is a normalized standard deviation of the model output (relative to observation) and the cosine of the azimuthal angle is given by the centered pattern correlation factor of model output with respect to observation [Wehner, 2013]. root mean square error relative to the observation whereas centered pattern correlation provides information of similarity in pattern between model and observation. In general, an agreement for the daily AM series is slightly better than the agreement in hourly AM series. The heat maps of skill score (Eqn. 7) across all durations are graphically presented in Figure 4 and Table S2. The skill score improves gradually from higher to lower temporal resolution (Table S2 -S6). The variability among models is small, and a moderate skill score (varies from more than 0.3 to 0.6) is noted for MM-Med climate ensemble across different durations (Table S1). Among individual model performance, RegCM4 shows the highest skill score over Oshawa whereas CRCM5 exhibits the lowest skill over Fergus Shand dam for 1-hour AM series. Likewise, for daily AM series, RegCM4 performed the best over London, whereas worst over Trenton. In many cases, we found that the single-model performance is better than or comparable to that of the multi-model climate ensembles. However, as noted in the previous studies [Kumar et al., 2014;Min et al., 2014;Martre et al., 2015;Ganguli and Ganguly, 2016], the main advantage of using a multi-model ensemble is not the vast improvement over the best performing single model but the consistently better performance of multimodel considering all aspects of projection.
Finally, we chose the method of bias correction based on the skill score. We further evaluate the fidelity of multi-model climate ensembles (MM-Med and associated bounds based on the method of bias correction) relative to observation using probability density function [PDF; Perkins et al., 2007]. Figures 5 and 6 show PDFs of hourly and daily AM series relative to baseline observations. For sub-daily AM series (2-to 12-hour), the density plots are shown in the supplement (Figures S6 -S8). The PDFs are created using a smoothed empirical distribution with Gaussian kernel density function. We found a reasonably well agreement between climate model ensembles and observed AM series; nonetheless, the multi-model RCMs perform better for the hourly extremes as compared to the daily. For hourly extreme, the best agreement is observed for Toronto, while the RCMs do not adequately capture a few outliers on far right tail of Hamilton. For daily extreme, MM-Med simulates reasonably well for Toronto and Stratford. However, for other sites, although multi-model bounds consistently simulate well, MM-Med ensemble exhibits higher peaks relative to the observations, indicating an overestimation of the frequency of extremes. Figure 7 shows time series plots of observed versus bias corrected AM series for all eight locations. Figure 7 indicates improvement in the performance of RCMs in simulating peaks than that of before bias correction . Next, we analyze the presence of trends in the AM series at different durations to detect the signature of nonstationarity in the time series. For this, first, we analyze the presence of a trend in historical AM series and compare it with CORDEX RCMs. The ability of RCMs in simulating historical trends provides the information regarding the credibility of models in simulating anthropogenic climate change in a future scenario. Tables S7-S11 summarizes results of trend analysis in the historical and MM-Median CORDEX RCMs. In general, we found agreement in the signs of the trend between observation and model, although exceptions exist in many cases. For example, at hourly and sub-hourly rainfall extremes in Toronto and Windsor. In both of these locations, while the observed AM series shows a decreasing trend, the model exhibits an upward trend. Further, in all cases, models over/under-estimate the magnitude of the trend. However, neither in observations nor models, the magnitude of trend estimate is zero. Further, to know the extent of warming we discuss the changes in historical air temperature (Table S12) in the daily AM time series. Except for Windsor, we found an upward trend in historical air temperature at all station locations. During baseline period (1970 -2010), Stratford shows the highest warming trend (at the rate of 0.46/decade), followed by Toronto (0.37/decade). However, none of these trends are statistically significant. At Windsor, a decrease in temperature extreme is at the rate of -0.035/decade. Since shape parameter directly influences the nature of simulated extremes [Ragulina and Reitan, 2017], we compare the shape parameters of observed versus RCMs for both stationary and nonstationary GEV models. Figures 8-9 and Table S13 present the shape parameters simulated by MM-Med climate ensembles relative to baseline observation. We found RCMs simulates shape parameters reasonably well, although variability is more for nonstationary GEV models . The shape parameter shows spatial and temporal variations , however, no specific pattern is identified. In most of the cases, the shape parameters vary in the range of -0.30 to +0.30 [Table S13]. For daily precipitation extreme, the spatial average of observed shape parameters across the stations is  = 0.13 with highest in Toronto recent study [Ragulina and Reitan, 2017], which estimated a global average of shape parameter around 0.14, with a 95% confidence interval between 0.13 and 0.15 using daily gridded precipitation extremes from the -Daily database (version 2.60). Further, the authors showed for Southern Ontario the estimates of shape parameter varies from 0.13 -0.16 [Fig. 5 in Ragulina and Reitan, 2017], which is little different in our case, given the fact that we are dealing with extreme precipitation and the ground-based observations are subjected to biases resulting from multiple sources [New et al., 2001;Adam and Lettenmaier, 2003]. Next, we investigate the variation of shape parameter in nonstationary models . The shape parameters for daily observed AM series varies between = 0.41) respectively. On an average, the shape parameter simulated by the nonstationary model is higher than that of the stationary model, indicating the former adequately simulates the presence of heavy upper tail (or higher likelihood of extreme magnitude events) [Markose and Alentorn, 2011;Halmstad et al., 2013]. Next, we compare the simulated IDF curves of observed versus CORDEX RCMs. Figures 10 -12 illustrate a design storm estimates for T = 2-, 5-, 10-, 25-, and 50-year return periods at three selected station locations: Toronto International Airport, Oshawa WPCP, and Windsor Airport. Weak trends also have considerable impacts on results of exceedance probability [Porporato and Ridolfi, 1998]. Therefore, we modeled nonstationarity considering time as a covariate in the location parameters of the GEV distribution. The uncertainty in design storm is shown using the vertical boxplots. Tables S14 -S16 show parameters of the fitted GEV distributions for stationary and nonstationary models. Table S17 presents a statistical performance of stationary versus nonstationary GEV models for the durations of 1-and 6-hour using Akaike Information Criteria with small sample correction [AICc; Akaike, 1974;Hurvich and Tsai, 1990]. The skill of the individual model, as measured by the AICc statistics does not exhibit any consistent patterns either based on the locations or the storm durations. In general, we found CORDEX RCMs able to simulate observed design storm intensity reasonably well for all three locations (Tables S14 -S16, and Figures 10 -12). This is in contrast with the results of trend estimates in Tables S7 -S11, where we noted the discrepancy in nature (sign) of trends simulated by the CORDEX RCMs. Further, Tables S14 -S16 show the fitted GEV parameters for MM-Med CORDEX ensemble are close to that of the observed AM series.
As noted from the plots, the uncertainty in design storm estimates is higher for the lower storm durations (d = 1-, and 2-hour) and higher return periods (T = 25-, and 50-year). For Toronto and Oshawa, at durations, d = 1-to 6-hour, the simulated rainfall intensity with higher return periods tend to underestimate the observed ones. However, for Windsor Airport, we find overestimation of design intensity by the models at d = 1-hour, whereas underestimation at d = 6-hour. Secondly, for Oshawa, we note an increase in observed rainfall intensity in nonstationary models reference to stationary one at short durations of 1-, and 2-hour and return periods, T = 10-year and beyond. This may be attributed to consistent upward trend in AM rainfall intensity across all durations in Oshawa. Further, we find at shorter durations, in nonstationary models, the length of top whisker in the boxplot tend to be longer than the bottom one, indicating the model adequately describes the asymmetric nature of extremes. Our findings are consistent with an earlier study [Cheng and AghaKouchak, 2014], in which authors reported the difference between nonstationary and stationary rainfall intensities are more prominent at shorter durations and higher return periods, and the differences gradually diminish at longer durations. For detailed exploratory analysis, we present relative changes in rainfall intensity for 1 in 10-year (i.e., T = 10-year) and 1 in 50-year events respectively . While the top panel in the figure shows stationary models, the bottom one shows the nonstationary ones, corresponding to design rainfall intensity at T = 10-year and 50-year return periods, which can be themselves are categorized as less and more intense events. At lower return period, T = 10-year, the relative difference between observation and model is close to zero, while the differences tend to increase at higher return period, T = 50-year. Few exceptions exist, for example, London International Airport shows large negative bias for 24-hour rainfall extreme. Next, we present ensemble median and spread of RCMs in simulating extreme rainfall intensity relative to observation (Figures 14-15). For stationary models and 1-hour extreme (Top panel; Overall, our analysis suggests CORDEX RCMs show higher confidence in simulating design storms at smaller return periods while lower at return periods of T = 25-year and beyond. In this section, we discuss changes in shape parameter of simulated extreme in future relative to the baseline period. 16) and nonstationary models in the near-future, the shape parameter shows a decrease in value, across most of the durations relative to the baseline period. Table S18 shows during 2050s, the shape parameters tend to be negative across all durations irrespective of the nature of simulations. A negative shape parameter indicates the GEV distribution has a bounded upper tail (i.e., an identified upper limit to those extreme events or a Weibull-type distribution), implying a tendency towards heavy lower tail. For stationary GEV models, the values of shape parameter are most negative at 12-hour duration, and ranged from -0.57 (at Trenton) followed by -0.49 (at Toronto) to -0.06 (at Hamilton) with spatial standard deviation of 0.29. At same duration, for nonstationary GEV models,  values ranged from -0.66 (at Trenton) followed by -0.65 (at Toronto) to 0.33 (at Stratford) with spatial standard deviation 0.09. These values in shape parameter are more negative than those conventionally considered as physically reasonable in hydrology literature [Martins et al., 2000]. While  is less than -0.33, the AM distribution has an infinite third moment, and when  is less than -0.5, the distribution has infinite variance as suggested by an earlier study [Morrison and Smith, 2002]. A very negative value of the shape parameter suggests that the distribution of AM precipitation tends to have a very heavy tail. However, we found at shorter durations (1 to 6-hour), the  values ranged between -0.33 and 0.23, and well within the physically consistent limit. Table S19 lists the performance of stationary and nonstationary GEV models using AICc statistics for the selected (i.e., 1-and 6-hour) storm durations. As observed from the To further evaluate how the future intensities are sensitive towards the choice of frequency analysis (i.e., stationary/ nonstationary), the boxplots of difference between future and historical precipitation extremes for 10-year and 50-year events are presented in Figures 20 and 21. The absolute changes in future design storm intensity relative to observations are shown for three different cases: stationary (future) versus stationary (historical); nonstationary (future) versus stationary (historical); and nonstationary (future) versus nonstationary (historical). For 1 in 10year events, except Stratford, at all locations, the quantile boxes are above zero for 6-and 24-hour durations. For Stratford, at 24-hour duration, the intensity of rainfall remains relatively unchanged over time as indicated by the location of quantile boxes, which is close to zero irrespective of the choice of the return period estimates. For hourly extreme, the increase is observed for Hamilton, Oshawa and Trenton. In contrast, for 1 in 50-year events, except a few stations, for most of the stations, the quantile boxes are below zero for all durations, indicating a decrease in precipitation intensity in the future period. However, we find an increase in rainfall intensity at 2-hour and 1hour extremes for Windsor and Oshawa respectively. For 24-hour precipitation extreme, the intensity of extreme precipitation remains relatively unchanged for 1 in 50-year events . Further, for both events, although the differences in rainfall intensity are apparent for durations up to 2-hour they follow the same trend irrespective of the choice of frequency analysis. However, the magnitude of differences fade away at durations of 6-hour and beyond. To quantify the magnitude of changes in projected intensity relative to the baseline period, we present heat maps of the difference in rainfall intensity of the 2050s versus 1990s for 10-, 25-and 50-year return periods (Figures 22-23 in contrast to nonstationarity assumption (Figure S11; bottom panel), in which magnitude of changes are deemed to be significant.
For 50-year return period , assuming nonstationarity, we found a statistically significant decrease in return period for sub-hourly rainfall extremes in Toronto, London, Trenton and Fergus Shand. However, the changes are not significant assuming stationarity in return period estimate. Considering stationarity, the percentage decrease in rainfall intensity ranges between 2 -37%, while it varies from 2 to 46% across different locations and durations. Taken together, the following broader insights emerge for future precipitation extremes in Southern Ontario: (i) For 1 in 10-year events, an increase in rainfall intensity is observed across most of the stations and durations, irrespective of the choice of frequency analysis. (ii) For more extreme events (i.e., 1 in 25-year and 1 in 50-year), a decrease in rainfall intensity is noted; however, the magnitude of changes in storm intensity are deemed to be statistically significant considering non-stationarity assumption in contrast to the stationary assumptions. Extreme precipitation poses significant risks for cities across Southern Ontario, Canada due to their dense population, valuable and geographically concentrated property, complex and interdependent infrastructure networks [Henstra and Thistlethwaite, 2017]. We characterize extreme rainfall frequency over Southern Ontario using IDF statistics with a suite of highresolution new generation of coordinated regional climate modeling experiments participating in the NA-CORDEX domain with an elevated global greenhouse forcing. Since the study area comprises urban catchments and most of the municipal infrastructures are designed based on design storm intensity of 50-year or less, the IDF statistics are compared for 2-, 5-, 10-, 25-, and 50-year return periods. The RCM performance in simulating extreme rainfall statistics in the present day climate is validated using ground-based AM observational records. The temporally downscaled NA-CORDEX RCM ensembles are able to simulate the observed IDF statistics reasonably well as demonstrated by various performance metrics. Further, we modeled design rainfall intensity using extreme value distributions under stationarity and nonstationarity assumptions with a time-varying location parameter. The following conclusions emerge from the present analysis:  The ground-based short-duration historical (1970 -2010) precipitation extremes often exhibit a statistically significant increase (or decrease) trends for the selected locations across Southern Ontario. On the other hand, extreme daily temperature anomaly shows an upward trend for most of the locations, although the apparent trend is not statistically significant across the stations. A robust bias-correction methodology is employed to correct systematic biases in the RCM outputs for different time scales and across individual station locations. However, the skill of bias correction does not follow any specific trend either based on locations or the durations. Statistically significant increase in rainfall intensity is observed for 10-year events over most of the station locations. However, opposite trend is noted for events with longer return periods. The paper presents a proof of principal results to compare changes in extreme precipitation intensity in the near-future using three regional climate models with the future atmosphere represented by the RCP8.5 scenario. While Canada's national climate change assessment indicates an increase in average surface air temperature and precipitation volume across several regions nationally [Bush et al., 2014]; still most of the ground-based individual station observations do not exhibit any significant trend. Further, typically a limited number of stations show a statistically significant nonstationary trend at a sub-daily and daily time scales [Shephard et al., 2014]. Recently, using station-based AM precipitation [Shephard et al., 2014] have shown that nationwide less than 5% of the stations exhibits significant increasing (or decreasing) trends. Based on their findings the authors recommend that in Canada, the traditional IDF design based on stationarity assumption of precipitation extreme has not been violated. Further, in projected scenario, we find an increase in rainfall intensity at 6-and 24-hour extremes for 10-year events over majority of stations. Our results are consistent with [Barbero et al., 2017], where authors reported except in winter months, there is a significant increase in daily AM precipitation trends as compared to the hourly extremes after analyzing a large number of ground-based observational records across the United States for the present day (1950-2011) condition. Their results indicate that at a station level, the trends in daily precipitation extremes are better detected than that of the hourly extremes in a changing climate. They attributed this to the limited spatial extent of the most extreme events and the inability of sparse station network to accurately measure such events [Barbero et al., 2017]. Some of the caveats of the study include: The heavy computational requirement involving successive spatiotemporal downscaling of extreme precipitation at individual station level has limited us to use only three RCMs. However, uncertainty in climate change projection results from both model responses and internal climate variability [Hawkins andSutton, 2009, 2011]. A series of previous studies [Deser et al., 2012[Deser et al., , 2014;;Wettstein and Deser, 2014;Kay et al., 2015] have examined the effect of intrinsic climate variability and anthropogenic climate change on projected climate using a 40-member climate model ensembles available from National Center for Atmospheric Research (NCAR)'s CESM (LENS). These studies explored various aspects of uncertainty resulting from intrinsic climate variability and found substantial climate uncertainty on a global and regional scale. The role of internal climate variability is considerably higher at a regional scale and in a shorter lead time (in a near-term planning horizon) [Figure 4 in Hawkins and Sutton, 2009;Hawkins and Sutton 2012]. Therefore, using three-member RCMs may not be sufficient to capture the full spectrum of climate variability in a near-term planning period. Next, hourly and sub-daily extreme precipitation, required for engineering design, are often produced by convective events. However, global and regional climate models are not able to simulate such events well because of models' limited spatial and temporal resolution, and convection is not explicitly resolved in these models [Lenderink and Fowler, 2017;Zhang et al., 2017]. Although high resolution convective permitting models may provide a reliable representation of local storm dynamics [Kendon et al., 2014], these models still suffer from inherent uncertainties in simulating extremes [Hagos et al., 2014;Singh and O'Gorman, 2014]. Secondly, the cascade of uncertainty may also stem from temporal downscaling of RCM output.
Third, we assume that the selected bias-correction method at historical period is time independent and same will holds good for correcting bias in the future period. The Cascade-based temporal downscaling method in this study is adopted from [Olsson, 1995[Olsson, , 31 1998a]]. The method implements a micro-canonical (exact conservation of mass in each cascade branching) 32 approach, in which daily rainfall is disaggregated using a uniformly distributed generator, dependent on 33 rainfall intensity and position of the rain sequence. The method was later developed and tested in different 34 climatic conditions [Guntner et al., 2001;Rana et al., 2013]. In the disaggregation process, each time 35 interval (box) at a given resolution (for example 1 day) is split into two half of the original length (1/2 day). 36 The procedure is continued as a cascade until the desired time resolution is reached, i.e., to 14 day, then to 37 1/8 of a day and so on. Each step is termed as a cascade step, with cascade step 0 as the longest time-period 38 with only one box (i.e., a day). The distribution of the volume between two sub-intervals (or smaller boxes) 39 is computed by multiplication with the cascade weights ( , Depending on the range of resolution involve, P01 either be assumed as resolution independent or 44 parameterized as a scaling law: Pr 01 (r) = c 1 r c 2 where c1 and c2 are constants. The distribution of , ir W is 45 termed as cascade generator, assumed to follow 1-parameter beta distribution [Olsson, 2012]. Following 46 [Olsson, 1998b], the probabilities P, the probability distribution of cascade generator are assumed to be 47 volume is large then it is more likely that both halves of the subintervals contribute to nonzero volume than 53 if the volume is small. Following [Olsson, 1998b] [Olsson, 1998b, n.d.]. 58 The objectives of employing cascade-based disaggregation tool are as follows: (i) to infill missing values 60 in the station-based AM precipitation time series (ii) to disaggregate NA-CORDEX regional climate model 61 output available at daily time step to an hourly time scale. First, the cascade parameters are calibrated using 62 good quality observed historical dataset available at an hourly resolution. The high-resolution historical 63 (1970 -2010) data are collected from (TRCA). At least 35-64 years or more available dataset is used to calibrate the model. Table S1 lists the length of the dataset used 65 for model calibration. The parameters are calibrated using 5-cascade steps, which implies, five successive 66 "halving" from one day to generate 45-minute (2700 seconds) data. After calibration, Monte Carlo 67 simulation is performed to gradually fine-grain the data and generate realizations at desired Rs. 68 Since Rs in this cases is not directly achieved by exact resolution doubling from original resolution, Rl, the 69 downscaled with bias-corrected precipitation for Toronto Airport. Following a similar procedure, we 79 infilled the missing gaps in the AM series of eight station locations across Southern ON. Fig. S3 shows 80 observed (red) versus disaggregated infilled (blue) hourly AM time series. As observed from the figure, the 81 infilled time series preserve the trend in the data reasonably well, without exhibiting any signature of 82 outliers (i.e., extremely high or low values). 83 84 Likewise, to disaggregate the daily NA-CORDEX output to hourly time step, we assume that dynamically 85 downscaled model output is statistically similar to that of station observed precipitation. For this, we 86 compare the distribution of daily precipitation extremes of observed versus NA-CORDEX RCMs (Fig. S4). 87 Fig. S4 presents boxplots of the anomaly of daily maximum precipitation for the observed and regional 88 climate model ensembles. As observed from the figure, the RCMs able to reproduce the observed daily 89 precipitation satisfactorily, although some bias exist. Among individual model performance, in general, 90 CanRCM4 shows maximum variability (or spread) in simulating precipitation extremes than that of the 91 other two RCMs. To disaggregate the daily RCM data to an hourly time step, we use the calibrated model 92 parameters derived in the historical simulation to downscale precipitation in the baseline and future periods. Then two bias correction schemes are compared (Kernel-and GEV-based) Gringorten's plotting position formula.