Assessment of Climate Change Impacts on Extreme High and Low Flows: An Improved Bottom-Up Approach

: A quantitative assessment of the likelihood of all possible future states is lacking in both the traditional top-down and the alternative bottom-up approaches to the assessment of climate change impacts. The issue is tackled herein by generating a large number of representative climate projections using weather generators calibrated with the outputs of regional climate models. A case study was performed on the South Nation River Watershed located in Eastern Ontario, Canada, using climate projections generated by four climate models and forced with medium- to high-emission scenarios (RCP4.5 and RCP8.5) for the future 30-year period (2071–2100). These raw projections were corrected using two downscaling techniques. Large ensembles of future series were created by perturbing downscaled data with a stochastic weather generator, then used as inputs to a hydrological model that was calibrated using observed data. Risk indices calculated with the simulated streamﬂow data were converted into probability distributions using Kernel Density Estimations. The results are dimensional joint probability distributions of risk-relevant indices that provide estimates of the likelihood of unwanted events under a given watershed conﬁguration and management policy. The proposed approach o ﬀ ers a more complete vision of the impacts of climate change and opens the door to a more objective assessment of adaptation strategies. that higher air temperatures lead to reduced summer low-ﬂows and less intense high ﬂows. The ensemble results indicate similar trends in both indices—fewer droughts and extreme ﬂooding events.


Introduction
Climate change has forced the research community to revisit assumptions and theories used for the design, planning, and management of water resource systems [1]. Hydrological processes depend to a great degree on the climate regime [2][3][4][5][6][7]. Perturbing the climate regime inevitably results in a perturbed hydrological regime, which in turn affects water-related risk and water resource system performance. Numerous researchers and practitioners have been working intensely to develop methodologies to identify future climate change impacts and viable adaptations strategies [8][9][10].
The most frequently used approach in climate impact assessments is the top-down approach, which is constrained by the availability of Global Circulation Models (GCM) and Regional Climate Models (RCM). Typically, it considers possible future changes in climatic conditions based on predetermined scenarios that are parameterized in large-scale models. This is achieved by downscaling a few distinct projections to be able to determine the impacts on the local hydrological system. However, a critical limitation of this approach is that it ignores plausible risks by not covering all possible future conditions despite using multi-model and multi-scenario projections. For example, while climate

Materials and Methods
One of the major limitations of the top-down approach is that the number of scenarios is too low to fully explore the climate-related risk space. The problem can be partially mitigated by cloning projections generated by corrected climate models and by slightly perturbing them through a stochastic weather generator to obtain a larger set of scenarios that cover the same statistical space as the synthetic time series.
Here is a brief explanation of the 'ensemble generation process' used herein for risk discovery. We firstly used downscaling methods to correct biased (i.e., raw) regional climate model data. Secondly, using corrected (i.e., downscaled) climate data, the stochastic weather generator (MulGETS) was utilized to generate 30-years of daily precipitation and minimum and maximum temperature data. The second step is repeated to generate a large number of climate projections: we run MulGETS 250 with each corrected RCM data. Thirdly, the calibrated SWAT model was forced with these climate projections individually to generate "ensemble" of streamflow projections. Subsequently, the newly-generated ensemble was utilized to determine the risk space that is largely overlooked by top-down methods. This was done to consider the inherent randomness of hydroclimatic variables. Figure 1 presents a schematic representation of all datasets used in this study and the relationships between different time series in the proposed framework. The subsequent sections describe the aforementioned process in details.

Study Area
The South Nation Watershed (SNW) of about 4000 km 2 in Eastern Ontario was used as a showcase. It is characterized by multigauge climatic data (from St. Albert, Russell, Morrisburg, and South Mountains metrological stations) and downstream daily discharge data (collected at Plantagenet gauging station (ID: 02LB005)) ( Figure 2). Hydrogeological investigations have indicated that an impermeable overburden watershed and a lack of streambed gradients in parts of the South Nation River pose a prominent flooding risk [26]. The average maximum and minimum air temperatures are 11.5 and 1.2 °C, respectively, with a mean annual precipitation of around 1000 mm. A more detailed description of the watershed has been presented previously [27].

