Assessment of Future Changes in Intensity-Duration-Frequency Curves for Southern Ontario using North American (NA)-CORDEX Models with Nonstationary Methods

1 Assessment of Future Changes in Intensity -Duration-Frequency Curves f or Southern Ontario using North American (NA)-CORDEX Models with Nonstationary Methods 1 McMaster Water Resourc es and Hydrologic Modeling Group, Department of Civil Engineering, McMaster University, 1280 Main Street West, Hamilton ON L8S 4L7, Canada Abstract The evaluation of possible climate chan ge consequence on extreme rainfall has significant implications for the de sign of e ngineering structure and socioeconomic resources dev elopment. While many studies have assessed the impact of climate change on design rainfall using global and r egional climate model (RCM) predictions, to date, there has been no comprehensive co mparison or evaluation of intensity-duration-frequency ( IDF) statistics at regional scale, considering both stationar y versus nonstationary models for the future climate. To underst and how extreme precipitation ma y respond to future IDF curves, w e used an ensemble of three RCMs participating in the North-American (NA)-CORDEX domain over eight rainfall stations across Southern Ontario, one o f the most dens ely popula ted and m ajor economic region in Can ada. The IDF relationships are de rived from multi-model RCM simulations and compare d with the station- based observations. We modeled precipitation e xtremes, at different durations using extreme value distributions considering parameters that are either stationary or nonstationar y, as a linear function of time. Our results showed that ex treme precipitation intensit y driven b y future climate forcing shows a significant increase in intensity for 10 - year events in 2050s (2030-2070 ) re l ative to 1970 - 2010 ba seline p eriod across most of the locations. How ever, for longer return periods, a n opposite trend is noted. Surprisingly, in term of design stor ms, no significant differences were found when comparing stationary and nonstationary IDF estimation methods for the future (2050s) for the larger return periods. Th e findings, which are specific to re gional precipitati on ex tremes, suggest no immediate reason for alarm, but the need for progressive updatin g of the design standards in light of global warming. 1. Introduction Extreme precipitation events have significant so cietal consequences a nd important im plications for urban and rural development, public infrastructure, agriculture and hum an health [ Rosenzweig et al. , 2001; Mailhot and Duchesne , 2009; Rosenberg et al. , 2010; Hatfield et al. , 2011; Murray 2 and Ebi , 2012]. Historical climate records indicate increasing trend in annual precipitation over Canada since 1950s, often statisticall y 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 ] . Furthe rmore, observed lon g-term trends (1951-2005) in daily annual maxim um (AM) precipitation over Northern Hemisphere land areas are attributed to anthropoge nic influences, a nd the magnitude of increase is consistent with an increase in global mean sur face temperature [ Zhang et al. , 2013] . Overall, the mean precipitation in Canada has incr eased b y ab out 12% (5 to 35% in Southern Canada ) during the period 1950 – 1998, although annual snowfall ha s been significantl y decreased over Southern region since 1950s [ Ligeti et al. , 2006]. Daily and sub-dail y precipitation are extensively used in engineering desi gn and infrastructural planning. The design of hydraulic infrastructure, in general, rel ies on short -duration extreme precipitation statis tics obtained from T- year return periods at different durat ions, which are in turn represented in the fo rm of intensity-duration-frequency (IDF) curves [ Bonnin et al. , 2006]. Although there is clear evidence that climate is changing, the traditional IDF anal yses assume stationarity of the climate, means trend in extremes will not vary si gnificantly over time . The concept of nonst ationarity in wa ter resou rces [ Milly et al. , 2008, 2015], poses a key question to the scientific communit y wh ether it is significant eno ugh at the scale r elevant to desi gn, operation, and maintenance of h ydraulic infrastructures . Further, it is an open questi on to stakeholders and decision makers how the changing climate inform ation will be translated in to the development of future IDFs. In Canada, t he nationwide historical I DF information is produced based on short-duration ground- based rainfall records and archived at Environment Canada (EC) Engineering database [ EC , 2014] . At regional and local scales, one of the fundam ental challenges is the sc ale mismatch between climate information provided b y the climate models and the impact relevant metrics [ McNie , 2007; Bonnin , 2013] . The global climate models (GCMs), which are typically d esigned to simulate large- scale atmospheric processes (large s cale models; horizontal grid spacing o f 150 -300 km), do not directly simulate all extreme rainfall processes, such as convection [ Zhang et al. , 2017] . To capture fine-scale regional infor mation for impact, ad aptation and vulnerabilit y (IAV) studies, different 3 downscaling methods, such as statisti cal [ Hayhoe et al. , 2004, 2008; Dibike and Coulibal y 2006 ; Stoner et al. , 2013] and dynamical [ Giorgi and Mearns , 1991; Giorgi , 2006] methods have been developed. The s tatistical method is less computationally dem anding and is based on finding statistical relationships between a set of predictors and predictands [ Wilby et al. , 1998]. How ever, its performance largel y relies on good quality long -term observ ational rec ords and assumes that the observed statisti cal r elationship would remain unchanged in the future [ Wilby and Wigley , 1997; Wilby et al. , 1998]. Furthermore, statisti cal downscaling approaches are generall y applied to aggregate (seasonal or monthl y time step) rather than dail y time scales. At a daily time scale, the perf ect prognosis ass umption of statistica l downsc aling makes them quite susceptible to GCM biases [ Ines and Hansen , 2006; Maraun , 2016]. On the other hand, the dynamical downscaling approach uses nested mo deling in which GCM boundary conditions drive regional climate models (RCMs) and involve explicit solving of process-ba sed physical dy n amics of the s ystem. Although dynamical downscaling is computationally ex pensive, it allows simulation of topographic complexity and fine -scale atmospheric dynamics and is not limited to either locations or length of records. H ence dynamically downscaled climate variables, in general, res pond 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 re solution typically by the factor of five - ten times as compared to the host GCMs, allows modelers to widely use the downscaled climate info rmation in impact assessment studies [ Gao et al. , 2012; Mailhot et al. , 2012; Glotter et al. , 2014; Ganguli and Ganguly , 2016a; Yang et al. , 2016; Y oo and Galewsky , 2016; Li et al. , 2017]. However, a common problem associated with the use of RCM outputs directly in h ydrologic impact assessment is that the simulated precipitation exhibits significant biases [ Teutschbein and Seibert , 2010, 2012]. RCM biases stem from the intrinsic R CM phy sics or from the errors in the latera l boundary conditions provided by the reanaly s is data or by a GCM [ Schoetter et al. , 2012] . Therefore, bias correction is often required to match RCM simulated rainfall to observations for impact assessme nt studies. The wide use of dynamical ly downscaled climate information for impact assessment has accelerated the growth of computational resources. Currently, man y large collaborative projects are generating downscaled climate information, which is archived and publicl y available. One such 4 effort includes the North American Re gional Climate Change Assessment Program (NARCCAP) since 2006, which generated high resolution (50 km) climate projections for the United States, Canada and Northern Mexico [ Mearns et al. , 20 12]. Two of such programs are available for Europe – the P rediction of Regional Scenarios and Uncertainties for D e fining European Climate Change Risks and Effects (PRUDENCE) [ Christensen and Christensen , 2007] and ENSEMBL ES [ Van der Linden and Mitchell , 2009]. Over North America, NARCCAP offers a suite of climate models consist of six different RCMs driven by fo ur different GCMs (http://www.narccap.ucar.edu/data/model-info.html), which helps in estimating climate uncertainties in multi-model ensembles on a regional scales. However, the NARCC AP simulations are av ailable onl y at a li mited 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 Coupled Model Intercomparison Project Phase III (CM I P3). In this context, th e latest generation global climate models that are part of the fifth phase of Climate Model Intercomparison P roject (CMIP5 ; Taylor et al. , 2012), involve over 20 modeling groups and more than 40 climate models , offers im proved horizontal resolutions and ph ysical processes in man y models over its earlier version (CMIP3) [ Sillmann et al. , 2013a, 2013b; Kumar et al. , 2014] . The Coordinated Regional D ownscaling Experiment (CORDEX) is a program started b y th e W orld Climate R esearch Program (WCRP) to develop fine-scale regional climate projections within the CMIP5 modeling framework in the I P CC AR5 [Inter governmental Panel of Climate Change F ifth Assessment Report Phase 5; Stocker et al. , 2013] ti meline. The CORDEX N orth American domain involves several universities and partner inst itutions to simulate high-resolution regional climate information (horizontal resolution of ~ 50 km), driven by several global climate models ( https:/ /na- cordex.org/ ) and archived at Earth S ystem Grid (ESG ) data portal for re gional impact assessment studies. So far, most efforts to d evelop futur e IDF curves are limited to individual city o r regional scale assuming stationa rity in precipitation extremes [ Coulibaly and Shi , 2005; Prodanovic and Simonovic , 2007; Hassanzadeh et al. , 2013; Rana et al. , 2013; Rodríguez 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 5 climate change in the IDF curves at the city of Auburn, Alabama using NARCC AP RCMs . Although their ana lysis suggests, future precipitation pattern tends towards less intense ra inf all for short duration events, fo r longer durations (> 4 -hour) the results are not consistent across the models. Wang et al. [2014] developed proje cted IDF curves for the city of Toronto using dynamically downscaled precipitation data from th e regional climate modeli ng s ystem (PREC IS) driven by a Hadle y Cent er Coupled Mod el version 3 forced with SRES A1B emission scenario . Kuo et al. [2015] assessed climate change im pact on future IDFs over the cit y o f Edmonton in central Alb erta using MM5 (the Penns ylvania State Universit y/National Center for Atmospheric Research numerical model) reg ional climate mode ls with boundary conditions from four different global climate models (CGCM3, ECHAM5, CCSM3, and M IROC3.2) forced with A2, A1 B, and B1 emission scenarios from CMIP3. Their analysis suggests overall; design storm intensities would g r adually increase fr om 2020 to 2080s for s hort dura ti on AM intensities of less than 1-hour and re turn period less than 25-y e ar. Simonovic e t al. [2016] a ssessed 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 y ear 2100) was based on the second generation Canadian Earth System Mod el (CanESM2) regional climate model and the multi-model ensemble median of 24 global cli mate models. The CanESM2 -based sim ulations indicated a reduction in A M intensities in central regions of Canada and increased in other reg ions, whereas ensemble median GCMs simulated lower value o f design storm i n tensity than the RCM . At the regional scale, using third generation Canadian Regional Climate Model (CRCM) Mailhot et al. , [2007] performed regional frequency analysis for grid box es over Southern Quebec region. Comparison of regional frequency estimates of control (1961 -1990) vers us future (2041-2070) climatic condition reveals return periods corresponding to 2- and 6- hour e vents would be approximately halved in the projected scenario, whereas the y will decrease b y a third for 12 - and 24-hour events. In an an other study, Mailhot et al. , [2012] compared r eturn periods of sub-dail y, daily, 3- , and 5-day AM precipitation intensity for the future (2041-2070) v ersus historical (1968- 2000) periods using NARCCAP simul ations over fourte en climatic regions of Canada and part of Northern United States. Their anal ysis suggested for 2- to 25-year return p eriods, inland regions, especially Southern Ontario, the Prairies, and Southern Quebec would experience the lar gest increase in AM intensities while the coastal regions, the least. Switzman et al. , [2017] characterized variability among an ensemble of future IDF curves g enerated using a combination of 5- different 6 regional and global clima te models, 3- climate change scenarios, and 2- do wnscaling 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, Toronto Pearson International Ai rport), with some ensemble members showin g increase and decrease in rainfall intensit y within the ensemble members with no specific trend. On the other h and, ove r Windsor re gions (W indsor Airport, Chatham WP CP, and Harrow) , a statistically significant increase in storm intensity is noted for the future. However, there is an evidence of va riability in the magnitude of change among ens emble members. Re cently, [ Agilan and Umamahesh , 2016, 2017], compared nonstationar y IDFs for the future (2015 -2056 and 2057- 2098) versus hist orical (1961 – 1990) period using an ensemble of 24 global climate model simu lations over the cit y o f Hyderabad, India. The results of the study indicate that the estimated return periods using nonstationary IDFs with the trend as a covariate can reasonabl y cope the ef fect of climate change for at least next fifty years. From existing li terature on future IDF developm ent, it may be conjectured that to date more scientific efforts are bein g towa rds on characterizing the uncertainty in climate change projections than on developing robust ada ptation strat egies fro m a range of plausible climate outcomes [ Wilby and Dessai , 2010]. Secondly, till date, no stud y 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 c hange scenarios into decision -relevant m etrics through appropriate bias correction schemes taking into account loca ti ons, durations and ensemble of climate model runs. Effective management of climate -induced risks requires robust characterization of the pr obability o f 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 especiall y v ulnerable to ex treme climatic events given their larger population concentration, property, a nd a ging infrastructure than any other part of Canada [ Kling et al. , 2003; Ligeti et al. , 2006; Hayhoe et al. , 2010 ; Henstra and Thistlethw aite , 2017]. Based on regional climate model experiments, two of the recent literature has indica ted an increase in frequency and 7 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 c limate models pa rticipating 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 approp riate bias cor rection methods based on various skill measures for translating RCM information to hydrometeorologic al impact assessment.  To evaluate the impact of climate change on IDF curves b y comparing design storms of historical versus projected climate assumin g both s tationary and nonstationa ry condition s. Thus we p erform ed a thr ee-way comparison betw een historical and projected IDFs for impact assessment: stationary ( future ) versus stat ionary (historical), nonstationary (future) versus stationary (historical) and nonstati onary (future) versus nonstati onary (historical) . To addre ss these research objectives, we ana lyze the trend and frequency of AM pre cipitation 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 pa rticular event to occur is once in every T- year) and develop stationary and nonstationary IDF rel ations for the historical and projected time period s. The t ypical planning ho rizon of 30 ~ 40 years is used b y climatologist to filter out n atural vari ability as well as a tradeoff betw een dat a len gth required to perform nonstationary anal ysis [ Palutikof et al. , 1999; WMO , 2009]. In this study, estimates of 50- year return period of AM hourly and dail y int ensity are considered as the representative of rare and extreme events. This return period len gth and the durations are in agreement with previous studies [ Palut ikof et al. , 1999; Kharin et al. , 2007; Mailhot et al. , 2012; Wehner , 2013]. Since RCM simulations are available at a dail y time scales, the RCM model output is temporally downscaled using a multiplicative cascade-based disaggregation technique [ Olsson , 1995, 1998; Güntner et al. , 2001; Rana et al. , 2013]. The downscaled climate information is c orrected b y employing the qu antile-based parametric and non - parametric bias correction methodologies [ Panofsky and Brier , 1965]. First, th e statisti cal characteristics of ground-based extreme rainfall events for the historical period, 1970 – 2010 (centered around 1990), are compared wit h ex treme rainfall statis tics from an ensemble of regional 8 climate models. Then, co mparative anal y ses 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 wa ter resources, re liability 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; W atts et al. , 2015], we consider RCM simu lations fr om 2 030 – 2070 as projec ted planning period. We use RCM outputs to project future precipitation for the 2050s considering the hi gh emission scenario, Represent ative Concentration Pathways 8.5 (RCP8.5 hereafter ; Meinshausen et al. , [2011]) to maximize adaptation and miti gation strategies. We also explore the effect of bias-correction taking into account the duration and location of obs ervation on model performance for historical simulation and projected scenario. A flowchart of thi s stud y is shown in Figure S1. The organization of the paper is as follows: Section 2 describes study area, ob servational data, regional climate models, the bias correction method a nd evaluation statis tics ; Section 3 presents historical evaluation at sub -daily and dail y scales at different return periods; Section 4 dis cusses relative ch anges in precipitation intensities under st ationary and nonstationary c onditions. The summary and discussion follow in Section 5. 2. Data and Methods 2.1 Study Region The study is conducted a cross eight rain gauge locations over Southe rn Ontario, specifically the Southwest-northeast transect, in the southernmos t region of Can ada. The S outhern Ontario is surrounded by th re e Great Lakes, Erie, Huron, and Ontario ( Jien and Gough , 2013 ). The stud y sites include eight rain gauge locations of the Windsor -Kingston corridor ( in the order from southwest to northeast ): Windsor Airport, L ondon International Airport, Stratford wastewater treatment plant (WWTP), Shand Dam in Fergus o n the Grand River, Hamilton Airport, Toronto International Airport, Oshawa Water Poll ution Control Plant (WPCP), Trenton Airport (Figure S2). The Digital Elevation Model ( DEM) o f the stud y area indicates a shallow slope with a maximum altitude of 67 0 m above mean sea level (MS L). The proximit y to Great Lakes and topographic effect, especially in areas to the lee of Lakes Erie, Lake Ontario, and the Georgian Bay significantly modifies the climate in the region. Convective storms and thunderstorms 9 primarily modul ate t he summer rainfall, but fall rainfall is dominated by reduced convective activity and increased lake effect precipitation [ Lapen and Hayhoe , 2003; Baldwin et al. , 2011]. 2.2 Station-based Observational Data Since a rec ent observation -based stud y has indi cated a steady increase in the g lobal w arming trend from the 1970s onwards [ Rahmstorf et al., 2017], we selected the b aseline period for the current analysis from 1970 to 2010. S tation-based historical (1970 – 2010) AM precipitation time series at durations ( d = 1-, 2-, 6-, 12-, 24- hour) are obtain ed from Environment Canada En gineering Database ( http://climate.weather.gc.ca/prods_servs/engineering_e.html ) fo r the eight rain gauge locations over Southern Ontario. The available da ta is thoroughly quality controlled [ Shephard et al. , 2014] and have be en previousl y anal yzed for the assessment of national ex treme rainfall trends [ Burn and Taleghani , 20 13; Shephard et al. , 2014; Simonovic et al. , 2016; Switzman e t al. , 2017] . The extent of missing values in the sub-daily AM rainfall ti me series r anges between 2 and 20% (average ~ 13%) with th e least being in Toronto ( only AM rainfall in the year 2005 is missing) and the highest are in Hamilton (1970, and 2004 – 2010 are missi ng) and Stratford (1973, 1999, 2005 – 2010) respectively. W e obtained dail y and hourly rainfall records and dail y maximum air temperature data from the EC Historical datab ase (htt p://climate.weather.gc.ca/historical_data/) and Toronto Region Con servation Authorit y (TRC A ; https://trca.ca/). W e infilled missing values and updated the extreme pre cipitation records till 2010 b y successively disaggregating daily rainfall values to hourl y and sub-hourl y ti me step s using multiplicative random cascade (MRC) - based disaggregation tool [ Olsson , 1995, 1998; Güntner et al. , 2001]. The details of the disaggregation procedure of historical precipitation time series are in the supplementary section. 2.3 Regional Climate Model We used archived prediction data from three RCMs available at NA -C ORDEX domain. The specific models include, fourth generation of the C anadian Regional Climate Model (CanRCM4) driven b y the second generation of the Canadian E arth S ystem Model (CanESM2); fifth generation of the Canadian Regional Climate Model (CRCM5) driven b y C anESM2; and Regional Climate Model version 4 (RegCM4) nested in Hadley Center Global Environmental Model, version 2 (HadGEM2- ES ) global climate model. The choice of RCMs are bas ed on their extensive use of the current and pr evious versions over North Am erica for high -resolution multi-decadal climate 10 change simulations [ Ashfaq et al. , 2010; Šeparović 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° horiz ontal resolution at a dail y time step for the historical and proj ected scenario, ex cept CanRCM4-CanESM2, for which historical simulations are available at an hourl y ti me step . Use of common horizontal resolution re moves large sources of va riability betwe en RCMs and helps us to evaluate how the difference s in th e confi guration of the RCMs influence simulation of extremes [ Whan and Zwiers , 2016 ]. Both CORDEX generation Canadian R CMs (CanRCM4 ; [ von Salzen et al. , 2013; Scinocca et al. , 2016] and CRCM5; [ Martynov et al., 2013] ) share the common dy namical core, however, diffe r in n esting strat egy emplo ye d, 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° latit ude/longitude resolution using bilinear interpolation scheme from the climate data operators [ CDO , 2017]. Further, we consider mul ti-ensemble approach (m ulti-model median and associated bounds, defined b y minimum and max imum simulation of AM series ) to take into account the internal variability of the climate s y stem, which is particularly suitable for near -t erm impact assessment [ Hawkins and Sutton , 2009; Meehl et al. , 2009; Hawkins and Sutton , 2011]. While the ensemble median (multi -model median, hereafter MM-Med) val ues represent the most likely case, the ensemble mi nima (multi -model minimum, hereafter MM-Min) and maxima (multi- model maximum, hereafter MM-Max) are considered as the best and worst case scenarios, which indicate the sp read of clim ate model [ Ganguli and Ganguly , 2016b; Vousdoukas et al. , 2017]. Grid-based RCM simulations are downloade d at the four nea rest neighbor values of station-based observation, and a dista nce weighted average r emapping was employed [ CDO , 2017]. Since RCMs, on average simulates a small amount of precipitation on re gular time steps , as a thr eshold for discrimination between wet and dr y da ys a value of 0.1 mm is chosen [ Wehner , 2013 ; Frei et al. , 2003]. 2.4 Bias Correction For bias correction of RCM output, we emplo y quantile mapping (QM) for hist orical (1970 – 2010) RCM simulations. For applying bias cor rection to climate model simulated rainfall, different options are available. For example, bias cor rection c an be applied to entire time series, whic h then 11 can be used to extract AM series for the IDF de velopment. On the other hand, si nce onl y AM rainfall is required for the IDF development, it is also possible to correct A M precipitation. [ Li et al. , 2017] found that bias correcting AM rainfall based on empirical distribution followed by frequency analysis y ields design storm closest to the observations. Hence, we employ quantile mapping bias correction on RCM simulated historical AM time serie s. In quantile mappin g, a quantile of th e present day simulated dist ribution is replaced b y the same quantile of the historical observed distribution [ Maraun , 2016]. Given a pr ecipitation time series x , the method is formulated as [ Ines and Hansen , 2006; Maraun , 2016; Vu et al. , 2017]:   1 m c o c m c m c x F F x         (1) Where, x is the bias-corrected precipitation,   F  is the CDF and   1 F   is its inverse of either the observations (‘o’) or model (‘m’) for the historical training period or current climate (‘c’) condition. Since RCMs simulate too many wet da ys (the ‘drizz le e ffect’), th e QM is automa tically able to adjust the number of wet d a y s [ Gutowski Jr et al. , 2003; Hay and Clark , 2003]. Based on distributional choices, QM can be both para metric and nonparametric. However, for high quantiles, where samplin g noise is high, nonparametric QM ma y produce noisy results and applies random correction [ Maraun , 2016]. Hence, we evaluate the performance of both parametric [based on Generalized Extreme Value (GEV)] and nonparametric [based on Kernel Density Estimate (KDE)] distributions for QM of hist orical AM series. For parametric QM, the GEV distribution with three parameters (location, scale and shape) are first fitted to observed and current climate model AM series. Then CDF of the GEV distribution is fitted to the current model data, which is again mapped to the CDF of the observed data. An inverse CDF transformation is then employed to bias correct AM series of the modelled time series [Eqn. 1]. For KDE -based bias correction, nonparametric kernel density function is fitted to the AM series. Mathematicall y, KDE   ˆ fx is [ McGinnis et al. , 2015] :     1 1 ˆ n hi i f x K x x n    (2) Where h K is the kernel function at a bandwidth ‘ h ’. Following a previous study [ McGinnis et al. , 2015] We employ Gaussian kerne l function and Silverman’s rule of thumb for ‘ h ’ calculation. 12 However, the basic assu mption of QM is that the future distribution properties (such as variance and skew) remain similar to the reference period, and only mean changes. However, o wing 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 e t 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:     11 m p m p o c m p m p m c m p m p x x F F x F F x                     (3) However, Eq. (3) can no t guarantee that the corrected AM series will have positi ve value (for example, Multim odel 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 equir atio QM to the AM series using following mathematical expression:     1 1 o c m p m p m p m p m c m p m p F F x x x F F x                (4) 2.5 Trend Test W e used non-parametric Mann-Kendall test with correction for ties and autocorrelation [ Hamed and Rao , 1998; Reddy and Ganguli , 2013] to detect sig nature of monotonic trends in the observed AM hist orical 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 wa rming we also present results of trend estimates for AM temperature anomaly records of historical observation s. Significance of tr ends are detected at 10% si gnificance level ( i.e. , p - value < 0.10). 2.6 Extreme Value Analysis: Stationary and Nonstationary G EV Models The Generalized Extreme Value (GEV) distribution is a combination of three different distributions, viz. , Fréchet ( 0   ), W eibull ( 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] 13   1 ex p 1 0 , , z Gz                                (5) Where,   m ax , 0 yy   with ‘+’ si gn indicates positive part of the argument,  is a location parameter,  is a scale parameter and the shape parameter,  determines the s ymmetry and heaviness of the tail [ El Adlouni et al. , 2008]. For nonstationary GEV model, we inco rporate ti me- varying covariates into GEV location (GEV I), to describe trends as a li near function of time (in years) , i.e. ,   10 tt     . Since we limit our planning horiz on to 40 y ears, 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 Ba yesian inference is performed combined with Differential Evaluation Markov Chain (DE-MC) Monte Carlo (MC) sim ulation as sugg ested b y [ Martins et al. , 2000; Che ng et al. , 2014b; Cheng and AghaKouchak , 2014]. For the AM time series, the p arameters 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 stationar y GEV model, the design storm int ensity, p q associated with T- year return period is obtained using following expression [ Coles and Tawn , 1996]     1 ln 1 0 p qp              (6) The computation of nonstationary design storm intensity is similar to the stationary model except the inclusion of time -varying loc ation parameter [ Cheng et al. , 2014b; Cheng and AghaKouchak , 2014]. We perform the calculations following Cheng and AghaKouchak , [2014] using an MATLAB-based softwa re package, Nonstationary Extreme Value Analysis (NEVA), Version 2.0]. 2.6 Historical Model Agreement and Relative Changes 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. 14 2. Taylor diagrams [ Taylor , 2001] of AM series to assess pattern e rror ov er eight locations. To better assess models’ performance, besides Ta ylor diagrams, we employ a skill score based on the model ’ s spread and correlation re ference to the observation [ Taylor , 2001]:       2 0 41 ˆˆ 11 ff R S R     (7) 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. 3. Densit y plots of observed versus multi -model ensemble climate models at different durations. 4. Box plot of the shape parameters for observed versus MM -Med ensemble to evaluate the credibility of the models’ to simulate the fat-tailed behavior of extremes. 5. Intensit y ( I ) v ersus Duration ( D ) plots of obs erved versus MM -Med ensemble at r eturn periods, T = 2-, 5-, 10-, 25-, and 50- year for stationar y and nonst ationary methods of analysis for thre e significant stations, viz., Toronto I ntern ational Airport, Oshawa WPCP , and Windsor Airport. 6. Vertical box plots of the relative difference between observed and MM -Med ensemble design storm for all eight locations for the duratio ns 1 - and 24-hour and T = 10-, and 50- year return p eriods. Th e relative changes are ob tained throu gh the difference between observed and MM -Med ensemble desi gn storm estimates normalized b y the baseline values. 7. Intensit y ( I ) ve rsus Frequency ( F ) plots of observe d versus multi-model climate ense mbles at return periods, T = 2-, 5-, 10-, 25-, and 50-year for stationar y and nonstationar y GEV models for all eight locations. 8. The percentage of change between future (calcul ated f rom MM -Med ensemble run) and historical period can be attributed to chang es in precipitation extremes. Further, the (statistically) significant differe nces of future versus ba seline design intensity is computed using standard z-statistics, which is given by [ Mikkelsen et al. , 2005; Madsen et al. , 2009]            󰇡 󰇥     󰇦     󰇢 (8) 15 Where    is the T -year intensit y obtained from the MM -Med climate model ensemble from a stationary/nonstationary model,    describe s the same but with the pre sent-d ay condition. The denominator indicates predictive uncertainty;      and      are the estimated variance obtai ned 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 norm al distribution, i.e., z = ± 1.96 and ± 1.64 correspond to 5 % and 10% si gnificance levels. The null h ypothesis of the test assumes that T -year event estimate in the future planning period is significantly differe nt from the baseline value. 3. RCM Evaluation of Historical Simulation 3.1 Climatological Statistics Figure 1 compares ti me series plot s of observed versus climate ensembles before bias -correction for 1- and 24-hour AM series at eight observation sites. Although models are in -phase with observations, in few cases (for example, London International Airport and Windsor Airport), multi-model ensemble underestimates obs erved AM series . Hence, we bias correct ed AM series by th e pro cedure d escribed in Section 2.4, before performing frequency a nalysis. As a first ste p, we obtain ed statistics from bias-corrected AM series and compared it with observed AM to comprehend if the model simulations are credible enough to project f uture climate change. Therefore, we graphically depict ed bias corrected models’ perform ance using Ta ylor dia grams [ Taylor , 2001] . A Tay lor diagram provides comprehensive statis tical summ ary of how well models agree with observation in the form of a radial plot where distance from the origin is a normaliz ed standard deviation of the model output (relative to observation) and the cosine of the azimuthal angle is given b y the centered pattern correlation factor of model output with respect to observation [ Wehner , 2013]. Figures 2 and 3 present Taylor diagrams of max imum hourly and dail y precipitation extremes for the eight rain gauge loca ti ons across Southern Ontario. The bias- corrected ( i.e., both GEV and KDE -based) thr ee RCMs, their multi -model ensembles, and observations are represented by different s ymbols (as in le gend ), and the different bias correction schemes are shown with different colours. The statisti cs of observed AM is shown on the x-ax is as a refer ence point . The perfect sim ulation would be plotted at a unit distance from the origin at an angle of 0°. The radial distance from model’s p oint in the Taylor diagram is proportional to its root mean square error relative to the observation whereas centered pattern correlation provides 16 information of similarity in pattern between model and observation. In general, an a greement for the daily AM series is slightly better than the agreement in hourly AM series. The h eat maps of skill score (Eqn. 7) across all durations are graphically presented in Fi gure 4 and Table S 2. The skill score im proves gradually from hi gher to lower temporal resolution (Table S2 – S6 ). Th e variability among models is small, and a moderate skill score (varies from more than 0.3 to 0.6) i s noted for MM -Med clim ate ensemble across diff erent durations (Table S1). Among indi vidual model performance, RegCM4 shows the highest skill score over Oshawa w hereas CRCM5 ex hibits the lowest skil l over Fergus Shand dam for 1 -hour AM series. Likewise, for dail y AM series, RegCM4 performed the best over London, whereas worst over Trenton. In many cases, we found that the single-model per formance 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 e t al. , 2014; Martre et al. , 2015; Ganguli and Ganguly , 2016], the mai n advantag e of usin g a multi-model ensemble is not the vast improve ment over the best performi ng single model but the consistently better performance of multimodel considering a ll aspects of projection. Finally, we chose the method of bias correction based on the skill score. W e further evaluate the fidelity of multi -model climate ensembles (M M-Med and associated bounds based on the method of bias correction) relative to observation usi ng probability densit y function [PDF ; Perkins et al. , 2007]. Figures 5 and 6 show PDFs of hourly and dail y AM series relative to bas eline observations . For sub-daily A M series (2- to 12- hour), the densit y plots are shown in the supplement (Figures S6 – S8). The P DFs are created usin g 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 ag reement 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 reasonabl y w ell for Toronto and Stratford. However, for other sit es, although multi-model bounds consistently si mulate well, MM-Med ensemble exhibits higher peaks relative to the observations, indi cating an overestimation of th e frequen cy of ex tremes. Figure 7 shows time se ries 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 (Figure 1). 17 3.2 Analysis of Trend Next, we analyze the pr esence of trends in the AM series at different durations to detect the signature of nonstationarity in the time serie s. For this, first, we analyze the presence of a tre nd in historical AM series and compare it with CORDEX RCM s. The abilit y o f RCM s in simulating historical trends provides the information regarding the credibilit y of models in simulating anthropogenic climate change in a future s cenario. 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 th e trend bet ween observation and m odel, although exceptions exist in many cases. For example, at hourly and sub -hourly rainfall extremes in Toronto and Windsor. In both of the se 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 w e di scuss the changes in hist orical air temperature (Table S12) in the dail y AM time se ries. Except for Windsor, we found an upward trend in historical air temperature at all station locations. During baseline period (1970 – 2 010), Stratford shows the highest warmin g trend (at the rate of 0.46/decade), followed b y Toronto (0.37/decade). However, none of these trends are statistically si gnificant. At Windsor, a decrease in temper ature ex treme is at the rate of - 0.035/decade. 3.3 Simulation of GEV Shape Parameter Since shape pa rameter directly influences the nature of simulated e xtremes [ Ragulina and Reitan , 2017], we compare the shape parameters of observed ve rsus RC Ms for both stationary and nonstationary GEV models. Fi gures 8-9 and Ta ble S13 present the shape parameters sim ulated by MM -Med climate ensembles relative to baseline observation. W e found RC Ms simulates shape parameters reasonabl y well, although variabilit y is more f or nonstationary GEV models (Figure 8). The shape parameter shows spatial and tempo ral variations (Figure 8), 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 dail y precipitation extreme, the spatial average o f observed shape parameters across the stations is  = 0.13 with highest in Toronto ( Med  = ~ 0.25 ; Med  indicates shape parameter obtained from the 50 th percentile of DE-MC simulat ed samples) and least in Fergus 18 Shand dam (median Med  = 0.05). On the other hand, t he spatial average simulated by MM -Med CORDEX RC M slightly overestimates the obs erved mean,  = 0.15 with highest being in Stratford ( Med  = ~ 0.2 9) and the least in London ( Med  = 0. 04). Our results are consi stent with a recent stud y [ 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 Global Historical Climatology Network -Dail y database (v ersion 2.60). Further, the author s 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 p arameter in nonstationar y models (Figure 9). The shape parameters for daily obs erved AM s eries var ies between Med  = 0.06 (for Fergus Shand Dam) and Med  = 0.31 (for Stratfor d) with an average estimate of  = 0.19. On the other hand, for CORDEX RCM, we fou nd models sl ightl y underestimates the spatial average shape parameter;  = 0.16, which varies between Med  = 0.31 in Stratford and Med  = 0.05 in London. Evaluation of shape par ameter across different ti me steps sug gests that the most positive shape pa rameter valu es occur for 12-hour storm duration in Stratford (in stationary simulation, Med  = 0.34) and Toronto (in nonstationar y simulation, Med  = 0.41) r espectively. On an average, the shape parameter simulated by th e nonstationary model is higher than th at of the stationary model, indicating the former adequately simulates the p resence o f heavy upper tail (or hi gher likelihood of extreme magnitude events) [ Markose and Alentorn , 2011; Halmstad et al. , 2013]. 3.4 IDF Curves: Stationary versus Nonstationary Simulation Next, we c ompare the simulated I DF 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 hav e considerable impacts on results of exceedance probabilit y [ Porporato and 19 Ridolfi , 1998]. Therefore, we model ed nonstationarity considering time as a covariate in the location parameters of the GEV dist ribution. The uncertainty in design storm is shown using the vertical box plots. Tables S14 – S16 show parameters of the fitted GEV dist ributions for stationar y and nonstationary mod els. Table S 17 pre sents a statistica l perf orman ce of stationary versus nonstationary GEV mod els for the durations of 1- and 6-hour using Akaike Information Criteria with small sample correction [ AIC c ; Akaike , 1974; Hurvich and Tsai , 1990] . The skill of the individual model, as measure d by the AIC c statistics does not exhibit any consistent patterns either based on the locations or the storm durations. In general, we found CO RDEX RCMs able to sim ulate observed design storm intensity reasonably well for all three locations (Tables S 14 – S1 6, and Fi gures 10 – 12). This is in contrast with the results of trend estimates in Tables S7 – S11, where we noted the discrepancy in nature (si gn) of trends simulated b y the CORDEX RCMs. Further, T ables S14 – S 16 show the fitted GEV parameters for MM-Med CORDEX ensemble are close to that of the observed AM series. As note d from the plot s, the unc ertainty in design storm estimates is higher fo r the lowe r storm dura ti ons ( d = 1-, and 2 -hour) and hi gher return periods ( T = 25-, and 50- year) . For Toronto and Oshawa, at durations, d = 1 - to 6-hour, t he sim ulated rainfall intensity with higher return periods tend to underestimate the observ ed ones. How ever, for Windsor Airport, we find overe stimation of design intensity b y the models at d = 1-hour, whereas underestimation at d = 6-hour. Secondly, for Oshawa, we note an increase in observed rainfall intensit y in nonstationary models reference to stationary one at short durations of 1 -, and 2-hour and return periods, T = 10- year and be yond. This may be attributed to consistent upward trend in AM rainfall intensity across all durations in Oshawa. Further, we fin d 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 as ymmetric natur e of extremes. Our findings are consistent with an ea rlier stud y [ Cheng and AghaKouchak , 2014], in which authors r eported the diff erence b etween nonstationary and stationary r ainfall intensities are more prominent at short er durations and higher return periods, and the differences gradually diminish at long er durations. For detailed exploratory analysis, we present relative changes in ra infall intensity for 1 in 10-year ( i.e., T = 10-year) and 1 in 50 -year events respectively (Figure 13). W hile the top panel in the 20 figure (Figure 13- a) shows stationary models , the bottom one (Figure 13-b) shows the nonstationary ones, corresponding to design rainfall intensit y at T = 10-year and 50- year return periods, which c an be themselves a re 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 t end to increase at hi gher return period, T = 50- year. Few ex ceptions exist, for example, L ondon International Airport shows large n egative bias for 24 -hour rainfall ex treme. 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 ; Figure 14), in g en eral, the models simulate observed rainfall intensity reasonabl y well for T = 2 to 10-year return periods, in which observed intensity is we ll within the bound s of simulated intensity. However, the variabilit y is larger for T = 25-year and be yond. Exceptions exist for London, in which observed intensity lies outsi de the simulated bounds, even in the case of shorter return periods, although d ifferences are small ( Top panel; Figure 14). For nonstationary sim ulation ( Bottom panel ; Figure 14), the un certainty or spread of ensemble members are hi gher than that of the stationary mod els. However, except London Interna tional Airport and Stratford WWTP, in all cases, th e observed ex tremes lie within the bounds of simulated extremes a t all return periods. For a 24-hour extreme (Figure 15), ex cept in a few locations and higher return periods ( T = 25-year and be yond), the multi -model climate ensembles reasonably simulates the observed intensity. 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. 4. Projection of Future IDF 4.1 Projection of GEV Shape Param eter In this section, we discuss changes in shape parameter of sim ulated extreme in future relative to the baseline period. Figures 16 – 17 and Table S1 8 compares the shape parameters of the GEV distribution in the baseline (1970 – 2010) versus projected (2030 – 2070) scenarios . Barring few exceptions, in general, for both stationary (Figure 16) and nonstationar y (Figure 17 ) models in the near-future, the shape pa rameter shows a decrease in value, across most of the durations relative to the baseline period. Table S 18 shows during 2 050s, 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 21 events or a Weibull -type distribution), impl ying a tend ency to wards 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 b y -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) foll owed by -0.65 (at Toronto) to 0.33 (at Stratfo rd) with spatial standard deviation 0.09. These values in shape parameter are more ne gative than those conventionall y considered as physicall y re asonable 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 infini te variance as suggested by an earlier stud y [ Morrison and Smith , 2002] . A very negative value of the sha pe parameter sugges ts that the distribution of AM pre cipitation tends to have a very heav y tail. However, we found at short er durations (1 to 6-hour), the  values ranged between -0.33 and 0.23, and well within the physically consistent limit. 4.2 Future Changes in IDF Statistics Table S 19 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 Table S19, for 1-hour storm duration, stationary GEV model performs the best in six out of eight sites. However, for longer storm duration (in thi s case, 6 -hour), the spa tial pattern is not c onsistent. Fig ures 18 and 19 present projected intensity versus duration relations for 1 in 10- year and 1 in 50 -year events for nonstationary simulations. The associated IDF curves under stationary assumption are presented in Figures S9 – S10. The uncertainty in rainfall intensities at different dura tions are shown using boxplots. We find that for 1 in 10-year event (F igures 18 and S9), the MM-Min projected Intensity versus duration ( ID ) curves are close to that of the historical c urves, and often superimposing over the later for some stations (for example, Toronto, Oshawa, and London). On the other hand, th e MM -Max yields the upper bound of ID relationships. For 1 in 50 - ye ar events (Figures 19 and S 10) , the nature of ID curves follows similar trends, except for the Toronto and Stratford, in which we note widening of uncertainty envelop of RC M ensembles (MM-Min for S tratford and MM -Max for Toronto) in the projected scenario. 22 To further evaluate how the future intensitie s are sensitive towards the choi ce of f requency anal y sis ( i.e. , stationary/ nonstationary), the box plots of difference between future and hist orical precipitation extremes for 10-year and 50 -year events are presented in Figures 20 and 21. The absolute changes in future desi gn storm intensity r elative to observ ations are shown for three different cases: stationary (future) versus stationary (historical) ; nonstat ionary (future) versus stationary (historical) ; and nonstationary (future) versus nonstationary (historical) . For 1 in 10- year events, except Stratford, a t all loca tions, the quantile boxes are above zero for 6- and 24-hour durations. For Stratford, at 24- hour duration, the intensity of rainfa ll re mai ns relatively unchanged over time as indicated b y the location of quantile boxes, which is close to zero irrespective o f the choice o f the return peri od estimates . For hou rl y 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 ar e below z ero for all durations, indi cating a decrease in precipitation intensity in the future p eriod. However, we find an increase in r ainfall intensit y at 2 -hou r and 1- hour extremes for Windsor and Oshawa respect ively. For 24-hour precipitation extreme, the intensity of extreme pre cipitation remains relatively unchanged for 1 in 50-year events (Figure 19). Further, for both events, although the differences in rainfall intensit y are apparent for dur ations up to 2-hour th e y follow t he same trend irrespective of the c hoice of fr equency anal ysi s. However , the magnitude of differences fade away at durations of 6-hour and beyond. To quantify the magnitude of changes in projected intensit y relative to the baseline period, w e 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, and S 10). While the left panel of the figures show the z- statistics indicating statistical significance o f the change, the right panel s ummarizes the median changes in projected rainfall intensities for the future reference to the baseline period is expressed in ter ms of pe rcentages. A positive va lue indicates an increase relative to baseline period, whe reas negative ones indicat es a decrease. For 10- year return pe riod (Figure 22), under stationarit y assumption (Fi gure 22 ; top panel ), the statis tical ly signific ant increase in design storm is obse rved for sub-hourl y rainfall e xtremes in Hamilton and Oshawa and hourl y precipitation extreme in Toronto International Ai rport. On th e other hand, assuming non-stationarit y ( Figure 20; bottom panel ), a significant increase in rainfall intensit y is observed for 6 -hour extremes in Toronto, Oshawa and Windsor and 24-hour extreme in Toronto and Trenton respectivel y. A significant 23 decrease in intensity is observed for London. In particular, the projected increase ran ges from 1 ~ 40% and 1 ~ 38 % respectively, assumin g stationa rity and non-stationarit y methods of frequency analysis. For 25- year ret urn period (Figure S 11) , the performance is mixed, with a few stations exhibit an in crease (for instance, 6- and 24 -hour rainfall extreme in To ronto, 6 -hour extreme in Oshawa) while the others show a decrease ( i.e. , 12-hour rainfall extreme in Toronto and Trenton, 1- and 2-hour hour extremes in London and Fergus Shand dam re spectively). However, considering stationar it y assumption, none of these changes are signific ant (Figure S11 ; top panel ) in contrast to nonstation arity assumption (Figure S11 ; bottom panel ), in which magnitude of changes are d eemed to be significant. For 50- year return period (Figur e 2 3), assumin g non- stationarity, we found a statistically significant d ecrease in return p eriod for sub -hourl y rainfall extremes in Toronto, London, Trenton and Fergus Shand. However, the changes are not significant assuming stationarit y in r eturn period estimate. Considering stationarit y, th e percentage decrease in rainfall intensit y ranges between 2 – 37%, w hile 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- ye ar events, an increase in rainfall intensity is obse rved a cross most of the stations and durations, irrespective of the choice of frequency an alysis. (ii) For more extreme events ( i.e. , 1 in 25-year and 1 in 50-year), a d ecrease in rainfall intensity is not ed ; however, the magnitude of changes in storm intensit y are deemed to be statistically signifi cant considering non -stationarity assumption in contrast to the stationary assumptions. 5. Discussion and Conclusions Extreme precipitation poses significant risks for cities across Southern Ontari o, Canada due to their dense popul ation, valuable and geographically concentrated p roperty, complex and interdependent inf rastructure networks [ Henstra and T histlethwaite , 2017]. We characterize extreme rainfall frequency over Souther n Onta rio using I DF statistics with a suite of high- resolution new generation of coordinated regional climate modelin g experiments participating in the NA -CORDEX doma in with an elevated global greenhouse forcin g. Since the study area comprises urban catchm ents and most of the municipal infrastructures are desi gned based on design storm intensit y of 50- ye ar or less, the IDF statistics are compared for 2-, 5-, 10-, 25-, and 50-year return periods. The RC M performance in simulating extreme rainfall statisti cs in the 24 present da y climate is validated using ground-based AM observational records. The temporall y downscaled NA-CORDEX RCM ensembles are able to simulate the observed IDF statisti cs reasonably well as d emonstrated by various p erformance metrics. Fu rther, we modeled design rainfall intensity using extreme value distributions under stationarity and nonst ationarity assumptions with a time-varying location parameter. The following conclusions emerge from the present analysis:  The ground-based short-duration historical (1970 – 2010) precipitation e xtremes often exhibit a statis tically significant increase (or decrease) trends for the selected locations across Southern Ontario. On the other hand, extreme dail y temperature anomaly shows an upward trend for most of the locations , although the apparent trend is not statisti cally significant across the stations.  A robust bias -correction methodolog y is emplo yed to correct s ystematic bias es in the RC M 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 si gnificant increase in rainfall intensity is observed for 10- year events over most of the station locations. However, opposite tr end is noted for events with longer return periods. The pape r presents a pr oof of principal r esults to compare changes in ex treme precipitation intensity in the near-future using three regional climate models with t he future atmosphere represented b y 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 e t al. , 2014] ; still most of the g rou nd-based individual station observations do not exhibit a ny significant trend. Further, t ypically a limited number of statio ns show a statisticall y significant nonstationary trend a t a sub-dail y and dail y time sc ales [ Shephard et al. , 2014]. Recently, usin g station-based AM precipitation [ Shephard et al. , 2014] have shown that nationwide less than 5% of the stations exhibi ts significant increasing (or d ecreasing) trends. Based on their findings the aut hors recommend that in C anada, the traditional IDF design based on stationarity assumption of precipitation extreme has not been violated. F urther, i n projected scenario, we find an increase in rainfall intensity at 6- and 24-hour extremes for 10- year events 25 over majority of stations. Our results ar e consistent with [ Barbero et al ., 2017], where authors reported except in winter mont hs, there is a significant increase in dail y AM precipitation trends as compared to the hourly extremes after a nalyzing a large number of ground-based observational records across the United S tates for the pres ent da y (1950 -2011) condition. Their results indicat e that at a station level, the trends in daily precipitation extremes are better detected than that of the hourly extremes in a c hanging climate. They attrib uted this to the limi ted spatial extent of the most extreme events and the inabilit y o f sparse station network to accurately measure su ch ev ents [ Barbero et al ., 2017]. Some of the caveats of the study include: The heavy computational requirement involving successive spatiotempor al downscaling of ex treme precipitation at individual station level has limited us to use only three RCMs. How ever, un certaint y in climate cha nge projection results from both model responses and internal climate variabilit y [ Hawkins and Sutton , 2009, 2011]. A series of previous studies [ Deser et al. , 2012, 2014; Wettstein and Deser , 2014; Kay et al. , 2015] have examined the effect of intrinsic climate va riability and anthropogenic climate change on projected climate using a 40-member climate model ensembles available from National Center for Atmospheric R esearch ( NCAR)’s CESM Large Ensemble Co mmunit y Project (LENS). These studies ex plored various aspe cts of uncertaint y resulting from intrinsic climate variability and found substantial climate uncertaint y on a global and re gional scale. The role of internal climate variability is considerably higher at a regional sc ale 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- me mber RCMs may not be sufficient to capture the full spectrum of climate variability in a near-term plannin g period. Next, hourl y and sub-dail y extreme precipitation, required for engineering design, are often produced b y convective events. However, global and r egional climate models are not able to sim ulate such events well because of models’ li mited 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 u ncertainties in sim ulating extremes [ Hagos et al. , 2014; Singh and O’Gorman , 2014]. 26 Secondly, the cascade of uncertaint y ma y also stem from temporal downscaling of RCM output . Third, we a ssume that the selected bias-correction method a t historical period is time independe nt and same will holds good for correcting bias in the future period. If global warming is solely responsible for changes i n e xtreme pre cipitation trends, we might expect to see continued increase in extreme precipitation intensit y . On the other hand, if large-scale cli mate variability and/or anthropogenic influences are responsible, the extremes (for instance, J uly, 2013, extreme storm event caused insured losses of over $850 million [ Canadian Underwrit er , 2013] ) ma y re-occur in the future leading to wide scale damages to life and property. As a final caveat, inferring extreme precipitation fr equency s olely b ased on ground -based obse rvation has considerable uncertainty, at least as much as in the e stimat ed projected trends in the future emission scenario. Although we examined sources of uncertaint y in model projection b y comp aring the mul ti -model spread of future IDF curves (Figures 20 – 21; S 9 – S10) in R CM ensembles i n the future, there is a need to explore a complete suite of RCM outputs [ Mote et al. , 2016] archived at the NA -CORDEX domain in developing probability-based IDF cu rves cons idering a ran ge of plausible scenarios fo r risk- based resilient water management [ Hall and Borgomeo , 2013]. Acknowledgement The annual maximum rainfall data used in this study is downloaded from Environment Canada website: ftp://ftp.tor.ec.gc.ca/Pub/Engineering_Climate_Dataset/IDF/ . Ho urly and dail y rainfall data are obtained from T oronto and Region Conservation Authorit y (TRC A; https:/ /trca.ca/) and Environment Canada Historical Climate Da ta Archive ( http://climate.weather.gc.ca/ ). Th e authors would li ke to acknowledge financial support from the Natural Science and Engineering Research Council (N SERC) of Canada and the NSERC Cana dian FloodNet (Grant Number: NETGP 451456). The first author of the manuscript would like to thank Dr. Jonas Olsson of S wedish Meteorological and Hydrological Institute (SMH I), Norrköping, Sweden for sharin g M ATLAB - based random cascade disaggregation tools and implementation details through email. The nonstationary GEV anal yses were performed using MAT LAB-based toolb ox NEVA, available at the Universit y of C alifornia, I rvine website: http://amir.eng.uci.edu/neva.php (as acc essed on May 2016). 27 References Adam, J. C., and D. P. Lettenmaier (2003), Ad justment of global g ridded precip itation fo r systematic bias, J. Geophys. R es. Atmospher es , 108 (D9). Agilan, V., and N. V. Umam ahesh (2016), Is the covar iate based non- st ationary rainfall I DF curve capable of encom passi ng future ra infall changes?, J. Hydrol. , 541 , P art B , 1441 – 1 455, doi:10.1016/j. jhydrol.2016. 08.052. Agilan, V., and N. V. Umam ahesh ( 2017), What are the best covar iates for dev el oping non- st ationary rainfall Intensity- Duration -Frequency r elationship?, A dv. Water Resour. , 1 01 , 11 – 22, doi:10.1016/j.adv watres.2016.12.016. Akaike, H. (1974), A new look at the statistical m odel identification, IEEE Tra ns. Autom. Con t rol , 19 (6), 716 – 723. Canadian Underwr i ter (2013), Prel iminary insured loss estim at e for July G TA floods at $850 m illi on, setting new reco rd for Ontario, Avai lable at: http://www.cana dianunder writer.ca/insurance /preliminary-insured-l oss-estimate- for-july-gta- floods-at-850-milli on-setting-new- record-for-ontar io-1002529855 /, retrieved on: May , 2017. Ashfaq, M., L. C. Bow li ng, K. Che rkauer, J. S. Pa l, and N . S. Diffenbaug h (2010), I nfl uence of climate model biases and da i ly-scale tempera ture and precip itation events on hydrolog ical im pacts assessment: A case st udy of the United States, J. G eophys. Res. Atmo spheres 198 4 – 2012 , 115 (D14). Baldwin, D. J., J. R. Deslog es, and L. E. Band (2011), 2 Phy si cal Geography of O ntario, Ecol. Manag. Terr. Landsc. Pat terns Process. For. Lan dsc. Ont. , 12. Bao, J., S. C. Sherwoo d, L. V. Alexander, and J. P. Evans (2017), Future increases in ex treme precipitation ex ceed observed sc aling rates, Nat. Cl im. Change , 7 (2), 128 – 132. Barbero, R., H. J. Fowler, G . Lenderink, and S. B l enkinsop (2017), I s the intensification of precip itation extremes with glob al warm ing better detected at hou rly than daily res olutions?, Geophys. Res. Lett. , 44 (2), 2016GL07191 7, doi:10.1002/2016GL071917. Bonnin, G. (2013 ), Precipitati on Frequency Estimates for the Nation and Extremes-A Perspective., Bonnin, G. M., D. M artin, B . Lin, T. Parzy bok, M. Yekta, and D. Riley (2006), P recipitation-frequency atlas of the United S tates, NO AA Atlas , 14 (2). Bourne, L. S., and J. Simmons (2 003), New fault l i nes? Recent trends i n the Canadian urba n system and their implicatio ns for planning and public policy , Can. J. Urban Res. , 12 (1) , 22. Burn, D. H., and A . Taleghani (2013), Estim ates of chang es in design rain fall values for Can ada, Hydrol. Process. , 27 (11), 1590 – 1599. Bush, E. J., J. W. Lode r, T. S. James, L. D. Mo rtsch, and S. J . Cohen (2014), An overview of Canad a’s changing clim ate, Can. Chang. Cl i m. Sect. Perspe ct. Impacts A dapt. FJ W arren Lemmen , 23 – 64. Chandra, R., U. Saha , and P. P. Mu jumdar (2015), Model and param eter uncertain ty in ID F relationships under clim ate change, Adv. W ater Resour. , 79 , 127 – 139. Cheng, L., and A. A ghaKouchak (2014), Nonsta tionary precipitat ion intensity- duration-frequency curv es for infrastructure design in a changing clim at e, Sci. Rep. , 4 . Cheng, L., A. A ghaKouchak , E. Gilleland, and R. W. K atz (2014a), Non- stationary extreme v al ue analysis in a ch anging clim ate, Clim. Change , 127 (2), 353 – 369. Christensen, J. H., and O. B. Christensen (2007 ), A summ ar y of the PRUDENC E model projec tions of changes in European climate by the end of this c entury, Clim. Change , 81 , 7 – 30. Coles, S., J. Bawa, L. T renner, and P. D orazio (2001) , An introduction to statistic al modeling of extreme values , Springer. 28 Coles, S. G., and J. A. Tawn (1996 ), A Bayesian an al ysis of extrem e rai nfall dat a, Appl. Stat. , 463 – 478. Coulibaly, P., and X. Shi (2005), I dentification of the effect of clim ate change on future desig n standards of drainage infras tructure in O ntario, Rep. Prep. Mc Master Univ. Fund ing Minist. Transp. Ont. , 82 . WMO (2009), G uidelines o n analysis of extremes in a chang i ng climate in s upport of info r med decisions for adaptation, W orld M eteorol. Organ. Deng, Z ., X. Qiu, J. Liu, N. Madras, X. Wang , and H. Z hu ( 2016), Trend in frequency of ex t reme precipitation ev ents over O ntario from ensem bl es of multipl e GCMs, Clim. Dyn. , 46 (9 – 10), 2909 – 2921, doi:10. 1007/s00382- 015 -2740-9. Deser, C., A. S. Phi llips, R. A. Tom as , Y. M. Okum ura, M . A. Alexand er, A. Capo tondi, J. D. Scot t, Y.- O. Kwon, and M. Oh ba (2012), EN SO and Pac ific decadal var iability in the Comm unity Climate System Model v ersion 4, J. Clim. , 25 (8), 2622 – 2651. Deser, C., A. S. Phi llips, M. A . Alexander, and B. V . Sm oliak (2014), Projecting N orth American c limate over the next 50 y ears: Uncertain ty due to inte rnal variability , J. Clim. , 27 (6), 2271 – 2296. Dibike, Y. B., Coulibaly, P. (2006). Te mporal neural networks for downscalin g climate variability and extremes. Special Is sue of Neural N et works , 19(2), 135- 144. Diffenbaugh, N. S. e t al. (2017), Q uantifying the influence o f global warm ing on unprecedente d extreme climate events, P roc. Natl. Acad. Sci. , 20 1618082, do i:10.1073/pnas.16180 82114. El Adlouni, S., B. B obée, and T. B. M. J. Ouarda ( 2008), On the t ails of extrem e event distributions in hydrology , J. Hydrol. , 355 (1 – 4), 16 – 33, doi:10.101 6/j.jhy drol.2008.02.011. Elshorbagy , A., A. Naz emi, and S. Alam (2015), Analyzing th e v ariations in intensity – duration – f requenc y (IDF) curves in t he city of S askatoon under climate change , Citeseer . Environm ent Canada (2014 ), Documentation on Env i ronment Canada Ra infall I ntensity- DurationFrequen cy (IDF) Tables and Gra phs. V2.30. Available at: ftp://ftp.tor.ec.gc. ca/Pub/Engine ering_Climate_Dat aset/IDF/ . Feser, F., B. Rock el, H. von Stor ch, J. Winter feldt, and M. Z ahn ( 2011), Reg i onal climate m odel s add value to global m odel data: a review and se l ected exam ples, Bull. Am. Meteoro l . Soc. , 92 (9), 1181 – 1192. Fischer, E. M., and R . Knut ti (2016), Observ ed heavy prec ipitation incr ease confirm s theor y and early models, Nat. Clim. C hange , 6 (11), 986 – 991. Ganguli, P., and A. R . Gang uly (2016a), Robustness o f Meteorological Droughts in Dynam ically Downscaled Clim ate Si mulations, JAW RA J. Am. Water Resour. Assoc. , 52 (1), 138 – 167, doi:10.1111/1752- 1688.12374. Ganguli, P., and A. R . Gang uly (2016b), Robustn ess of Meteo rological Droug ht s in Dynam ic al ly Downscaled Clim ate Si mulations, JAW RA J. Am. Water Resour. Assoc. , 52 (1), 138 – 167, doi:10.1111/1752- 1688.12374. Ganguly, A . R., D. Kum ar, P. Ganguli, G. Short, and J. Klausn er (2015), Cl imate Adapta tion Inform ati cs: Water Stress on Pow er Product ion, Comput. Sci. Eng. , 17 (6), 53 – 60. Gao, Y., J. S. Fu, J. B. D rak e, Y. Liu, and J.-F. Lamarque (2012), Pro j ected chang es of extreme weather events in the east ern United States based on a hig h resolution cl imate modeling system, Environ. Res. Lett. , 7 (4), 044 025, doi:10.1088 /1748-9326/7/4/0 44025. Gilbert, R. O. (198 7), Statistical me thods for envir onmental po l lution monitoring , John Wiley & Sons. Giorgi, F. (2006 ), Regional clim at e modeling : Status and perspec tives, in Jour nal de Phys ique IV (Proceedings) , v ol. 139, pp. 101 – 118, EDP sciences. Giorgi, F., and L. O . Mearn s (1991), Approaches to t he simulation o f regional clim ate change: a review, Rev. Geophys. , 29 (2), 191 – 216. 29 Glotter, M., J. El liott, D. McI nerney, N. Best, I. Foster, and E. J. Moyer (2014), Evaluating t he utility of dynamical downsca ling in agricultural im pacts project ions, Proc. Natl. Acad. Sci. , 111 (24), 8776 – 8781. Güntner, A., J. Olsson, A. C alver, and B. Gannon ( 2001), Cascad e-based disagg regation of continuous rainfall time se ries: the influence of c limate, Hydrol. E art h Syst. Sci. Di scuss. , 5 (2), 145 – 1 64. Gutmann, E. D ., R. M. Rasm uss en, C. Liu, K. Ik eda, D. J . Gochis, M. P. Clark, J. Dudhia, and G. Thompson (2012 ), A comparison o f statistical and dy nam i cal downscaling of w i nter precipitation over complex terr ain, J. Clim. , 25 (1) , 262 – 281. Gutowski Jr, W. J., S. G. Deck er , R. A. Donavon, Z . Pan, R. W. Arritt, and E. S. T akle (2003), Temporal – spatia l scales of observed a nd simulated precipitation in cen t ral US cl im ate, J. Clim. , 16 (22), 3841 – 3847. Hagos, S., Z . Feng, C. D . Burleyson, K.-S. S. Lim , C. N. Long , D. Wu, and G. Th ompson (2014), Evaluation of conv ection-perm itting model sim ul ations of cloud popu l ations asso ciated with the Madden- Julian Oscillation using data collec ted during the AMIE/DY NAMO field campaign, J. Geophys. Res. A tmospheres , 119 (21 ), 2014JD022 143, doi:10.1002 /2014JD02214 3. Hall, J., and E. Borg omeo (2013), Risk-based principle s for defining and manag ing water security, Philos. Transact. A Ma th. Phys. Eng. Sci. , 371 (2002), doi:1 0.1098/rsta.2012. 0407. Hall, J. W., D. G rey, D. Garrick , F. Fung, C. Brow n, S. J. Dadson, and C. W. Sadoff (2014), Cop i ng with the curse of freshwa t er variabili ty, Science , 346 (6208), 429 – 430. Hay, L. E., and M. P. C l ark (2003 ), Use of stati stically and dy namically downscaled atm ospheric model output for hydro logic sim ulations in three m ount ainous basins in the western U nited States, J. Hydrol. , 282 (1), 56 – 75. Halmstad, A., M. R. Najafi, and H. Moradkhani (2 013), Analysis of p r ecipitation e xtremes with the assessment of reg ional climate m odel s over the W i llamette Riv er Basin, USA , Hydrol. Process. , 27 (18), 2579 – 2590. Hamed, K. H., and A. R . Rao (1998), A modified Mann-Kendall trend test for autocor related data, J. Hydrol. , 204 (1 – 4), 182 – 196. Hassanzadeh, E., A . Nazemi, and A. E l shorbagy (2013), Quant ile-based downscali ng of precipitat ion using genetic prog ramming: application to I DF curves in Saskatoon, J. Hydrol. Eng. , 19 (5), 943 – 955. Hatfield, J. L., K. J. Boote, B . A. Kimball, L. H. Z iska, R. C. Izaurralde, D. Or t, A. M. Thom son, a nd D. Wolfe (2011), Cl imate impacts on agricultu re: implications for c rop production, Agron. J. , 103 (2), 351 – 370. Hawkins, E., and R. Su t ton (2009 ), The potential to nar row uncertainty in regional climate predic tions, Bull. Am. Meteo rol. Soc. , 90 (8), 1095 – 1107. Hawkins, E., and R. Su t ton (2011 ), The potentia l to narrow uncertain ty in pro jections of regiona l precipitation ch ange, Clim. Dyn. , 37 (1 – 2) , 407 – 418. Hawkins, E., and R. Su t ton (2012 ), Time of emergence of climate sig nals, Geophys. Res. L ett. , 39 (1), L01702, doi:10.1029 /2011G L050087. Hayhoe, K. et al. ( 2004), Em issi ons pathways, c limate change, and i mpacts on C alifornia, Proc. Natl. Acad. Sci. U. S. A. , 101 (34), 12422 – 12 427. Hayhoe, K., C. Wak e, B. Anderson, X.-Z. Liang, E. Maurer, J . Zhu, J. Bradbu ry, A. DeGaetano, A . M. Stoner, and D. Wue bbles (2008), R egional cl imate chang e pr ojections for th e Northeas t USA, Mitig. Adapt. Strateg. G lob. Change , 13 (5 – 6), 425 – 436. 30 Hayhoe, K., J. VanD orn, T. C roley, N. Schleg al, and D. Wuebb l es (2010), R egional Clim ate Change Projections for C hicago an d the US Great Lak es , J. Gt. Lakes Res. , 36 (sp2), 7 – 21, doi:10.1016/j. jglr.2010.03. 012. Henstra, D., and J . Thistlethwaite (2017), C limate cha nge, floods and mun icipal risk sharing in Canad a , Institute on Muni cipal Finance an d Governance, Univ ersity of Toronto. Hurvich, C. M., and C.- L. Tsai (1990), The im pac t of m odel selection on infe rence in linear reg ression, Am. Stat. , 44 (3), 2 14 – 217. Ines, A. V. M., and J . W. Hansen (2006), Bias correction of dai ly GCM rainfa ll for crop sim ulat ion studies, Agric. For. M et eorol. , 138 (1 – 4 ), 44 – 53, doi:10.1016/j.agrf ormet.2006. 03.009. Jalbert, J., A.- C. Favre, C. Bélisle, and J. -F. Angers (2017), A sp atiotemporal model for extreme precipitation s imulated by a climate model, w ith an app lication to asse ssing changes in return levels over Nort h America, J. R. Stat. S oc. Ser. C App l. Stat. Jien, J. Y., and W. A. G ough (2013 ) , The influenc e of A tlantic hurrican es on Southern On tario’s precipitation ex tremes, Theor. App l . Climatol. , 114 (1 – 2), 55 – 60. Kay, J. E. et al. ( 2015), The C ommunity Earth System Mode l (CESM) large en semble project: A community resource for study ing climate change in the presence of intern al clim ate variability, Bull. Am. Meteo rol. Soc. , 96 (8), 1333 – 1349. Kendon, E. J., N. M. Rob erts, H. J. Fowler, M. J. Roberts, S. C. C han, and C. A. S eni or (2014), H eavier summer downpours w i th clim ate change revealed by weather forecast resolu tion model, Nat. Clim. Change , 4 (7), 5 70 – 576, doi:10.1 038/nclim ate2258. Kharin, V. V., F. W. Z wiers, X. Zhang, and G. C. H egerl (2007), Cha nges in temperature and precipitation ex tremes in the I PCC ensemble of global coupled m odel simulations, J. Cli m. , 20 (8), 1419 – 1444. Kling, G. W. et a l. (2003), C onfronting clim ate change in the Great Lakes region : impacts on our communities and ecosyste ms, Union Concerned Sc i. Camb. Mas s. Ecol. Soc. Am . Wash. DC , 92. Kumar, D., E. Kodra, and A. R. Ganguly ( 2014), Regional and s easonal in tercomparison o f CMIP3 and CMIP5 clim at e model ensem bl es for temperature and precipi tation, Clim. Dyn. , 43 (9 – 10), 2491 – 2518, doi:10.1007 /s00382-014- 2070-3. Kunkel, K. E. (200 3) , North Am erica n Trends in Ex trem e Prec ipitation, Nat. Hazards , 29 (2), 291 – 3 05, doi:10.1023/A:1023694 115864. Kunkel, K. E., K. An dsager, and D. R. Easterling (199 9), Long-term t rends in extrem e precipitation events over the c ontermino us United States and Canada, J. Clim. , 12 ( 8), 2515 – 2527. Kuo, C.-C., T. Y. Gan, and M. Gizaw (2015), Potenti al impact of cl imate chang e on intensity durat ion frequency curv es of central Alberta, Clim. Change , 13 0 (2), 115 – 129. Lapen, D. R., and H. N . Hay hoe ( 2003), Spatial Ana lysis of Season al and Annu al Temperature and Precipitation Norm als in Southern Ontario, Can ada, J. Gt. Lakes Res. , 29 (4), 529 – 544, doi:10.1016/S0380- 1330(03)70457 -2. Lenderink, G., and H. J. Fowler (2017), Hydroclimate: Understanding rainfall extremes, Nat. Clim. Change , 7 (6), 391 – 393, do i:10.1038/nclim at e3305. Lehmann, J., D. C oumou, and K. F rieler (2015), I ncreased record-breaking precipitation ev ent s under global warm i ng, Clim. Change , 132 (4), 501 – 5 15, doi:10.1007/s10584-015- 1434 - y. Li, H., J. Sheffield, and E. F. Wood (2010), Bias cor r ection of month l y precipitation and t emperature fields from Intergovernm ental Panel on Clim at e Change AR 4 model s using equidistant qu antile matching, J. Geophy s. Res. Atmosphe re s , 115 (D10). 31 Li, J., J. Evans, F. Johns on, and A. S harm a (2017), A com pari son of methods for estimating cli mate change impact on design rainfa ll using a hig h-resolution RCM, J. Hydrol. , 547 , 413 – 427. Ligeti, E., I. Wie ditz, and J. Penney ( 2006), A scan of clima t e change impac ts on Toront o , Clean Air Partnership. Lima, C. H., H .-H. Kwon, and J.-Y. Kim (2016), A B ayesian beta di stribution m odel for estim at ing rainfall ID F curves in a cha nging cli mate, J. Hydrol. , 5 40 , 744 – 756. Madsen, H., K. Arn bj erg-Nielsen, and P. S. Mi kkelsen (2009), Update o f regional intensity – duration – frequency curv es in Denmark: tendency towa r ds increased s torm intensities, Atmospheric Res. , 92 (3), 343 – 349. Mailhot, A., and S. D uchesne (20 09), Design cri teria of urban d rainage infrast ructures unde r climate change, J. Water Re sour. P lan. Manag. , 136 (2), 201 – 208. Mailhot, A., S. Duc hesne, D . Caya, and G. Talb ot (2007), Ass essment of futu r e change in int ensity – duration – frequency ( IDF) curv es for Southern Q uebec using the Canadian R egional Cl imate Model (CRCM), J. Hy drol. , 347 (1 ), 197 – 210. Mailhot, A., I . Beauregard, G. Talbot, D. Caya, and S. Biner (2012), Fu ture chang es in intense precipitation ov er Canada assessed from multi-model NARCCA P ensemble simulat ions, Int. J. Climatol. , 32 (8), 115 1 – 1163, doi:10.1002/ joc.2343. Maraun, D. (2016 ), Bias Corre cting Clim ate Change Sim ulat ions - a Critical Review, C urr. Clim. Chang e Rep. , 2 (4), 211 – 220, doi : 10.1007/s40641- 016-0050-x. Markose, S., and A. A lentorn (2011), The gener alized extrem e value distribution, im pl ied tail index, and option pricing, J. Der i v. , 18 (3), 35 – 60. Martins, E. S., J. R. S tedinger, an d others (2000), G eneralized maximum- likelihood generalized ex treme- value quantile es t imators for hy drologic data, Water Re sour. Res. , 36 (3), 737 – 7 44. Martre, P. et al. (2 015), Multim odel ensembles of wh eat growth : many models are better than one, Glob. Change Biol. , 21 (2 ), 911 – 925, doi:10.1 111/ gcb.12768. Martynov, A ., R. Laprise, L. S ushama, K. Wing er, L. Šeparov ić, and B. Dug as (2013), Re analysis -driven climate sim ulation over COR DEX North Am eri ca domain using th e Canadian Regional Clim at e Model, version 5 : model perfo rmance evalua tion, Clim. Dyn. , 41 (11 – 12), 297 3 – 3005. McGinnis, S., D. N ychka, and L. O . Mearns (2015), A new distribu tion m apping technique for cl imate model bias correc tion, in Machine learn i ng and data m ining approaches t o climat e science , pp. 91 – 99, Springer. McNie, E. C. (2007), Reconciling the supply of scie ntific inform ati on with us er demands: an analysis of the problem and review of the lite rature, Environ. Sc i. Policy , 10 (1), 17 – 38 . Mearns, L. O. et al. (2012), The North Am er ican regional clim at e change asses sment program : overview of phase I results, B ull. Am. Me teorol. Soc. , 93 (9 ), 1337 – 1362. Meehl, G. A. et a l. (2009), Decadal predict ion: can it b e skillful?, Bull. Am. Meteo rol. Soc. , 90 (10), 1467 – 1485. Meinshausen, M. e t al. (2011), T he RCP g reenhouse gas concen t rations and their extensions fr om 1765 to 2300, Clim. Chang e , 109 (1 – 2), 213 . Mikkelsen, P. S., H . Madsen, K. A rnbjerg-N ielsen, D. Rosbjerg , and P. Harremoë s ( 2005), Select ion of regional historica l rainfall t ime series as in put to urban dra inage simulatio ns at ung auged locations, Atmosph eric Res. , 77 (1 – 4 ), 4 – 17, doi:10.1016/j.atm osr es.2004.10.016. Mills, E. (2005), I nsur ance in a C limate of Change, Sci ence , 309 (5737), 1040 – 1 044, doi:10.1126/science.1 112121. 32 Milly, P. C. D., J . Betancourt, M. Fa l kenm ark, R. M. Hirsch, Z . W. Kundzewicz, D . P. Lettenmaier, and R. J. Stouffer (2008 ) , Stationarity Is Dead: Whithe r Water Man agem ent ?, Science , 319 (5863), 573 – 574, doi:10.11 26/science.1151 915. Milly, P. C. D., J . Betancourt, M. Fa l kenm ark, R. M. Hirsch, Z . W. Kundzewicz, D. P. Le ttenmaier, R. J. Stouffer, M. D. De ttinger, a nd V. Krysanov a (2015), On C ritiques of “Stationarity is Dead: Whither Water M anagem ent?,” Water Resour. Res . , 51 (9), 7785 – 7789, doi:10.1002/2015WR017 408. Min, Y.-M., V. N. Kryjov, and S. M. Oh (2014), Assessm ent of APCC m ultimodel ensem ble prediction in seasonal clim ate f orecasting : Retrospective (198 3 – 2003) and real-tim e forecasts (200 8 – 2013), J. Geophys. Res. Atmos pheres , 119 (21), 2014JD02223 0, doi:10.1002/2014JD022230. Mirhosseini, G., P. S r ivastav a, and L. Stefanova (2013), The im pact of clim ate chang e on rainfall Intensity – Duration – Fre quency (I DF) curves in Alab am a, Reg. Environ. Change , 13 (1) , 25 – 33, doi:10.1007/s10113- 012-0375-5. Morrison, J. E., and J. A. Sm it h (2002), Stochastic m odeling of flood peaks using t he generalized ex t reme value distribution, W ater Resour. Res. , 38 (12), 1305, d oi :10.1029/2001WR0005 02. Mote, P. W., M. R. A l len, R. G . Jones, S. Li, R. Me r a, D. E. Rupp, A . Salahuddin , and D. Vickers (2016), Superensemble reg i onal clim at e modeling for the wes tern United Sta tes, Bull. Am. Meteor ol. Soc. , 97 (2), 203 – 21 5. Murray, V., and K. L. Ebi (2012), IPCC special repor t on managing t he risks of e xtreme events and disasters to advance climate change adaptation (SREX) , BMJ Publishing Group Ltd. Naz, B. S., S.- C. Kao, M. A shfaq, D. Rastog i , R. Mei, and L. C. B owling (2016 ), Regional hy drologic response to clim ate change in the conterm inous United States using hig h-resolution hydroclim at e simulations, Glob. P lanet. Change , 143 , 100 – 117. New, M., M. Todd, M. Hulm e, a nd P. Jones (2001), P recipitation m easurements and trends in the twentieth century , Int. J. Climatol. , 21 (15 ), 1889 – 1922. Olsson, J. (1995 ), Limits and cha r acteristics o f the multifrac t al behaviour of a high- resol ution rainfall time series, Nonlinea r Process. Geo phys. , 2 (1), 23 – 29. Olsson, J. (1998 ), Evaluation of a sca ling cascade model for tem poral rain-fall disagg r egation, Hydrol. Earth Syst. Sci. Di scuss. , 2 (1), 19 – 3 0. Pfahl, S., P. A. O’Gorman, and E. M. Fischer (2017), Understanding the regiona l pattern of pr ojected future changes in extrem e pr ecipitation, Nat. Clim. Ch ange , 7 (6), 4 23 – 427, doi:10. 1038/nclim ate3287. Palutikof, J. P., B . B. Brabson, D . H. Lister, and S . T. Adcock (1999), A review of m ethods to calculate extreme wind spe eds, Meteorol. Appl. , 6 (2), 119 – 132, doi:10.1017/S135048279 9001103. Panofsky, H. A ., and G. W. Brier (1 965), Some appl ication of statist ics to meteoro l ogy, The Pennsylvania State University , University Park, PA, 224. Partridge, M., M. R. Olfert, and A. Alasia (2007), C anadian cities as reg ional engines o f growth: agglom eration and amenitie s, Can. J. Econ. Can. D économique , 40 (1 ), 39 – 68. Perkins, S. E., A. J . Pitman, N. J. Holb r ook, and J. Mc Aneney (2007 ), Evaluation of the AR4 Cl imate Models’ Simulat ed Daily Maxim um Tem perature, Minim um Tem perature, and Precipitatio n over Australia Using Probability Density Functions, J. Clim. , 20 (17), 43 56 – 4376, doi:10.1175/JCLI 4253.1. Porporato, A., and L. Ridolfi (1998), Influence of weak t rends on exceedance probability , Stoch. Hydrol. Hydraul. , 12 (1), 1 – 14. Prein, A. F., R. M. Rasmussen, K. Ik eda, C. Liu, M. P. Clark, and G. J. Holland (2016), The future intensification of hourly pr ecipitation ext remes, Nat. Clim. Chang e . 33 Prodanovic, P., an d S. P. Sim ono vic (2007), Development o f rainfall intensity duration freque ncy curves for the City of Lo ndon under the changing climat e , Department of Civ il and Env i ronmental Engineering, Th e Universit y of Western Onta rio. Ragulina, G., and T . Reitan (2017 ), Generalized e xtrem e value shape param eter and its natu re for extrem e precipitation us ing long time ser ies and the Bay esi an approach, Hydrol. Sci. J. , 62 (6), 863 – 879, doi:10.1080/02626667. 2016.1260134. Rana, A., L. Bengts son, J. O lsson, and V. Joth iprakash (2013 ) , Developm ent of IDF-curves for tropical india by random cascade model ing, Hydrol. Earth Sy st. Sci. Discuss. , 10 (4 ), 4709 – 4738. Rahmstorf, S., G. F oster, and N . Cahill (2017), G lobal tem per ature evolution : recent trends a nd some pitfalls, Environ. R es. Lett. , 12 (5), 054001, doi:1 0.1088/1748-9326/a a6825. Reddy, M. J., and P. G anguli (2013), Spatio- t emporal analysis and d erivation of copula-based intens ity – area – frequency curves for droug hts in western Ra j asthan (I ndia), Stoch. Env iron. Res. Ri sk Assess. , 27 (8), 1975 – 1 989. Rodríguez, R., X. Nav arro, M. C. Casas, J. Ribalay gua, B. Russo, L. Po uget, and A. Redaño (2014), Influence of c limate change on I DF curves for the metropolitan ar ea of Barce lona (Spain), Int. J. Climatol. , 34 (3), 643 – 654. Rosenberg, E. A ., P. W. Key s, D. B. Booth, D. Har tley, J. Burkey , A. C. Steinemann, and D. P. Lettenmaier (2010 ), Precipitation extremes and the impacts of clim ate change on storm water infrastructure in Washingto n State, Clim. Change , 102 (1), 319 – 349. Rosenzweig, C ., A. Ig lesias, X. B. Yang, P. R. Eps tein, and E. Ch i vian (2001), C limate change and extreme weather ev ents; implication s for food p roduction, plant dis eases, and pe sts, Glob. Change Hum. Hea lth , 2 (2), 90 – 104. von Salzen, K. et al. (2013) , The Canadian fourth g eneration atmospheric glob al clim ate model (CanAM4). Part I : r epresentation o f physica l processes, Atmosphere-Ocean , 51 ( 1), 104 – 125. Schoetter, R., P. Ho f fmann, D . Rechid, and K. H. Sc hlünzen (2012), Evaluation a nd Bias Correc tion of Regional Clim ate Model R esults Using Mode l Evaluation M easures, J. Appl. Met eorol. Climato l. , 51 (9), 1670 – 1684, d oi:10.1175/JA MC-D-11-0161.1. Scinocca, J. F., V. V. Kharin, Y . Jiao, M. W. Qian, M. Lazare, L. Solheim , G. M. Flato, S. Biner, M. Desgagne, and B . Dugas (2016), Coo rdinated globa l and regional c limate m odel ing, J. Clim. , 29 (1), 17 – 35. Sen, P. K. (1968), Estimate s of the regression coefficie nt based on Kenda ll’s tau, J. Am. Stat. Assoc. , 63 (324), 1379 – 1389. Šeparović, L., A. A lexandru, R . Laprise, A. Ma rt ynov, L. Susham a, K. Winger, K. T ete, and M. V alin (2013), Present c limate and climate change ov er North Am er ica as simulated by t he fifth- generation Canad i an regional c limate model, Clim. Dy n. , 41 (11 – 12), 3167 – 3201. Shephard, M. W., E. Mekis, R. J . Morris, Y. Feng, X. Z hang, K. Kilcup, and R. Fleetwood (2014), Trends in Canadian Short‐Du ration Extreme Rainfa ll: Including an Intensity– Duration – Frequency Perspective, Atmos phere-Ocean , 52 (5 ) , 398 – 417, doi:1 0.1080/07055900.2014.969677. Sillmann, J., V. V . Kharin, X. Zhang , F. W. Zwiers, and D. Bronaug h (2013a), Cl imate extrem es i ndices in the CMIP5 m ultimodel ensemble: Part 1. Mo del evaluation in the presen t climate, J. Geophys. Res. Atmospheres , 118 (4 ) , 1716 – 1733. Sillmann, J., V. V . Kharin, F. W. Z wier s, X. Zhang, an d D. Bronaugh (2 013b), Clim ate extremes indices in the CMIP5 m ultimodel ensemble: Part 2. Fu t ure clim ate pr ojections, J. Geophys. Res. Atmospheres , 118 (6), 2473 – 2493. 34 Simonov ic, S. P., A. Schardong , and D. Sandink (2016), Mapping ext reme rainfa ll statistics for Canad a under clim ate change using updated I ntensity-Duration-Frequency curv es , J. Water Resou r. Plan. Manag. , 04016078. Singh, D., M. Ts iang, B. Rajaratnam , and N. S. Diffenbaugh (2013), Pr ecipitation ext remes over the continental United S t ates in a transien t, high-resolution , ensemble climate model e xperiment, J. Geophys. Res. A tmospheres , 118 (13 ), 7063 – 7086. Singh, M. S., and P. A . O’Gorm an ( 2014), Influence of microphysics on the scaling of precipitation extremes with tem perature, Geophys. Res. Lett. , 41 (16) , 6037 – 6044, do i: 10.1002/2 014GL061222. Srivastav, A. S. R ., and S. P. Simonovic (2014), Gene ralized Too l for Updating i ntensity-Duration- Frequency Curv es under Clim at e Change, Srivastav, R. K., A. Sch ar dong, and S. P. S imonov ic (2014), Equidistance Quan t ile Matching Me thod for Updating I DFCurves under C limate Change, W ater Resour. Manag. , 28 (9), 2 539 – 2562. Stocker, T. F., D . Qin, G. K. Pla ttner, M. Tignor, S. K. Allen, J. Boschung, A . Nauels, Y. Xia, B. Bex, and B. M. Midg ley (2013), IPCC , 2013: climate cha nge 2013: the physical scienc e basis. Contribution of wor king group I to the fifth assessmen t report of the intergovernm ent al panel o n climate change , Cam bridge University Press. Stoner, A. M., K. H ayhoe, X. Yang , and D. J. Wuebbles (2013), An asynchronou s regional reg ression model for statist ical downs caling of daily climate variable s, Int. J. Climatol. , 33 (11), 247 3 – 2494. Switzm an, H., T. Razavi, S. Traore, P. C oulibaly, D . Burn, J. H enderson, E. Faus to, and R. Ness ( 2017), Variability of Fut ure Extrem e Rainfall Statistics: A Com par ison of Multiple I DF Proj ections, J. Hydrol. Taylor, K. E. (2001 ), Summarizing m ultiple aspects of m odel perform ance in a single diagram , J. Geophys. Res. A tmospheres , 106 (D7 ) , 7183 – 7192. Taylor, K. E., R. J . Stouffer, and G. A . Meehl (2012), An overview of CMIP5 and the expe r iment design, Bull. Am. Meteo rol. Soc. , 93 (4), 485 – 498. Teutschbein, C., and J. Seibert (20 10), Regional climate m odel s for hydrolog ical im pact st udies at the catchment scale : a review of re cent modeling strategies, Geogr. Compa ss , 4 (7), 834 – 860. Teutschbein, C., and J. Seibert (20 12), Bias correc tion of reg i onal climate m odel simulations for hydrological clim at e-change im pact studies: Rev iew and evaluatio n of differen t m ethods, J. Hydrol. , 456 , 12 – 29. Van der Linden, P., a nd Jfb. Mitchell (2009), ENSE MBLES: Clim at e Change and i ts Im pacts: Summ ary of research and results from the ENSEMBLES project, Met Off. Hadl ey Cent. FitzRoy Road Exeter EX1 3PB UK , 1 60 . Vousdoukas, M. I ., L. Mentaschi, E. Voukouv alas, M. Verlaan, a nd L. Feyen (2017), Extrem e sea levels on the rise along Euro pe’s coasts, Earths Future . Vu, M. T., V. S. Rag havan, and S.-Y. Liong (2017), D er iving short- duration rainfall I DF curves from a regional clim ate model, Nat. H azards , 85 (3), 1877 – 189 1. Wang, J., and V. R. Kotamarthi ( 2015), High- resolution dynam i cally downscaled projections of precipitation in th e mid and late 21st century over North Am eri ca, Earths Future , 3 (7), 268 – 28 8. Wang, L., and W. Chen (2014), A CMIP5 m ul timodel projection of future tem perature, precipita tion, and climatologica l drought in C hina, Int. J. Climatol. , 34 (6), 2059 – 2078. Wang, X., G. Hu ang, and J. Liu (2 014), Projected increases in intensity and f requency of rainfall extremes throug h a r egional clim ate modeling approach, J. Geophys. R es . Atmospheres , 119 ( 23). Watts, G. et al. (2015), Clim ate change and water in the UK – past chang es and future prospects, Prog. Phys. Geogr. , 39 (1), 6 – 28. 35 Wehner, M. F. (2013) , Very extreme seasonal preci pitation in th e NARCCAP ensemble: m odel performance an d projection s, Clim. Dyn. , 40 (1 – 2 ), 59 – 80. Wettstein, J. J., an d C. Deser (2014), Internal v ar iability in projections of twen t y -first- century Arctic sea ice loss: Role o f the large-scale atm ospheric circulation, J. Clim. , 27 (2), 527 – 550. Whan, K., and F. Z wiers (2016), Evaluation of ex treme rainfall and temperature o ver North America in CanRCM4 and CRC M5, Clim. Dyn. , 46 (11 – 1 2), 3821 – 3843. Wilby, R. L., and S. De ss ai (2010 ), Robust adap tation to climate c hange, Weather , 65 (7), 180 – 185. Wilby, R. L., and T . M. L. Wigley (1997), Downs caling general circulation m odel output: a review of methods and lim i tations, Prog. Phys. Geogr. , 21 (4), 530 – 548. Wilby, R. L., T. M. L. Wi gley , D. Conway , P. D. Jones, B. C. H ewitson, J. Ma in, and D. S. Wilk s (1998), Statistical downsc aling of gener al circulation m odel output : a comparison of m ethods, Water Resour. Res. , 34 (11), 2995 – 3008. Yang, Z ., S. Taraphdar, T. Wang , L. R. Leung, and M. Grear (201 6), Uncerta inty and feas ibility of dynamical downsca ling for modeling tropical cy clones for storm surg e simulation, Nat. Hazards , 84 (2), 1161 – 1184, d oi:10.1007/s11069- 016-2482- y. Yoo, J., and J. Galewsky (2016), Dy namical downscaling of the western No r th Pacific from CCSM4 simulations dur ing the last g lac ial maxim u m and late 20th c entury using th e WRF m odel: model configuration and v al idation, Clim Past D i scuss , 2016 , 1 – 33, do i :10.5194/cp- 2015-170. Zadra, A., M. Roch, S. L aroche, and M. Ch arron (2003), The subg rid-scale orographic blo cking parametrization o f the GEM Model, Atmosp here-Ocean , 41 (2), 15 5 – 170. Zhang, X., H. Wan, F. W. Zwiers, G. C. Hegerl, and S . -K. Min (2013), A ttributing i ntensification of precipitation ex tremes to hum an i nfluence, Geophys. Res. Lett. , 40 (19), 5 252 – 5257, doi:10.1002/grl.5 1010. Zhang, X., F. W. Z wier s, G. Li, H. W an, and A. J. Cannon (2017), Com plexity in estimating past and future extrem e short-duration rainfal l, Nat. Geosci. , 10 (4), 255 – 259. 36 List of Main Figures Figure 1. Time series of (a) 1- ( top panel ) and ( b) 24- ( bottom panel ) hour annual maxim a precipitation before bias-correction in 8 selected stations across Southern Ontario. The b lac k lines indicate observed AM precipitation and the blue multi -m odel median NA-CORDEX simulated precipitation. The shaded bounds in blue indicate the m inimum and the maxim u m RCM simulated pr ec ipitation ext remes. Figure 2. Taylor diagrams to evaluate the skill of bias -correct ed 1-hour AM r egional climate mode l ensembles relative to observ ation during baseline (1970 – 2010) period. Sym bols represent the method of bias-correction (red: Kernel -based and black: GEV-based). The redial di stance is the spatial standar d deviation ( s  ) normalized by the observation s and is a m easure of spa t ial variat ion. s  contours f rom t he origin are shown in bl ack. Contours showing t he RMSE di fference between models and the observation are shown in green. The cosine of the azimuthal angle i s given by the centered pattern correlation ( r ) between the observation and t he RCM. Figure 3. Sam e as Figure 2 but f or 24-hour AM pre cipitation extrem e. Figure 4. Heat maps of the skill scores for bias-corrected AM precipitation relative to observ at ion for different durations during baseline (1970 – 2010) period. The top six rows show k er nel -based bias correction while the bo ttom six rows indicate GE V-based bias-correction. The dar k shade indicates higher skills, whereas the lighter ones show, the lower skill s core. Figure 5. PDF of observ ed versus multi-m odel median AM precipit ation at 1- hour duration in eig ht st ation locations. The m ini mum an d the maxim um bounds are shown us ing dotted blue and red lines respect ively. Figure 6. Sam e as Figure 5 but f or 24-hour precip i tation extrem e. Figure 7. Sam e as Figure 1 but a fter bias co rrection. Figure 8. T he distribution of shape parameters (Y -axes) of the simulated GEV di stribution at diffe rent durations (x - axes) for the pr esent day cl imate (1970-2010) assum ing st ationary GEV models. The box-plot in blue indicates the observed distribution and in red represent multi-model median NA-CORDEX RCM simulated distribution. The horizontal segm ent inside the box shows the median, and the whiskers above and below the box show the loc ations of the minimum and m aximum. T he span of the boxes represen ts the interquartile rang e or the variabi l ity in the distrib ution. Figure 9. Sam e as Figure 8 but f or nonstationary G EV models. Figure 10. Comprehensive comparison of nonstationary versus stationary ID F distributions of observed versus RCM sim ul ated precipitation extrem es f or differ ent return periods a t Toronto Internationa l Airport. 37 The boxplots in darker shade i ndicate observed extreme modeled using stat ionary (red) and nonstationary (blue) GEV models, whe reas the one in lighter shade s hows NA-CORDEX sim ul ated extreme (stationary : green and orange: nonstationary). T he vertical span of the boxes in t he boxplot represents uncertainty in the estimated de sign storm. The leg end applies to all fig ure panels. Figure 11. Sam e as Figure 10 but f or Oshawa WPC P. Figure 12. Sam e as Figure 11 but f or Windsor Ai rport. Figure 13. The relative difference in design storm estim ates of observed versus multi-model median desig n storm in all station l ocations at 1 - and 24-hour storm durations corresponding to 10 - and 50-year return periods during present day cl imate modelled us i ng (a) sta tionary (top pan el) and (b) nonstat ionary (bott om panel) GEV m odel s. Figure 14. Return periods versus design storm estimates of present day 1 -hour precipitation extreme modelled using (a) stationary and (b) nonstationary GEV models. Colored -boxes expr ess the inter- mo del variability represented by the m ult i-model m inimum and m axi mum of RCM sim ulations (best -worse case). The subplots are show n for the 8 stat ions across South ern Ontario. Figure 15. Sam e as Figure 14 but f or 24-hour precip itation extrem e. Figure 16. The distribution of shape parameters of the simulated GEV distribution at different durations for the future period (2030-2070) assuming stationary GEV models. The box -plot in blue indicates the observed distribution and in red represents multi-model median NA-CORDEX RCM si mulated distribution. Figure 17. Sam e as Figure 16 but f or nonstationary G EV models. Figure 18. Nonstationary pr esent-day (in blue) and Projected (i n red; 2030-2070) IDF curves for 10 -year return period across eight stations in Southern Ontario. The uncertainty in IDF simulation at different duration is ex pressed using boxplots. The minimum ( in green) and m aximum (in red) bounds in I DF curves are shown to expres s the best and w orst plausible s cenarios. Figure 19. Nonsta tionary present- day and Projected IDF curves for 50- year return per iod. Figure 20. Design storm intensity of future scenario versus present day climate for 10 -year return period. The variability in the differences in design storm intensity are shown using box -plots representing three different cases: Stationary ( 2030-2070; or 2050s) versus stationary (1970 – 2010; or 1990s) [in purple]; Nonstationary (2050s) versus stationary (1990s) [in blue]; Nonstationary (2050s) versus nonstationary (1990s) [in orang e]. 38 Figure 21. Sam e as Figure 20 but for 50-year retu rn period. Figure 22. Projected c hanges in rainfall intensity for 10 -year r eturn period in three different cases: stationary (2050s) versus stationary (1990s) [ top panel ]; Nonstationary (2050s) versus stationary (1990s) [ middle panel ]; Nonstatio nary (2050s) versus nonstationary (1990s) [ bottom panel ] . The comparative assessment i s performed bet ween projec ted storm intensity modeled using multi-model median NA- CORDEX RCM ensemble and ob served baseline intensity. The shades of the changes express high and low end, with a dark red indicating increase in storm intensity while dark blue show decrease in the intensity . Very sm all changes [ i.e., in and around z ero values] ar e mark ed with white. Figure 23. Sam e as in Figure 21 b ut for 50-year re turn period. 39 Figure 1. Ti me series of (a) 1- ( top panel ) an d (b) 24- ( bottom pane l ) hour annual m axima precipita tion befor e bias-correctio n in 8 selected stat ions across Southern O ntario. The bl ack lines ind icate observed AM precipitat ion and the blue multi -model median N A-CO RDEX sim ulat ed precipitation. The shaded bounds i n blue indicate the m inimum and the m axi mum RCM sim ulat ed precipita tion extremes. 40 Figure 2. Taylor diagram s t o evaluate the skill of bias-correct ed 1-hour AM r egional climate model ensembles relative to observation during basel ine (1970 – 2010) pe riod. Symbols r epresent th e m et hod of bias-correct ion (red: Kernel-based and black: GEV-based ). The r edial d istance i s the spatial standard deviation ( s  ) normalized by the observations and is a measure of spatial variation. s  contours from the origin are shown in black. Contours showing the RMSE difference be tween m odels and the observation are shown in green. The cosine of the azimutha l angle is given by the centered pattern correlation ( r ) betw een the observ at ion and the RCM. 41 Figure 3. Same as Fig ure 2 but for 2 4-hour AM pre cipitation extrem e. 42 Figure 4. Heat maps of the skill scores for bias-corrected AM precipitation relative to observation for different durations during baseline (1970 – 2010) period. T he top six rows show kernel-based bias correction while the bottom si x r ows indicate GEV -based bias- cor rection. The dark shade indicates higher sk ills, whereas the lighter ones show, the lower sk ill score. 43 Figure 5. PDF of observed versus multi-m odel median AM precipitation at 1-hour duration in eight station locations. The minimum and the maximum bounds are shown using do tt ed blue and re d lines respectively. 44 Figure 6. Same as Fig ure 5 but for 2 4 -hour precip itation extrem e. 45 Figure 7. Sam e as Figure 1 but a fter bias co rrection. 46 Figure 8. The distribution of shape parameters (Y -axes) of the simulated GEV distribution at different durations (x -axes ) for t he present day clim ate (1970-2010) assuming st ationary GEV models. The box-plot in blue indicates the observed distribution and in red represent multi-model median NA -CORDEX RCM simulated distribution. The horizontal segment inside the box shows t he median, and the whi skers above and below the box show the location s of the minim um and maximum. The span of the boxes represents the int erquartile rang e or the variabi l ity in the distri bution. 47 Figure 9. Sam e as Figure 8 but f or nonstationary G EV models 48 Figure 10. Com prehensive comparison of nonstatio nar y versus stationary IDF distributions of observed versus RCM simulated precipitation extremes for different return periods at T oronto International Airport. T he boxplots in darker shade indicate observ ed extrem e modeled using stationary (red) and nonstationa ry (blue) GEV models, whereas the one in lighter shade shows NA -CORD EX simulated extreme (stationary : green and orang e: nonstat ionary). The v ertical span of the bo xes in the boxp l ot repre sents uncertainty in t he estimated des ign sto rm. The legen d applie s to all figure panels . 49 Figure 11. Sam e as Figure 10 but f or Oshawa WPC P. 50 Figure 12. Sam e as Figure 11 but f or Windsor Ai rport. 51 Figure 13. The relative difference in design storm estimates of observed versus multi-model median design storm in all station l ocations at 1 - and 24-hour st orm dur ations corresponding to 10 - and 50 -year return periods during present day climate modelled using (a) stationary ( top panel) and (b) nonstationary ( bottom panel) G EV m odel s. 52 Figure 14. Return periods versus design stor m estimates of present day 1 -hour precipitation extrem e modelled using (a) stationary and (b) nonstationary GEV models. Colored-boxes express the i nter-model variability represented by the multi-model minimum and maximum of RCM simulation s (best-worse case ). The subplo t s are shown for the 8 sta tions across Southe rn Ontario. 53 Figure 15. Same as F i gure 14 but for 24- hour precipitation extrem e. 54 Figure 16. The distribution of shape parameters o f the sim ulated GE V d i stribution at differ ent durations for the fu ture period (2030 -2070) assum i ng stationary GEV models. The box-plot in blue indicates the observed distribution and in red represents multi-model median NA-CORDEX RCM simulated distrib ution. 55 Figure 17. Same as F i gure 16 but for nonstationary GEV m odels. 56 Figure 18. Nonstationary present-day (in blue) and Projected (in red; 2030-2070) IDF curves for 10-year return period across eight stations in Southern Ontario. T he uncertainty in IDF si mulation at different duration is expressed using boxpl ots. T he minimum (in green) and maximum (in red) bounds i n IDF curv es are shown to express the best and wo rst plausible scenarios. 57 Figure 19. Nonsta tionary present- day and Projected IDF curves for 50- year return period. 58 Figure 20. Design storm i ntensity of future scenario versus present day clim ate f or 10 -year re turn period. The v ar iability i n t he differences in desig n storm intensity are shown using box-plots representing three different cases: Stationary (2030 -2070; or 2050s) versus stationary ( 1970 – 2010; or 1990s) [in purple] ; Nonstationary (2050s) versus st ationary (1990 s) [in blue] ; Nonstationary (2050s) v ersus nonstat ionary (19 90s) [in ora nge]. 59 Figure 21. Sam e as Figure 20 but f or 50-year retu rn period. 60 Figure 22. Projected changes in rainfall intensity for 1 0-year return period i n three differen t cases: stationary (2050s) versus stationary (1990s) [ top panel ]; Nonstationary (2050s) versus stationary (1990s) [ middle panel ]; Nonstat ionary (2050s) versus nonstationary (1990s) [ bottom panel ]. T he comparative assessm ent i s performed between projected st orm int ensity modeled using multi -model median NA-CORDEX RCM ensem ble and observed baseline intensity. The shades of t he changes express high and l ow end, with a dark red i ndicating increase in st orm intensity whil e dark blue show decrease in the intens ity. Very sm al l changes [ i.e., in and aro und zero values] are mark ed with white. 61 Figure 23. Sam e as in Figure 21 b ut for 50-year re turn period. 1 Assessment of Future Changes in Intensity -Duration-Frequency Curves f or Southern 2 Ontario using North American CORDEX Models w ith Nonstationary Methods 3 4 1 McMaster Water Res ources and Hydr ologic Modelin g Group, D epartment of Civil Engineering , 5 McMaster Universit y, 1280 Main Street West, Ha m ilton ON L 8S 4L7, Canada 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 SI 1. T emporal Do w nscaling o f Precipitation : Multiplicat i ve Rando m Cascade-based disa ggre gatio n 29 Tool 30 The Cascade-based temporal downscaling method in this study is adopted from [ Olsson , 1995, 31 1998a]. The m et hod i mplem ents a micro-canonical (exact conservation of mass in each c as cade branching) 32 approach, in which daily rai nfall is disaggregated using a uniformly distributed generator, dependent on 33 rainfall intensity and position of the rain sequen ce. The m ethod was later developed and t ested in differen t 34 climatic conditions [ Güntner et al. , 2001; Rana et al. , 2013] . In the disaggregation process, each time 35 interval (box) at a given re solution (for exam ple 1 day) is split into two ha l f of the original length (1/2 day). 36 The procedure is continued as a cascade until the desired t ime resolution is reached, i.e., to ¼ day, t hen to 37 1/8 of a day and so on. E ach st ep is termed as a cascade step, with cascade step 0 as t he longest tim e -period 38 with only one box ( i.e. , a day). The distribu tion of th e volume between two sub-interv als ( or smaller boxes) 39 is computed by multiplication with the cascade weights ( , ir W ), , 01 ir W  that assigns ,, i r i r WR  to the 40 first half of the period and 󰇡   , ir W 󰇢 , ir R  to the next half. In each branching two possibilit ies exist: (1) 41 W 1 = 0, W 1 = 1 (2) 0 < W 1 < 1. T he occurrence of (1) and ( 2) may be expressed in terms of probabilities, 42 P 01 =  󰇛    󰇜 or  󰇛    󰇜 = P(W 1 = 0 or W 1 = 1) and P xx =  󰇛    󰇜 = P(0 < W 1 < 1) = 1 - P 01 . 43 De pending on t he range of resolution involve, P 01 either be assumed as resolution independent or 44 parameterized as a scaling law:   󰇛  󰇜       where c 1 and c 2 are constants. The distribution of , ir W is 45 termed as cascade generator, assumed to follow 1 -param eter beta distribution [ Olsson , 2012]. Following 46 [ Olsson , 1998b], the probabilities P , the probability distribution of cascade generator are assumed to be 47 related to (1) position in rainfall sequence, and (2) rainfall volume. The wet boxes, with a rainfall volume 48 V > 0, can be characterized by their position in the rainfall series: (1) the starting box, box preceded by a 49 dry box ( V = 0) and succeeded by a wet box ( V > 0); (2) the enclosed box, box preceded and s ucceeded by 50 wet boxes; (3) the e ndi ng box, box preceded by a wet box and suc ceeded by a dry box, and (4) the i solated 51 box, box preceded and succeeded by dry boxes. On the other hand, based on volume dependence, if the 52 volum e is lar ge then it is m ore likely that both halves o f the subinterval s con t ribute to nonzero volume than 53 if the volume is small. Following [ Olsson , 1998b], a partition into three volum e classes ( v c = 1, 2, 3) was 54 used, separated by percentiles 33 rd and 67 th of the valu es at the cascade step. Next, the variation of  󰇛    󰇜 55 with volume is param et erised as,   / m c xx P a b v  , where a is the intercept a t c v = 0, m b is the mean slope 56 of l inear regression ob t ained from all cascade steps and c v is volum e clas s. The theor y and implem entation 57 issues of MRC- bas ed disaggreg at ion tool are discussed fur ther in [ Olsson , 1998b, n.d.]. 58 59 The objectives of employing cascade -based disagg regation tool are as follows: (i) to infill missing values 60 in t he station-based AM precipitation time series (ii ) t o disaggregate NA -CORD EX r egional climate m odel 61 output available at daily time step to an hourly time scale. First, the cascade parameters are calibrated using 62 good quality observed historic al dataset available at an hourly resolution. The high-resolution historical 63 (1970 – 2010) data are collected from Toronto Resources Conservati on Authority (TRCA). At l east 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, f ive successive 66 “halving” f r om one day to generate 45 -minute (27 00 seconds ) data. After calibration, Monte Ca rlo 67 simulation is performed to gradually fine-grain the data and generate realizations at desired resolution, R s . 68 Since R s in this cases is not directly achieved by exact resolut ion doubling from original resolution, R l , the 69 target resolutions are obtained by geometric i nterpola tion of the disagg regated model out put at t he f inal 70 time step. Then annua l maxim a value is obtained fr om hourly disagg regated time series f or each year. The 71 statistic, we use for the validation is by comparing maximum disaggregated rainfall time series re ference 72 to observed rainfall. Fig S2 (t op panel) show t he perform ance of disaggreg at ion tool for modeling hourly 73 precipitation extremes at Toronto Internation al Ai rport. The daily AM series is also shown (Fig S2; bo tt om 74 panel) to com pare the magnitude of downs caled rainfall volume relative to the observ ed dai ly extreme. T he 75 bars in red shows observed series and blue disaggregated AM series. We adjusted occasional overestim ation 76 of disaggreg at ed rainfall extreme (fo r instance, 1980), by a GEV-based bias correction. The bias cor rected 77 AM series is shown in black bars. A reasonable agreem ent is not ed between observed and temporally 78 downscaled with bias-corrected precipitation for Toronto Airpo rt. Following a similar procedure, we 79 infilled the missing gaps in the AM series of eight station locations across Southern ON. Fig. S 3 shows 80 observed (red) versus disag gregated infilled (blu e) hour ly AM tim e ser ies. As obse rved from the figure, t he 81 infilled time series preserve the trend in the data reasonably well, without exhibiting any signature of 82 outliers ( i.e. , extr emely high or low v al ues). 83 84 Likewise, to disag gregate the daily N A-CORDEX output to hourly time step, we assum e that dynamically 85 downscaled model output is statistical ly similar t o that of station observed precipitation. For t his, we 86 compare t he distribution of daily prec ipitation 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 ens embles. As observed from the figure, the RCMs able to reproduce the observed daily 89 precipitation satisfactorily , al though some bi as exist. Among individual model performance, i n general, 90 CanRCM4 shows maximum variability (or spread) i n simulating precipitation extremes than that of the 91 other two RCM s. To disaggreg ate the daily RCM d at a to an hourly time step, we use the calibrated m odel 92 parameters deriv ed in the historical simulation t o downscale precipitation in the baseline and f uture periods . 93 Then two bias correction sc hemes are compared (Kernel- and GEV-based) to match station observed 94 precipitation. The performance of bias-corrected RCMs are assessed using various skill metrics as discussed 95 in the main m anus cript. 96 97 Table S1. Detai ls of available hou rly data set for Casca de-based disaggreg ation model calibration 98 Station nam e Start Year End Year Number of y ears Toronto I nternational Airpo r t 1961 2010 50 Hamilton Airpor t 1970 2002 40 Oshawa WPCP 1970 2005 36 Windsor Airport 1960 2006 47 London Internationa l Airport 1960 2001 42 Trenton Airport 1965 2000 36 Stratford WP TP 1966 2005 40 Fergus Shand Dam 1961 2006 46 99 100 101 102 103 104 105 106 107 108 109 Table S2. Skill scores of the Kerne l and GEV-based bias correction relative to observations for hourl y AM precipitation 1 Locations Kernel-based GEV-based CanRCM4 CRCM5 RegCM4 MMin MMed MMax CanRCM4 CRCM5 RegCM4 MMin MMed MMax Shand 0.502 0.327 0.474 0.434 0.337 0.449 0.499 0.317 0.489 0.455 0.342 0.457 Hamilton 0.420 0.631 0.453 0.398 0.416 0.641 0.424 0.626 0.449 0.429 0.405 0.636 London 0.412 0.426 0.433 0.404 0.331 0.448 0.394 0.431 0.420 0.391 0.331 0.438 Oshawa 0.394 0.413 0.612 0.405 0.535 0.505 0.399 0.426 0.609 0.409 0.541 0.504 Stratford 0.423 0.431 0.451 0.424 0.376 0.464 0.460 0.453 0.454 0.462 0.373 0.458 Toronto 0.382 0.455 0.505 0.395 0.540 0.424 0.373 0.467 0.507 0.385 0.539 0.423 Trenton 0.544 0.440 0.436 0.544 0.492 0.383 0.542 0.441 0.441 0.541 0.481 0.392 Windsor 0.599 0.559 0.425 0.561 0.569 0.507 0.604 0.569 0.427 0.555 0.581 0.541 2 3 4 5 Table S3. Skill scores of the Kerne l and GEV-based bias correction relative to observations for 2-h ourly AM precipitation 6 Locations Kernel-based GEV-based CanRCM4 CRCM5 RegCM4 MMin MMed MMax CanRCM4 CRCM5 RegCM4 MMin MMed MMax Shand 0.555 0.407 0.493 0.435 0.405 0.570 0.552 0.404 0.489 0.435 0.398 0.561 Hamilton 0.464 0.485 0.325 0.352 0.346 0.440 0.443 0.579 0.410 0.425 0.413 0.546 London 0.411 0.391 0.488 0.438 0.360 0.411 0.403 0.398 0.489 0.432 0.359 0.416 Oshawa 0.445 0.484 0.607 0.515 0.540 0.558 0.445 0.498 0.610 0.517 0.540 0.557 Stratford 0.391 0.530 0.390 0.402 0.488 0.394 0.457 0.531 0.364 0.397 0.479 0.256 Toronto 0.386 0.456 0.489 0.415 0.498 0.395 0.379 0.464 0.486 0.455 0.486 0.397 Trenton 0.542 0.506 0.422 0.572 0.578 0.383 0.551 0.509 0.423 0.567 0.577 0.381 Windsor 0.592 0.557 0.490 0.564 0.607 0.540 0.596 0.561 0.493 0.564 0.605 0.543 7 Table S4. Skill scores of the Kerne l and GEV-based bias correction relative to observations for 6-h ourly AM precipitation 8 Locations Kernel-based GEV-based CanRCM4 CRCM5 RegCM4 MMin MMed MMax CanRCM4 CRCM5 RegCM4 MMin MMed MMax Shand 0.477 0.494 0.453 0.456 0.410 0.521 0.484 0.485 0.481 0.451 0.414 0.531 Hamilton 0.399 0.505 0.413 0.386 0.294 0.442 0.404 0.564 0.413 0.389 0.411 0.443 London 0.470 0.414 0.517 0.412 0.491 0.464 0.474 0.416 0.539 0.429 0.507 0.464 Oshawa 0.407 0.438 0.710 0.537 0.477 0.542 0.408 0.458 0.711 0.535 0.482 0.550 Stratford 0.414 0.674 0.399 0.400 0.373 0.549 0.389 0.623 0.387 0.418 0.407 0.530 Toronto 0.433 0.456 0.447 0.326 0.350 0.513 0.439 0.485 0.490 0.359 0.382 0.509 Trenton 0.558 0.552 0.380 0.645 0.546 0.365 0.562 0.553 0.380 0.641 0.556 0.368 Windsor 0.584 0.512 0.557 0.582 0.633 0.538 0.568 0.510 0.558 0.591 0.619 0.511 9 10 Table S5. Skill scores of the Kerne l and GEV-based bias correction relative to observations for 12- hourly AM precipitation 11 Locations Kernel-based GEV-based CanRCM4 CRCM5 RegCM4 MMin MMed MMax CanRCM4 CRCM5 RegCM4 MMin MMed MMax Shand 0.515 0.483 0.433 0.401 0.479 0.502 0.537 0.484 0.442 0.408 0.483 0.533 Hamilton 0.382 0.525 0.459 0.442 0.430 0.439 0.382 0.522 0.462 0.442 0.433 0.433 London 0.573 0.478 0.672 0.629 0.642 0.596 0.582 0.483 0.682 0.641 0.633 0.609 Oshawa 0.436 0.363 0.535 0.443 0.395 0.494 0.423 0.415 0.606 0.530 0.478 0.487 Stratford 0.421 0.654 0.467 0.385 0.469 0.490 0.364 0.657 0.464 0.384 0.479 0.480 Toronto 0.439 0.441 0.450 0.345 0.384 0.490 0.400 0.457 0.282 0.313 0.417 0.452 Trenton 0.582 0.516 0.373 0.633 0.443 0.468 0.584 0.503 0.379 0.592 0.449 0.469 Windsor 0.479 0.520 0.628 0.451 0.557 0.575 0.379 0.513 0.619 0.459 0.540 0.458 12 13 14 15 16 Table S6. Skill scores of the Kerne l and GEV-based bias correction relative to observations for daily AM precipitation 17 Locations Kernel-based GEV-based CanRCM4 CRCM5 RegCM4 MMin MMed MMax CanRCM4 CRCM5 RegCM4 MMin MMed MMax Shand 0.577 0.551 0.534 0.409 0.591 0.647 0.585 0.559 0.535 0.412 0.598 0.648 Hamilton 0.593 0.480 0.476 0.586 0.462 0.538 0.590 0.452 0.481 0.573 0.470 0.539 London 0.618 0.438 0.635 0.557 0.523 0.630 0.593 0.419 0.653 0.349 0.548 0.606 Oshawa 0.530 0.349 0.566 0.418 0.396 0.524 0.518 0.370 0.584 0.438 0.454 0.547 Stratford 0.559 0.643 0.433 0.489 0.490 0.602 0.578 0.610 0.441 0.469 0.492 0.626 Toronto 0.526 0.471 0.472 0.371 0.493 0.511 0.507 0.462 0.470 0.410 0.477 0.498 Trenton 0.526 0.550 0.303 0.551 0.492 0.338 0.502 0.540 0.306 0.506 0.474 0.338 Windsor 0.448 0.544 0.634 0.443 0.548 0.542 0.443 0.544 0.636 0.407 0.548 0.535 18 19 Table S7. Mann-Kendall trend statistics and Theil-Sen slope per decade of observed versus MMed CORDEX RCM for hourly AM 20 precipitation 21 Stations Observed MMed RCM MK Statistics b MK Statistics b Hamilton -0.048 (0.962) -0.019 -0.416 (0.678) -0.689 London -1.011 (0.312) -1.293 -0.607 (0.544) -0.823 Oshawa 0.315 (0.753) 0.153 0.876 (0.381) 0.940 Shand 0.146 (0.884) 0.143 0.124 (0.902) 0.131 Stratford 0.247 (0.805) 0.261 1.820 (0.069 * ) 2.164 Toronto - 1.775 (0.076 * ) -2.567 1.135 (0.256) 1.540 Trenton 1.011 (0.312) 1.199 1.079 (0.281) 1.085 Windsor -0.809 (0.419) -1.017 0.180 (0.857) 0.287 Value is bracket indicated p-value of the test statistics. Letters in bold indicates statistically significant trend at 10 % sign ificance level. 22 23 24 25 Table S8. Mann-Kendall trend statistics and Theil-Sen slope per decade of observed versus MMed CORDEX RCM for 2-hourly AM 26 precipitation 27 Stations Observed MMed RCM MK Statistics b MK Statistics b Hamilton -0.048 (0.962) -0.034 -0.595 (0.552) -0.753 London -0.887 (0.375) -1.109 -0.966 (0.334) -1.181 Oshawa 0.539 (0.590) 0.578 1.292 (0.196) 1.481 Shand 0.640 (0.522) 1.080 1.269 (0.204) 2.129 Stratford 1.855 (0.064*) 0.961 1.719 (0.086*) 3.864 Toronto -1.955 (0.051*) -2.484 2.202 (0.028*) 2.871 Trenton 1.178 (0.239) 1.579 0.135 (0.893) 0.262 Windsor -0.112 (0.911) -0.270 2.101 (0.036*) 3.445 28 29 Table S9. Mann-Kendall trend statistics and Theil-Sen slope per decade of observed versus MMed CORDEX RCM for 6-hourly AM 30 precipitation 31 Stations Observed MMed RCM MK Statistics b MK Statistics b Hamilton -0.573 (0.567) -0.585 -0.124 (0.902) -0.114 London -0.622 (0.534) -0.652 0.887 (0.375) 1.957 Oshawa 0.697 (0.486) 0.988 1.143 (0.253) 1.382 Shand 0.371 (0.711) 0.935 0.244 (0.808) 0.207 Stratford -0.708 (0.479) -0.927 0.708 (0.479) 1.282 Toronto -1.393(0.164) -2.022 1.910 (0.056*) 2.505 Trenton -0.622(0.534) -0.652 2.483 (0.013*) 2.707 Windsor 0.607(0.544) 1.016 1.720 (0.085*) 2.560 32 33 34 35 Table S10. Mann-Kendall trend statistics and Theil-Sen slope per decade of observed versus MMed CORDEX RCM for 12-hourly 36 AM precipitation 37 Stations Observed MMed RCM MK Statistics b MK Statistics b Hamilton -0.584 (0.559) -1.11 0.191 (0.849) 0.312 London -0.371 (0.711) -0.58 1.101 (0.271) 1.851 Oshawa 1.370 (0.171) 1.47 -0.247 (0.805) -0.540 Shand -0.687 (0.492) -0.92 0.382 (0.703) 0.873 Stratford -1.022 (0.307) -1.74 1.674 (0.094*) 4.113 Toronto -0.539 (0.590) -0.99 0.663 (0.507) 1.154 Trenton -0.371 (0.711) -0.58 0.748 (0.455) 1.453 Windsor 0.247 (0.805) 0.69 2.917 (0.004*) 2.227 38 39 Table S11. Mann-Kendall trend statistics and Theil-Sen slope per decade of observed versus MMed CORDEX RCM for daily AM 40 precipitation 41 Stations Observed MMed RCM MK Statistics b MK Statistics b Hamilton -1.13 (0.257) -2.48 -0.326 (0.745) -0.515 London -0.831 (0.406) -1.48 1.359 (0.174) 1.754 Oshawa 2.572 (0.010*) 4.188 -1.134 (0.257) -1.746 Shand -1.574 (0.116) -1.899 0.326 (0.744) 0.665 Stratford -1.461 (0.144) -2.228 1.032 (0.302) 1.965 Toronto -0.438 (0.661) -0.638 -0.775 (0.438) -1.410 Trenton 2.741 (0.006*) 5.018 1.606 (0.108) 2.799 Windsor 0.899 (0.369) 1.460 0.573 (0.567) 0.961 42 43 44 45 Table S12. Mann-Kendall trend statistics and Theil-Sen slope per decade of historical daily AM temperature anomaly 46 Stations Observed MK Statistics b Hamilton 1.33 (0.18) 0.254 London 1.34 (0.18) 0.308 Oshawa 0.22 (0.82) 0.070 Shand 0.25 (0.80) 0.00 Stratford 1.48 (0.14) 0.463 Toronto 1.16 (0.25) 0.370 Trenton 0.60 (0.55) 0.142 Windsor -0.35 (0.73) -0.035 47 Table S13. Estimates of the median shape parameter of observed vs modeled MM-Med RCMs for baseline (1970 – 2010) simulation 48 GEV-T ype Sites Observed MM - Med 1-hour 2-hour 6-hour 12 -hour 24 -hour 1-hour 2-hour 6-hour 12 -hour 24 -hour Stationary Hamilton 0.166 0.140 0.278 0.221 0.184 0.169 0.109 0.150 0.225 0.175 London 0.029 0.060 0.128 0.077 0.177 0.053 0.047 0.123 0.060 0.037 Oshawa -0.066 0.113 0.078 0.235 0.156 -0.046 0.122 0.186 0.268 0.211 Shand 0.159 0.121 -0.073 0.024 0.049 0.172 0.262 0.152 0.053 0.055 Stratford 0.170 0.158 0.276 0.337 0.052 0.074 0.159 0.407 0.275 0.287 Toronto -0.102 0.039 0.149 0.294 0.247 -0.079 0.026 0.036 0.419 0.158 Trenton 0.142 0.172 0.188 0.178 0.064 0.146 0.130 0.183 0.238 0.151 Windsor 0.111 -0.118 0.102 0.198 0.100 0.165 -0.068 -0.032 0.144 0.190 Nonstationary Hamilton 0.184 0.046 0.404 0.190 0.243 0.222 0.075 0.216 0.260 0.152 London -0.060 -0.031 0.150 0.068 0.291 0.127 0.105 0.224 -0.022 0.048 Oshawa -0.052 0.333 0.008 0.276 0.152 -0.361 0.157 0.148 0.376 0.247 Shand 0.144 0.283 -0.102 0.024 0.061 0.155 0.182 0.202 0.054 0.073 Stratford -0.037 0.270 0.351 0.314 0.311 0.157 0.196 0.499 0.369 0.308 Toronto -0.044 0.049 0.101 0.408 0.281 -0.180 0.063 0.066 0.435 0.146 Trenton 0.074 0.279 0.268 0.230 0.080 0.114 0.248 0.166 0.333 0.118 Windsor 0.116 -0.101 0.109 0.193 0.072 0.196 -0.055 -0.016 0.132 0.228 49 Table S14. GEV shape parameters fo r stationa ry and nonstationar y models of obs erved ve rsus MM-Med RCM simulation during 50 baseline period for Toronto International Airpor t 51 * GEV paramete rs are estimated using Baye sian Inference. The reported values are med ian of DE -MC sampled parameters. Bracketed values indicate 5 th and 52 95 th percentile of simu lated GEV parameters ob tained from DE -MC simulatio ns 53 54 55 Time Slice Data Model Location par ameter Scale param eter Shape param eter 0  1    1-hr Observati on GEV sta 21.38 (19.07, 23.14) 7.59 (6.32,9.1) -0.10 (- 0.24,0.08) GEV Nonsta 27.28 (23.51, 31.26) -0.27 (- 0.41, -0.14) 7.46 (6.38,8.9 7) -0.04 (- 0.20,0.10) MM -Med CORDEX GEV sta 21.04 (18.85, 22.87) 7.73 (6.64, 9. 88) -0.07 (- 0.22, 0.10) GEV Nonsta 19.51 (15.03, 23.54) 0.10 (-0.07, 0.2 9) 8.0 (6.76, 9.8 9) -0.18 (- 0.26, -0.12) 2-hr Observati on GEV sta 11.82 (10.75, 12.96) 4.23 (3.46, 5. 26) 0.03 (-0.08, 0.2 5) GEV Nonsta 13.65 (11.53, 15.71) -0.08 (- 0.17, -0.006) 4.11 (3.39, 5. 05) 0.05 (0.03, 0. 07) MM -Med CORDEX GEV sta 11.83 (10.63, 12.92) 4.39 (3.70, 5. 05) -1.53 e- 5 ( -0.14, 0.20) GEV Nonsta 9.47 (7.18, 11 .66) 0.12 (0.02, 0. 21) 4.08 (3.39, 4. 93) 0.04 (0.003, 0 .096) 6-hr Observati on GEV sta 5.08 (4.6, 5.4 3) 1.55 (1.25, 1. 89) 0.18 (0.005, 0 .43) GEV Nonsta 5.58 (5.19, 6. 01) -0.02 (- 0.04, -0.002) 1.55 (1.29, 1. 92) 0.10 (0.008, 0 .21) MM -Med CORDEX GEV sta 3.95 (3.57, 4. 34) 1.48 (1.43, 1. 53) 0.036 (-0.12, 0. 27) GEV Nonsta 3.20 (2.54, 3. 93) 0.033 (0.0002 , 0.061) 1.37 (1.12, 1. 68) 0.06 (-0.01, 0.1 4) 12-hr Observati on GEV sta 2.88 (2.71, 3. 06) 0.65 (0.52, 0. 83) 0.28 (0.11, 0. 54) GEV Nonsta 2.87 (2.52, 3. 25) 0.0001 (-0.015, 0.014) 0.67 (0.52, 0. 87) 0.40 (0.32, 0. 47) MM -Med CORDEX GEV sta 2.82 (2.68, 3. 04) 0.62 (0.46, 0. 82) 0.42 (0.26, 0. 66) GEV Nonsta 2.90 (2.73, 3. 04) -0.001 (- 0.009, 0.006) 0.64 (0.49, 0. 81) 0.43 (0.25, 0. 67) 24-hr Observati on GEV sta 1.68 (1.59, 1. 79) 0.37 (0.31, 0. 46) 0.25 (0.08, 0. 43) GEV Nonsta 1.81 (1.63, 1. 98) -0.006 (- 0.013, 0.001) 0.36 (0.28, 0. 46) 0.28 (0.26, 0. 31) MM -Med CORDEX GEV sta 1.47 (1.39, 1. 52) 0.40 (0.35, 0. 49) 0.16 (0.005, 0 .34) GEV Nonsta 1.58 (1.49, 1. 67) -0.005 (- 0.012, 0.001) 0.42 (0.3 5, 0.49) 0.15 (0.003, 0 .35) Table S15. GEV shape parameters fo r stationa ry and nonstationar y mo dels of observed versus MM -Med RCM simulation during 56 baseline period for Oshawa WPCP 57 58 59 60 Time Slice Data Model Location par ameter Scale param eter Shape param eter 0  1    1-hr Observati on GEV sta 17.39 (15.27, 19.03) 6.69 (5.89, 7. 93) -0.06 (- 0.23, 0.22) GEV Nonsta 17.06 (13.61, 20.07) 0.02 (-0.10, 0.1 5) 6.66 (5.58, 8. 26) -0.05 (- 0.06, -0.045) MM -Med CORDEX GEV sta 17.44 (15.61, 19.41) 6.60 (5.58, 7. 98) -0.046 (- 0.24, 0.12) GEV Nonsta 14.46 (10.24, 18.06) 0.17 (0.05, 0. 32) 8.66 (7.80, 9. 98) -0.36 (- 0.41, -0.29) 2-hr Observati on GEV sta 10.70 (9.54, 1 1.69) 3.43 (2.71, 4. 55) 0.11 (-0.08, 0.3 3) GEV Nonsta 10.02 (8.67, 1 1.25) 0.01 (-0.03, 0.0 6) 3.32 (2.67, 4. 32) 0.33 (0.30, 0. 35) MM -Med CORDEX GEV sta 10.75 (9.91, 1 1.76) 3.23 (2.58, 4. 22) 0.12 (-0.07, 0.2 6) GEV Nonsta 9.07 (7.54, 10 .64) 0.07 (0.01, 0. 14) 3.11 (2.54, 3. 93) 0.15 (0.003, 0 .27) 6-hr Observati on GEV sta 4.81 (4.55, 5. 18) 1.37 (1.11, 1. 69) 0.078 (-0.06, 0. 19) GEV Nonsta 4.33 (3.93, 4. 75) 0.025 (0.009, 0.04) 1.45 (1.22, 1. 77) 0.007 (-0.11, 0. 15) MM -Med CORDEX GEV sta 3.53 (3.25, 3. 89) 1.08 (0.88, 1. 37) 0.19 (0.006, 0 .33) GEV Nonsta 3.06 (2.36, 3. 68) 0.022 (-0.005, 0.05) 1.10 (0.89, 1. 39) 0.15 (0.11, 0. 19) 12-hr Observati on GEV sta 2.82 (2.62, 3. 01) 0.77 (0.64, 0. 93) 0.23 (0.07, 0. 43) GEV Nonsta 2.62 (2.33, 2. 83) 0.007 (-0.002, 0.02) 0.69 (0.58, 0. 85) 0.28 (0.16, 0. 41) MM -Med CORDEX GEV sta 2.81 (2.58, 3. 05) 0.75 (0.59, 0. 94) 0.27 (0.08, 0. 51) GEV Nonsta 2.72 (2.35, 3. 08) 0.003 (-0.011, 0.016) 0.72 (0.57, 0. 94) 0.37 (0.36, 0. 38) 24-hr Observati on GEV sta 1.64 (1.48, 1. 77) 0.52 (0.41, 0. 64) 0.15 (0.006, 0 .28) GEV Nonsta 1.54 (1.36, 1. 71) 0.005 (-0.002, 0.011) 0.50 (0.43, 0. 62) 0.15 (0.07, 0. 20) MM -Med CORDEX GEV sta 1.43 (1.34, 1. 50) 0.37 (0.31, 0. 46) 0.21 (0.06, 0. 38) GEV Nonsta 1.46 (1.28, 1. 64) -0.001 (- 0.007, 0.004) 0.38 (0.32 , 0.45) 0.25 (0.10, 0. 39) Table S16. GEV shape parameters fo r stationa ry and nonstationar y mo dels of observed versus MM -Med RCM simulation during 61 baseline period for Wndsor Airport 62 63 64 65 Time Slice Data Model Location par ameter Scale param eter Shape param eter 0  1    1-hr Observati on GEV sta 23.08 (21.18, 25.32) 7.73 (6.38, 9. 60) 0.11 (-0.06, 0.3 4) GEV Nonsta 26.35 (22.69, 30.5) -0.15 (- 0.32, -0.03) 7.67 (6.47, 9. 60) 0.11 (0.01, 0. 22) MM -Med CORDEX GEV sta 22.48 (20.54, 24.82) 8.10 (6.6, 10. 07) 0.16 (-0.02, 0.3 1) GEV Nonsta 21.12 (17.91, 24.2) 0.06 (-0.05, 0.2 0) 7.67 (6.14, 9. 61) 0.20 (0.16, 0. 23) 2-hr Observati on GEV sta 14.38 (13.07, 15.83) 5.15 (4.33, 6. 11) -0.11 (- 0.29,0.042) GEV Nonsta 15.15 (12.54, 17.62) -0.04 (- 0.13, 0.05) 5.23 (4.46, 6. 27) -0.10 (- 0.13, -0.07) MM -Med CORDEX GEV sta 14.32 (13.05, 15.60) 5.12 (4.31, 6. 37) -0.07 (- 0.23, 0.13) GEV Nonsta 10.39 (7.51, 1 3.16) 0.19 (0.074, 0 .30) 4.72 (4.06, 5. 66) -0.05 (- 0.09, -0.017) 6-hr Observati on GEV sta 5.82 (5.35, 6. 30) 1.78 (1.46, 2. 28) 0.10 (-0.068, 0. 29) GEV Nonsta 5.72 (4.81, 6. 54) 0.004 (-0.034, 0.044) 1.78 (1.46, 2. 19) 0.10 (0.06, 0. 15) MM -Med CORDEX GEV sta 4.73 (4.26, 5. 12) 1.71 (1.46, 2. 04) -0.03 (- 0.16, 0.17) GEV Nonsta 4.06 (3.09, 4. 84) 0.035 (-0.004, 0.077) 1.65 (1.42, 1. 98) -0.016 (- 0.04, 0.004) 12-hr Observati on GEV sta 3.20 (2.96, 3. 45) 0.90 (0.75, 1. 13) 0.19 (0.0008, 0.37) GEV Nonsta 3.19 (2.75, 3. 58) 0.0009 (-0.016, 0.018) 0.88 (0.71, 1. 08) 0.19 (0.17, 0. 21) MM -Med CORDEX GEV sta 3.22 (2.94, 3. 48) 0.96 (0.76, 1. 22) 0.14 (-0.02, 0.3 0) GEV Nonsta 2.99 (2.61, 3. 33) 0.013 (-0.005, 0.03) 0.98 (0.80, 1. 20) 0.13 (-0.007, 0. 32) 24-hr Observati on GEV sta 1.90 (1.77, 2. 08) 0.57 (0.46, 0. 68) 0.10 (-0.12, 0.3 3) GEV Nonsta 1.69 (1.63, 1. 75) 0.008 (0.004, 0.014) 0.53 (0.43, 0. 69) 0.07 (-0.13, 0.3 1) MM -Med CORDEX GEV sta 1.60 (1.48, 1. 74) 0.42 (0.33, 0. 57) 0.19 (0.02, 0. 45) GEV Nonsta 1.41 (1.23, 1. 63) 0.009 (0.007, 0.016) 0.40 (0.32 , 0.52) 0.23 (0.12, 0. 32) Table S17. The AICc values of stationary and nonstationary GEV models at 1- and 6-hour duration during baseline (1970-2010) period 66 Duration (hour) Stations Observation MM -Med NA -CORDEX RCMs GEV sta GEV Nonsta GEV sta GEV Nonsta 1- Hamilton -273.68 -273.67 -290.91 -285.85 London -252.00 -251.15 -286.09 -283.10 Oshawa -268.21 -272.28 -244.83 -238.16 Shand -302.63 -305.21 -275.15 -273.87 Stratford -266.32 -256.17 -260.26 -277.12 Toronto -275.29 -270.06 -275.08 -269.83 Trenton -276.34 -278.01 -298.49 -295.16 Windsor -268.32 -268.59 -285.05 -284.94 6- Hamilton -285.46 -295.95 -245.22 -246.73 London -306.55 -298.91 -241.78 -245.62 Oshawa -288.58 -277.32 -300.48 -288.05 Shand -296.63 -278.46 -277.43 -278.89 Stratford -262.79 -269.36 -284.68 -292.52 Toronto -284.60 -277.39 -284.48 -295.38 Trenton -252.40 -264.43 -283.42 -289.46 Windsor -287.92 -299.88 -289.94 -289.89 * The AICc values are calcu lated between d esign storm estimate from GE V mod el and emp irical d istribution , fitted with Gringorten’ s plotting 67 position formula. The expression for the perfo rmance statistics, co rrected AI C is given as,       2 1 1 c A IC m A IC m m n m      , in 68 which n is the n umber o f observations, m denotes the number of fitted mod el parameters. The best distribu tion is marked as bold. In this case, 69 the best distribution is selected from the min imum AICc statistics. 70 71 72 73 74 75 Table S18. Estimates of the median shape parameter of observed baseline (1970 – 2010) vs MM -Med RCM modeled projected (2030 76 – 2070) scenario 77 GEV-T ype Sites Observed MM - Med 1-hour 2-hour 6-hour 12 -hour 24 -hour 1-hour 2-hour 6-hour 12 -hour 24 -hour Stationary Hamilton 0.166 0.140 0.278 0.221 0.184 -0.075 0.125 0.180 -0.059 0.044 London 0.029 0.060 0.128 0.077 0.177 -0.012 -0.040 0.024 0.021 0.188 Oshawa -0.066 0.113 0.078 0.235 0.156 0.230 -0.083 -0.078 -0.191 -0.087 Shand 0.159 0.121 -0.073 0.024 0.049 -0.058 0.007 -0.331 -0.052 0.102 Stratford 0.170 0.158 0.276 0.337 0.052 -0.092 -0.094 0.087 0.273 0.221 Toronto -0.102 0.039 0.149 0.294 0.247 -0.123 -0.063 0.098 -0.490 -0.032 Trenton 0.142 0.172 0.188 0.178 0.064 -0.260 -0.067 -0.123 -0.567 -0.278 Windsor 0.111 -0.118 0.102 0.198 0.100 0.022 -0.097 -0.225 0.108 -0.112 Nonstationary Hamilton 0.184 0.046 0.404 0.190 0.243 0.017 0.073 0.157 -0.157 0.067 London -0.060 -0.031 0.150 0.068 0.291 -0.036 -0.029 -0.012 0.014 0.230 Oshawa -0.052 0.333 0.008 0.276 0.152 0.210 -0.294 -0.079 -0.220 -0.075 Shand 0.144 0.283 -0.102 0.024 0.061 -0.049 -0.056 -0.326 -0.073 0.030 Stratford -0.037 0.270 0.351 0.314 0.311 -0.190 -0.068 0.084 0.335 0.209 Toronto -0.044 0.049 0.101 0.408 0.281 0.056 -0.036 0.121 -0.653 -0.079 Trenton 0.074 0.279 0.268 0.230 0.080 -0.303 -0.063 -0.152 -0.659 -0.381 Windsor 0.116 -0.101 0.109 0.193 0.072 0.021 -0.012 -0.255 0.131 -0.112 78 79 80 81 82 83 84 85 86 Table S19. The AICc values of stationary and nonstationary GEV models at 1- and 6-hour duration during projected (2030-2070) period 87 Duration (hour) Stations GEV sta GEV Nonsta 1- Hamilton -278.51 -189.54 London -277.65 -268.92 Oshawa -290.30 -281.95 Fergus Shand -288.61 -279.86 Stratford -242.03 -240.96 Toronto -271.06 -263.22 Trenton -263.46 -284.98 Windsor -220.39 -221.36 6- Hamilton -277.22 -280.78 London -281.53 -252.72 Oshawa -245.97 -280.53 Fergus Shand -270.07 -250.75 Stratford -308.54 -311.28 Toronto -289.57 -291.09 Trenton -259.90 -253.57 Windsor -299.06 -224.51 * The AICc values are calculated between de sign storm estimate from GEV model and empirical distribution , fitted with 88 Gringorten’s plotting position formula. 89 List of Supplementary Figure Captions 90 Figure S1. Flow cha r t of the com parison procedure. 91 Figure S2. Selec ted station loc ations in South er n Ontario. The Southe rn Ontario ( 41° - 44°N , 84° - 76°W) 92 is the southernmost region of Canada and is situated on a southwest -northeast transect, bounded by lakes 93 Huron, Erie, and O ntario. The eight locations on the map are ( from southwest to northeast corner ): Windsor 94 Airport, London International Airport, Stratford Wastewater Treatment Plant (WWTP), Fergus Shand Dam, 95 Hamilton Airport, Toronto International Airport, Oshawa Water Pollution Control Plan t (WPCP) and 96 Trenton Airport. Topography map is obtained from 90-m digital elevation model (DEM; SRTM-90m, 97 [ Jarvis et al. , 2008 ] indicates shal low slope with eleva t ion of maxim um 670 m above mean sea level. 98 Figure S3. Observed versus temporally downscaled hourly AM series ( top panel ). The daily AM series 99 ( bottom panel ). 100 Figure S4. Observ ed versus infilled hourly AM series. 101 Figure S5. Distribution of annual maximum preci pitation anomaly in observations and NA -CORDEX 102 simulated RCMs d uring baseline (1970-2010) period. 103 Figure S6. PDF of ob served versus multi-model median 2-hour AM precip i tation in eight station locations. 104 The minim um and the maxi mum bounds are shown using dott ed blue and red lines respec tively. 105 Figure S7. Sam e as Figure S6 but for 6-hour prec ipitation extr eme. 106 Fi gure S8. Same as Fig ure S7 but for 12-hour precip itation extrem e. 107 Figure S9. Present-day (in blue) versus Pr ojected (in red; 2030 -2070) ID F curv es for 10-year return period 108 across eight stations in Sou thern Ontario using stationa ry GEV models. The uncert ai nty in IDF sim ulation 109 at diffe rent duration is expressed using boxplots. The minimum ( in green) and maximum ( in red) bounds 110 in IDF curv es are shown to express the bes t and worst plaus ible scenar ios. 111 Figure S10. Sam e as Figure S7 bu t for 50-y ear return period. 112 Figure S11 . Projected chang es in rainfall intensity for 25 -year return period in three different cases: 113 stationary (2050s) versus stationary (1990s) [ top panel ]; Nonstationary (2050s) versus stationary (1990s) 114 [ middle panel ]; Nonstatio nary (2050s) versus nonstationary (1990s) [ bottom panel ] . The comparative 115 assessment i s performed bet ween projec ted storm intensity modeled using multi-model median NA- 116 CORDEX RCM ensemble and ob served baseline intensity. The shades of the ch anges express h igh and low 117 end, with dark red indicating increase in storm intensity while dark blue show decrease in t he intensity. 118 Very sm all changes [ i.e., in and around z ero values] ar e mark ed with white. 119 References 120 Güntner, A., J. Olsson, A. C alver, and B. Gannon ( 2001), Cascad e-based disagg regation of continuous 121 rainfall time se ries: the infl uence of clim ate, Hydrol. E art h Syst. Sci. Di scuss. , 5 (2), 145 – 164. 122 Jarvis, A., H. I . Reuter, A. Nelson, and E. Gu evara (2008), Hole-fill ed seamless SRT M data V4, 123 International Ce ntre for Tropical Agr i culture (CIAT) . 124 Olsson, J. (1995 ), Limits and cha r acteristics o f the multifrac t al behaviour of a high- resol ution rainfall 125 time series, Nonlinea r Process. Geo phys. , 2 (1), 23 – 29. 126 Olsson, J. (1998a), Ev al uation of a scaling casca de model for temporal rain-fall disag gregation, Hydrol. 127 Earth Syst. Sci. Di scuss. , 2 (1), 19 – 3 0. 128 Olsson, J. (1998b ), Evaluation o f a scaling cas cade m odel for temporal rain-fall disag gregation, Hydrol. 129 Earth Syst. Sci. Di scuss. , 2 (1), 19 – 3 0. 130 Olsson, J. (2012), Random Cascade Model: docum entation, In: Willem s e t al. Im pacts of Clim ate Change 131 on Rainfall Extr emes and U rban Drainag e Systems, I WA Publishing. 132 Rana, A., L. Bengts son, J. O lsson, and V. Joth iprakash (2013), Dev el opment of I DF-curves for tropical 133 india by random cascade model ing, Hydrol. Earth Sy st. Sci. Discuss. , 10 (4 ), 4709 – 4738. 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 Figure S1. Flow cha r t of the com par ison procedure Figure S 2. Selected station l ocations in Southe rn O ntario. T he Southern Ontario (4 1° - 44°N, 84° - 76°W) is the sou t hernmost region of Canada and is situated on a southwest-northeast transect, bounded by lakes Huron, Erie, and Ontario. T he ei ght locations on the map are ( from southwes t to northeast corner ): Windsor Airport, London International Airport, Stratford Wastewater Treatment Plant (WWTP), Fergus Shand Dam, Hamilton Airport, Toronto I nternational Airport, Oshawa Water Pollu tion Control Plan t (WPCP) and Trenton Airport. Topography map is obtained from 90- m digital elevat ion m odel (D EM; SR T M-90m, [ Jarvis et al. , 2008] indicates shallow sl ope with maxim u m elevation of 670 m above mean sea level. Figure S3. Observ ed versus temporally downscaled h our ly AM series ( top panel ). The da i ly AM ser ies ( bottom panel ). Figure S4. Observ ed versus infilled hourly AM series Figure S5. Distribution of annual maximum precip itation anomaly in observations and N A-CORDEX si mulated RCMs during baseline (1970-2010) period. Figure S6. PD F of ob served versus multi-model median 2-hour AM precip i tation in eight station locations. The minim u m and the maxim um bounds are shown using dotted blue and red lines resp ectively. Figure S7. Sam e as Figure S6 but f or 6-hour precip itation extrem e . Figure S8. Sam e as Figure S7 but f or 12-hour precip itation extrem e. Figure S9. Present-day (in blue) versus P rojected (in red; 2030-2070) ID F curves for 10- year return period across eight stations in Southern Ontari o using stationary GEV models. The uncerta inty in IDF simulation at different duration is expressed using boxplots. The minim um (in green ) and maximum (in red) bounds i n IDF curv es are shown to express the best and wo rst plausible scenarios. Figure S10. Sam e as Figure S7 bu t for 50-y ear return period. Figure S11. Projected chang es in rainfall intensity for 25 -year return period in three different cases: st ationary (2050s) versus stationary (1990s) [ top panel ]; Nonstationary (2050s) versus st ationary (1990s) [ middle panel ]; Nonstationary (2050s) versus nonstationary (1990s) [ bottom panel ]. The comparative assessment is perform ed bet ween projected storm intensity modeled using multi -model median NA-CORDEX RCM ens em ble and observed baseline i ntensity. The sh ades of the changes express hig h and low end, with dark red i ndicating increase in s torm intensity while dark blue show decrease in the intensity. Ve ry sm all changes [ i.e., in and around z ero values] are m ar ked with white.