Study Area
The South Nation Watershed (SNW) of about 4000 km 2 in Eastern Ontario was used as a showcase. It is characterized by multigauge climatic data (from St. Albert, Russell, Morrisburg, and South Mountains metrological stations) and downstream daily discharge data (collected at Plantagenet gauging station (ID: 02LB005)) ( Figure 2). Hydrogeological investigations have indicated that an impermeable overburden watershed and a lack of streambed gradients in parts of the South Nation River pose a prominent flooding risk [26]. The average maximum and minimum air temperatures are 11.5 and 1.2 • C, respectively, with a mean annual precipitation of around 1000 mm. A more detailed description of the watershed has been presented previously [27].

Study Area
The South Nation Watershed (SNW) of about 4000 km 2 in Eastern Ontario was used as a showcase. It is characterized by multigauge climatic data (from St. Albert, Russell, Morrisburg, and South Mountains metrological stations) and downstream daily discharge data (collected at Plantagenet gauging station (ID: 02LB005)) ( Figure 2). Hydrogeological investigations have indicated that an impermeable overburden watershed and a lack of streambed gradients in parts of the South Nation River pose a prominent flooding risk [26]. The average maximum and minimum air temperatures are 11.5 and 1.2 °C, respectively, with a mean annual precipitation of around 1000 mm. A more detailed description of the watershed has been presented previously [27].

Regional Climate Model Data
The availability of daily climatic variables at fine spatial scales is essential for hydrological impact studies [28]. The gridded daily data of precipitation, maximum near-surface air temperature, and minimum near-surface air temperature were extracted from the North American domain of Coordinated Regional Climate Downscaling Experiment (NA-CORDEX) project archive, driven by state-of-the-art GCMs as part of the the Coupled Model Inter-comparison Project Phase 5 (CMIP5) experiment ( Table 1). The CORDEX project was initiated by the World Climate Research Program (WCRP) to provide reliable predictions of regional future climate change. The four independent models were forced with medium-to high-emission scenarios RCP4.5 and RCP8.5 (resulting in 8 different projections). They are originally based on outputs from the CanESM2 and EC-EARTH Global Circulation Models, whose outputs are downscaled with three Regional Climate Models: CanRCM4, RCA4, and HIRHAM5. The horizontal spatial resolution was chosen to be NAM-44 (North American domain at 0.44 • horizontal grid resolution or approximately 50 km).

Downscaling Methods
Downscaling is a procedure undertaken to reconcile the large-scale data representativeness of observations and produce climate change data at a finer resolution. This study makes use of two statistical downscaling techniques-namely Change Factors (CF) and Quantile-Quantile transformation (QQ)-to produce more reliable station-level climatic information to be used later as inputs to a distributed impact model.
The linear correction method applied in this study, namely CF, requires three time series: observations (or synthetic climate series), historical, and future raw RCM data. The primary assumption in this method is that the relationship between historical and future projections is the same as the relation between historical RCM data and observed time series. However, the Quantile-Quantile method applies a different concept that is based on the distribution of raw simulations and observations. For a given variable, downscaling is conducted using a monthly time-step. Observed precipitation and temperature variables of 30 years (1981 to 2010) were used to define a reference climatology and to spatially downscale the coarse-scale climate data in the future (2071-2100).

Change Factor Method
The Change Factors (CF) method is a commonly used approach to linearly downscale a variable (e.g., precipitation) by calculating the difference between the control and future GCM simulations before scaling the baseline observations accordingly [29]. Yet, the conventional CF method is usually used a top-down approach and as such not always reliable. Nevertheless, it has been used excessively in recent climate change assessments [30][31][32][33][34][35][36][37]. The CF approach lacks the ability to correct future projections as it predicts the same variability for the future climate by keeping the same historical temporal structure, which is unrealistic [29,30,38]. One of its limitations is that it produces the exact same temporal structure of wet days (occurrences) in the future [39], which is by no means certain for a stochastic variable such as precipitation. Furthermore, extreme events that are vital in risk assessments cannot be adequately corrected by such a method when using only one observational dataset. To alleviate these disadvantages, we employed stochastic weather generators as they are capable of capturing climate variability. We explored the potential risks by producing a new set of weather data and applying factors of change obtained from raw daily RCM outputs between the reference period and the future.
Each climate station is perturbed using the simplest version of the CF method, that is, we applied the following deterministic transformation to downscaled climate model outputs: where PCP RCM, CF t is the perturbed RCM output, PCP RCM t is the original RCM output, and PCP OBS t denotes the observed data. Analogously for temperature, we use: where TMP RCM, CF t is the perturbed RCM output, TMP RCM t is the original RCM output, and TMP OBS t is the observed climate. The expectation is that after transformation, the means of PCP RCM, CF t and TMP RCM, CF t will each be closer to the respective means of PCP OBS t and TMP OBS t over the reference period.

Quantile-Quantile Transformation
Quantile-Quantile (QQ) transformation is an empirical routine applied to correct systematic errors in regional climate models [33,40,41]. The raw RCM outputs typically do not have the same distribution as the observations [42], and both distributions must be reconciled. It is more sophisticated than mean-based methods, and its procedure includes remapping the probability density function (PDF) of uncorrected RCM data onto the PDF of observations. Corrected RCM simulations, X CORR , of the future are produced by applying the following transformation to their cumulative distribution functions (F): where X RCM refers to the climate variable extracted from raw RCM data. This technique has been proven to be more accurate for reproducing a valid agreement between the corrected and observed PDFs than the linear method discussed earlier [29,[41][42][43]. However, one of the main drawbacks of the QQ technique is its inability to predict values beyond the original range of historic extremes [44], which may ultimately affect risk analysis.

Generating an Ensemble of Corrected-RCM-Like Realizations
The Change Factors and Quantile-Quantile downscaling methods lack the ability to produce more than one possible future state and underestimates the impact variability [45]. The variability issue is tackled in this paper by generating an ensemble of corrected RCM-like realizations, as illustrated in Figure 1, which consists of multiple runs (the number of realizations is 250) of a given stochastic weather generator fed with corrected future climate for a 30-year period. This produces a total of 2000 unique future climate projections (= 250 × 4 × 2) for each downscaling method (= 2000 × 2). Methods involving stochastic weather generators are generally capable of creating a limitless number of sequences of weather data with novel scenarios which helps with the uncertainty analysis [31]. The selection of an adequate stochastic weather generator is crucial to efficiently explore a wide range of climate risk scenarios for water resource systems [21].
This study made use of the Multi-site weather Generator of the École de Technologie Supérieure (MulGETS) [46], a multisite, multivariate weather generator that correctly reproduced historical climatic and hydrological characteristics for the SNW [27]. Wet-day precipitation sequences were reproduced using a multi-Gamma distribution (a combination of several gamma distributions). MulGETS can account for the coherency among multiple climate variables. Details of climatic data and MulGETS configuration have been previously presented in Alodah and Seidou [27]. The final result is a super-ensemble combining all ensembles containing a large number of scenarios (hereafter the "ensemble") that inherit the trends projected by climate models while covering the range of variability from the historical period. Given that the variability range of the historical period is matched, the assumption that all scenarios are equally probable will now be more defendable.

Hydrological Response
Ultimately, the development of such new future-climate scenarios, essentially driven by corrected GCM/RCM projections, enables us to discover the risks presented by changing hydrology using altered flows at a local scale. The Soil and Water Assessment Tool (SWAT-2012) was used to simulate hydrological variables and utilizing the semi-automated SUFI-2 optimization algorithm (Sequential Uncertainty Fitting ver. 2) to evaluate the most accurate simulation based on an uncertainty analysis routine [47]. SWAT is a semi-distributed, watershed-scale hydrological model that has been extensively utilized to address water quality and quantity issues [48]. Hydrologists, conservationists, and policy makers have extensively used SWAT to predict a variety of water-related issues and their environmental impacts [48][49][50][51][52][53][54][55]. Weather information is the main physically-based input that controls the transformation of precipitation into runoff in SWAT [47].
The SNW was partitioned into 31 different Hydrologic Response Units (HRUs), each with its own land, management, and soil characteristics. A detailed description of the SWAT model configuration and parameterization used in this study has been presented previously in Alodah and Seidou [27]. A set of different statistic measures was used to assess the goodness of fit of the calibrated model, including the Nash-Sutcliffe efficiency, the RMSE-observations standard deviation ratio, and the percent bias (Table 2). Hydrological impact models, distributed models in particular (e.g., SWAT), are particularly sensitive to small-scale climate variations that might be underestimated by large-scale models [56], and it could be reasonably argued that relying on a selected few models, either in a top-down or bottom-up approach, cannot be adequately justified. Table 2. Goodness-of-Fit metrics used for SWAT model performance evaluation.

Statistic Measure Formula 1
The Nash-Sutcliffe efficiency coefficient (NSE) The RMSE-observations standard deviation ratio (RSR) where O i stands for observed and P i for predicted values; O is the mean of the observed values.

Quantifying the Risk Spaces
Because the goal is to build on the bottom-up approach, the same framework as Brown et al. [57] was adopted. A climate state represented by a time series is summarized by subsets of climate and hydrological statistics that are relevant to the problem under investigation. These subsets can be calculated for any observational time series or any downscaled climate model outputs. A climate state will yield a level of risk and performance that is measured by a set of risk/performance indicators. These indicators are obtained by feeding SWAT with the perturbed time series data, each of which is unique. We placed a greater focus on hydrological risks, due to their direct impacts on the watershed systems. This includes analyzing extreme values based on statistical models and investigating the timing and intensity of peak spring flow.
2.6.1. Extreme Value Statistical Probability Models (AM and 7Q) The 3-parameter Generalized Extreme Value (GEV) and 2-parameter Weibull (WBL) distribution functions were fitted to the annual maximum streamflows (AM) and minimum extreme events (the lowest 7-day average flow, or 7Q), respectively, based on seven time-intervals (2, 5, 10, 20, 50, 100, and 500 years). The cumulative distribution functions (CDFs) of the GEV and WBL models are [58]: where ξ, µ, and σ are the shape, location, and scale parameters, respectively. Maximum likelihood (ML) estimators, as recommended by Das et al. [59], are used to compute the parameters of these statistical models based on 95% confidence intervals (α = 0.05).

Spring Flow Timing and Intensity
Another important phenomenon with significant impacts on the hydrology of many river systems in temperate regions is the timing of the annual snowmelt and changes to river-ice conditions. In the SNW, the most common flooding events are caused by snowmelt-driven runoff. Hence, future changes in the spring snowmelt and river-ice breakup, particularly at high latitudes or in large mountain regions, seem inevitable. These changes in river-ice processes are useful indicators of climate change because of their sensitivity to air temperature [60,61]. Warmer winters may promote dramatic alterations in discharge patterns and in the severity and timing of the snowmelt. A threshold amount of energy, including the mean daily temperature and soil heat flux, as well as the areal coverage of snow are critical in controlling the release of the stored water. In the SWAT model, the snowmelt starts when the daily maximum temperature (SMTMP) is above 0 • C and its uniformly released water is considered to be precipitation [48]. For the sake of simplicity, the peak spring flow is defined as the event of maximum daily flow in the spring season from January 1 to May 31.

Likelihood Estimation Using KDE
It is assumed that the synthetically generated time series will be a fair representation of the natural climate variability. Using all sample data, we employed a nonparametric density estimation function, or the kernel density estimator (KDE) to build a continuous probability density function (PDF) without making any a priori assumptions with regard to the underlying distribution. This approach involves exploring the extent of possible changes to hydroclimatic variables in the context of climate change. This includes analyzing the main characteristics of hydrological variables. The kernel density estimator of a sample of a d-variate random vector (x 1 , x 2 , . . . , x n ), drawn from an unknown distribution, is given by:f where K is a non-negative kernel smoothing function, and h the non-negative bandwidth. The bivariate kernel density estimation in this study was based on a diagonal bandwidth matrix and a Gaussian kernel [62].

SWAT Calibration and Validation
The results of the proposed method for the South Nation Watershed under changing climate conditions using a stochastic data-oriented approach are presented. Applying the two downscaling methods described above to different CORDEX data with corrected daily time series of precipitation, maximum and minimum temperature from four RCMs under the conditions of RCP4.5 and RCP8.5 emission scenarios were generated for every metrological station. The SWAT model was calibrated with 15 years of observed streamflow at a monthly resolution for the period from Jan. 1981 to Dec. 1995 and validated on the period from Jan. 1996 to Dec. 2005. Table 3 provides results of a set of statistics that show a satisfactory agreement between observed and simulated streamflows in both periods.

Time Series Generation for the Reference and Future Periods
The weather generator, MulGETS, was used to generate 250 30-year climate datasets for each climate scenario. Then, the SWAT model, calibrated on historical climate data, was forced with the corrected future climate information to investigate hydrological changes. This resulted in 4000 realizations of future daily streamflow sequences, each consisting of 30 years of predictions. The exact number of simulations is presented in Table 4. To account for uncertainties, the total of 480,000 years of predictions are summarized using extreme-streamflow indicators that describe the associated risks and allow for comparisons against observations.

Flood and Drought Indicators
Hydrologically-based extreme streamflow metrics simulated by SWAT, including the design flood and 7-day low flow based on selective return periods, were analyzed. These results were compared to the equivalent SWAT-estimated extreme indices on the historical period to quantify possible changes ( Figure 3). The isolines in Figure 3 depict the intensity of projections as the probability scale is given by the colour varying smoothly from yellow (higher intensity) to blue (less intensity). The analysis of both datasets showed that flooding events will be affected more than low flow extremes, presumably due to predicted changes in the duration of snow cover. Results of the ensemble (Figure 3e) imply significant decreases in the annual maximum flows and increasing summer minimum flows by mostly all data configurations (i.e., RCMs, RCPs, and downscaling methods). summer minimum flows by mostly all data configurations (i.e., RCMs, RCPs, and downscaling methods).

Spring Flow Variability
Increasing temperatures have a direct effect on the hydrology of the watershed by changing to the snowmelt-driven runoff that occurs in late winter/early spring. Shorter duration of snowpack and an increase in stream flow in the winter are projected, triggered by a warming winter, particularly in cold regions. This is predicted by most models and scenarios with some variations in timing and magnitude of the largest flow with more substantial warming occurring in the winter months. Figure  5 shows the disparity between the timing of the largest flow in the historical period and in the late-21st century. Between 1981 and 2010, the median date of the largest spring flow was March 12th (red dashed line in Figure 5b). The domination of most of the models and downscaling methods suggests that the peak spring flow will occur earlier, which is consistent with an earlier spring warming. As a result, the intensity of the annual spring flow at the outlet of the SNW is expected to significantly decrease by up to 60% (Figure 5a). The ensemble results suggest a 50% decrease in the quantity and the occurrence of the peak spring flow to be before the month of March.

Spring Flow Variability
Increasing temperatures have a direct effect on the hydrology of the watershed by changing to the snowmelt-driven runoff that occurs in late winter/early spring. Shorter duration of snowpack and an increase in stream flow in the winter are projected, triggered by a warming winter, particularly in cold regions. This is predicted by most models and scenarios with some variations in timing and magnitude of the largest flow with more substantial warming occurring in the winter months. Figure 5 shows the disparity between the timing of the largest flow in the historical period and in the late-21st century. Between 1981 and 2010, the median date of the largest spring flow was March 12th (red dashed line in Figure 5b). The domination of most of the models and downscaling methods suggests that the peak spring flow will occur earlier, which is consistent with an earlier spring warming. As a result, the intensity of the annual spring flow at the outlet of the SNW is expected to significantly decrease by up to 60% (Figure 5a). The ensemble results suggest a 50% decrease in the quantity and the occurrence of the peak spring flow to be before the month of March. Water 2019, 11, x FOR PEER REVIEW 12 of 18

Discussion
The main contribution of this paper is to represent the variability of future outcomes through ensembles of realizations of climate time series that are converted into ensemble of impacts. The spread of these ensembles is a representation of the uncertainty in future climate and impacts. The uncertainty comes from the fact that there is no consensus approximation of climate response to a doubling of atmospheric carbon dioxide (or climate sensitivity) although the initial estimation by National Academy of Sciences in 1979 [65] to be in the range of 1.5-4.5 °C was still adopted by the recent reports of the Intergovernmental Panel on Climate Change (IPCC) [66,67]. Andrews et al. [68] suggested that the response of global temperature might be greater and more uncertain than previously estimated, specifically at the upper end. The ensemble also covers the range of variability observed during the baseline period. The spread of the ensemble depends on the climate model, emission scenario and the choice of downscaling methods. The following sections discuss the impacts of these factors on the results.

Comparison of Downscaling Methods
By comparing downscaling methods in Figures 3 and 5, it becomes apparent that the QQ technique outperformed the mean-based CF by correcting the distributions of RCM data rather than a single statistic. The cloud of points is denser (yellow contours) where the probability is higher. QQ outcomes have a tendency to consistent clouds representing different climate models and scenarios. Indeed, the convergence in the projected QQ results across RCMs with different parameterization and climate sensitivity may suggest a more refined and skillful representation of future climate and hydrological systems, especially in the estimation of extreme daily precipitation [42]. However, in the context of risk analysis, Daniels et al. [69] stated that one cannot utterly exclude an outlier scenario as it may characterize some underlying physical processes that were missed by other models or scenarios.

Discussion
The main contribution of this paper is to represent the variability of future outcomes through ensembles of realizations of climate time series that are converted into ensemble of impacts. The spread of these ensembles is a representation of the uncertainty in future climate and impacts. The uncertainty comes from the fact that there is no consensus approximation of climate response to a doubling of atmospheric carbon dioxide (or climate sensitivity) although the initial estimation by National Academy of Sciences in 1979 [65] to be in the range of 1.5-4.5 • C was still adopted by the recent reports of the Intergovernmental Panel on Climate Change (IPCC) [66,67]. Andrews et al. [68] suggested that the response of global temperature might be greater and more uncertain than previously estimated, specifically at the upper end. The ensemble also covers the range of variability observed during the baseline period. The spread of the ensemble depends on the climate model, emission scenario and the choice of downscaling methods. The following sections discuss the impacts of these factors on the results.

Comparison of Downscaling Methods
By comparing downscaling methods in Figures 3 and 5, it becomes apparent that the QQ technique outperformed the mean-based CF by correcting the distributions of RCM data rather than a single statistic. The cloud of points is denser (yellow contours) where the probability is higher. QQ outcomes have a tendency to consistent clouds representing different climate models and scenarios. Indeed, the convergence in the projected QQ results across RCMs with different parameterization and climate sensitivity may suggest a more refined and skillful representation of future climate and hydrological systems, especially in the estimation of extreme daily precipitation [42]. However, in the context of risk analysis, Daniels et al. [69] stated that one cannot utterly exclude an outlier scenario as it may characterize some underlying physical processes that were missed by other models or scenarios.
The mean-based downscaling method used in this study provides more diverse results (Figure 3). Results, however, are scattered, suggesting that the choice of the downscaling method will greatly affect the estimated impacts and thus adaptation strategies. In the absence of evidence that one method is better than another, the safest approach is to use all available methods. Several researchers or practitioners use one downscaling method without justification, neglecting that their results will be tainted by that choice. Research is needed to assess the relative credibility of downscaling methods in order to choose the most reliable approach.

Sensitivity Domains Assessment
Based on the probabilistic assessment of climate change, the 7-day mean annual minima and annual maxima indices suggested positive changes in the future period compared to the reference period. While low streamflow is a product of a complex combination of several factors in the physical processes, a possible reason of such projected changes in low flow is attributed to enhanced evapotranspiration (due to increased air temperature) from the wet-land surfaces of the watershed in summer [70], which leads to increased precipitation frequency through more intense convective storms and consequently increasing flows during the summer months. These projected changes of volume in low flows would be welcome especially in terms of their impacts on aquatic habitat and water quality.
Higher values of future low flows compared to the baseline period due to increased rainfall intensity in North America were reported by Shrestha et al. [71] based on SWAT simulations forced by data from ten GCMs from CMIP5 archives in response to RCP4.5 and RCP8.5 scenarios in a river located in southwest Ohio, USA. Similar conclusions of extremely low flow conditions to occur less frequently in the future were reported by modeling studies applied to different watersheds across Europe [72,73] and Asia [74,75]. In fact, these projections are also in agreement with previously reported findings of upward trends in observed low flow values across Canada [76,77]. For example, Khaliq et al. [78] investigated observed values of 1-, 7-, 15-, and 30-day annual low-flow indicators in some Canadian rivers and found significantly increasing trends in southern Ontario. Similarly, Moore et al. [79] found a statistically large positive trend in streamflow during low flow season in southwest British Columbia. Moreover, Novotny and Stefan [80] found 7-day low flows to be increasing at a significant rate in three major river basins of Minnesota during both the summer and winter. In general, both lowand high-flows have been investigated in various regional studies, with no consistent conclusions in terms of how they will be altered by climate change [21]. Rather, it is believed that such extremes are very localized and dependent on the selected GCMs and hydrological models; thus results cannot be simply generalized.
On the other hand, spring peak flow is generally a product of average air temperature and snowmelt runoff. Most of the ensemble data suggest an early occurrence of peak spring flow. This result is consistent with the finding of Gunawardhana and Kazama [72], who projected an early occurrence of peak spring flow in Italy by approximately 12 days during the 2080s in comparison to the baseline period. Also, it is supported by observed and projected changes in spring streamflow-timing across parts of western North America [81]. An advance in the timing of spring peak flow could negatively affect water supply, ecosystem, and reservoirs' storage management [81]. However, such change in timing implies a reduction in snowmelt runoff severity (and hence the magnitude of the largest daily mean flows in spring), and the economic impacts will be apparent particularly in Canada [82]. These findings can likely explain the projected reduction in annual maxima values in a mainly snowmelt dominated river, where typically the highest flow occurs in early spring each year caused by snowmelt-generated runoff. These results are also consistent with previous works done in some Canadian watersheds that detected the same behaviour in the past such as Zhang et al. [83].
When multiple RCMs and downscaling approaches are used, the best hydrological estimate can be constructed based on weighting of combined results of rainfall-runoff modelling. By adopting the versatile methodology considered herein, it is believed that the risk of extremes and their consequent losses under an uncertain climate can be significantly reduced, as risk cannot be completely eliminated.
Ideally, further improvement to the super-ensemble can be achieved by considering more stochastic weather generators, climate-change models, and downscaling methods to ensure a thorough evaluation. To overcome the high computational requirements to construct the super-ensemble, further research is needed to determine the optimal length of the chain containing realizations needed to represent observed data in climatic and hydrological modelling.

Conclusions
This paper proposed a novel approach to associate a credibility measure to the climate and sensitivity measure in the bottom-up approach by combining synthetically generated climate series with downscaled climate model output to populate the climate and sensitivity spaces. The large number of projections allowed us to quantify future uncertainties and provided better coverage of the risk space based on probabilistic assessment of climate change. Furthermore, the likelihood measure provides a more practical assessment of the plausibility of future risks to serve infrastructure design and allow more confidence in water-related management decisions. Ideally, the generation of such ensembles should contain-to a feasible extent-as many uncertainty elements as possible. Given that the choice of climate models, downscaling techniques, and weather generator affects the results, a super-ensemble of future climate series was used to derive flooding and drought indices, while bearing in mind the uncertainties inherent in the extreme value modelling under a changing climate at regional scales. The results of the modelled risk analysis on the South Nation Watershed located in Ontario, Canada, indicate a marginal increase in the 7-day low flow, whereas design floods will noticeably weaken. A shorter and warmer winter is expected to result in an earlier disappearance of accumulated winter snow cover, an early onset of snowmelt, and consequently an earlier and less intensive peak spring flow. While an application was centered in one pilot watershed, the methodology delivers new insights into hydrological processes under changing climate conditions in general and will be of interest for Canada and beyond. The need is also manifest for examining a broader range of hydrologic indicators. Finally, the broader natural uncertainty posed by climate change and land use, and their implicit disruptions to various agro-socioeconomic and water sectors at an operational level, may constitute an interesting and worthwhile extension of this study.

Conflicts of Interest:
The authors declare no conflict of interest.