ASSESSING THE EFFECTS OF UNCERTAINTY AND CLIMATE CHANGE ON HYDROLOGICAL SIMULATIONS ACROSS A PERMAFROST GRADIENT IN NORTH-CENTRAL CANADA by Rajtantra Lilhare M.Tech., Indian Institute of Remote Sensing, India, 2013 B.Tech., Jawaharlal Nehru Agriculture University, India, 2011 DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN NATURAL RESOURCES AND ENVIRONMENTAL STUDIES UNIVERSITY OF NORTHERN BRITISH COLUMBIA September 2020 © Rajtantra Lilhare, 2020 ii ABSTRACT Hudson Bay, a vast inland sea in northern Canada, receives the highest average annual freshwater from the Nelson River system among all other contributing rivers. A rapidly changing climate and flow regulation from hydroelectric developments alter Nelson River streamflows timing and magnitude, affecting Hudson Bay’s physical, biological, and biogeochemical state. Despite recent developments and advances in climate datasets, hydrological models, and computational power, modelling the Hudson Bay system remains particularly challenging. Therefore, this dissertation addresses crucial research questions from the Hudson Bay System (BaySys) project by informing how climate change impacts variability and trends of freshwater-marine coupling in Hudson Bay. To that end, I present a comprehensive intercomparison of available climate datasets, their performance, and application within the macroscale Variable Infiltration Capacity (VIC) model, over the Lower Nelson River Basin (LNRB). This work aims to identify the VIC parameters sensitivity and uncertainty in water balance estimations and investigates future warming impacts on soil thermal regimes and hydrology in the LNRB. An intercomparison of six climate datasets and their equally weighted mean reveals generally consistent air temperature climatologies and trends (1981–2010) but with a prominent disagreement in annual precipitation trends with exceptional wetting trends in reanalysis products. VIC simulations forced by these datasets are utilized to examine parameter sensitivity and uncertainties due to input data and model parameters. Findings suggest that infiltration and prescribed soil depth parameters show prevailing seasonal and annual impacts, among other VIC parameters across the LNRB. Further, VIC simulations (1981–2070) reveal historical and possible future climate change impacts on cold regions iii hydrology and soil thermal conditions across the study domain. Results suggest that, in the projected climate, soil temperature warming induces increasing baseflows as future warming may intensify infiltration processes across the LNRB. This dissertation reports essential findings in the application of state-of-the-art climate data and the VIC model to explore potential changes in hydrology across the LNRB’s permafrost gradient with industrial relevance of future water management, hydroelectric generation, infrastructure development, operations, optimization, and implementation of adaptation measures for current and future developments. iv TABLE OF CONTENTS ABSTRACT ............................................................................................................................. ii TABLE OF CONTENTS .........................................................................................................iv LIST OF TABLES .................................................................................................................. vii LIST OF FIGURES ..................................................................................................................ix GLOSSARY ...........................................................................................................................xvi DEDICATION ..................................................................................................................... xviii ACKNOWLEDGEMENTS ....................................................................................................xix CHAPTER 1: INTRODUCTION ........................................................................................ 1 1.1. Dissertation Objectives and Research Questions ........................................................ 3 1.2. Background ................................................................................................................. 5 1.2.1. Climate Datasets .................................................................................................. 5 1.2.2. Uncertainties in Hydrological Modelling ............................................................ 6 1.2.3. Climate Change Impacts on Cold Regions Hydrology ....................................... 8 1.3. Dissertation Outline .................................................................................................. 10 CHAPTER 2: INTERCOMPARISON OF MULTIPLE HYDRO-CLIMATIC DATASETS ACROSS THE LOWER NELSON RIVER BASIN, MANITOBA, CANADA ... .................................................................................................................... 12 2.1. Abstract ..................................................................................................................... 12 2.2. Introduction ............................................................................................................... 14 2.3. Study Area: the Lower Nelson River Basin .............................................................. 16 2.4. Data and Methods ..................................................................................................... 19 2.4.1. Datasets .............................................................................................................. 19 2.4.2. Methods ............................................................................................................. 26 2.5. Results ....................................................................................................................... 28 2.5.1. Intercomparison of Gridded Climate Data with Station Observations .............. 28 2.5.2. Intercomparison of Precipitation and Mean Air Temperature Datasets ............ 32 v 2.5.3. 2.6. Intercomparison of Annual and Seasonal Trends (1981–2010) ........................ 42 Discussion and Conclusions ..................................................................................... 45 CHAPTER 3: SENSITIVITY ANALYSIS AND UNCERTAINTY ASSESSMENT IN WATER BUDGETS SIMULATED BY THE VARIABLE INFILTRATION CAPACITY MODEL FOR CANADIAN SUB-ARCTIC WATERSHEDS ............................................... 51 3.1. Abstract ..................................................................................................................... 51 3.2. Introduction ............................................................................................................... 53 3.3. Study Area ................................................................................................................ 56 3.4. Materials and Methods .............................................................................................. 58 3.4.1. Datasets .............................................................................................................. 58 3.4.2. The Variable Infiltration Capacity (VIC) Model ............................................... 60 3.4.3. Calibration and Evaluation ................................................................................ 62 3.4.4. Experimental Set-up and Analysis Approach .................................................... 64 3.5. Results ....................................................................................................................... 67 3.5.1. Intercomparison of the VIC Simulations ........................................................... 67 3.5.2. Uncertainty in the Water Budget Estimation ..................................................... 70 3.5.3. Model Parameter Sensitivity and Uncertainty ................................................... 75 3.5.4. Uncertainty Assessment of the VIC Model Parameters using OLH ................. 79 3.6. Discussion ................................................................................................................. 83 3.6.1. Input Data Uncertainty ...................................................................................... 83 3.6.2. Model Parameter Sensitivity and Uncertainty Assessment ............................... 85 3.7. Conclusion ................................................................................................................ 87 CHAPTER 4: WARMING SOIL TEMPERATURES AND INCREASING BASEFLOWS IN RESPONSE TO RECENT AND POTENTIAL FUTURE CLIMATE CHANGE ACROSS NORTHERN MANITOBA ................................................................... 89 4.1. Abstract ..................................................................................................................... 89 4.2. Introduction ............................................................................................................... 91 4.3. Study Area, Data, and Methodology......................................................................... 93 4.3.1. Study Area ......................................................................................................... 93 vi 4.3.2. Data and Hydrological Model ........................................................................... 95 4.3.3. Methods ........................................................................................................... 100 4.4. Results ..................................................................................................................... 102 4.4.1. VIC Performance ............................................................................................. 102 4.4.2. CMIP5 Precipitation and Air Temperature Climatology................................. 106 4.4.3. Projected Changes in VIC Simulated Variables .............................................. 109 4.4.4. Intermodel Uncertainties in the Projected Water Balance Terms ................... 112 4.4.5. Projected Changes in Soil Thermal Regime .................................................... 115 4.4.6. Long Term Projected Changes in the LNRB’s Hydrology and Soil Temperatures ................................................................................................................. 120 4.5. Discussion ............................................................................................................... 125 4.5.1. Projected Changes in the LNRB’s Hydrology................................................. 125 4.5.2. Warming Impacts on Soil Temperatures and Baseflows ................................. 128 4.6. Conclusions ............................................................................................................. 130 CHAPTER 5: CONCLUSIONS ...................................................................................... 132 5.1. Major Findings ........................................................................................................ 132 5.2. Recommendations ................................................................................................... 134 5.3. Limitations .............................................................................................................. 136 5.4. Implications of This Work ...................................................................................... 138 BIBLIOGRAPHY ................................................................................................................. 141 APPENDICES ....................................................................................................................... 164 Appendix A........................................................................................................................ 164 Appendix B ........................................................................................................................ 169 Appendix C ........................................................................................................................ 170 Appendix D ........................................................................................................................ 171 Appendix E ........................................................................................................................ 172 vii LIST OF TABLES Table 2.1 Main characteristics of the datasets used in this study. ........................................... 21 Table 2.2 List of selected weather stations for daily precipitation and mean air temperature, maintained by Environment and Climate Change Canada (ECCC) and Manitoba Hydro (MH), for the IDW interpolation (1981-2010). ....................................................................... 22 Table 2.3 Seasonal and annual precipitation and mean air temperature statistics for the domain-averaged ANUSPLIN, NARR, ERA-I, WFDEI, HydroGFD and ENSEMBLE datasets against four ECCC stations average values across the LNRB, 1981–2010. .............. 31 Table 3.1 VIC intercomparison experiments performed using different forcings (Lilhare et al., 2019). ................................................................................................................................. 60 Table 3.2 List of ten selected unregulated hydrometric stations, maintained by the Water Survey of Canada and Manitoba Hydro, for the VIC model calibration and evaluation with sub-watershed characteristics and mean annual discharge (Water Survey of Canada, 2016). 62 Table 3.3 Components of the simulated water budgets in the LNRB’s sub-watersheds, with average annual values for 1981–2010. The average annual precipitation (PCP) based on the mean of six forcing datasets, and other terms are the total runoff (TR), ET, and average soil moisture (SM), all based on the mean of VIC simulations from six different input forcing datasets. Standard deviation (SD) shows inter VIC simulations variation in the water balance estimations. .............................................................................................................................. 71 Table 3.4 Average annual precipitation (1981–2010) from six forcing datasets and standard deviation (SD) across the datasets over the LNRB’s sub-watersheds. .................................... 72 Table 4.1 CMIP5 simulations used in this study (MacDonald et al., 2018). ........................... 96 Table 4.2 List of ten selected unregulated hydrometric stations, maintained by the Water Survey of Canada and Manitoba Hydro, for the VIC model calibration and evaluation with sub-watershed characteristics and mean annual discharge (Water Survey of Canada, 2016). ............................................................................................................................................... 100 viii Table 4.3 The VIC model calibrated (1990s) values of the D1, D2, and D3 soil layer thicknesses for all sub-basins of the LNRB. .......................................................................... 103 Table 4.4 The MME mean projected seasonal and annual changes, 2030s-1990s and 2050s1990s, in VIC simulated soil temperatures at D1, D2, and D3 under RCP 8.5 across the study domain. .................................................................................................................................. 116 Table 4.5 The MME mean projected annual trend magnitudes (30 years) of VIC simulated D3 soil temperatures (°C 30-year-1) and baseflows (mm 30-year-1) under RCP 8.5 for the 1990s, 2030s, and 2050s. Trend plots are provided in Figure 4.18 (all significant). ............ 116 Table 4.6 The MME mean projected seasonal and annual trend magnitudes (all significant) (30-years) of VIC simulated soil temperatures at D1, D2, and D3 under RCP 8.5 across the study domain. Trend plots are provided in Figure 4.14. ........................................................ 119 Table 4.7 The MME mean projected seasonal and annual changes, 2030s-1990s and 2050s1990s, in VIC simulated baseflow under RCP 8.5 over all sub-basins across the study domain. .................................................................................................................................. 124 Table 4.8 The MME mean historic (1990s) and projected (2030s and 2050s) annual precipitation (P) partitioning in % within water balance components (ET, surface runoff (SR), and baseflow (BF)) under RCP 8.5 for all ten sub-basins............................................ 125 ix LIST OF FIGURES Figure 1.1 (a) Location of the Nelson River basin (NRB), Churchill River basin (CRB), and Lower Nelson River basin (LNRB). (b) Major rivers and sub-watersheds selected for the study; yellow and black triangles show the Water Survey of Canada hydrometric and Environment and Climate Change Canada (ECCC) climate stations, respectively, which are selected for the climate data analysis (Chapter 2); red diamonds denote current generating stations; the yellow circle shows a generating station under construction by Manitoba Hydro; the green box represents the Notigi Control Structure; and the red star indicates the Churchill River diversion. Base map source: World physical map (http://services.arcgisonline.com/arcgis/services) .................................................................... 3 Figure 2.1 Maps of (a) The Nelson River Basin (NRB), Churchill River Basin (CRB), and Lower Nelson River Basin (LNRB). (b) major rivers within the LNRB are labelled, black triangles show the selected ECCC stations for this study, red diamonds denote current generating stations, and the yellow circle shows a generating station under construction by Manitoba Hydro, a green box represents the Notigi Control Structure, and a red star indicates the Churchill River diversion. (c) domain elevation distribution and selected sub-watersheds (black line). .............................................................................................................................. 18 Figure 2.2 Map of selected weather stations (red triangles) for the IDW interpolation across the LNRB, 1981–2010. ............................................................................................................ 23 Figure 2.3 Comparison of the interpolated (IDW, 14 stations) and averaged (ECCC, four stations) data over the LNRB, 1981–2010. ............................................................................. 28 Figure 2.4 Area-averaged annual (a) total precipitation and (b) mean air temperature over the LNRB for the ANUSPLIN, NARR, ERA-I, WFDEI, and HydroGFD datasets compared to four ECCC stations average values across the basin, 1981–2010. .......................................... 30 Figure 2.5 Area-averaged climatology of (a) monthly total precipitation (solid lines, dotted lines for ENSEMBLE) and (b) monthly mean air temperature over the LNRB for the IDW, ANUSPLIN, NARR, ERA-I, WFDEI, HydroGFD, and ENSEMBLE datasets, 1981–2010. 33 x Figure 2.6 Spatial distribution of the annual total precipitation (mm month-1) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. ................................................................................................................................................. 35 Figure 2.7 Bias as measured against the ENSEMBLE total annual precipitation (mm month -1) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. ................................................................................................................ 36 Figure 2.8 Ensemble mean and standard deviation of (left) precipitation and (right) mean air temperature, 1981–2010. ......................................................................................................... 37 Figure 2.9 Spatial distribution of the mean annual air temperature (°C) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. .. 38 Figure 2.10 Bias as measured against the ENSEMBLE mean annual air temperature (°C) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. ................................................................................................................ 38 Figure 2.11 Area averaged seasonal (a) precipitation totals and (b) standard deviations (1981–2010) for all datasets over the ten sub-watersheds. Data are shown for each subwatershed and all seasons: DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. ............................................................................................. 40 Figure 2.12 Area averaged seasonal (a) air temperature means and (b) standard deviations (1981–2010) for all datasets over the ten sub-watersheds. Data are shown for each subwatershed and all seasons: DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. ............................................................................................. 41 Figure 2.13 Seasonal coefficient of variation (1981–2010) for all precipitation datasets over the ten sub-watersheds. Data are shown for each sub-watershed and all seasons: DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. ................ 42 Figure 2.14 Spatial trends of the annual total precipitation (mm decade-1) from different datasets, 1981–2010. The grid cells with significant trends (p < 0.05) are indicated by dots. 44 xi Figure 2.15 Seasonal and annual precipitation trends for all datasets and sub-watersheds, 1981–2010: dots (Annual); DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. Red colour symbols denote significant trends (p < 0.05). . 45 Figure 2.16 Spatial trends of the mean annual air temperature (°C decade-1) from different datasets, 1981–2010. Dots indicate grid cells with significant trends (p < 0.05). ................... 47 Figure 2.17 Seasonal and annual mean air temperature trends for all datasets and subwatersheds, 1981–2010: dots (Annual); DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. Red colour symbols denote significant trends (p < 0.05). ..................................................................................................................................... 49 Figure 3.1 Maps of the LNRB. (a) The Nelson River Basin (NRB), Churchill River Basin (CRB), and Lower Nelson River Basin (LNRB). (b) major rivers and sub-watersheds within the LNRB, yellow triangles show the hydrometric stations used in this study, white circles denote existing generating stations, and the yellow circle shows a future generating station (currently under construction) by Manitoba Hydro. A red star indicates the Churchill River diversion and the Digital Elevation Model (DEM) represents the VIC model domain at 0.10° resolution. ................................................................................................................................ 57 Figure 3.2 Schematic representation of the overall methodology adopted for the propagation, sensitivity, and uncertainty assessment in the VIC modelling over the LNRB. Coloured boxes indicate various objectives of this study. ................................................................................. 65 Figure 3.3 Boxplots for monthly calibration (a1-d1) and validation (a2-d2) performance metrics, NSE (a1-a2), KGE (b1-b2), r (simulated vs observed monthly streamflow) (p-value < 0.05 for all) (c1-c2) and PBIAS (d1-d2), for ten sub-watersheds within the LNRB based on IDW-VIC, ANUSPLIN-VIC, NARR-VIC, ERA-I-VIC, WFDEI-VIC, HydroGFD-VIC and Ensemble-VIC simulations. The black dots within each box show the mean, the red lines show the median, the vertical black dotted lines show a range of minimum and maximum values excluding outliers, and the red + signs show the outliers defined as the values greater than 1.5 times the interquartile range of each metrics. ............................................................ 68 xii Figure 3.4 Simulated and observed daily runoff (mm day-1) averaged over water years 1981– 2010 for the LNRB’s ten unregulated sub-basins. An external routing model is used to calculate runoff for the IDW-VIC, ANUSPLIN-VIC, NARR-VIC, ERA-I-VIC, WFDEI-VIC, HydroGFD-VIC and Ensemble-VIC simulations. Note that y-axis scales vary between panels. ...................................................................................................................................... 69 Figure 3.5 Spatial differences of seasonal total runoff (TR) (mm) for the LNRB’s ten unregulated sub-basins based on Ensemble-VIC (ENSEM) minus (1st row) IDW-VIC, (2nd row) ANUSPLIN-VIC, (3rd row) NARR-VIC, (4th row) ERA-I-VIC, (5th row) WFDEI-VIC, and (6th row) HydroGFD-VIC simulations, water years 1981–2010, for winter (DJF), spring (MAM), summer (JJA) and autumn (SON). ............................................................................ 73 Figure 3.6 Area-averaged VIC simulated seasonal water balance mean (mm) of precipitation (PCP, a1-a4), total runoff (TR, b1-b4), ET (c1-c4) and soil moisture (SM, d1-d4), represented by different columns, for the LNRB’s ten unregulated sub-basins, water years 1981–2010, for the winter (DJF, 1st row), spring (MAM, 2nd row), summer (JJA, 3rd row) and autumn (SON, 4th row) seasons. Red bars show the average of VIC simulations from six forcing datasets (except Ensemble), black error bars indicate inter-VIC simulation variation among IDWVIC, ANUSPLIN-VIC, NARR-VIC, ERA-I-VIC, WFDEI-VIC, and HydroGFD-VIC, while black dots represent the water balance estimation from Ensemble-VIC. ................................ 74 Figure 3.7 Ratio of factor sensitivity (%) of IVARS 50 for each parameter at annual scale over all sub-watersheds of the LNRB for the three model performance metrics (1981–2010) (a) KGE, (b) NSE, and (c) PBIAS. Ratio of factor sensitivity is estimated by normalizing IVARS 50 values in each case, so they add up to 100% for all parameters. ........................... 77 Figure 3.8 Annual streamflow sensitivity to parameter uncertainty for all ten LNRB subwatersheds. The green dots show streamflow associated with the control run (calibration), the red lines show the median, the vertical black dotted lines show a range of minimum and maximum values excluding outliers, and the red + signs show the outliers defined as the values greater than 1.5 times the interquartile range of annual streamflow. ........................... 79 xiii Figure 3.9 Average likelihood distribution of the VIC parameters using 600 samples generated after the OLH sampling over all ten sub-watersheds of the LNRB. Red bars show the maximum likelihood of parameter range to the model performance metric (KGE) for the LNRB. ...................................................................................................................................... 81 Figure 3.10 Streamflow prediction uncertainty associated with estimated parameters from the OLH. Top 10% (shown in blue colour) of OLH samples, based on KGE, used for the prediction of streamflow for all ten sub-watersheds, water years 1981–2010. Note that y-axis scales vary between panels. Shaded area (grey colour) shows the envelope of VIC runs from 600 OLH samples. ................................................................................................................... 82 Figure 4.1 Maps of the LNRB. (a) The Nelson River Basin, Churchill River Basin, and LNRB. (b) Major rivers and sub-watersheds within the LNRB, yellow triangles show the hydrometric stations used in this study, green circles denote current generating stations, the yellow circle shows a future generating station (currently under construction) by Manitoba Hydro and a red star indicates the Churchill River diversion. Basemap shows the permafrost distribution across the domain (Natural Resources Canada, 2010). Base map source: World light grey reference (http://services.arcgisonline.com/arcgis/services) .................................. 95 Figure 4.2 Schematic representation of the overall methodology adopted to quantify climate change impacts on the soil thermal regime and LNRB’s hydrology using the VIC model. ... 99 Figure 4.3 Spatial distribution of monthly calibration (1981–1985, 1995–1999) performance metrics, NSE, KGE, r (p < 0.05) and PBIAS, for ten selected sub-watersheds within the LNRB based on the HydroGFD-VIC simulations. ................................................................ 104 Figure 4.4 Spatial distribution of monthly validation (1986–1994) performance metrics, NSE, KGE, r (p < 0.05) and PBIAS, for ten selected sub-watersheds within the LNRB based on the HydroGFD-VIC simulations. ................................................................................................ 105 Figure 4.5 Simulated and observed daily runoff (mm day-1) averaged over water years 1981– 2010 for the LNRB’s ten unregulated sub-basins. An external routing model is used to calculate runoff for the HydroGFD-VIC simulations. Note that y-axis scales vary between panels. .................................................................................................................................... 106 xiv Figure 4.6 Area-averaged water year cycles of CMIP5 Multi-Model Ensemble (MME) mean daily precipitation (a, b) and air temperature (c, d) over the LNRB for the 1990s (black), 2030s (blue), and 2050s (red) (RCP 4.5, left and RCP 8.5, right). Shading corresponds to intermodel ranges................................................................................................................... 107 Figure 4.7 Projected changes (2050s-1990s; %) in the spatial distribution of precipitation during (a) winter, (b) spring, (c) summer, and (d) autumn for RCP 8.5................................ 108 Figure 4.8 Projected changes (2050s-1990s; °C) in the spatial distribution of air temperature during (a) winter, (b) spring, (c) summer, and (d) autumn for RCP 8.5................................ 109 Figure 4.9 Future changes (2050s-1990s; %) in the spatial distribution of mean annual (a) rainfall, (b) snowfall, and (c) total precipitation for RCP 8.5................................................ 110 Figure 4.10 VIC partitioned CMIP5 MME (1981-2070) mean total precipitation into (a) daily rainfall and (b) snowfall; CMIP5-VIC-simulated (c) mean ΔSWE (daily SWE rate) and (d) snowmelt, in (c) values greater than 0 represent snow accumulation, while those below 0 indicate snow ablation; VIC-simulated MME (e) mean daily runoff for the LNRB. All values are spatial averages over the basin under RCP 8.5 scenario. Units are mm day-1. Runoff units are an equivalent regional average rainfall rate rather than a discharge rate. ........................ 111 Figure 4.11 The VIC-CMIP5 MME mean daily runoff for the LNRB’s sub-basins. Black, blue, and red curves represent the daily climatology for the base period (1990s), 2050s RCP 4.5, and 2050s RCP 8.5, respectively. Shading corresponds to intermodel uncertainties. .... 112 Figure 4.12 The MME mean projected changes (%) in (a, b) ET, (c, d) surface runoff, and (e, f) streamflow at annual scale under RCP 4.5 and 8.5 scenarios for the 2030s and 2050s with respect to the base period (1981–2010). Bars show MME mean change, while error bars show intermodel variation using one standard deviation. ..................................................... 114 Figure 4.13 The VIC model simulated (red) vs observed (black) monthly mean soil temperatures (June 1994 to December 1998) at three depths (D1, D2, and D3) over one grid cell (latitude: 55.8°N, longitude: 98.7°W) within the LNRB. Black dashed line shows 0°C. ............................................................................................................................................... 117 xv Figure 4.14 The MME mean projected (a) winter (b) summer and (c) annual trends of VIC simulated soil temperatures at D1 (red), D2 (black), and D3 (blue) under the RCP 8.5 scenario over the LNRB sub-basins. The solid (dashed) trend lines show significant (insignificant) trends for the 1990s, 2030s, and 2050s and shading corresponds to intermodel uncertainties. Black dashed lines show 0°C. Trend magnitudes are provided in Table 4.6. . 118 Figure 4.15 The MME mean projected changes (reference period 1990s) (a, b) % for precipitation, (c, d) °C for air temperature, and (e, f) % for surface runoff in the 2030s and 2050s under the RCP 8.5 scenario over the LNRB. .............................................................. 121 Figure 4.16 The MME mean projected % changes (reference period 1990s) in (a, b) ET and (c, d) total soil moisture for the 2030s and 2050s under RCP 8.5 scenario over the LNRB. 121 Figure 4.17 The MME mean projected annual climatology (1990s) of VIC simulated (a) D3 soil temperature (°C) and (d) baseflow (mm) under RCP 8.5 scenario over the LNRB. The second and third rows show changes (reference period 1990s), (b, c) °C for soil temperature, and (e, f) % for baseflow, in the 2030s and 2050s, respectively. .......................................... 123 Figure 4.18 The MME mean projected annual trends of VIC simulated D3 soil temperature (red) and baseflow (black) under RCP 8.5 scenario over the LNRB sub-basins. The solid trend lines show significant trends for the 1990s, 2030s, and 2050s and shading corresponds to intermodel uncertainties. Trend magnitudes are provided in Table 4.5. ........................... 124 Figure 4.19 Schematic representation of the projected hydrology in warmer and wetter climate over the LNRB. Positive and negative signs show indicate the direction of changes, respectively in the water balance terms and thermodynamic states relative to the 1990s. .... 127 xvi GLOSSARY ANUSPLIN Australian National University spline interpolation BaySys Hudson Bay System CaPA Canadian Precipitation Analysis CMIP5 Coupled Model Intercomparison Project 5 CRB Churchill River Basin CRCM Canadian Regional Climate Model CRD Churchill River Diversion DEM Digital Elevation Model ECCC Environment and Climate Change Canada ECMWF European Centre for Medium-Range Weather Forecasts ERA-I European Reanalysis-Interim ET Evapotranspiration GLUE Generalized Likelihood Uncertainty Estimation GPCC Global Precipitation Climatology Centre HydroGFD Hydrological Global Forcing Data IAHS International Association of Hydrological Sciences IDW Inverse Distance Weighted IPCC Intergovernmental Panel on Climate Change IVARS Integrated Variogram Across a Range of Scales KGE Kling–Gupta Efficiency LAI Leaf Area Index LNRB Lower Nelson River Basin LWR Lake Winnipeg Regulation MH Manitoba Hydro MME Multi-Model Ensemble MOCOM-UA University of Arizona Multi-Objective COMplex Evolution Algorithm NARR North American Regional Reanalysis NCAR National Center for Atmospheric Research NCEP National Centers for Environmental Prediction NEMO Nucleus for European Modelling of the Ocean xvii NLDAS North American Land Data Assimilation System NRB Nelson River Basin NRCan Natural Resources Canada NSE Nash–Sutcliffe Efficiency NSERC Natural Sciences and Engineering Research Council of Canada OLH Orthogonal Latin Hypercube PBIAS Percent Bias PCP Precipitation PET Potential Evapotranspiration RCP Representative Concentration Pathway RMSE Root Mean Square Error ROS Rain-On-Snow SA Sensitivity Analysis SD Standard Deviation SM Soil Moisture SOC Soil Organic Carbon SRTM Shuttle Radar Topography Mission SWE Snow Water Equivalent TR Total Runoff VARS Variogram Analysis of Response Surfaces VIC Variable Infiltration Capacity WATCH European Union Water and Global Change WFDEI WATCH Forcing Data ERA-Interim xviii DEDICATION This work is respectfully dedicated to my father (Shri Devram Lilhare) and mother (Smt. Sunita Lilhare), who passed away during the initial and final stages of this dissertation, respectively. They are and will be a real motivation, and without their encouragement, I would have never thought of carrying out this doctoral research. Love you, Maa and Papa… I miss you so much! xix ACKNOWLEDGEMENTS This doctoral research would not have materialized without the guidance and support of many people. First and foremost, I would like to extend sincere thanks and deep regards to my supervisor Dr. Stephen J. Déry, for providing help and support in all possible ways during my course tenure and research work. His broad knowledge, suggestions, experiences, and professionalism have contributed greatly to the successful completion of this project. I want to thank my co-supervisor Dr. Tricia A. Stadnyk (University of Calgary) for being with me in all difficulties and for her constant help, suggestions, and logistical support during my site visits. I also thank the other members of my supervisory committee, Drs. Phil Owens (UNBC) and John Xiaogang Shi (University of Glasgow, UK) for their helpful comments on this dissertation. Special thanks to Dr. Francis Zwiers (Pacific Climate Impacts Consortium, University of Victoria) for his review and comments on this dissertation. I acknowledge financial and in-kind support for this research that were provided by Manitoba Hydro, ArcticNet, and the Natural Sciences and Engineering Research Council of Canada (NSERC) through the BaySys project. I thank Ms. Kristina A. Koenig, Messrs. Mark Gervais, Philip Slota, Mike Vieira, and Shane Wruth (Manitoba Hydro) for their logistical support, useful discussions, and insights into this project. I acknowledge the constructive discussions on climate data, model sensitivity, and uncertainty that I have entertained with Messrs. Scott Pokorny and Andrew Tefs (University of Manitoba). Observed hydrometric records and meteorological stations data were generously provided by Manitoba Hydro and Water Survey of Canada. I would like to thank Dr. Siraj Ul Islam (UNBC) for an in-kind contribution to the project and assistance with the VIC hydrological model. I extend my special thanks to Drs. Marco Braun (OURANOS, Québec) and Matt MacDonald (Ontario Power Generation) for graciously constructing and providing downscaled and bias-corrected future climate datasets. Both the Taylor & Francis and Wiley publishers are thanked for permitting the reproduction herein of papers published in one of their journals. Special thanks are due to the Northern Hydrometeorology Group (NHG) member and friend Dr. Aseem Raj Sharma for answering promptly the many questions, both personal and professional, I have asked him. I must thank all NHG members for their constant support and fruitful discussions during the group meetings. A huge thanks, love, and care to my loving wife, Pratiksha. Her positive attitude, love, constant care, and support are always with me. Her support and encouragement helped me a lot in overcoming various tensions and burdens during my research. I would not be withstanding this research without the help of my UNBC friends Ankush Barad, Bunu Sharma, Dhruv Desai, Giridhar Krishnan, Jeremy Morris, Mehul Solanki, Tabrez Alam, only to name a few, were also most helpful in my doctoral dissertation. I also wish to acknowledge Furquana Khan for her excellent editing skills and perspicacity. Finally, I would like to acknowledge my elder brother Mr. Swatantra Lilhare, family, relatives, and friends for their continued support during my long years as a student. 1 CHAPTER 1: INTRODUCTION Hudson Bay is a large inland sea in northern Canada that receives 94 km3 year-1 average annual freshwater from the Nelson River alone (Déry, Mlynowski, Hernández-Henríquez, & Straneo, 2011; Déry, Stieglitz, McKenna, & Wood, 2005), which is the highest among all contributing rivers of the Hudson Bay drainage system. Therefore, it is imperative to examine changing conditions within the Hudson Bay system since it is a vast inland sea that offers numerous opportunities for marine shipping and natural resource exploration and extraction. This dissertation is a part of Hudson Bay System (BaySys), a Natural Sciences and Engineering Research Council (NSERC) of Canada Collaborative Research and Development project with Manitoba Hydro as its main industrial partner. The overall purpose of the BaySys project is to examine the contributions of climate change and hydroelectric regulation to the variability and trends of freshwater-marine coupling in the Hudson Bay system. The specific objective of the BaySys project is to separate climate change and regulation effects to understand their relative contributions and impacts on Hudson Bay’s physical, biological, and biogeochemical conditions. To address the overarching goal, various sub-themes are assigned to the following teams: 1) marine and climate system, 2) freshwater systems, 3) marine ecosystem, 4) carbon cycling, 5) contaminants, and 6) Nucleus for European Modelling of the Ocean (NEMO) modelling. This Ph.D. dissertation comes under the freshwater systems team of the BaySys project. The overarching goal of the team is “To investigate the role of freshwater timing and magnitude on contemporary and future projections of freshwater-marine coupling in Hudson Bay as a means of understanding the relative contributions of regulation and climate change to the system” (Source: BaySys proposal). The team provides freshwater export outcomes to 2 other teams, allowing them to evaluate the impacts of regulation and climate change on the physical, biological, and biogeochemical processes in Hudson Bay. Moreover, results from the freshwater systems team will be beneficial to stakeholders, regulators, and industry by providing valuable information that can be used to help guide and inform the operation and development of the Nelson-Churchill hydroelectric system into the future. Future projections of net changes in runoff under climate change and flow regulation are anticipated to be highly uncertain due to the typical nonlinear characteristic behaviour of the hydrologic system (Blöschl & Zehe, 2005; Koutsoyiannis, 2010). Our study area, the Lower Nelson River Basin (LNRB), is highly susceptible to those changing scenarios (Figure 1.1). Therefore, to achieve the overall goal, our team incorporates multiple hydrological modelling approaches for robust projections of future runoff that quantify the uncertainty associated with these projections (Chen, Brissette, & Leconte, 2011; Chen, Brissette, Poulin, & Leconte, 2011; Kumar, Singh, Jena, Chatterjee, & Mishra, 2015). This partnership between researchers and industry involves detailed, physically based, multi-model streamflow projections at continental (~1 × 106 km2) and regional (~1 × 105 km2) scales. Several hydrologic models are selected, and their implementation is carried out under both naturalized and regulated conditions using multiple climate projections. I implement one of the hydrologic models, the Variable Infiltration Capacity (VIC) model (Liang, Lettenmaier, Wood, & Burges, 1994; Liang, Wood, & Lettenmaier, 1996), under naturalized conditions and simulations are performed at a regional scale (i.e., LNRB). This dissertation provides essential inputs to other BaySys team members that eventually lead to the regulation and climate change impact assessments on physical, biological, and biogeochemical processes of the Hudson Bay system. 3 Figure 1.1 (a) Location of the Nelson River basin (NRB), Churchill River basin (CRB), and Lower Nelson River basin (LNRB). (b) Major rivers and sub-watersheds selected for the study; yellow and black triangles show the Water Survey of Canada hydrometric and Environment and Climate Change Canada (ECCC) climate stations, respectively, which are selected for the climate data analysis (Chapter 2); red diamonds denote current generating stations; the yellow circle shows a generating station under construction by Manitoba Hydro; the green box represents the Notigi Control Structure; and the red star indicates the Churchill River diversion. Base map source: World physical map (http://services.arcgisonline.com/arcgis/services) 1.1. Dissertation Objectives and Research Questions Critical research gaps related to available gridded climate datasets, their performance, and utilization in the VIC model, VIC parameters sensitivity and uncertainty analyses, and the VIC model implementation for future climate change analysis exist over sub-arctic watersheds. In this context, I aim to critically analyse available gridded climate data against observations and improve our understanding on the usability of these datasets in the VIC hydrological model, parameter sensitivity, and uncertainty analysis for water balance 4 estimation over the LNRB including ten of its sub-watersheds. Future warming and extreme weather events associated with climate change impact the hydrology, snow, rain, and soil thermal regime that primarily govern the hydro-climatology of this region. However, our knowledge of the impacts of climate change on the soil thermal regime, and subsequent alteration of the sub-arctic’s hydrology, including the LNRB, remains limited (Andresen et al., 2020; Wang et al., 2019; Westermann et al., 2016). Thus, this dissertation addresses three main objectives as follows: 1. to intercompare hydro-climatic datasets, their variability, and trends across the LNRB and its ten sub-watersheds (Chapter 2); 2. to perform uncertainty propagation from multiple climate datasets in the VIC model outputs, model parameter sensitivity, and uncertainty in VIC simulated water budgets (Chapter 3); and 3. to quantify climate change impacts on the LNRB’s hydrology and soil thermal regime (Chapter 4). The seasonal variability, upstream flows, and the amount of water released from the LNRB’s snowpacks provide vital information for water management and hydropower generation. Assessing how the LNRB’s water resources may change in the future supports land and water resource managers, planners, and government to make informed decisions and requires scientific attention. Therefore, this dissertation addresses the following research questions in detail for our study domain: 1. What are the suitable datasets for the hydro-climatic study of the LNRB? How do historic (1981–2010) climatologies and trends vary across the available climate datasets? 5 2. What is the contribution of uncertainty propagated through multiple climate datasets in water balance estimations using the VIC model? What are the most sensitive model parameters in water balance modelling and parameter uncertainty to streamflow generation across the LNRB sub-basins? 3. What are the impacts of a changing climate on the LNRB’s water balance? How do projected climate changes affect soil thermal regimes and baseflows across the domain? 1.2. Background 1.2.1. Climate Datasets Recent high-resolution regional and global climate datasets that incorporate various observational sources and sometimes numerical modelling provide improved inputs for hydro-climatic modelling studies. It remains challenging, however, to obtain appropriate and credible climate data for relatively small and remote areas with sparse climate observing stations. Therefore, understanding the spatial and temporal accuracies of such gridded datasets at watershed and sub-watershed levels is essential for modelling hydrologic responses to climate change studies. Previous evaluations among precipitation and temperature products have focused mainly on the intercomparison of (1) satellite-derived products (Skok et al., 2016; Turk, Arkin, Sapiano, & Ebert, 2008), (2) climate model simulations (Mearns et al., 2013), and (3) reanalysis datasets (Bosilovich, Chen, Robertson, & Adler, 2008; Kim, Kim, Boo, Shim, & Kim, 2019). However, in situ measurements (i.e., precipitation) may contain errors arising from wind effects (Kochendorfer et al., 2017), undercatch (Mekonnen, Matula, Doležal, & Fišák, 2015), evaporation (Leeper & Kochendorfer, 2015), human errors and instrumental problems (Dahri et al., 2018). Steps 6 toward recognizing limitations in various precipitation and temperature observations were undertaken in many previous studies (Eum, Dibike, Prowse, & Bonsal, 2014; Maggioni, Sapiano, Adler, Tian, & Huffman, 2014; Shen, Xiong, Wang, & Xie, 2010). Numerical climate and weather prediction models including atmosphere-ocean General Circulation Models (GCMs) and Regional Climate Models (RCMs) offer another potential source of precipitation and air temperature estimates. Climate variables from these models, however, produce systematic biases due to inadequate model conceptualization, discretization, spatial averaging within grid cells and remain relatively coarse in resolution (Teutschbein & Seibert, 2010; Xu, Widén, & Halldin, 2005). Access to reliable and accurate atmospheric data, especially precipitation and air temperature, remains essential to understand the climate system and hydrological processes. Reliable precipitation and air temperature measurements provide constructive information to climatologists, meteorologists, hydrologists and other decision-makers in various applications. Prior to the application of precipitation and air temperature from different available products, intercomparison of these variables is essential to ensure their reliability, particularly for specific spatial and temporal domains (Kim et al., 2019; Tanarhte, Hadjinicolaou, & Lelieveld, 2012). 1.2.2. Uncertainties in Hydrological Modelling The input forcing analysis is an essential part for any hydrological modelling related studies; along with this, many studies have their main interest either on hydrological model structure (Jiang et al., 2007; Poulin, Brissette, Leconte, Arsenault, & Malo, 2011; Velázquez et al., 2013; Wilby & Harris, 2006) or on calibration parameters (Bennett, Werner, & Schnorbus, 2012; Teutschbein, Wetterhall, & Seibert, 2011). Renard, Kavetski, Kuczera, 7 Thyer, & Franks (2010) identified the total predictive uncertainty in hydrological models and found that model inputs and structural components significantly affect predicted flows. Moreover, Arsenault & Brissette (2014) estimated parameter set uncertainty in hydrologic modelling over Québec and concluded that model implementation and flow estimates are particularly sensitive to parameter set selection. This shows that a hydrological model can have multiple equivalent local optima within a model parameters range (Poulin et al., 2011) and several different parameter sets may reflect the same “optimal” efficiency measures during the optimization process (Beven, 2006). Thus, Islam & Déry (2017) analysed parameter uncertainties involved in the VIC model calibration process and modelling uncertainties related to the selection of the calibration and validation time periods. They found that choice of an initial parameter range during the optimization process plays a crucial role and different climatic conditions during calibration and validation processes result in biases within the model setup. Biemans et al. (2009) assessed the uncertainty in discharge calculations associated with various precipitation datasets for 294 river basins worldwide and found substantial differences, especially in mountainous, arctic, and small basins, between seven global gridded precipitation datasets. Although recent developments and advances have been achieved in hydrological modelling and computational power (Fatichi et al., 2016; Singh, 2018), efficiently addressing the uncertainties in hydrological simulations remains a critical challenge (Liu & Gupta, 2007). There is a growing need for sensitivity and uncertainty assessments associated mainly with the model parameterization and input forcing datasets to achieve the hydrological model’s optimal performance for decision-making. While there may be other uncertainties (e.g., model structure, calibration, etc.), this doctoral research focuses primarily on the 8 sensitivity of hydrological model parameters to model outputs, and uncertainties in model simulations due to model parameters and input forcing datasets, which are perhaps the most significant source of uncertainty (Zhang, Li, Huang, Wang, & Cheng, 2016). 1.2.3. Climate Change Impacts on Cold Regions Hydrology Climate change profoundly impacts terrestrial freshwater export to the coastal ocean (Haddeland, Skaugen, & Lettenmaier, 2006; Schewe et al., 2014). There is strong evidence that climate change impacts are altering the water cycle and influencing water availability and demand (Haddeland et al., 2006; Liu, Tian, Hu, & Sivapalan, 2014). Climate change impact assessment on the water cycle has been extensively studied using physically based hydrologic and statistical models at the global, continental, sub-continental, watershed and local scales (Arora & Boer, 2001; Elsner et al., 2010; Haddeland et al., 2014; Islam, Déry, & Werner, 2017; Mishra & Lilhare, 2016). In Canada, several studies have examined climate change impacts on Churchill-Nelson River water yields; for example, Westmacott and Burn (1997) concluded that the magnitude of hydrologic events decreased over time while snowmelt runoff events occurred earlier, and the timing of a hydrologic event was highly sensitive to temperature change. Moreover, Sushama et al. (2006) examined projected climate change signals using simulations from two versions of the Canadian Regional Climate Model (CRCM) and found an increase in the average annual precipitation and a significant decrease in snow cover over the Churchill-Nelson River basin. Further, Poitras et al. (2011) reported projected increases in 2041–2070 winter streamflows for several northern Canadian watersheds and found significant changes in extreme flows (low and high) for high-latitude basins of western Canada including the Churchill-Nelson River watershed. 9 Permafrost comprises about a quarter of the exposed land in the Northern Hemisphere (Zhang, Barry, Knowles, Heginbottom, & Brown, 1999). Changes in permafrost distribution and summer active layer depth would have significant impacts on hydrology, soil organic carbon (SOC), vegetation distribution, and infrastructure at high latitudes (ACIA, 2004; Anisimov et al., 2001; Chapin et al., 1992; Nelson, 2003). For example, permafrost disappearance or deepening summer active layer depth could alter terrain and hydrologic conditions (Hinzman et al., 2005), affect growth and distribution of vegetation (Jorgenson, Racine, Walters, & Osterkamp, 2001), and increase SOC decomposition and enhance CO2 and CH4 emissions from the soil to the atmosphere (Goulden et al., 1998; Oechel et al., 1993). Climate models project large increases in mean surface temperature (~8°C) across present-day permafrost areas of the Canadian landmass by the end of the 21 st century under the Representative Concentration Pathway (RCP) 8.5 scenario (Koven, Riley, & Stern, 2013). While this enhanced warming affects permafrost temperatures and conditions (Chadburn et al., 2017; Guo & Wang, 2016; Slater & Lawrence, 2013), it is challenging to project permafrost extent reductions from climate models due to inadequate representation of soil properties, including ice content, and uncertainties in understanding the response of deep permafrost. Simulations from a process-based permafrost model driven by six GCMgenerated climate scenarios, considering deeper permafrost, project that the area underlain by permafrost in Canada will decline by approximately 16-20% by 2090, relative to a 1990s baseline (Zhang, Chen, & Riseborough, 2008a). Further simulations also show that permafrost thaw would continue during the late 21 st century, even if air temperatures stabilize by mid-century (Zhang, Chen, & Riseborough, 2008b). 10 Overall, high latitudes are more susceptible to climate change and air temperature has increased at a higher rate than the global mean over that region during the 20th century (Serreze et al., 2000). Most GCMs project the same pattern will continue in the 21st century in high latitude regions (McCarthy, Canziani, Leary, Dokken, & White, 2001). Therefore, quantifying changes in soil thermal regime under a warming climate is important for assessing its impacts on baseflows, water and energy fluxes, SOC, ecosystems, and northern communities and infrastructure in arctic and sub-arctic regions such as the LNRB, and for including their feedbacks in GCMs for projecting future climate change. Climate change impacts on the LNRB’s water resources require a systematic analysis of the GCM-driven hydrological model simulations that has not been achieved so far. This study, therefore, investigates increasing air temperatures and changing precipitation impacts on the LNRB’s water balance, snowmelt, soil temperatures, baseflows, and its contribution to runoff generation. 1.3. Dissertation Outline To accomplish the above-mentioned objectives, this dissertation is organised into five chapters. Except for Chapters 1 and 5, all chapters independently address one or more objectives. Chapters 2 and 3 are written as stand-alone journal articles. Chapters 2 and 3 of this dissertation have been published in Atmosphere-Ocean (Lilhare, Déry, Pokorny, Stadnyk, & Koenig, 2019) and Hydrological Processes (Lilhare, Pokorny, Déry, Stadnyk, & Koenig, 2020), respectively, while Chapter 4 is written as a manuscript for future submission to a journal. Each of these three chapters is presented as a stand-alone effort with its own introduction, study area, methodology, results, discussion, and conclusion sections. There are slight modifications in these chapters that include a combined list of references from all 11 chapters into one bibliography, renumbering and positioning of the figures and tables, english and layout reformatting. Chapter 2 (Lilhare et al., 2019) covers a detailed intercomparison of multiple hydroclimatic datasets, available at global and regional scales, across the LNRB. Moreover, it includes the 1981–2010 variability and trend analyses of climatic variables such as precipitation and air temperature on annual and seasonal bases over the LNRB’s sub-basins. In Chapter 3 (Lilhare et al., 2020), I quantify the uncertainty propagated from available forcing datasets in their surface water balance estimations, using the VIC model, over the LNRB during 1981–2010. Further, this chapter presents the VIC model parameter sensitivity in simulating seasonal and annual streamflow and assesses uncertainty envelopes in streamflow generation across the LNRB sub-basins. Chapter 4 describes the impacts of climate change using output from GCMs on the LNRB’s hydrology and its major water balance components. Furthermore, it explores projected climate impacts on the LNRB’s soil thermal regime and water balance by the end of the 2070s under the RCP 4.5 and 8.5 scenarios. Chapter 5 concludes the dissertation by synthesizing the previous chapters, provides recommendations for future research, and outlines some limitations and the broader implications of this study. 12 CHAPTER 2: INTERCOMPARISON OF MULTIPLE HYDRO-CLIMATIC DATASETS ACROSS THE LOWER NELSON RIVER BASIN, MANITOBA, CANADA Publication details: This chapter has been published in Atmosphere-Ocean. Lilhare, R., Déry, S. J., Pokorny, S., Stadnyk, T. A., & Koenig, K. A. (2019). Intercomparison of multiple hydro-climatic datasets across the Lower Nelson River Basin, Manitoba, Canada, Atmosphere-Ocean, 57(4), 262–278. https://doi.org/10.1080/07055900.2019.1638226 2.1. Abstract This study evaluates the 1981–2010 spatio-temporal differences in six available climate datasets (daily total precipitation and mean air temperature) over the Lower Nelson River Basin (LNRB) in ten of its sub-watersheds at seasonal and annual timescales. We find that the Australian National University spline interpolation (ANUSPLIN), and Inverse Distance Weighted (IDW) interpolated observations from 14 Environment and Climate Change Canada (ECCC) meteorological stations show dry biases, whereas reanalysis products tend to overestimate precipitation across most of the basin. All datasets exhibit prominent disagreement in precipitation trends whereby the European Reanalysis-Interim (ERA-I) and European Union Water and Global Change (WATCH) Forcing Data ERA-Interim (WFDEI) show exceptional wetting trends, while the IDW and ANUSPLIN manifest drying trends. Mean air temperature trends generally agree across most of the datasets; however, the North American Regional Reanalysis (NARR) and IDW show stronger warming relative to other datasets. Overall, analyses of the different climate datasets and their ensemble reveal that the choice of input dataset plays a crucial role in the accurate estimation of historical climatic conditions, particularly when assessing trends, for the LNRB. Using the ensemble has the 13 distinct advantage of preserving the unique strengths of all datasets and affords the opportunity to estimate the uncertainty for hydrologic modelling and climate change impact studies. 14 2.2. Introduction Access to reliable and accurate atmospheric data, especially precipitation and air temperature, remains essential to understand the climate system and hydrological processes. These variables form vital elements of the water and energy cycles and are key for driving land surface models. Reliable precipitation and air temperature measurements provide constructive information to climatologists, meteorologists, hydrologists and other decisionmakers in various applications such as in agricultural and environmental research (e.g., Hively, Gérard-Marchant, & Steenhuis, 2006; Zhang, Sun, Singh, & Chen, 2012), climate and/or land use change studies (e.g., Cuo et al., 2011; Dore, 2005; Huisman et al., 2009), hydrological and water resources planning (e.g., Hong, Adler, Huffman, & Pierce, 2010; Middelkoop et al., 2001), and mitigation of natural hazards (e.g., Blenkinsop & Fowler, 2007; Kay, Davies, Bell, & Jones, 2009; Taubenböck et al., 2011). Prior to the application of precipitation and air temperature from different available products, intercomparison of these variables is an essential prerequisite to ensure their reliability, particularly for specific spatial and temporal domains. In recent decades, multiple global and regional datasets have been generated using different forcing products, such as climate model simulations and interpolations of remotely sensed and/or in situ observations (Dee et al., 2011; Mesinger et al., 2006). These datasets tend to systematically agree over the major temporal trends and spatial distribution of climate variables (i.e., precipitation and air temperature), but they often show notable differences at regional scales (Adler, Kidd, Petty, Morissey, & Goodman, 2001; Costa & Foley, 1998). Precipitation and air temperature stations provide direct physical readings with relatively accurate measurements at specific points. Moreover, these station measurements are often 15 spatially interpolated and thus are unable to capture the actual spatial variability of precipitation and temperature fields due to sparse station networks. Satellite-based precipitation estimates provide good spatial coverage but have inaccuracies resulting primarily from instrumental, temporal sampling, and algorithm errors (Gebremichael, Krajewski, Morrissey, Huffman, & Adler, 2005; Nijssen & Lettenmaier, 2004). Reanalyses provide an alternative source of precipitation and air temperature by assimilating all available data (meteorological stations, aircraft, satellite, etc.); however, the accuracy of these products depends on the choice of physical parameterizations and specific analysis-forecast systems (Betts, Ball, & Viterbo, 2003). Several gridded datasets for precipitation and air temperature – based on available observations, post-processing techniques and sometimes climate modelling – are available for hydro-climatic studies over the Canadian landmass (Hopkinson et al., 2011; Mesinger et al., 2006). Since observational data are incorporated to derive the gridded datasets, they may also contain measurement errors and missing records. These missing values translate into the data interpolation and add to the overall uncertainty in the resulting gridded products (Eum et al., 2014; Horton, Schaefli, Mezghani, Hingray, & Musy, 2006; Kay et al., 2009). Understanding the observed climate trends from different gridded datasets for a river basin is essential to assess regional changes and gauge data consistency and reliability. Several datarelated difficulties arise when attempting to analyse Canadian climate trends at regional and national scales. For example, its vast land mass and high latitude result in a wide range of climates and relatively large spatial climate variability (Phillips, 1990). Given the discrepancies in available gridded datasets, errors during data development, and the importance of data intercomparison prior to hydro-climatic application, our paper 16 presents a detailed intercomparison of variability and trends calculated from available global and regional air temperature and precipitation datasets over the Lower Nelson River Basin (LNRB). Intercomparisons are quantified for ten sub-basins of the LNRB on seasonal and annual bases. Further, no studies have explored the spatio-temporal differences and historic trends from different gridded datasets for the LNRB. Thus, this paper seeks to fill an important knowledge gap in the understanding of historic trends and climate datasets application for future hydrologic modelling studies over the LNRB. 2.3. Study Area: the Lower Nelson River Basin The Nelson River Basin (NRB) is one of the largest river systems in Canada (third largest by area and volumetric discharge to the coastal ocean) that drains water mainly from the interior of Canada, cutting through the Canadian Shield of northern Manitoba before emptying into Hudson Bay (Figure 2.1a) (Newbury & Malaher, 1973). The Churchill River system covers the north-western part of the NRB and is considered here since it was joined to the Nelson River by a diversion in 1976. The entire Nelson-Churchill River Basin extends geographically between ~45.5°N to 59.5°N, and ~90.0°W to 117.5°W. In this study, the downstream segment of the Nelson River system fed by the Churchill River basin and Lake Winnipeg constitute the LNRB (Figure 2.1). The LNRB spans an area of ~90,500 km2 and collects all water from the drainage area upstream of the Nelson River (~970,000 km 2) before discharging into Hudson Bay. In the LNRB, the main stem river (Nelson) and its largest tributary – the Burntwood, whose downstream segment carries diverted flows from the Churchill River – have less seasonal flow variability due to streamflow regulation and a large drainage area. Most of the LNRB has gentle slopes, with common channelized lakes moderating flow variability. Wetlands abound within the LNRB and store significant 17 volumes of water, cover large areas and moderate streamflow responses to rainfall and snowmelt events. Shallow soils and permafrost limit infiltration, groundwater storage and groundwater flows. To increase its hydroelectric capacity, Manitoba Hydro manages flows in the LNRB with two major sources of streamflow regulation: the Churchill River Diversion (CRD) and Lake Winnipeg Regulation (LWR) (Figure 2.1b). On the LNRB’s north-western boundary, Manitoba Hydro operates the CRD. In 1977, a portion of the Churchill River Basin was diverted into the LNRB and regulated at Notigi Lake by the Notigi Control Structure on the Rat River. In 1972, Manitoba Hydro started the LWR project, which is key to hydropower development on the Nelson River system. Presently, Manitoba Hydro operates six hydroelectric generating stations while one station remains under construction (Keeyask) (Figure 2.1b) within the LNRB (Déry, Stadnyk, MacDonald, Koenig, & Guay, 2018). 18 Figure 2.1 Maps of (a) The Nelson River Basin (NRB), Churchill River Basin (CRB), and Lower Nelson River Basin (LNRB). (b) major rivers within the LNRB are labelled, black triangles show the selected ECCC stations for this study, red diamonds denote current generating stations, and the yellow circle shows a generating station under construction by Manitoba Hydro, a green box represents the Notigi Control Structure, and a red star indicates the Churchill River diversion. (c) domain elevation distribution and selected sub-watersheds (black line). The LNRB experiences a sub-arctic continental climate characterized by moderate precipitation and humidity, cool summers, and cold winters. The snow-free season remains brief, generally beginning in May and ending in October. Most of the precipitation that occurs during the summer months falls as rain, accounting nearly two-thirds of the annual total precipitation. The precipitation peaks in July, the warmest month of the year with an average daily air temperature of 16.2°C. The average annual precipitation over the LNRB totals ~500 mm while evapotranspiration in the region attains ~300–350 mm annually with 19 more annual surface water evaporation (~450 mm) (Environment and Climate Change Canada, 2016; Smith, Delavau, & Stadnyk, 2015). 2.4. Data and Methods 2.4.1. Datasets The observation-based ANUSPLIN, NARR, ERA-I, WFDEI, and HydroGFD datasets are available for hydro-climatological studies over the LNRB (Table 2.1).These datasets were selected because they are derived using advanced interpolation (for ANUSPLIN) and data assimilation (for NARR, ERA-I, WFDEI, and HydroGFD) techniques, are open source, cover the entire study period (1981–2010) and domain, and have continuous records at the required temporal and spatial resolutions for hydrological modelling. However, the data assimilation capacity of global and regional reanalysis products over different parts of the world does not extend to Canada, where weather station coverage is much lower and not assimilated into many real-time reanalysis updates, especially in northern Canada (Bukovsky & Karoly, 2007; Essou, Sabarly, Lucas-Picher, Brissette, & Poulin, 2016; Langlois et al., 2009). To compare these products against observations, a gridded dataset (IDW) from 14 Environment and Climate Change Canada (ECCC) meteorological stations is constructed across the LNRB using a squared IDW interpolation technique (Table 2.2 and Figure 2.2). Other observation-based datasets such as the Canadian Precipitation Analysis (CaPA) are also available for an intercomparison but CaPA contains records from 2002 onwards and so remains unsuitable for the present study (Fortin et al., 2018; Mahfouf, Brasnett, & Gagnon, 2007). 20 The observation-based daily gridded dataset, ANUSPLIN, developed by Natural Resources Canada (NRCan) is available for the Canadian landmass south of 60°N at 10 km resolution (Hopkinson et al., 2011; McKenney et al., 2011; Natural Resources Canada, 2014). Preliminary analysis shows the updated version (1950–2015) of ANUSPLIN exhibits a consistent dry bias in precipitation over the LNRB; therefore we retained the earlier version (1950–2011) of this dataset (Hopkinson et al., 2011). This dataset uses a trivariate thin-plate smoothing spline technique and includes daily total precipitation (mm), minimum and maximum air temperature (°C) at 10 km spatial resolution based on 7514 meteorological stations (1950–2011) over the Canadian landmass south of 60°N (Eum et al., 2014; Sharma & Déry, 2016). 21 Table 2.1 Main characteristics of the datasets used in this study. Description Temporal and Spatial Resolution Geographical Coverage ECCC The Environment and Climate Change Canada meteorological stations (Figure 2.1a): Norway House A (506B047), Flin Flon A (5050960), Gillam A (5061001), Thompson A (5062922) Daily (1981–2010), Point data LNRB, Land only IDW Inverse Distance Weighted interpolated observations from 14 ECCC meteorological stations (Table 2.2 and Figure 2.2) (Gemmer, Becker, & Jiang, 2004; Shepard, 1968) Daily (1981–2010), 0.10° LNRB, Land only ANUSPLIN The Australian National University spline interpolation (Hopkinson et al., 2011; Natural Resources Canada, 2014) Daily (1950–2011), 0.10° Canada, Land only NARR North American Regional Reanalysis (Mesinger et al., 2006) 3-hourly and daily (1979 to near present), 0.30° North America, Land and ocean ERA-I European Reanalysis-Interim (Dee et al., 2011) 3, 6-hourly and daily (1979 to near present), 0.125° Global, Land and ocean WFDEI European Union Water and Global Change (WATCH) Forcing Data ERA-Interim (Weedon et al., 2014) 3, 6-hourly and daily (1979 to near present), 0.50° Global, Land and ocean HydroGFD Hydrological Global Forcing Data (Berg et al., 2018) 3, 6-hourly and daily (1979 to near present), 0.50° Global, Land only ENSEMBLE Average of above mentioned six gridded datasets Daily (1981–2010), 0.10° LNRB Dataset (Reference(s)) 22 Table 2.2 List of selected weather stations for daily precipitation and mean air temperature, maintained by Environment and Climate Change Canada (ECCC) and Manitoba Hydro (MH), for the IDW interpolation (1981-2010). Stations (Operating agency) Climate IDs Latitude Longitude Elevation (°N) (°W) (m) Data availability Cross Lake Jenpeg (ECCC) 5060623 54.54 98.03 218.8 1972–present Flin Flon A (ECCC) 5050960 54.68 101.68 303.9 1968–2017 Gillam A (ECCC) 5061001 56.36 94.71 145.1 1970–2014 Grand Rapids Hydro (MH) 5031111 53.16 99.28 222.5 1966–present Island Lake A (ECCC) 5061376 53.85 94.65 235.6 1970–2015 Kelsey Dam CS (MH) 5061422 56.04 96.51 182.9 2000–present Kelsey Hydro (MH) 5061423 56.03 96.53 182.9 1983–1995 Limestone GS (MH) 5061567 56.50 94.18 88.4 1985–1992 Norway House A (ECCC) 506B047, 53.95 97.85 223.7 1973–2010 Snow Lake (ECCC) 5062706 54.88 100.03 295.7 1983–1998 South Indian Lake (ECCC) 5062734 56.78 98.97 259.1 1976–1989 The Pas A (ECCC) 5052880 53.97 101.10 270.4 1943–2014 Thompson A (ECCC) 5062922 55.80 97.86 224.0 1967–2014 Wabowden (ECCC)1 5063041, 5063044 54.92 98.63 232.9 1982–2008 5062045 The NARR product is a high resolution, regional extension of the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) global reanalysis data (Kalnay et al., 1996; Kistler et al., 2001). It is developed at 32 km spatial and 3-hourly temporal resolution by utilizing a version of the Eta Model and its 3D variational data assimilation system (Mesinger et al., 2006) for the North American continent, available from 1979 to present. 1 Records spliced from multiple stations 23 The ERA-I dataset is a global reanalysis product from the European Centre for MediumRange Weather Forecasts (ECMWF). Originally developed at 0.8° spatial resolution, data are also available for download at different spatial (0.125° to 3° grids) and 3-hourly temporal resolutions for January 1979 through near real-time (Dee et al., 2011). In this application, the ERA-I products are downloaded at 0.125° (~13 km) spatial resolution. The product combines observations with a prior estimate of the atmospheric state generated by a global forecast model in a statistically optimal way. The ERA-I datasets have been evaluated and widely used in a variety of studies related to pan-Arctic hydro-climatology (Betts & Beljaars, 2017; Bromwich et al., 2016; Simmons et al., 2014; Simmons & Poli, 2015). Figure 2.2 Map of selected weather stations (red triangles) for the IDW interpolation across the LNRB, 1981– 2010. 24 The WFDEI relies on a method developed by the European Union’s WATCH project (http://www.eu-watch.org) and incorporates in situ observations in the reanalysis (Weedon et al., 2014). The WFDEI derived from the ERA-I (Dee et al., 2011) and was improved by an elevation correction for numerous variables. Further, to retain the monthly statistics similar to in situ observations of the Global Precipitation Climatology Centre (GPCC) (Schneider, Fuchs, Meyer-Christoffer, & Rudolf, 2008), an undercatch correction is adopted whereby the daily variability of the reanalysis product is conserved (Weedon et al., 2014). We obtained the WFDEI-GPCC precipitation and mean air temperature data at ~55 km spatial and daily temporal scale for this study. The recently developed HydroGFD dataset combines different reanalysis datasets to produce near real-time, updated hydrologic forcing of precipitation and air temperature (Berg, Donnelly, & Gustafsson, 2018). The HydroGFD resembles the established WFDEI method but uses updated climatological observations and interim products to produce nearreal time estimates of precipitation and air temperature at 3- and 6-hourly temporal and 0.50° spatial resolutions (Berg et al., 2018). The IDW dataset of daily precipitation and mean air temperature derives from 14 ECCC meteorological stations using the WATFLOOD interpolation tool (Kouwen, 1988) (Table 2.2 and Figure 2.2). These observation stations are spatially interpolated by applying the IDW interpolation method, and gridded datasets are procured at 0.10° (~10 km) horizontal resolution for the LNRB. The WATFLOOD model is optimized through a calibration procedure over the LNRB with specified parameter values for the radius of influence (241.2 km), smoothing distance (52.8 km), and precipitation and air temperature vertical lapse rates (0 mm m-1 and -0.004°C m-1, respectively). 25 The NARR, ERA-I and HydroGFD daily precipitation and mean air temperature are obtained from the sum and average, respectively, of 3-hourly values for a 24-hour period. To simplify the datasets intercomparison and to provide consistent hydrological model inputs for future work over the LNRB, the NARR, ERA-I, WFDEI and HydroGFD are then regridded to 10 km (~0.10°) spatial resolution using bilinear interpolation. The NARR (32 km) curvilinear grid and the ERA-I, WFDEI, and HydroGFD Gaussian grids are interpolated from coarser resolution to higher resolution (10 km). No elevation correction during the interpolation from coarser to finer spatial resolutions is performed as elevations vary no more than ±10% in the study area; hence regridding of the NARR (32 km), ERA-I (~13 km), WFDEI (~55 km) and HydroGFD (~55 km) datasets to 10 km spatial resolution results in negligible elevation-dependent uncertainty. LNRB grid cells exhibit almost no difference in orography and none of these datasets contain discontinuous flux fields (precipitation and air temperature); therefore, atmospheric variables (i.e., average air temperature) and basin elevation remain nearly identical at both spatial resolutions. Spatially regridded datasets (IDW, ANUSPLIN, NARR, ERA-I, WFDEI and HydroGFD) at daily temporal and 10 km spatial resolutions are assembled to produce an ensemble mean dataset (i.e., ENSEMBLE hereafter) from 1981 to 2010. For the multi-product ensemble dataset, daily precipitation and mean air temperature are derived from the equally-weighted average of all six gridded products. We assign equal weights to each dataset and so ignore prior knowledge of their modelling capacity. This is one of the most commonly used methods, where it assumes that the equally weighted ensemble provides the best estimates of contemporary and future conditions, since each model is equally likely to represent the truth (Krishnamurti et al., 1999; Suh et al., 2012). The concept of a multi-product ensemble has 26 been used widely over global and regional domains to examine uncertainty in forcings and provide input to land surface models under historical and projected future climate conditions (Fowler, Ekström, Blenkinsop, & Smith, 2007; Fowler & Kilsby, 2007; Mishra & Lilhare, 2016; Wang, Bohn, Mahanama, Koster, & Lettenmaier, 2009). 2.4.2. Methods Observational data for an intercomparison are obtained from the average of four ECCC meteorological stations (ECCC hereafter) within the LNRB (Table 2.1 and Figure 2.1). The ECCC data with adjusted, continuous and homogenous climatic records over 1981–2010 represent domain-averaged observations of mean annual air temperature and total precipitation (Figure 2.3). They may not represent the ‘best-guess’ of precipitation and air temperature but an additional comparison of ECCC data provides acceptable results against interpolated (i.e., IDW) values (Figure 2.3). We examine the Root Mean Square Error (RMSE), Percent Bias (PBIAS) and bias (air temperature only) for the long-term seasonal and annual precipitation and air temperature error statistics. Analyses are performed for all four seasons: winter (DJF), spring (MAM), summer (JJA) and autumn (SON). Daily total precipitation (mean air temperature) from all datasets are averaged and aggregated to annual and seasonal totals (means) for each grid cell. Temporal and spatial analyses of the mean aggregated gridded data are then performed. As the first step of temporal analysis, the mean of all grid cells for each year and its seasons is calculated for all datasets; for spatial analyses, calculations are performed at the grid scale. The long-term seasonal variability is measured by the standard deviation (SD) for precipitation and mean air temperature. 27 For trend analyses, precipitation and air temperature datasets are verified for autocorrelation and pre-whitened if any serial correlation exists using the method of Yue, Pilon, Phinney, & Cavadias (2002). Further, we used the Mann-Kendall, non-parametric test to estimate trends in the total seasonal and annual precipitation and mean air temperature (Kendall, 1975; Mann, 1945; Wilks, 2011). The Mann-Kendall trend test follows the relative ranking of the data rather than model fitting, and is resistant to outliers, non-normal and dominant data. This test remains robust for trend detection in similar hydro-climatic studies (Burford, Déry, & Holmes, 2009; Déry & Brown, 2007; Déry & Wood, 2005b; Gan & Kwong, 1992; Gocic & Trajkovic, 2013; Modarres & Sarhadi, 2009). The trend magnitudes are calculated using the Theil-Sen trend estimate and the statistical significance is determined at the 95% confidence level (Mondal, Kundu, & Mukhopadhyay, 2012; Sen, 1968). 28 Figure 2.3 Comparison of the interpolated (IDW, 14 stations) and averaged (ECCC, four stations) data over the LNRB, 1981–2010. 2.5. Results 2.5.1. Intercomparison of Gridded Climate Data with Station Observations To examine the consistency and pattern of gridded datasets against the ECCC, each dataset is spatially averaged over the LNRB from 1981 to 2010 (Figure 2.4). Overall, annual precipitation from the ERA-I and WFDEI exceeds that from the ANUSPLIN, ECCC, NARR and HydroGFD datasets across the entire study period. ANUSPLIN consistently underestimates annual total precipitation, whereas HydroGFD and ENSEMBLE show better 29 agreement with the observations for most years. Differences in annual total precipitation from six different datasets increase in recent years, mainly from 2004 to 2010; however, the ENSEMBLE shows better agreement with the ECCC dataset. These emerging differences (post 2003) are likely because of the Canadian precipitation observations not being assimilated into most of the gridded products as of 2004 (Mesinger et al., 2006). Long-term annual total precipitation for the HydroGFD dataset shows less positive PBIAS (0.5%) and RMSE (25.4 mm) relative to other datasets, while ERA-I and WFDEI show high RMSE and PBIAS due to systematic overestimation of precipitation (Table 2.3). The ANUSPLIN data exhibit a dry bias (-8.2%) in annual precipitation amongst other datasets. Long-term seasonal analyses reveal ANUSPLIN underestimates precipitation during all seasons, apart from winter; whereas HydroGFD and ENSEMBLE better represent seasonality with lower RMSE with ranges from 6.3 to 15.9 mm compared to ECCC stations (Table 2.3). The NARR precipitation shows less PBIAS (-0.5% to 10.2%) during winter, summer, autumn and at an annual timescale while it overestimates spring precipitation (17.1%). The ERA-I and WFDEI substantially overestimate seasonal precipitation, which range from 4.3% to 18.6% during spring, summer and autumn. Apart from precipitation differences, the NARR dataset does not show agreement on mean annual air temperature variability with other datasets (Figure 2.4). NARR exhibits ~1.0°C deviation in annual air temperature and high RMSE (1.5°C) over the LNRB relative to the ECCC observations whereas the ANUSPLIN and ERA-I show better agreement with the lowest RMSE (0.2°C) among all other datasets (Table 2.3). ANUSPLIN and HydroGFD are ~0.5°C colder than the observations, a negative bias (-0.5°C to -0.1°C) that persists throughout the study period (Figure 2.4 and Table 2.3). The seasonal analysis reveals colder 30 mean air temperatures from the ANUSPLIN, ERA-I and HydroGFD, which ranges from -0.5°C to -0.1°C, with similar inter-annual variability during spring, summer and autumn. NARR is warmer across all seasons, which affects the ENSEMBLE mean air temperature (i.e., highest (lowest) positive biases, 1.7°C (1.1°C) in winter (spring)). In general, the ANUSPLIN has lower biases and RMSEs than other datasets for mean seasonal air temperature, while ERA-I and HydroGFD fall statistically between ANUSPLIN and WFDEI for spring and summer. Moreover, NARR shows larger RMSEs than its counterparts and has a strong positive bias in mean seasonal and annual air temperature over the LNRB. The NARR biases may translate into the ENSEMBLE dataset that exhibits positive seasonal and annual biases in air temperature. Figure 2.4 Area-averaged annual (a) total precipitation and (b) mean air temperature over the LNRB for the ANUSPLIN, NARR, ERA-I, WFDEI, and HydroGFD datasets compared to four ECCC stations average values across the basin, 1981–2010. 31 Table 2.3 Seasonal and annual precipitation and mean air temperature statistics for the domain-averaged ANUSPLIN, NARR, ERA-I, WFDEI, HydroGFD and ENSEMBLE datasets against four ECCC stations average values across the LNRB, 1981–2010. Precipitation (1981–2010) PBIAS (%) RMSE (mm) Errors Datasets Winter Spring Summer Autumn Annual ANUSPLIN 6.0 14.6 28.1 15.3 47.2 NARR 12.6 21.6 34.1 17.1 43.5 ERA-I 8.6 14.5 36.6 17.8 53.8 WFDEI 21.8 19.1 19.0 24.7 71.8 HydroGFD 9.7 9.8 15.9 14.0 25.4 ENSEMBLE 6.9 6.3 13.3 11.0 21.2 ANUSPLIN 3.9 -12.8 -10.6 -6.4 -8.2 NARR 10.2 17.1 -4.5 -0.5 2.2 ERA-I -4.7 12.4 9.9 7.4 7.9 WFDEI 33.9 18.6 4.3 15.9 13.5 HydroGFD -2.9 -7.8 1.5 6.2 0.5 ENSEMBLE 7.9 4.2 -0.2 3.9 2.6 Mean air temperature (1981–2010) BIAS (oC) RMSE (oC) Errors Datasets Winter Spring Summer Autumn Annual ANUSPLIN 0.4 0.5 0.1 0.3 0.2 NARR 1.9 1.2 1.6 1.6 1.5 ERA-I 0.8 0.4 0.3 0.5 0.2 WFDEI 0.3 0.7 0.5 0.4 0.4 HydroGFD 0.3 0.6 0.1 0.3 0.3 ENSEMBLE 0.7 0.4 0.6 0.6 0.5 ANUSPLIN -0.2 -0.5 -0.1 -0.1 -0.2 NARR 1.7 1.1 1.6 1.5 1.5 ERA-I 0.8 -0.2 -0.2 -0.2 0.0 WFDEI 0.1 0.6 0.4 0.3 0.3 HydroGFD -0.2 -0.5 -0.1 -0.1 -0.2 ENSEMBLE 0.7 0.3 0.5 0.5 0.5 32 2.5.2. Intercomparison of Precipitation and Mean Air Temperature Datasets 2.5.2.1. Basin Average Intercomparison of Datasets The domain-averaged monthly total precipitation magnitudes vary substantially among datasets, with greater inter-dataset differences over the LNRB from March to September (Figure 2.5a). ANUSPLIN consistently underestimates precipitation throughout the study period relative to the IDW and NARR datasets, with 5 to 10 mm month-1 differences, especially in summer. For peak spring and summer precipitation, the range of inter-dataset spread varies from 70 to 90 mm month-1 as overestimated by the ERA-I and WFDEI datasets, respectively, during the study period. On the other hand, the 1981–2010 observed seasonal climate normals for total precipitation compare favourably with the NARR dataset (Environment and Climate Change Canada, 2016). ERA-I and WFDEI, however, overestimate summer precipitation while winter, spring and autumn precipitation correspond more favourably with climate normals. The dry bias in ANUSPLIN precipitation arises possibly from the thin plate smoothing spline surface fitting technique used in its preparation, a feature reported in previous studies (Islam & Déry, 2017; Milewska, Hopkinson, & Niitsoo, 2005; O’Neil, Prowse, Bonsal, & Dibike, 2017). In the reanalysis products, HydroGFD shows the best agreement followed by NARR with IDW, while other products (ERA-I and WFDEI) reveal considerable differences in precipitation, which may have been induced by the climate model used to assimilate and generate these products. However, HydroGFD shows an improvement over the WFDEI and ERA-I datasets when compared to IDW, in agreement with previous studies (Berg et al., 2018). 33 The monthly mean air temperature from all datasets falls below 0.0°C from October to April but rises above 0.0°C in early spring over the LNRB domain (Figure 2.5b). While the inter-datasets seasonal variability of air temperature is quite similar, IDW and NARR monthly mean air temperatures are ~1.5°C warmer compared to the remaining datasets. Seasonal climate normals for monthly mean air temperature of ANUSPLIN, ERA-I, WFDEI and HydroGFD datasets match relatively well, except during spring where NARR overestimates mean air temperature when compared with the observed records (Manitoba Hydro, 2015). Figure 2.5 Area-averaged climatology of (a) monthly total precipitation (solid lines, dotted lines for ENSEMBLE) and (b) monthly mean air temperature over the LNRB for the IDW, ANUSPLIN, NARR, ERA-I, WFDEI, HydroGFD, and ENSEMBLE datasets, 1981–2010. 34 2.5.2.2. Spatial Characteristics (1981–2010) Due to the large number of datasets, and for the sake of clarity, we only compare annual means. A more detailed intercomparison is carried out for the seasonal and annual standard deviations and trends in precipitation and air temperature at the sub-watershed scale in the following sections. To obtain the overall biases among datasets and their spatial distribution, we calculate the spread among the precipitation and air temperature datasets as represented by the ensemble standard deviation. 2.5.2.2.1. Annual Total Precipitation Annual total precipitation patterns across the LNRB exhibit notable differences from relatively wet conditions in WFDEI to relatively dry ones in ANUSPLIN (Figure 2.6); however, the IDW, NARR, and HydroGFD precipitation data are in general agreement with each other. Other noticeable features (wet and dry) in the rainfall distribution are present over the south-eastern part of the basin and Hudson Bay coast with different precipitation magnitudes in all datasets. 35 Figure 2.6 Spatial distribution of the annual total precipitation (mm month-1) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. For the bias calculations, the NARR and HydroGFD datasets deviate only marginally and thus influence the ensemble mean relatively strongly (Figure 2.7). This also applies to the IDW dataset, though the bias is slightly larger in the central part of the basin. Generally, ERA-I and WFDEI exhibit wetter conditions compared to other datasets for most of the regions, whereas ANUSPLIN appears to be unrealistically dry over the LNRB. The NARR and HydroGFD show improvements in modelled precipitation resulting in less differences compared to the ENSEMBLE dataset. The differences are attributed to the number of observation stations used in the IDW and ANUSPLIN and the climatological observations, interpolation and data assimilation procedures used to develop the NARR, ERA-I, WFDEI and HydroGFD. HydroGFD precipitation has low biases in all regions, with dry (wet) biases over the north-eastern part of the basin. The WATCH method with updated climatological observations applied in HydroGFD, appears to have improved precipitation records for the 36 LNRB. NARR is generally slightly wetter in most parts of the basin, which shows an improvement over global reanalysis products (i.e., ERA-I and WFDEI). IDW and ANUSPLIN exhibit dry conditions over the eastern LNRB compared to the other datasets, which may be due to an interpolation technique associated with the poor gauge coverage in the underlying observation stations. The ensemble spread is found to be large (small) over the north-eastern and south-eastern (north-western and middle) sections of the basin (Figure 2.8). This high deviation in precipitation may result from the scarcity of observation stations that were incorporated during development of the climate datasets. The corresponding ensemble mean exceeds 43.0 mm month-1, which is relatively high with low spatial variation in the higher elevations and southern part of the basin. Figure 2.7 Bias as measured against the ENSEMBLE total annual precipitation (mm month -1) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. 37 Figure 2.8 Ensemble mean and standard deviation of (left) precipitation and (right) mean air temperature, 1981– 2010. 2.5.2.2.2. Mean Annual Air Temperature The ANUSPLIN, ERA-I, WFDEI and HydroGFD datasets capture mean air temperature patterns (e.g., the longitudinal gradient between the north-eastern and the south-western parts of the domain) quite similarly, which ranges from 0.4°C to -2.0°C (Figure 2.9). The IDW and NARR datasets exhibit similar spatial patterns, ranging from 2.0°C at higher elevations to -2.0°C near Hudson Bay in mean annual air temperature although both show overall warm biases relative to the ensemble mean (Figure 2.10). Other datasets, including global reanalysis products, show overall cold biases over the domain that may be related to relatively large influences of data assimilation and numerical modelling. Moreover, coarser resolutions of the reanalysis products impact biases due to the proximity of Hudson Bay. Large ensemble spreads (0.9-1.1°C) over the northern part of the basin exist while relatively small deviations (0.5-0.7°C) emerge in the south (Figure 2.8). 38 Figure 2.9 Spatial distribution of the mean annual air temperature (°C) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. Figure 2.10 Bias as measured against the ENSEMBLE mean annual air temperature (°C) for the (a) IDW, (b) ANUSPLIN, (c) NARR, (d) ERA-I, (e) WFDEI, and (f) HydroGFD datasets, 1981–2010. 39 2.5.2.3. Seasonal Means, Totals and Standard Deviations (1981–2010) In general, the precipitation totals and mean air temperature values from all datasets approach each other for all sub-watersheds and seasons (Figures 2.11 and 2.12). The precipitation totals show relatively large differences between datasets and sub-basins (Figure 2.11). Further, summer and spring values show more variability among the datasets than the other seasons. In general, the ERA-I and WFDEI datasets appear to overestimate the seasonal precipitation in all sub-basins, and ANUSPLIN clearly underestimates the precipitation in all seasons and sub-watersheds. Over the southern and north-eastern sub-basins (Gunisao, Limestone and Weir), all datasets show high variability in seasonal precipitation. The temperature means are similar, except for the IDW and NARR, which have a warm bias in all seasons, particularly during summer and winter in all sub-watersheds (Figure 2.12). 40 Figure 2.11 Area averaged seasonal (a) precipitation totals and (b) standard deviations (1981–2010) for all datasets over the ten sub-watersheds. Data are shown for each sub-watershed and all seasons: DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. Although each precipitation dataset is based on its own selected network of stations, with a probable substantial overlap, the IDW and ANUSPLIN include a higher density of stations than the other datasets, notably in the middle portion of the basin (McKenney et al., 2011). The comparison of ANUSPLIN with higher density network datasets over different Canadian watersheds, however, has shown that ANUSPLIN underestimates higher precipitation values (Eum et al., 2014; Wong, Razavi, Bonsal, Wheater, & Asong, 2017). Over the domain, reanalysis products indicate high precipitation values, which lead to a high standard deviation over all sub-watersheds (Figure 2.11). The latter is larger in winter and spring because of the high coefficient of variation in precipitation over the Gunisao sub-watershed (Figure 2.13). 41 The standard deviation over most of the sub-regions is similar for all seasons, except for the Gunisao and Weir, which are the ones with the highest variability in precipitation among all datasets (Figures 2.11 and 2.13). The air temperature standard deviation from all datasets, except for NARR, is higher in winter over all sub-watersheds (Figure 2.12). Other seasonal spreads, except autumn for NARR, in temperature are similar in all the datasets over all subwatersheds. Thus, apart from NARR, every dataset captures similar seasonal air temperature variation. Moreover, the IDW dataset shows high spread over the Weir sub-watershed for all seasons, which might be related to that region’s sparse observation station density. Figure 2.12 Area averaged seasonal (a) air temperature means and (b) standard deviations (1981–2010) for all datasets over the ten sub-watersheds. Data are shown for each sub-watershed and all seasons: DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. 42 Figure 2.13 Seasonal coefficient of variation (1981–2010) for all precipitation datasets over the ten subwatersheds. Data are shown for each sub-watershed and all seasons: DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. 2.5.3. Intercomparison of Annual and Seasonal Trends (1981–2010) 2.5.3.1. Precipitation The 1981–2010 maximum (minimum) annual precipitation trends among all datasets are 30.0 (-30.0) mm decade-1 over the LNRB (Figure 2.14). Trend magnitudes from the IDW and ANUSPLIN show similar decreasing trend magnitudes (5.0-20.0 mm decade-1) whereas reanalysis products such as NARR and ERA-I show significant increasing trends (20.0-30.0 mm decade-1) with disagreement over most of the basin. The WFDEI and HydroGFD show better agreement in the north but opposite trend magnitudes (±20.0 mm decade-1) in the south. Overall, the ENSEMBLE dataset does not show any significant trend over the LNRB; however, its precipitation increases by 5.0-25.0 mm decade-1 over most of the basin. Sub-watershed annual trend signals are statistically insignificant for all regions, excluding the Gunisao, where ERA-I shows significant increasing trend (>30.0 mm decade-1) (Figure 2.15). All datasets show increasing trends in annual precipitation (>10.0 mm decade -1) for all sub-watersheds except ANUSPLIN and IDW, which show similar negative trends. The 43 Gunisao sub-watershed shows different trend magnitudes in each dataset with slight increasing trends from WFDEI and ENSEMBLE. The WFDEI and HydroGFD datasets produce similar trend magnitudes in most of the sub-watersheds, followed by the ENSEMBLE dataset, which indicates the substantial influence these two datasets have on the magnitude and pattern of the ENSEMBLE precipitation. This may also be due to the IDW, ANUSPLIN, ERA-I, and NARR offsetting each other so that the ENSEMBLE values approach the HydroGFD and WFDEI in all sub-basins except the Gunisao. Significant increasing trends (~25.0 mm decade-1) arise in NARR during summer for all sub-watersheds whereas ERA-I manifests similar significant trends (25.0-30.0 mm decade-1) across the Footprint, Gunisao and Grass sub-basins (Figure 2.15). Winter precipitation exhibits no trend across all datasets while autumn precipitation decreases by < 8.0 mm decade-1 for the IDW, ANUSPLIN, HydroGFD and ENSEMBLE datasets in most sub-basins with high drying tendencies (>15.0 mm decade-1) in HydroGFD over the Gunisao sub-watershed. 2.5.3.2. Mean Air Temperature The 1981–2010 maximum and minimum spatial annual warming rate among all datasets, except IDW, are 0.6°C decade−1 and 0.1°C decade−1, respectively over the LNRB (Figure 2.16). Trend magnitudes from the IDW, ANUSPLIN and ERA-I show similar significant warming trends (>0.5°C decade−1) approaching those in NARR, WFDEI and HydroGFD albeit with insignificant increasing trends over the eastern LNRB. The IDW temperature trend shows disagreement (i.e., a cooling trend of 0.2°C decade−1), among all other datasets across the upstream part of the basin. The WFDEI and HydroGFD do not show any significant warming with similar magnitudes across the basin, which may affect the ENSEMBLE air temperature trends with some outliers at the lower elevations. Overall, the 44 ENSEMBLE air temperature does not demonstrate any significant warming across most of the LNRB; however, it shows warming of 0.1-0.6°C decade−1 across the basin. Figure 2.14 Spatial trends of the annual total precipitation (mm decade -1) from different datasets, 1981–2010. The grid cells with significant trends (p < 0.05) are indicated by dots. Almost all datasets show insignificant rising annual air temperatures over all LNRB sub-watersheds, apart from NARR for which we find significant warming trends (>0.4°C decade-1) in the Burntwood and Footprint sub-basins (Figure 2.17). The ANUSPLIN, NARR, ERA-I, WFDEI, HydroGFD, and ENSEMBLE manifest similar annual trend magnitudes for all regions with some exceptions: the NARR trends (>0.3°C decade -1) over all sub-basins are much greater than the other datasets and seem unrealistic. Consequently, trends analyses for 45 the LNRB using NARR should be interpreted with caution. IDW is the only dataset that shows minimum annual temperature warming over the LNRB’s sub-watersheds, which are in fact negative in the Gunisao and Grass. All datasets in autumn show insignificant warming trends except NARR and ERA-I, which show significant and strong warming trends (0.8°C decade-1) over all sub-watersheds (Figure 2.17). NARR exhibits strong warming trends in winter (>0.8°C decade-1); however other datasets (except ANUSPLIN) show more modest winter warming trends over all sub-watersheds. Summer and spring show insignificant negative trends whereas most of the datasets show no trend during winter over each subbasin. Overall, air temperature trend analysis shows insignificant temperature warming during winter and autumn, and less warming in spring and summer over most of the basin. Figure 2.15 Seasonal and annual precipitation trends for all datasets and sub-watersheds, 1981–2010: dots (Annual); DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. Red colour symbols denote significant trends (p < 0.05). 2.6. Discussion and Conclusions An intercomparison of the 1981–2010 ERA-I and WFDEI precipitation showed wetter conditions as compared to four other datasets across the LNRB including ten of its subbasins. The NARR and HydroGFD showed similar interannual variability and magnitudes for 46 most of the sub-watersheds and across the LNRB. The ANUSPLIN showed considerably dry seasonal and annual precipitation for the entire domain and sub-watersheds. The discrepancy in all precipitation datasets was notably greater for southern and low elevation regions of the LNRB and its sub-watersheds (and for all seasons except winter) than those at the middle and higher latitude sub-basins. These results demonstrate that there is greater uncertainty in precipitation over the downstream portion of the basin. This is possibly due to the combined effects of sparse observation stations (Figure 2.2), the coarse spatial resolution of reanalysis products, and the possible undercatch of snowfall and rainfall in the observational data for the downstream and south-eastern part of the LNRB. Moreover, we found much larger differences between ANUSPLIN and reanalysis datasets during 2004 to 2010 than that in previous years. This may be due to the non-assimilation of precipitation observations in reanalysis products over Canada as of 2004 (Bukovsky & Karoly, 2007; Mesinger et al., 2006). Furthermore, NARR precipitation exhibits a structural breakpoint in January 2004 over Canada’s Athabasca watershed, further highlighting this issue (Eum et al., 2014). Nonetheless, previous studies found NARR to provide reliable climate input data relative to the global reanalysis product (NCEP/NCAR Global Reanalysis-1) over northern Manitoba (Choi, Kim, Rasmussen, & Moore, 2009). Precipitation data are more problematic compared to temperature data over Canada due to cumulative daily problems such as evaporation, wetting loss and trace measurements (Mekis & Hogg, 1999). Therefore, it is necessary to conduct a data homogeneity test for different available datasets; moreover, reanalysis products such as NARR should be tested over different Canadian regions before use and limited to the qualified period from 1979 to 2003 for hydro-climatic studies (Choi et al., 2009; Eum et al., 2014; Mesinger et al., 2006). 47 Figure 2.16 Spatial trends of the mean annual air temperature (°C decade -1) from different datasets, 1981–2010. Dots indicate grid cells with significant trends (p < 0.05). The intercomparison showed that NARR and IDW provided warmer seasonal and annual air temperatures than all other datasets over the northern LNRB and its sub-watersheds. The IDW was slightly warmer than the NARR air temperature in the upper basin, and somewhat similar across the rest of the basin. Intercomparison with the four ECCC stations showed that NARR exhibited warmer annual air temperature than other datasets and, therefore, resulted in larger RMSEs with strong positive bias for air temperature over the basin. The ANUSPLIN, ERA-I, WFDEI, and HydroGFD showed better agreement for seasonal and annual air 48 temperatures over the entire basin, although the WFDEI was warm over the upstream part of the LNRB. Moreover, the WFDEI and ERA-I air temperatures showed less bias as compared to other datasets for most of the basin. The intercomparison of different air temperature datasets suggests that the global reanalysis products (ERA-I and WFDEI) compared more favourably against observations than the regional reanalysis (NARR) and stations interpolated (IDW) datasets over the study domain. Consequently, there is a need of an improved air temperature monitoring system to increase the reliability of the observed gridded datasets. However, the ANUSPLIN and IDW datasets have been derived from similar temperature observations within the LNRB but showed different spatial patterns, which exhibit interpolation errors and outlier effects in the resulting datasets. On the other hand, the NARR data showed strong disagreement in air temperatures with other global reanalysis products for the majority of the LNRB. This quantifies discrepancies in data assimilation and numerical modelling techniques during regional and global reanalysis product development. Overall, this analysis showed that HydroGFD may be an optimal option in selecting reanalysis products for both daily precipitation and air temperature for the LNRB. In the observational datasets, IDW and ANUSPLIN, both have major issues with either precipitation or air temperature. Moreover, because of indirect assimilation in NARR through latent heat profiles over Canada during 2004–2010, this study found noticeable differences between selected datasets (Essou et al., 2016; Mesinger et al., 2006). Therefore, combining all gridded datasets and generating a hybrid climate product (i.e., ENSEMBLE) can provide a reliable long-term climate dataset for further hydro-climatic study over the LNRB (Christensen & Lettenmaier, 2007; Eisner et al., 2017; Fowler et al., 2007; Wang et al., 2009). 49 Figure 2.17 Seasonal and annual mean air temperature trends for all datasets and sub-watersheds, 1981–2010: dots (Annual); DJF (Winter) stars; MAM (Spring) diamonds; JJA (Summer) triangles; SON (Autumn) pluses. Red colour symbols denote significant trends (p < 0.05). Trend analyses for the relatively small area could be influenced by different spatial resolutions (e.g., global and regional climate datasets) and the time period used for the study. For instance, precipitation records indicated below average precipitation in many regions since the mid-1980s; therefore, the trend analysis showed small increasing trends in eastern Canada for 1948–1992 (Boden, Kaiser, Sepanski, & Stoss, 1994; Environment Canada, 1995). Quantitative analyses of the influence of different spatial resolutions, seasonal and annual variations on the climate trends are difficult. Therefore, we investigated climate variability and trends from available datasets (global and regional) at both seasonal and annual timescales. Given many previous findings on Canadian climate trends are supported by our analysis where they indicate insignificant increasing trends, which are not artificial (e.g., changes in observation sites, locations or programs), in precipitation and air temperature over northern Manitoba (LNRB) (e.g., Environment Canada, 1995; Skinner & Gullett, 1993; Zhang et al., 2012). Using multiple datasets for trend analysis from different sources confirms data quality and the outcome to be reliable and robust. Perhaps the most interesting finding based on our analysis was the disagreement in historical precipitation 50 trends from different gridded datasets, which involved wet and warm patterns across the basin. The eastern sub-watersheds are affected by abnormally warm and dry conditions in summer and autumn despite the general increase in precipitation throughout the study period. The cause of different spatial and temporal trends from available gridded datasets may be associated with increasing atmospheric greenhouse gases or natural climate variability and cannot be addressed by a study of this nature. In this study, our primary goal was to obtain the data that have better comparable historical climatic trends with observations. There is evidence, nevertheless, suggesting that a certain degree of agreement exists between observed trends in Canadian climate and those estimated by the WFDEI, HydroGFD and ENSEMBLE precipitation, and air temperature from all datasets except NARR and ERA-I (Environment Canada, 1995; Skinner & Gullett, 1993; Zhang et al., 2012). One of the greatest challenges presently facing the hydro-climatic and hydrological modelling community remains the reliability of available gridded data at different locations, including the distribution and amounts of regional and continental scale precipitation. Most reanalysis products predict a concomitant increase in regional and global precipitation as mean air temperatures rise during the historical time period (1981–2010). This study suggests that available gridded datasets have issues, mainly at a seasonal time scale, with either precipitation or air temperature. The NARR precipitation, WFDEI air temperature, and HydroGFD precipitation and air temperature more closely resemble observations. Moreover, it is useful to have all datasets to quantify the uncertainty in precipitation and air temperature across the basin and its sub-watersheds. Confidence in the historical and future hydrologic cycle simulations could be improved if ENSEMBLE precipitation and air temperature are used to drive land surface models. 51 CHAPTER 3: SENSITIVITY ANALYSIS AND UNCERTAINTY ASSESSMENT IN WATER BUDGETS SIMULATED BY THE VARIABLE INFILTRATION CAPACITY MODEL FOR CANADIAN SUB-ARCTIC WATERSHEDS Publication details: This chapter has been published in Hydrological Processes. Lilhare, R., Pokorny, S., Déry, S. J., Stadnyk, T. A., & Koenig, K. A. (2020). Sensitivity analysis and uncertainty assessment in water budgets simulated by the Variable Infiltration Capacity model for Canadian sub-arctic watersheds. Hydrological Processes, 34(9), 2057– 2075. https://doi.org/10.1002/hyp.13711 3.1. Abstract In this study, we evaluate uncertainties propagated through different climate datasets in seasonal and annual hydrological simulations over ten sub-arctic watersheds of northern Manitoba, Canada, using the Variable Infiltration Capacity (VIC) model. Further, we perform a comprehensive sensitivity and uncertainty analysis of the VIC model using a robust and state-of-the-art approach. The VIC model simulations utilize the recently developed Variogram Analysis of Response Surfaces (VARS) technique that requires in this application more than 6000 model simulations for a 30-year (1981–2010) study period. The method seeks parameter sensitivity, identifies influential parameters, and showcases streamflow sensitivity to parameter uncertainty at seasonal and annual time scales. Results suggest that the Ensemble-VIC simulations match observed streamflow closest whereas global reanalysis products yield high flows (0.5-3.0 mm day-1) against observations and an overestimation (1060%) in seasonal and annual water balance terms. VIC parameters exhibit seasonal importance in VARS, and the choice of input data and performance metrics substantially affect sensitivity analysis. Uncertainty propagation due to input forcing selection in each 52 water balance term (i.e., total runoff, soil moisture, and evapotranspiration) is examined separately to show both time and space dimensionality in available forcing data at seasonal and annual time scales. Reliable input forcing, the most influential model parameters, and the uncertainty envelope in streamflow prediction, are presented for the VIC model. These results, along with some specific recommendations, are expected to assist the broader VIC modelling community, and other users of VARS, and land surface schemes, to enhance their modelling applications. 53 3.2. Introduction Numerical modelling of a river basin remains essential for both climate and ecological studies as it provides vital information about the hydrological cycle and water availability for societies and ecosystems. Although recent developments and advances have been achieved in hydrological modelling and computational power, addressing efficiently the uncertainties in hydrological simulation remains a critical challenge (Liu & Gupta, 2007). There is a growing need for sensitivity and uncertainty assessments associated mainly with the model and input forcing datasets to achieve the hydrological model’s optimal performance for decision making. Input climate forcing for numerical modelling, primarily precipitation and air temperature, are essential for accurate streamflow simulations and water balance calculations (Eum et al., 2014; Fekete, Vörösmarty, Roads, & Willmott, 2004; Reed et al., 2004; Tobin, Nicotina, Parlange, Berne, & Rinaldo, 2011). For cold regions, these input forcing alter the phase and magnitude of modelled variables and cascade through all hydrological processes during numerical simulations, impacting the reliability of model output (Anderson et al., 2008; Tapiador et al., 2012; Wagener & Gupta, 2005). In Canada, however, numerous studies have also used multiple forcing datasets to assess the performance of hydrological simulations. For example, Sabarly et al. (2016) used four reanalysis products to evaluate the terrestrial branch of the water cycle over Québec, Canada with acceptable results for the period 1979–2008. The question of which forcing dataset is the most suitable and accurate to drive hydrological models remains elusive and inconclusive. Steps toward answering that question were undertaken by Pavelsky & Smith (2006) who concluded that observations covered the trends significantly better than two reanalysis products when they assessed the quality of four global precipitation datasets against the discharge observations from 198 pan- 54 Arctic rivers. The bias and uncertainty in global hydrological modelling due to input datasets and associated over- or underestimations in modelled streamflows in several river basins have also been identified in previous studies (e.g., Döll, Kaspar, & Lehner, 2003; Gerten, Schaphoff, Haberlandt, Lucht, & Sitch, 2004; Nijssen, Schnur, & Lettenmaier, 2001). While there may be other uncertainties (e.g., model structure, calibration, soil type, land use, etc.), this paper focuses primarily on the uncertainties due to model parameters and input forcing datasets, which are perhaps the most significant source of uncertainty for any hydrological modelling study (Zhang et al., 2016). In practice, many (from tens to hundreds) parameters in most hydrological models lead to dimensionality issues where parameter estimation becomes mostly nonlinear and a high dimensional problem. Numerous optimization algorithms are available to address these problems (e.g., Abebe, Ogden, & Pradhan, 2010; Aster, Borchers, & Thurber, 2013; Beven & Binley, 1992; Duan, Sorooshian, & Gupta, 1992; Hill & Tiedeman, 2007; Vrugt, Diks, Gupta, Bouten, & Verstraten, 2005; Vrugt, Gupta, Bouten, & Sorooshian, 2003) but it is not often feasible or necessary to include all these parameters in the calibration and Sensitivity Analysis (SA) process to obtain efficient optimization and sensitive parameters, respectively. For instance, over-parameterization is another well-known problem in land surface modelling (Van Griensven et al., 2006). At present, various SA methods (e.g., qualitative or quantitative, local or global, and screening or refined methods) are used widely in different fields, such as complex engineering systems, physics, and social sciences (Frey & Patil, 2002; Iman & Helton, 1988). Given the extensive range of SA methods available, users should have a clear understanding of the methods that are appropriate for a specific application. In general, the Variable Infiltration Capacity (VIC) hydrological model 55 incorporates many parameters (some with physical significance and some statistical), which are used to calibrate the model by various methods. In some cases, parameters with physical significance may be adjusted interactively during calibration. Some parameters may have less influence on model output such that they could be easily ignored. One of the objectives of this study is thus to explore the sensitivity of VIC calibration parameters to reduce the dimensionality issue in model optimization at different time scales and to establish their interannual importance in the calibration and model performance. In this study, we quantify the uncertainty propagated from available forcing datasets in their surface water balance estimations over the Lower Nelson River Basin (LNRB) in northern Manitoba, Canada. To achieve this goal, seven input forcing datasets that are intercompared in our companion paper are ingested into the VIC model over the LNRB (Lilhare et al., 2019). These datasets are used in various other studies over different Canadian regions (Boucher & Best, 2010; Islam & Déry, 2017; Sauchyn, Vanstone, & Perez-Valdivia, 2011; Seager et al., 2014; Woo & Thorne, 2006). To our knowledge, this is perhaps the first comprehensive study that collectively utilizes available gridded datasets in hydrological modelling, establishes the most suitable datasets, minimizes the input data uncertainty by evaluating the best performing product, and then propagates input and parameter uncertainty through the model output. Moreover, we consider not only the total uncertainty (i.e., total runoff) but also the apportioned uncertainty in runoff generating processes such as precipitation, evapotranspiration (ET hereafter), and soil moisture at annual and seasonal time scales. The main objectives of this study are to: (i) examine uncertainty propagated through various input forcing datasets in the VIC model; (ii) identify parameter sensitivity of 56 the VIC model to streamflow; and (iii) assess streamflow sensitivity to parameter uncertainty in the VIC model over the LNRB. 3.3. Study Area In this study, the lower Nelson River, which is the downstream segment of the Nelson River system, is selected for the VIC modelling, sensitivity, and uncertainty analyses (Figure 3.1). The LNRB spans an area of ~90,500 km2 and collects all water from the drainage area upstream of the Nelson River (~970,000 km2) before discharging into Hudson Bay. In the LNRB, the main stem river (Nelson) and its largest tributary – the Burntwood, whose downstream segment carries diverted flows from the Churchill River – have less seasonal flow variability due to streamflow regulation and a large drainage area. Most of the LNRB has gentle slopes, with common channelized lakes moderating flow variability. Wetlands abound within the LNRB and store significant volumes of water, cover large areas and moderate streamflow responses to rainfall and snowmelt events. Shallow soils and permafrost limit infiltration, groundwater storage and groundwater flows. To increase its hydroelectric capacity, Manitoba Hydro manages flows in the LNRB with two major sources of streamflow regulation: the Churchill River Diversion and Lake Winnipeg Regulation. The LNRB experiences a sub-arctic continental climate characterized by moderate precipitation and humidity, cool summers, and cold winters. The snow-free season remains brief, generally beginning in May and ending in October. Most of the precipitation that occurs during the summer months falls as rain, accounting nearly two-thirds of the total annual precipitation. The most expansive land cover class in the LNRB is temperate or subpolar needleleaf forest covering ~33% of its total area with secondary classes being mixed forests (19%) and temperate or sub-polar shrublands (9%) (North American Land Change 57 Monitoring System, 2010). Wetlands (bogs and fens, 21%) and open surface water (13%) also prevail in the region. The entire region exhibits low relief with a maximum elevation of 390 m.a.s.l. and average basin slope of 0.037%. Permafrost abounds in the LNRB with the downstream, northeastern portion underlain by continuous (between 90% to 100%) and extensive discontinuous (between 50% to 90%) permafrost (approximately 0.8% and 9% of the LNRB, respectively) while sporadic discontinuous (between 10% to 50%) and isolated permafrost spans ~68% and 16% of the LNRB’s total area, respectively (Natural Resources Canada, 2010). Figure 3.1 Maps of the LNRB. (a) The Nelson River Basin (NRB), Churchill River Basin (CRB), and Lower Nelson River Basin (LNRB). (b) major rivers and sub-watersheds within the LNRB, yellow triangles show the hydrometric stations used in this study, white circles denote existing generating stations, and the yellow circle shows a future generating station (currently under construction) by Manitoba Hydro. A red star indicates the Churchill River diversion and the Digital Elevation Model (DEM) represents the VIC model domain at 0.10° resolution. 58 3.4. Materials and Methods 3.4.1. Datasets Soil parameters for the VIC model are sourced from the multi-institution North American Land Data Assimilation System (NLDAS) project at 0.50° resolution (Cosby, Hornberger, Clapp, & Ginn, 1984). These soil parameters are then aggregated to the VIC model resolution (0.10°) following Mao & Cherkauer (2009). Frost-related parameters (e.g., bubbling pressure) are extracted from Miller & White (1998), or set to default values (Mao & Cherkauer, 2009). Land cover data are obtained from the Natural Resources Canada’s (NRCan) GeoGratis - Land Cover, circa 2000-Vector (LCC2000-V) product and vegetation parameters are estimated by following Sheffield & Wood (2007). All landcover classes are mapped into standard VIC model vegetation classes and the Leaf Area Index (LAI) for each vegetation class in every grid cell is estimated (Myneni, Ramakrishna, Nemani, & Running, 1997). Rooting depths are obtained from Maurer et al. (2002), while other vegetation parameters are taken from Nijssen et al. (2001). We obtain various gridded forcing datasets for further analysis: the Australian National University spline interpolation (ANUSPLIN), North American Regional Reanalysis (NARR), European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERAInterim), European Union Water and Global Change (WATCH) Forcing Data ERA-Interim (WFDEI), and Hydrological Global Forcing Data (HydroGFD). As well, an Inverse Distance Weighted (IDW) dataset constructed from 14 Environment and Climate Change Canada (ECCC) meteorological stations across the LNRB using a squared IDW interpolation technique is also used (see Table 3.1 for more details). These datasets are assembled to produce the Ensemble dataset from 1981 to 2010. Our companion paper (Lilhare et al., 2019) 59 and Appendix A provide a comprehensive intercomparison and additional details of these datasets. The NARR, ERA-I and HydroGFD daily precipitation and wind speed are obtained from the sum and average of 3-hourly values for a 24-hour period, respectively. To obtain daily maximum and minimum air temperature (Tmax and Tmin) for these products, we extract the maximum and minimum value for one day from the 3-hourly NARR, ERA-I and HydroGFD air temperature products. Daily wind speed is not available for the ANUSPLIN and IDW forcing datasets. The observed wind speeds, both upper air and near-surface values, are assimilated in the NARR reanalysis product and they show satisfactory correspondence with ECCC observations (Hundecha, St-Hilaire, Ouarda, Adlouni, & Gachon, 2008). Therefore, we use NARR wind speeds to run VIC in combination with the ANUSPLIN and IDW datasets for input forcing. For the Ensemble, daily precipitation, Tmax and Tmin are derived from the equally weighted average of all six gridded products, while the daily wind speed ensemble is calculated from four reanalysis products (NARR, ERA-I, WFDEI and HydroGFD) as the other two datasets (IDW and ANUSPLIN) do not have such records. The equally weighted ensemble approach has been used previously over global and regional domains to evaluate changes in water balance components under historical and projected future climate conditions (Fowler et al., 2007; Fowler & Kilsby, 2007; Mishra & Lilhare, 2016; Wang et al., 2009). 60 Table 3.1 VIC intercomparison experiments performed using different forcings (Lilhare et al., 2019). VIC model input forcing datasets Description IDW Inverse Distance Weighted interpolated observations from 14 ECCC meteorological stations (Gemmer et al., 2004; Shepard, 1968) ANUSPLIN The Australian National University spline interpolation (Hopkinson et al., 2011; Natural Resources Canada, 2014) NARR North American Regional Reanalysis (Mesinger et al., 2006) ERA-I European Reanalysis-Interim (Dee et al., 2011) WFDEI European Union Water and Global Change (WATCH) Forcing Data ERAInterim (Weedon et al., 2014) HydroGFD Hydrological Global Forcing Data (Berg et al., 2018) Ensemble Average of above mentioned six gridded datasets VIC configuration Domain: 53°−58°N, 91°−103°W Resolution: 0.10°×0.10° Time step: Daily Soil Layers: 3 Vertical elevation band: No Natural lakes and frozen ground: On Calibration period: 1981–1985 (dry/cool) and 1995–1999 (wet/warm) Validation period: 1986–1994 (average) Overall simulation period: 1981–2010 3.4.2. The Variable Infiltration Capacity (VIC) Model In this study, the VIC (version 4.2.d) model (Liang et al., 1994, 1996) with more recent modifications is used to simulate daily streamflow in full water and energy balance mode (Bowling & Lettenmaier, 2010; Bowling et al., 2003; Cherkauer & Lettenmaier, 1999) (Table 3.1). The VIC model grid cells over the LNRB comprise 41 rows and 90 columns with a 5° range of latitudes (53°-58°N) and a 12° range of longitudes (103°-91°W). The VIC model uses three soil layers, five soil thermal nodes (the default value) and a constant bottom boundary temperature at a damping depth of 10 m for our study region (Williams & Gold, 1976). The LNRB’s tiles are characterized by soil and vegetation fractions, which are partitioned proportionally within a grid cell. For cold regions hydrology, VIC follows the 61 U.S. Army Corps of Engineers’ empirical snow albedo decay curve (USACE, 1956), the total precipitation is distributed based on the 0.10° grid cells, and the air temperature is adjusted to resolve the precipitation type with a 0°C threshold to discriminate rainfall/snowfall. The default single elevation band is used whereby VIC assumes that each grid cell is flat and takes the mean grid elevation into account for simulations over the LNRB. A finite difference algorithm for frozen soil, which tracks soil ice content and represents permafrost, is implemented into the VIC model to improve its modelling abilities (Cherkauer & Lettenmaier, 1999, 2003). The frozen soil algorithm solves heat fluxes through the soil column using a heat transfer equation (Cherkauer & Lettenmaier, 1999). This algorithm supersedes the original soil thermal flux equations (Liang, Wood, & Lettenmaier, 1999) in favour of a more robust numerical technique (Cherkauer & Lettenmaier, 1999) that simulates soil temperatures at five thermal nodes through the soil column. Natural lakes and wetlands are considered in the model implementation; however, anthropogenic structures (i.e., dams, reservoirs) and flow regulation are not incorporated in the VIC model. The VIC model lake and wetlands algorithm represents the effects of hydrologically disconnected lakes and wetlands by creating its land class that can be added to the grid cell mosaic, in addition to the vegetation and bare soil land classes (Bowling & Lettenmaier, 2010). It does not represent riparian systems that receive water from overbank flow. The energy balance of open water in VIC builds on the work of Hostetler (1991), Hostetler & Bartlein (1990), and Patterson & Hamblin (1988), while that of the exposed wetland vegetation follows Cherkauer & Lettenmaier (1999). Ten of the lower Nelson River’s unregulated tributaries (including the unregulated, upstream portion of the Burntwood River) are selected for model calibration, evaluation, and subsequent analyses (Table 3.2). The routing network and other essential inputs for the routing model (e.g., flow direction, fraction, and mask) are created at 10 km 62 resolution for the entire LNRB using the 30 m Shuttle Radar Topography Mission (SRTM) digital elevation model (United States Geological Survey, 2013). Table 3.2 List of ten selected unregulated hydrometric stations, maintained by the Water Survey of Canada and Manitoba Hydro, for the VIC model calibration and evaluation with sub-watershed characteristics and mean annual discharge (Water Survey of Canada, 2016). Station name (Gauge ID) Latitude Longitude (°N) (°W) Mean subwatershed elevation (m) Gauged drainage area (km2) Mean annual discharge (m3 s-1) Burntwood River above Leaf Rapids (05TE002) 55.49 99.22 302.4 5,810 22.9 Footprint River above Footprint Lake (05TF002) 55.93 98.88 273.8 643 3.2 Grass River above Standing Stone Falls (05TD001) 55.74 97.01 265.0 15,400 64.6 Gunisao River at Jam Rapids (05UA003) 53.82 97.77 260.9 4,610 18.0 Kettle River near Gillam (05UF004) 56.34 94.69 164.7 1,090 13.2 Limestone River near Bird (05UG001) 56.51 94.21 173.6 3,270 21.5 Odei River near Thompson (05TG003) 55.99 97.35 253.5 6,110 34.3 Sapochi River near Nelson House (05TG006) 55.90 98.49 259.1 391 2.2 Taylor River near Thompson (05TG002) 55.48 98.19 236.2 886 5.1 Weir River above the mouth (05UH002) 57.02 93.45 125.8 2,190 15.6 3.4.3. Calibration and Evaluation For VIC model calibration, an optimization process using the University of Arizona Multi-Objective COMplex Evolution Algorithm (MOCOM-UA) minimizes the difference 63 between observed and simulated monthly streamflow at unregulated hydrometric gauge locations within the LNRB (Shi, Wood, & Lettenmaier, 2008; Yapo, Gupta, & Sorooshian, 1998). Here, the training parameter set used in the sensitivity and calibration processes comprises six soil parameters: binf (infiltration parameter that controls the amount of water infiltrating into the soil with values ranging from 0 to 0.4, in fractions), Dsmax (the maximum velocity of baseflow for each grid cell ranging from 0 to 30 mm day-1), Ws (the fraction of maximum soil moisture where nonlinear baseflow occurs ranging from 0 to 1), D2 and D3 (thickness of the second and third soil layers, which affects the soil moisture storage capacity, ranging from 0.3 to 1.5 m), and Ds (fraction of the Dsmax parameter at which nonlinear baseflow occurs ranging from 0 to 1). The Nash–Sutcliffe efficiency (NSE) (Nash & Sutcliffe, 1970), Kling–Gupta efficiency (KGE) (Gupta, Kling, Yilmaz, & Martinez, 2009), and Pearson’s correlation (r) coefficients (for simulated vs observed monthly streamflows) in addition to Percent Bias (PBIAS) provide metrics to summarize model performance. Separate calibration using each forcing dataset is applied to all ten sub-basins within the LNRB to determine the most optimized parameters against the observed streamflow. We use a split-sample approach to span the variety of relatively dry/wet/warm/cool years. Years 1981–1985 (dry/cool) and 1995–1999 (wet/warm) are used for calibration, and 1986–1994 (average) forms the validation period (Table 3.1) (Lilhare et al., 2019). The MOCOM-UA optimizer searches a group of VIC input parameters using the population method; it attempts to maximize the NSE coefficient between observed and simulated streamflow at each iteration. At each trial, the multi-objective vector consisting of VIC parameters is determined, and the population is ordered by the Pareto rank of Goldberg (1989). In the MOCOM-UA optimization process, the user defines the training parameter set and these parameters are selected based on the calibration experience from previous studies 64 (Islam & Déry, 2017; Kang, Gao, Shi, Islam, & Déry, 2016; Nijssen, Lettenmaier, Liang, Wetzel, & Wood, 1997; Shi et al., 2008). 3.4.4. Experimental Set-up and Analysis Approach A series of different VIC model setups is conceived to: (i) compare the VIC model’s response when forced by different gridded datasets (with each simulation referred to as a given “dataset-VIC” hereafter), (ii) evaluate the uncertainties propagated in the water budget estimation using different input forcings, (iii) assess VIC parameter sensitivity using the Variogram Analysis of Response Surfaces (VARS) at seasonal and annual time scales, and (iv) gauge streamflow sensitivity to the VIC model parameter uncertainty (Figure 3.2). The sensitivity, parameter sampling, and uncertainty methodology are discussed in the following sub-sections. Moreover, the VIC simulations driven by each forcing dataset from 1981–1985 are used to generate a VIC initial state parameter file, to allow model spin-up time for five years. This diminishes simulation uncertainty in the calibration, validation, and water balance estimation during the study period. Intercomparison of the seven meteorological datasets from our companion paper suggests that the Ensemble dataset provides more robust historical meteorological forcing (Lilhare et al., 2019); therefore, the VIC model forced by the Ensemble dataset (i.e., Ensemble-VIC) is used as a reference calibration simulation to investigate the propagated uncertainties in water balance estimation from different input forcing datasets. 65 Figure 3.2 Schematic representation of the overall methodology adopted for the propagation, sensitivity, and uncertainty assessment in the VIC modelling over the LNRB. Coloured boxes indicate various objectives of this study. 3.4.4.1. Sensitivity Analysis VARS, a model parameters SA approach, is applied to the VIC model (Razavi & Gupta, 2016a, 2016b). The SA approach reduces the number of parameters that numerical models require to consider in the optimization process. Moreover, the setup is useful for highdimensional optimization problems and can reduce the parametrization uncertainty (Razavi, Sheikholeslami, Gupta, & Haghnegahdar, 2019). We utilize six VIC model parameters in the “star-based” sampling strategy (STAR-VARS), to incorporate VARS in the VIC model and subsequent uncertainty assessment (Razavi & Gupta, 2016b). Parameter selection is based on the optimal number of VARS simulations and SA conducted by various VIC model users using different SA methods (Demaria, Nijssen, & Wagener, 2007; Kavetski, Kuczera, & Franks, 2006; Liang & Guo, 2003; Liang, Xie, & Huang, 2003). Specifically, we evaluate the sensitivity of the Kling-Gupta criterion (Gupta et al., 2009) (which measures goodness-of-fit between simulated and observed streamflows) to variations in the six VIC model parameters across their feasible ranges. VARS determines parameter reliability through a bootstrapping process and ranks them based on similar parameter occurrence and relative sensitivity 66 (Razavi & Gupta, 2016a). The SA is performed at seasonal and annual time scales and if a given model parameter suffers with identifiability issues then it varies temporally in relative sensitivity and reliability. We use 35 star centres (i.e., 1925 VIC model runs for each subwatershed) and 0.10 variogram resolution to generate efficient and robust estimates of the VARS sensitivity ranking (Razavi & Gupta, 2016b). 3.4.4.2. VIC Parameter Sampling and Uncertainty Analysis Parametric uncertainty is assessed by utilizing the Ensemble input forcing and Generalized Likelihood Uncertainty Estimation (GLUE) methodology using 1925 STARVARS samples (Beven & Binley, 1992). Model parameters are sampled from uniform prior distributions and behavioural parameter sets and then used to generate parameter likelihood distributions. The pseudo-likelihood function of KGE is used to assess model performance. The less subjective selection criteria are a common practice in the literature thus we use a behavioural parameter set, which subjectively meets the desired performance criteria (Li & Xu, 2014; Shafii, Tolson, & Matott, 2015; Stedinger, Vogel, Lee, & Batchelder, 2008). These methods fail to account for output uncertainty; therefore, we use a simple method of selecting the top 10% of the model simulations. The STAR-VARS generates directional variograms in the dimension of each parameter. This implies that once a parameter’s directional variogram is sampled for a star centre, it is held constant until being varied for the next star centre; this creates a high density of sampling at one parameter value per star centre. To determine sufficient sampling towards reasonably well-defined parameter uncertainty, we perform a visual inspection of the parameter distribution. If the most sensitive parameters, determined by VARS, show notable deviation from their uninformed priors by visual inspection, then we assume sufficient sampling of the parameter. Further, if the GLUE method, using 1925 67 STAR-VARS samples, fails to accommodate reasonable likelihood distributions, then we additionally perform Orthogonal Latin Hypercube (OLH) sampling. The OLH offers uniform sampling and generates an additional 600 VIC parameter samples (Gan et al., 2014). 3.5. Results 3.5.1. Intercomparison of the VIC Simulations The NSE (KGE) average scores during calibration and validation are much higher for the NARR-VIC, Ensemble-VIC, ANUSPLIN-VIC, and HydroGFD-VIC (NARR-VIC, Ensemble-VIC, HydroGFD-VIC, and WFDEI-VIC) compared to other simulations (Figure 3.3). Despite low (< 0.5) NSE and KGE scores from the IDW-VIC, ANUSPLIN-VIC, ERAI-VIC, and WFDEI-VIC simulations, the correlation coefficients remain substantially high for all sub-basins. The negative (positive) biases from IDW-VIC, ANUSPLIN-VIC and HydroGFD-VIC (ERA-I-VIC and WFDEI-VIC) contribute to the lower NSE and KGE coefficients, whereas the timing of seasonal flows resembles the observed flows in the IDWVIC and ANUSPLIN-VIC. The ERA-I-VIC and WFDEI-VIC simulations reveal strong positive biases for most of the sub-watersheds due to their wet biases in the precipitation datasets (Lilhare et al., 2019); however, they show acceptable NSE and KGE coefficients (≥ 0.5) for most of the sub-basins. 68 Figure 3.3 Boxplots for monthly calibration (a1-d1) and validation (a2-d2) performance metrics, NSE (a1-a2), KGE (b1-b2), r (simulated vs observed monthly streamflow) (p-value < 0.05 for all) (c1-c2) and PBIAS (d1d2), for ten sub-watersheds within the LNRB based on IDW-VIC, ANUSPLIN-VIC, NARR-VIC, ERA-I-VIC, WFDEI-VIC, HydroGFD-VIC and Ensemble-VIC simulations. The black dots within each box show the mean, the red lines show the median, the vertical black dotted lines show a range of minimum and maximum values excluding outliers, and the red + signs show the outliers defined as the values greater than 1.5 times the interquartile range of each metrics. 69 Figure 3.4 Simulated and observed daily runoff (mm day-1) averaged over water years 1981–2010 for the LNRB’s ten unregulated sub-basins. An external routing model is used to calculate runoff for the IDW-VIC, ANUSPLIN-VIC, NARR-VIC, ERA-I-VIC, WFDEI-VIC, HydroGFD-VIC and Ensemble-VIC simulations. Note that y-axis scales vary between panels. Comparison of simulated daily runoff against the observed hydrometric records reveals satisfactory model performances from the NARR-VIC and Ensemble-VIC, while the IDW-VIC and ANUSPLIN-VIC runoff is considerably low for all sub-watersheds (Figure 3.4). ANUSPLIN-VIC and IDW-VIC runoff shows substantial disagreement with the 70 observed hydrograph, especially in the Kettle, Limestone, Odei, Sapochi, and Weir subbasins, owing to the dry bias and undercatch issue in the precipitation data, respectively. The ERA-I-VIC and WFDEI-VIC simulations overestimate summer and autumn runoff and capture reasonably well winter and spring flows for all sub-watersheds. Simulated flows for the Burntwood, Footprint and Taylor sub-watersheds from all VIC simulations are comparable in magnitude with observations, but the timing is considerably shifted (~20 days), yielding more spring runoff and earlier decline of summer recession flows. The NARR air temperature is warmer among all other datasets during winter, spring and autumn. This improves the snowmelt-driven runoff in the NARR-VIC simulation, causing a better representation of simulated flows for these seasons over each sub-watershed. The EnsembleVIC and NARR-VIC simulations exhibit satisfactory hydrographs with ≥ 0.5 NSE and KGE scores in most of the sub-basins (Figures 3.3 and 3.4). 3.5.2. Uncertainty in the Water Budget Estimation The average annual precipitation and VIC simulated water budgets, which are important factors driving changes in total runoff, from all input forcing and their standard deviations (SD) are estimated to find the uncertainty in annual water budgets over the LNRB’s sub-basins (Tables 3.3 and 3.4). For 1981–2010, the Gunisao sub-watershed shows high average annual inter-dataset variability (52.7 mm year-1) in precipitation that results in 61.5, 50.0, and 88.8 mm year-1 SD in the total runoff, ET, and average soil moisture, respectively. The Gunisao (southern outlier in the study area, where air temperatures are much warmer than other sub-basins) generates the lowest total runoff despite the highest annual precipitation. This situation arises through the compensating effect of the highest ET values in the LNRB and appreciable precipitation variability that contributes to overall 71 uncertainty for this sub-basin. A somewhat similar pattern arises in the Grass where a decrease in precipitation uncertainty yields less deviation in total runoff; for example, the Grass sub-watershed exhibits 26.7 mm year -1 deviation in precipitation, which results in 19.1 mm year-1 deviation in total runoff, smallest among all sub-basins. The smaller Sapochi, Footprint and Taylor (gauged area < 900 km2) sub-basins manifest similar inter-dataset errors (> 29 mm year-1) for annual precipitation. Further, relatively larger sub-watersheds (i.e., Gunisao and Odei) show significant differences in the SD, which reveal higher spatial variability from different forcing datasets. Table 3.3 Components of the simulated water budgets in the LNRB’s sub-watersheds, with average annual values for 1981–2010. The average annual precipitation (PCP) based on the mean of six forcing datasets, and other terms are the total runoff (TR), ET, and average soil moisture (SM), all based on the mean of VIC simulations from six different input forcing datasets. Standard deviation (SD) shows inter VIC simulations variation in the water balance estimations. Sub-basins PCP (mm) TR (mm) ET (mm) SM (mm) Mean SD Mean SD Mean SD Mean SD Burntwood 502.8 29.6 97.1 23.8 408.7 17.5 77.4 14.0 Footprint 521.0 35.5 109.6 31.3 408.3 34.6 167.1 78.3 Grass 508.2 26.7 92.8 19.1 412.0 21.4 168.8 55.6 Gunisao 546.6 52.7 93.3 61.5 451.8 50.0 150.5 88.8 Kettle 519.6 47.2 148.6 46.9 373.4 23.3 83.5 15.2 Limestone 511.6 44.9 132.2 48.0 380.1 23.4 93.1 24.2 Odei 525.3 35.3 148.2 44.6 379.8 29.9 89.7 14.7 Sapochi 524.5 34.3 109.1 28.0 417.7 36.7 94.4 18.8 Taylor 522.2 29.4 137.6 32.6 385.9 27.0 91.5 13.4 Weir 508.3 44.7 129.6 45.8 380.6 21.9 91.3 25.2 Mean 519.0 38.0 119.8 38.1 399.8 28.6 110.7 34.8 Area-averaged seasonal total runoff (TR) shows higher uncertainty for relatively large sub-watersheds (e.g. Gunisao, Kettle, Limestone, Odei, and Weir), especially in spring and summer (Figures 3.5 and 3.6 b1-b4). The Ensemble-VIC TR (black dots) matches closely the 72 average of the remaining VIC simulations (red bars, referred to as multidata-VIC hereafter) for all sub-basins (Figure 3.6). The Ensemble-VIC captures multidata-VIC spring TR closely in eight out of ten sub-watersheds whereas two others show underestimation (Figure 3.6 b2). This underestimation persists in summer, which could arise from the extension of calibrated parameters to the entire study period (Figure 3.6 b2-b3). This approach may be unable to represent long term daily and seasonal runoff for these sub-basins. Inter-seasonal air temperature analysis shows that due to extreme minimum air temperature in winter, simulated multidata and Ensemble-VIC TR over each sub-watershed are low and result in less uncertainty between simulations. The simulated error increases in early spring and persists until late autumn, consistent with seasonal precipitation for all sub-watersheds. For annual TR estimates, the Gunisao, Kettle, Limestone, and Weir sub-watersheds reveal high inter-simulation error whereas relatively smaller sub-basins show less deviation and better TR estimates from Ensemble-VIC (Appendix B). Table 3.4 Average annual precipitation (1981–2010) from six forcing datasets and standard deviation (SD) across the datasets over the LNRB’s sub-watersheds. Average annual PCP (mm) IDW ANUSPLIN NARR ERA-I WFDEI HydroGFD SD (mm) Burntwood 494.9 450.4 512.2 529.4 530.6 499.6 29.6 Footprint 522.9 470.6 516.5 539.4 575.7 501.1 35.5 Grass 529.5 479.2 522.1 615.3 605.0 528.4 26.7 Gunisao 504.7 461.7 511.3 528.6 538.7 504.1 52.7 Kettle 482.9 473.1 516.8 554.1 595.8 495.2 47.2 Limestone 480.4 465.1 513.9 537.2 586.7 486.4 44.9 Odei 516.5 478.3 518.8 541.7 584.2 512.5 35.3 Sapochi 524.2 475.1 517.8 536.6 580.3 513.0 34.3 Taylor 526.6 477.8 521.4 534.2 566.4 507.0 29.4 Weir 473.5 457.8 509.1 529.2 583.6 496.6 44.7 Mean 505.6 468.9 516.0 544.6 574.7 504.4 38.0 Sub-basins 73 Figure 3.5 Spatial differences of seasonal total runoff (TR) (mm) for the LNRB’s ten unregulated sub-basins based on Ensemble-VIC (ENSEM) minus (1 st row) IDW-VIC, (2 nd row) ANUSPLIN-VIC, (3 rd row) NARRVIC, (4th row) ERA-I-VIC, (5th row) WFDEI-VIC, and (6 th row) HydroGFD-VIC simulations, water years 1981–2010, for winter (DJF), spring (MAM), summer (JJA) and autumn (SON). The Gunisao and Kettle sub-basins attain 451.8 and 373.4 mm average annual ET, respectively, which are the maximum and minimum values among all other sub-basins (Table 3.3). Regional ET maps from Natural Resources Canada (2016) show 350 to 450 mm average annual ET over the LNRB that satisfy the average ET estimate from VIC for the study period (1981–2010). Due to cold air temperatures in winter, Ensemble-VIC ET is lower (< 3 mm) and corresponds well with the average value (red bars) of other VIC simulations for all sub-watersheds (Figure 3.6 c1–c4). It increases through spring (~100 mm) and peaks in summer (~250 mm) with 35 mm multidata-VIC simulation error, which can be attributed to a substantial rise in air temperature and precipitation. The multidata-VIC SD shows 74 identical values in autumn that essentially reveals less regional variability in ET estimates (~60 mm) from all forcing datasets over the LNRB’s sub-basins. Depleted soil moisture conditions induce basin water limitations that yield uncertainty in ET estimates (Appendix C); for example, the largest sub-watersheds (Grass and Gunisao) within the LNRB show higher uncertainty in ET estimates. As the Ensemble dataset assigns equal weight and each dataset is equally likely to represent the truth, then the Ensemble-VIC simulation indeed better represents the winter, spring, and autumn ET with overestimation in summer for all sub-basins (Figure 3.6 c1–c4 and Appendix C). For annual ET, the Gunisao and Sapochi subbasins show high variability within VIC simulations, but other sub-watersheds have less inter-simulation error (Appendix B). Figure 3.6 Area-averaged VIC simulated seasonal water balance mean (mm) of precipitation (PCP, a1-a4), total runoff (TR, b1-b4), ET (c1-c4) and soil moisture (SM, d1-d4), represented by different columns, for the LNRB’s ten unregulated sub-basins, water years 1981–2010, for the winter (DJF, 1 st row), spring (MAM, 2nd row), summer (JJA, 3 rd row) and autumn (SON, 4 th row) seasons. Red bars show the average of VIC simulations from six forcing datasets (except Ensemble), black error bars indicate inter-VIC simulation variation among IDW-VIC, ANUSPLIN-VIC, NARR-VIC, ERA-I-VIC, WFDEI-VIC, and HydroGFD-VIC, while black dots represent the water balance estimation from Ensemble-VIC. 75 The LNRB’s water balances vary within each sub-basin depending on the magnitude and timing of precipitation and air temperatures obtained from different input forcing datasets (Tables 3.3, 3.4 and Appendix C). For instance, long-term ERA-I-VIC and WFDEIVIC simulations show higher mean TR that results in higher soil moisture (SM) levels for each sub-basin (Figure 3.5 and Appendix D). However, Ensemble-VIC estimates nearly similar seasonal SM conditions as calculated by the average of the six other VIC simulations for most sub-basins (Figure 3.6 d1-d4). Among all other seasons, the highest SM occurs in spring followed by summer and autumn due to seasonal transitions and snowmelt runoff, which are more evident in relatively large sub-watersheds (Burntwood, Gunisao, Grass, Limestone, and Weir) (Figure 3.6 d1-d4). These increased SM values for spring, summer, and autumn show concomitant effects on runoff during these seasons. Furthermore, the Footprint sub-watershed is smaller relative to others; however, it shows considerable interdataset variation (~90 mm) in SM for all seasons. Moreover, eight out of ten sub-basins demonstrate substantial multi-datasets uncertainty in SM for all seasons but mean seasonal SM is well captured by the Ensemble-VIC for these sub-watersheds. However, large uncertainty in SM also suggests excess water availability that goes either into the TR or ET during these seasons. The highest annual SM arises in the Grass, Footprint, and Gunisao subbasins with substantial inter-datasets variation whereas other sub-watersheds show less error in SM simulations with nearly identical annual values (Appendix B and Appendix C). 3.5.3. Model Parameter Sensitivity and Uncertainty Figure 3.7 (a-c) shows results from VARS using values from the Integrated Variogram Across a Range of Scale (IVARS) between 0 and 50% of the parameter ranges (IVARS 50), as suggested in the VARS-Tool manual for a single global sensitivity metric (Razavi & 76 Gupta, 2016a, 2016b). The hydrologically active depth (D2) of soil for movement and storage of water is by far the most important parameter that contributes 24–66% of the sensitivity across all sub-watersheds. The next most important parameter is binf, which accounts for approximately 8–48% of KGE sensitivity across all sub-watersheds. Together, these two parameters contribute to nearly 40–88% of the KGE sensitivity. In the Grass River sub-watershed, Dsmax, which is the maximum velocity of baseflow for each grid cell, also becomes important (~30%) in controlling the amount of runoff generated at the sub-basin outlet. Note that physically inter-linked parameters (D2 and Dsmax) together have almost the same sensitivity ratio in the Grass River. Ds (fraction of the Dsmax parameter) is the third most important parameter, and Ws, the fraction of maximum soil moisture, is also among the more influential parameters in most of the sub-watersheds. Seasonal sensitivity of model parameters changes substantially; for example, in winter Ds and Dsmax, which control baseflow, become the most sensitive parameters (>25%) over all sub-basins whereas in spring and summer D2 still plays a dominant role in computing sub-surface flow (Appendix E). Autumn shows Ds max as the most sensitive parameter since most of the water comes from baseflow during this season. For NSE, the D2 parameter becomes dominant by a large margin in six out of ten subbasins, responsible for 28-70% of the model sensitivity in these sub-watersheds (Figure 3.7b). This is not the case for the other four sub-basins (Kettle, Limestone, Odei, and Weir) where binf remains the most influential factor controlling predictions of low flows. For the Footprint and Grass, Ws is also influential (~17%). Ws emerges as the third most important parameter in most of the sub-basins. Seasonal sensitivity of model parameters changes substantially; for example, in winter Ds and Dsmax, which control baseflow, become the most 77 sensitive parameters over all sub-basins whereas in spring and summer, binf and D2, respectively, play dominant roles in establishing streamflow (Appendix E). The binf parameter shows >45% ratios of factor sensitivity in spring for most of the sub-watersheds that reflects excess water availability for infiltration during snowmelt seasons. Autumn shows Dsmax and D2 as the most sensitive parameters since they are responsible in generating seasonal peak flows. Figure 3.7 Ratio of factor sensitivity (%) of IVARS 50 for each parameter at annual scale over all subwatersheds of the LNRB for the three model performance metrics (1981–2010) (a) KGE, (b) NSE, and (c) PBIAS. Ratio of factor sensitivity is estimated by normalizing IVARS 50 values in each case, so they add up to 100% for all parameters. 78 The total flow volume measured by PBIAS shows D2 as the most influential parameter that determines maximum water storage in soils and thus streamflow (Figure 3.7c). The ratio of sensitivity for this parameter exceeds 50% for nine sub-basins and nears 80% in three of them. Unlike the KGE and NSE cases, the baseflow parameters (Ds and Ds max) are not important for PBIAS because they have no effect on the total flow volume. Overall, for PBIAS, binf and D2 are influential parameters followed by Dsmax. This is due to the influence of all these parameters in controlling surface and subsurface water storages. In all sub-basins the depth of the second soil layer is more important than the third soil layer. This is perhaps because (a) the third layer is much thicker than the other two layers and (b) the second layer has a larger control on infiltration and ET. Seasonal sensitivity of model parameters for PBIAS is similar to KGE and NSE as the sensitive parameters (i.e., D2, binf, Ds, and Dsmax) are responsible for streamflow magnitudes and inter-linked with each other (Appendix E). Minimum and maximum values for streamflow over the LNRB’s sub-watersheds are 74274 mm year-1 (Sapochi) and 495-955 mm year-1 (Kettle) (Figure 3.8). Greater range is observed for streamflow in relatively small sub-basins such as the Weir, Kettle, and Footprint, as the VARS provides parameter samples within a broad range over small areas. This may also occur due to the combined overall uncertainty from other mass fluxes (i.e., ET, soil moisture, etc.). This analysis shows a range of potential variability that the mean annual streamflow has, which is intended to provide the reader with an estimation of the model uncertainty. 79 Figure 3.8 Annual streamflow sensitivity to parameter uncertainty for all ten LNRB sub-watersheds. The green dots show streamflow associated with the control run (calibration), the red lines show the median, the vertical black dotted lines show a range of minimum and maximum values excluding outliers, and the red + signs show the outliers defined as the values greater than 1.5 times the interquartile range of annual streamflow. 3.5.4. Uncertainty Assessment of the VIC Model Parameters using OLH In this section, we demonstrate applicability and performance of the OLH to identify and estimate VIC model parameters and their associated uncertainty bounds. Input forcing data and model structure are held constant in this analysis, so that the entire uncertainty in streamflow simulation may be attributed to VIC parameters. Uniform distributions are obtained on the parameter ranges from OLH and the behavioural parameter sets are used to generate parameter likelihood distributions. As stated earlier, we obtained 600 samples by using six major VIC model parameters in the OLH; therefore, parameter likelihood distributions derive from 600 VIC simulations (Figure 3.9). These distributions illustrate two points; first, the likelihood distribution of b inf, which is also the most sensitive parameter among others, nears the upper boundary of the predefined parameter range. This can be an 80 indication of a higher binf value that is important for these sub-basins. Second, the likelihood distribution for D2 captures only a small space of the pre-defined range, whereas D3 covers almost the entire parameter range. However, the hydrograph uncertainty bounds, which come from the top 10% of OLH runs, associated with these parameter ranges do not cover the expected number of observed streamflow values (dark blue region in Figure 3.10). This can be argued as a problem of over-conditioning the selected relationships between observed and modelled output (Bermúdez et al., 2017). The Footprint and Weir sub-watersheds have the widest uncertainty envelopes whereas the Burntwood, Limestone, Odei, Sapochi, and Taylor sub-watersheds show relatively narrow uncertainty bounds from the OLH simulation. This may depend on watershed area as parameter variation among a broad range of values over small sub-basins (e.g. Footprint, Weir) yields greater streamflow uncertainty than for relatively larger ones (e.g. Burntwood, Limestone, Odei, Sapochi, and Taylor). Even though the other 90% prediction uncertainty range (light grey region in Figure 3.10) captures all observations, it remains quite wide compared to observations and reveals a notable uncertainty range (0.5-1 mm day-1) in the model parameters as other conditions are static in this analysis. 81 Figure 3.9 Average likelihood distribution of the VIC parameters using 600 samples generated after the OLH sampling over all ten sub-watersheds of the LNRB. Red bars show the maximum likelihood of parameter range to the model performance metric (KGE) for the LNRB. 82 Figure 3.10 Streamflow prediction uncertainty associated with estimated parameters from the OLH. Top 10% (shown in blue colour) of OLH samples, based on KGE, used for the prediction of streamflow for all ten subwatersheds, water years 1981–2010. Note that y-axis scales vary between panels. Shaded area (grey colour) shows the envelope of VIC runs from 600 OLH samples. 83 3.6. Discussion 3.6.1. Input Data Uncertainty The underestimation in flows from the IDW-VIC and ANUSPLIN-VIC simulations reflect the precipitation undercatch and dry bias in these datasets over the LNRB (Lilhare et al., 2019). As the model resolution and other configuration (i.e., soil type, land use, etc.) are similar for all VIC simulations, different values of model performance metrics exhibit uncertainty associated only with input forcing datasets. These simulations show substantial disagreement in the runoff with observed hydrographs, especially in the Kettle, Limestone, Odei, Sapochi, and Weir sub-basins, owing to the dry bias and undercatch issues in the precipitation data. Consistent with our previous findings, the wet (warm) ERA-I and WFDEI precipitation (mean air temperature) over the LNRB in spring, summer and autumn induce more surface runoff and snowmelt that overestimate simulated flows (Figure 3.4) (Lilhare et al., 2019). Moreover, shifts in the hydrographs may be associated with warmer air temperatures over these sub-basins that cause earlier snowmelt runoff. Such variation in the simulated runoff, especially during the snowmelt period (April-July), is either due to the uncertain amount and timing of precipitation or air temperature in input forcing datasets (Lilhare et al., 2019). Other source of uncertainties in water budget estimation may be the dry (wet) bias in precipitation that results in poor calibration where model parameters cannot achieve optimal values due to less (excess) water availability. Consequently, these precipitation uncertainties among all sub-watersheds translate to a minimum (maximum) 14.0 (88.8) mm year -1 error in the water balance estimates. These results correspond well with those of Fekete et al. (2010) who found that the uncertainty in precipitation forcing translates to at least the same, or 84 typically more substantial, uncertainty in runoff and related water balance terms. The simulated TR uncertainty is higher in spring and summer than fall and winter, which is mainly due to the more substantial seasonal variation in inter-datasets precipitation and air temperature. However, there remains much uncertainty in air temperature records over the LNRB from the different forcing datasets. This uncertainty can be translated into interseasonal water balance estimation through runoff and ET processes, which are more sensitive to air temperature. The sub-basins that are more susceptible to ET loss with relatively more surface water coverage, i.e., the Grass and Gunisao will, therefore, have higher uncertainty propagation from input data. This represents both time and space dimensionality to the uncertainty and plays a critical role in climate change studies where changes in runoff are most important. Given that the domain average air temperature and precipitation differ across VIC forcings, the choice of VIC calibration and validation periods may induce uncertainty in water balance simulations. Therefore, calibration using different forcing data over 5-10 years may generate biases in simulated water balance conditions, as only a wet or dry period may be captured. Intercomparison of precipitation partitioning across various land surface models showed that specific representation and parameterizations for water balance components (i.e., ET and TR) were not consistent across models (Andresen et al., 2019). However, some models maintained similar runoff and precipitation ratios throughout the simulation; in contrast, VIC showed shifts from a runoff-dominated system to an ET-dominated system over permafrost regions in the Northern Hemisphere north of 45°N (Andresen et al., 2019). 85 3.6.2. Model Parameter Sensitivity and Uncertainty Assessment The sensitivity of model outputs to selected parameters is justified given the formulations of the variable infiltration and baseflow generation curve that form the foundation of the VIC architecture (Liang et al., 1994, 1996) and as these parameters are traditionally applied in model calibration (i.e., Elsner et al., 2010). As reported in previous studies, sensitivity to these parameters hold in both current and future climate scenarios (Bennett et al., 2018; Christensen & Lettenmaier, 2007; Demaria et al., 2007). A previous effort used various objective functions and found binf and D2 were the most sensitive parameters followed by the drainage parameter among ten VIC parameters across four American river basins of different hydro-climates (Demaria et al., 2007). Xie & Yuan (2006) manually varied four VIC parameters from ±10% to ±25% to perform SA over 12 watersheds in France, finding baseflow and soil depth as the most sensitive parameters. These studies either used manual analysis methods or limited objective functions at annual time scale to examine the parameter sensitivity. In contrast, we apply a more robust, automated, and efficient approach at seasonal and annual time scales to determine the parameter sensitivity and its seasonal importance over ten sub-arctic watersheds. Moreover, the SA method takes multiple objective functions into account that provide more robust estimates of parameter sensitivity. The high values of IVARS 50 for D2 are caused partly by the interaction of this parameter with other model parameters (e.g., soil profile and root depth) and its relatively large range of values. The Ds and Dsmax parameters influence baseflow, which have a higher impact on low flow predictions. Therefore, these parameters become important for NSE in some sub-basins (Grass and Footprint) as they control the timing of low flows. Moreover, in 86 PBIAS, high values of D2 indicate its considerable interaction on model responses as D2 characterizes the seasonal soil moisture behaviour but by no means binf as being perhaps also an important parameter over the LNRB. Overall, for PBIAS, binf and D2 are influential parameters followed by Dsmax. This is due to the influence of all these parameters in controlling surface and subsurface water storages. In all sub-basins D2 is more important than D3. This is perhaps because of (a) the third layer is much thicker than the other two layers and (b) the second layer has a larger control on infiltration and ET. There is general agreement between the NSE and KGE sensitivity experiments over most sub-watersheds, particularly in identifying the most influential parameters. Nonetheless, parameter sensitivity depends on metric choice and varies significantly according to model performance metrics. For example, the D2 parameter is quite important for KGE and PBIAS in most of the rivers but has slightly less impact on NSE. This is because D2 controls baseflow and thus the timing of flows, which is important particularly for peak flows represented by NSE. However, flow timing is not important when assessing total flow volume represented by PBIAS. Similarly, binf is less important for KGE and PBIAS but more influential for NSE over most sub-basins through its control of available infiltration capacity, thereby influencing peak flows and soil water volumes. Moreover, seasonality and wet or dry years may yield different SA results, which should be noted as a cautionary tale for using SA as a pre-calibration methodology. Consequently, to better understand the dominant controls on model behavior, multiple criteria should be considered. These results reinforce the well‐known conclusion that for most effective SA results, one should select SA criteria in alignment with the final goals of the modelling application (e.g., flood forecasting, drought analysis, or water balance assessment). Regardless of the 87 metric choice, often a limited number of parameters control most of the model response variations. This has important implications such as minimizing the dimensionality of the optimization process (i.e., calibration) through emphasis on a few influential parameters to generate reliable results. Even if fixed values for these influential parameters cannot be prescribed, any available information including observational data may reduce parameter ranges during calibration. This is generally true for all parameters and greatly increases the identifiability of our modelling application, which is often overlooked. Moreover, this also fits with the International Association of Hydrological Sciences’ (IAHS’) 23 unsolved problems in hydrology initiatives focused on understanding process changes, which control changing runoff response (Blöschl et al., 2019). Moreover, using these SA results, one can focus on specific model parameters and their value ranges, thus diminishing computational burdens, by fixing the value of non‐essential parameters. 3.7. Conclusion This exercise provides valuable new insights into the internal functioning of models and allows the provision of impactful recommendations for improving development and application of the VIC model. In this respect, we found that daily precipitation is more important than air temperature for annual and seasonal water balance estimates. The choice of model performance metric significantly affects the sensitivity assessment. Therefore, to obtain in-depth understanding of model behaviour, SA using multiple criteria should be adopted, which capture distinct characteristics of the model response. SA results can be used more effectively when aligned with the final goals of the model application (e.g., flood forecasting and drought monitoring). SA results depend on various factors such as hydro-climatic conditions, model configuration, input forcing, 88 landcover classes, initial state, vegetation parameters, etc., and these can have a large impact on model behaviour. We considered a full range of parameters that can influence their ratio of factor sensitivity if the range changes in other applications. SA can identify aspects of the model internal functioning that are counter‐intuitive and thus assist modellers to diagnose possible model deficiencies and make recommendations for end users. The calibration process identified a set of influential parameters that assists VIC users in reducing prediction uncertainty by providing a more robust, accurate and less computationally intensive calibration effort. Overall, parameters for the second soil layer depth and variable infiltration curve dominate the control of streamflow prediction in VIC followed by the Ds and Ds max parameters. The VIC community may prioritize these parameters during model calibration for similar physical and climatic environments. While this study focused on VARS sensitivity and OLH uncertainty analysis, a multi-criteria SA approach under various conditions may lead to improved understanding of model structure, reductions in prediction uncertainty, and more efficient parameter calibration. Potential future work could investigate the effects of initial or boundary conditions and/or other model variables such as soil moisture or ET in model sensitivity assessments. 89 CHAPTER 4: WARMING SOIL TEMPERATURES AND INCREASING BASEFLOWS IN RESPONSE TO RECENT AND POTENTIAL FUTURE CLIMATE CHANGE ACROSS NORTHERN MANITOBA 4.1. Abstract This study investigates climate change impacts on the hydrology and soil thermal regime of ten sub-arctic watersheds (northern Manitoba, Canada) using the Variable Infiltration Capacity (VIC) model. We utilize statistically downscaled and bias-corrected forcing datasets based on 17 climate models from phase 5 of the Coupled Model Intercomparison Project (CMIP5) to run the VIC model for three 30-year periods: a historical baseline (1981–2010), and future projections (2021–2050: 2030s and 2041–2070: 2050s), under representative concentration pathways (RCPs) 4.5 and 8.5. The CMIP5 Multi-Model Ensemble (MME) mean based VIC simulations show 15-20% increases and 10% decreases in the projected annual precipitation and snowfall, respectively over the southern part and >20% rainfall increase over the northern part of the domain by the 2050s. Snow accumulation is projected to decline across all sub-basins, particularly in the lower latitudes. Projected uncertainties in major water balance terms (evapotranspiration, surface runoff, and streamflow) are more substantial in the wetland and lake-dominated Grass and Gunisao watersheds than their eight counterparts. Future warming increases soil temperatures >2.5°C by the 2050s and results in 40-50% more baseflow. Further analyses of soil temperature trends at three different depths show the most pronounced warming in the top soil layer (1.6°C 30-year-1 in the 2050s), whereas baseflows increase significantly by 14.0% and 26.7% during the 2030s and 2050s, respectively. These results provide crucial information on the 90 potential future impacts of warming soil temperatures on the hydrology of sub-arctic watersheds in north-central Canada and similar hydro-climatic regimes. 91 4.2. Introduction Cold regions and snow-dominated river basins are particularly sensitive to warming air temperatures, as this can decrease seasonal and long-term water availability that otherwise provides a consistent and reliable source of streamflow during spring and summer melt periods (Aygün, Kinnard, & Campeau, 2019; Barnett, Adam, & Lettenmaier, 2005). These changes lead to increased streamflow interannual variability that has been observed across Canadian watersheds in recent decades (Déry, Hernández‐Henríquez, Burford, & Wood, 2009; Déry, Hernández-Henríquez, Owens, Parkes, & Petticrew, 2012). It is unclear, however, whether future climate changes will further impact the water availability and interannual variability in flows for cold regions, as previous studies focused solely on changes in mean annual and seasonal flows (e.g., Gelfan et al., 2017; MacDonald et al., 2018; Shrestha, Schnorbus, Werner, & Berland, 2012). Results from an ensemble of climate models project that runoff in eastern Canada and north-central Manitoba will increase 20– 30% relative to 1900-1970 by the middle of the 21 st century (Milly, Dunne, & Vecchia, 2005). At a regional scale, trends toward declining flows have also been reported for rivers in the Canadian Rockies (Rood, Samuelson, Weber, & Wywrot, 2005), the Hudson Bay basin (Déry & Wood, 2005a), and the Canadian Prairies (Westmacott & Burn, 1997; Yulianti & Burn, 1998). Permafrost covers about a quarter of the exposed land in the Northern Hemisphere (Zhang, Barry, Knowles, Heginbottom, & Brown, 1999) and, as defined by ground temperature, is sensitive to atmospheric climate changes. Changes in permafrost distribution and summer active layer depth will have significant impacts on hydrology, soil organic carbon (SOC), vegetation distribution, and infrastructure at high latitudes (ACIA, 2004; 92 Anisimov et al., 2001; Chapin et al., 1992; Nelson, 2003). For example, permafrost disappearance or deepening summer active layer depth could alter terrain and hydrologic conditions (Hinzman et al., 2005; St. Jacques & Sauchyn, 2009), modify the growth and distribution of vegetation (Jorgenson et al., 2001), increase SOC decomposition, and enhance carbon dioxide and methane emissions from the soil to the atmosphere (Goulden et al., 1998; Oechel et al., 1993). Climate models project substantial increases in mean surface temperature (~8ºC) across present-day permafrost areas of the Canadian landmass by the end of the 21st century under the RCP 8.5 scenario (Koven et al., 2013) (see Chapter 3, Section 3.3.3). This amplified warming affects permafrost temperatures, distribution, and conditions (Chadburn et al., 2017; Guo & Wang, 2016; Slater & Lawrence, 2013). Simulations from a process-based permafrost model driven by six general circulation model (GCM)-generated climate scenarios project that the area underlain by permafrost in Canada will decline by approximately 16-20% by 2090, relative to a 1990s baseline (Zhang, Chen, & Riseborough, 2008). Smith and Burgess (1998, 2004) categorized the present state and recent trends in ground temperatures, assessed the physical impacts of permafrost thaw and its sensitivity to climate warming in Canada. About half of Canada's landmass is classified as permafrost regions (Kettles, Tarnocai, & Bauke, 1997), which is about a third of the total permafrost area in the Northern Hemisphere (Zhang et al., 1999). In our study domain, the Lower Nelson River Basin (LNRB), sporadic discontinuous permafrost (between 10% to 50%) makes up ~68% of the total area (Figure 4.1) (Natural Resources Canada, 2010). The extensive discontinuous (i.e., between 50% to 90%), continuous (i.e., between 90% to 100%), and isolated patches of permafrost cover approximately 9%, 0.8%, and 16% of the total basin area, respectively. 93 To better quantify the impacts of climate warming on the LNRB’s hydrology and soil thermal conditions, we apply a distributed high-resolution hydrologic model over the LNRB (Lilhare et al., 2020) and utilize future climate projections from the fifth Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2011) to examine future water balance and soil thermal regimes. While the historical hydrology of the study area has been reported in previous studies, the current literature lacks detailed information on historical and future changes in soil temperatures over sub-arctic watersheds. Moreover, a quantitative evaluation of the CMIP5 projected climate changes and systematic analysis of GCM-driven hydrological model simulations have not been reported over the study domain. This study, therefore, investigates the impacts of increasing air temperatures coupled with changing precipitation on the LNRB’s water balance and the sensitivity of rainfall and snowfall to a changing climate. To better understand climate change impacts on the LNRB’s hydrology, and to provide a context for any detected changes in annual baseflow due to future warming, this paper has the second goal of investigating future changes in soil thermal regimes relative to the historical period. This will fill an important research gap on climate change impacts over sub-arctic regions and its contribution to changes in the overall hydrology of the study area. 4.3. Study Area, Data, and Methodology 4.3.1. Study Area In this study, we selected a downstream segment of the Nelson River system (i.e., the LNRB) fed by the Churchill River basin and Lake Winnipeg (Figure 4.1a). The LNRB covers ~90,500 km2 in area, collects all water from the upstream drainage area (~970,000 km2), and discharges into Hudson Bay. In the LNRB, the main stem river (Nelson) has less 94 seasonal flow variability due to regulation and large drainage areas. The LNRB has gentle slopes, with common channelized lakes and wetlands altering flow variability. Wetlands within the LNRB store significant volumes of water, cover large areas, and moderate streamflow responses to rainfall and snowmelt events. Permafrost and shallow soils limit groundwater flows, storage, and infiltration (Figure 4.1). A major part of the LNRB lies in the Canadian Shield, a recently glaciated area with shallow soil depth, leaving the precambrian igneous and metamorphic rocks of the Canadian Shield near the surface (Centre for Land and Biological Resources Research, 1996). The total exposed rock percentage over the LNRB is about 8% – this covers most of the north-western and southern part of the basin. The entire LNRB is dominated by 45% organic and 48% mineral soils and contains both silt and clay soil types (Centre for Land and Biological Resources Research, 1996). Due to the erosion from glaciers in the Holocene, organic content is low in LNRB soils, which are dominated by clays and poor drainage capacity. To increase the LNRB’s hydroelectric capacity, Manitoba Hydro manages flows in the watershed with two primary sources of streamflow regulation: the Lake Winnipeg Regulation and Churchill River Diversion (Figure 4.1b). The study domain experiences a sub-arctic continental climate with moderate precipitation and humidity, cold winters, and cool summers. The snow-free season spans May to October and most of the precipitation occurs during summer and falls as rain, contributing nearly two-thirds of the annual total precipitation. The precipitation peaks in July, which is the warmest month of the year with 16.2°C average daily air temperature. The average annual precipitation totals ~500 mm, while evapotranspiration (ET) over the LNRB attains ~300–350 mm annually with more annual surface water evaporation (Environment and Climate Change Canada, 2016). 95 Figure 4.1 Maps of the LNRB. (a) The Nelson River Basin, Churchill River Basin, and LNRB. (b) Major rivers and sub-watersheds within the LNRB, yellow triangles show the hydrometric stations used in this study, green circles denote current generating stations, the yellow circle shows a future generating station (currently under construction) by Manitoba Hydro and a red star indicates the Churchill River diversion. Basemap shows the permafrost distribution across the domain (Natural Resources Canada, 2010). Base map source: World light grey reference (http://services.arcgisonline.com/arcgis/services) 4.3.2. Data and Hydrological Model 4.3.2.1. Climate Models and Observational Data We use the CMIP5 data under the Representative Concentration Pathways (RCP) 4.5 and 8.5 scenarios (Taylor, Stouffer, & Meehl, 2011; van Vuuren et al., 2011). Data from CMIP5 simulations assembled, pre-processed and analysed by MacDonald et al. (2018) are used to force the VIC model for climate change analyses and cover ~90% of the full ensemble’s uncertainty range for precipitation and air temperature over the entire Hudson 96 Bay domain (Table 4.1). Note however that the simulation from the IPSL-CM5A-LR GCMs used by MacDonald et al. (2018) exhibits spurious trends and spatial patterns relative to other models, especially over the Weir River sub-basin and is excluded from analysis. This yields a total of 17 simulations (10 for RCP 4.5 and seven for RCP 8.5) for climate change analyses in this study. The bias correction of daily precipitation, maximum and minimum air temperature is based on the daily translation quantile-mapping approach by Mpelasoka and Chiew (2009) (see details in Stadnyk et al. (2019) for the climate scenarios preparation). Further, the Hydrological Global Forcing Data (HydroGFD hereafter) is used as a reference in the bias correction and forced to drive the VIC model over the LNRB (Berg et al., 2018; Lilhare et al., 2020). Table 4.1 CMIP5 simulations used in this study (MacDonald et al., 2018). Model Centre name RCPs ACCESS1-0 Australian Community Climate ACCESS1-3 and Earth System Simulator CMCC-CM 4.5 CMCC-CMC Euro-Mediterranean Center on Climate Change, Bologna, Italy CanESM2 Canadian Centre for Climate Modelling and Analysis 4.5 CRNM-CM5 National Center for Meteorological Research, France 4.5 GFDL-CM3 National Oceanic and Atmospheric Administration's Geophysical Fluid Dynamic Laboratory, USA inmcm4 Russian Institute for Numerical Mathematics MIROC5 University of Tokyo Center for Climate System Research, National Institute for Environmental Studies, Japan, Japan Agency for Marine-Earth Science and Technology Frontier Research Center for Global Change MIROC-ESMCHEM MIROC-ESM 4.5, 8.5 8.5 4.5 4.5 4.5, 8.5 8.5 8.5 MRI-CGCM3 Meteorological Research Institute, Japan 4.5, 8.5 NorESM1-M Norwegian Climate Center 4.5 97 4.3.2.2. Variable Infiltration Capacity Model and Simulation Strategy We use the semi-distributed Variable Infiltration Capacity (VIC) model (version 4.2.d) (Liang et al., 1994, 1996) to simulate hydrological processes in the LNRB. Our companion paper (Lilhare et al., 2020) provides detailed information on the VIC model implementation over the LNRB. The VIC model has been widely used for climate change research over various river basins (e.g., Wu et al., 2007; Adam et al., 2009; Cuo et al., 2009; Hidalgo et al., 2009; Elsner et al., 2010; Kang et al., 2016; Islam and Déry 2017; Sharma et al., 2018). Simulations of soil thermal fluxes are based on temperatures at the surface, at the bottom of the first layer, and at the bottom of the soil column whereby the surface temperature is solved explicitly as a part of the surface energy balance (Liang et al. 1999). Thermal fluxes are solved numerically via an explicit finite difference approximation of a general equation for heat flux through the soil column. Initial soil profile temperatures are estimated from the current air temperature. The lower boundary temperature is set to the average annual air temperature, and thermal conductivity and volumetric heat capacity of the soil layers are calculated at each time step from the revised soil water and ice contents. VIC conserves surface water and energy balances for regional-scale watersheds such as the LNRB and takes sub-grid variability into account by dividing each grid cell into elevation bands, land use, and land cover types (Cherkauer, Bowling, & Lettenmaier, 2003; Nijssen et al., 2001). Snow, energy, and water balances are defined for each grid cell separately (Gao et al. 2010). The VIC model is coupled with a standalone routing scheme that routes runoff from grid cells using a predefined channel network (Lohmann, Nolte- 98 Holube, & Raschke, 1996; Lohmann, Raschke, Nijssen, & Lettenmaier, 1998). Consequently, routed flow is extracted at user-defined outlet points of specific sub-basins. After the CMIP5 bias-corrected and downscaled data collection, the resulting fields are then regridded to 10 km (~0.10°) spatial resolution to match the VIC model grid using bilinear interpolation (Lilhare et al., 2019, 2020). In addition to daily precipitation and air temperatures, VIC also requires daily wind speeds. We obtain daily wind speeds from all CMIP5 models at a 10 m height above ground; however, these data are not downscaled nor bias-corrected, and the original data are simply regridded to the VIC model resolution (0.10° × 0.10°) (Lilhare et al., 2020). The VIC model monthly calibration (1981–1985, 1995–1999) and validation (1986– 1994) are conducted using the HydroGFD data (Lilhare et al., 2020) applied to ten unregulated rivers within the LNRB (Table 4.2). The VIC model runs in full energy and water balance mode at a daily time step, and three soil layers for each grid cell (Figure 4.2). In this application, the VIC model utilizes a no flux bottom boundary with the finite difference soil thermal solution (Cherkauer & Lettenmaier, 1999). The VIC setup is adopted from our companion paper (Chapter 3) where VIC was calibrated using an optimization process based on the Nash–Sutcliffe efficiency (NSE) coefficient (Lilhare et al., 2020; Nash & Sutcliffe, 1970). We use the NSE, Kling–Gupta efficiency (KGE) (Gupta et al., 2009), and Pearson’s correlation (r) coefficients (for simulated vs observed monthly streamflows) in addition to Percent Bias (PBIAS) to summarize model performance during calibration and evaluation. Additionally, we utilize monthly mean soil temperature observations from Sutton et al. (1999) (data available from June 1994 to December 1998 within the LNRB) to evaluate the VIC simulated soil temperatures. The soil temperature point observation (55.8°N 98.7°W 99 near the Burntwood River sub-basin) is measured at various depths (0.1 m, 0.2 m, 0.5 m, and 1.0 m). We extract soil temperature observations at 0.2 m, 0.5 m, and 1.0 m to compare corresponding VIC simulated soil temperatures, for a grid cell, at D1, D2, and D3, respectively, which range from 0.3 to 1.5 m (Table 4.3). High monthly NSE values for most of the LNRB’s sub-basins in the calibration and validation periods reveal the reliable application of the VIC model (Lilhare et al., 2020). The calibrated VIC model runs at a daily time step from 1981 to 2070 using each of the 17 CMIP5 GCM simulations (using historical, RCP 4.5 and 8.5 forcings) at 0.10° horizontal resolution. We initiate each simulation by running VIC for five years using the 1981–1985 meteorological forcings to allow model spin-up. The 1990s (1981 to 2010), 2030s (2021 to 2050) and 2050s (2041 to 2070) and changes relative to historical (1990s) conditions are examined in detail. Figure 4.2 Schematic representation of the overall methodology adopted to quantify climate change impacts on the soil thermal regime and LNRB’s hydrology using the VIC model. 100 Table 4.2 List of ten selected unregulated hydrometric stations, maintained by the Water Survey of Canada and Manitoba Hydro, for the VIC model calibration and evaluation with sub-watershed characteristics and mean annual discharge (Water Survey of Canada, 2016). Station name (Gauge ID) Latitude Longitude (°N) (°W) Mean subwatershed elevation (m) Drainage area (km2) Mean annual discharge (m3 s-1) Burntwood River above Leaf Rapids (05TE002) 55.49 99.22 302.4 5,810 22.9 Footprint River above Footprint Lake (05TF002) 55.93 98.88 273.8 643 3.2 Grass River above Standing Stone Falls (05TD001) 55.74 97.01 265.0 15,400 64.6 Gunisao River at Jam Rapids (05UA003) 53.82 97.77 260.9 4,610 18.0 Kettle River near Gillam (05UF004) 56.34 94.69 164.7 1,090 13.2 Limestone River near Bird (05UG001) 56.51 94.21 173.6 3,270 21.5 Odei River near Thompson (05TG003) 55.99 97.35 253.5 6,110 34.3 Sapochi River near Nelson House (05TG006) 55.90 98.49 259.1 391 2.2 Taylor River near Thompson (05TG002) 55.48 98.19 236.2 886 5.1 Weir River above the mouth (05UH002) 57.02 93.45 125.8 2,190 15.6 4.3.3. Methods A standalone model routes streamflow to the sub-basin outlet, which is then converted into precipitation equivalent units by dividing the corresponding sub-basin area. To aid our analyses of projected changes in streamflow, we evaluate several other variables such as total precipitation, air temperature, snowfall, rainfall, and VIC-simulated snow water 101 equivalent (SWE), and snowmelt. In the current implementation, VIC classifies total precipitation as 100% rainfall for temperature above 1.5°C, 100% snowfall for temperature below -0.5°C, and as a mixture of both for temperatures between these two thresholds. We perform analyses of the VIC variables using water years, defined as 1 October to 30 September of the following calendar year. We consider each bias corrected climate model output as a valid approximation to the real climate (represented by HydroGFD); therefore, a multi-model ensemble (MME hereafter) approach is suitable for the projected climate change analysis (Christensen & Lettenmaier, 2007; Fowler et al., 2007; Mishra & Lilhare, 2016). Spatial patterns of the MME outputs (i.e., precipitation and air temperature) are reasonably similar across the basin under RCP 4.5 relative to RCP 8.5, with lower magnitudes; therefore, we present projected change analyses mainly for the high emission (RCP 8.5) scenario. Here, results from the MME are summarized using four statistics: the temporal mean and standard deviation of each model and the multi-model mean and intermodel variation using standard deviations. For each 30-year analysis period (i.e., the 1990s, 2030s, and 2050s) and a daily variable x and model i, we calculate the mean climatology ( ̅ i) and interannual standard deviation ( ̅i). Further, the MME means of these quantities, ̅ and ̅, are calculated, while the intermodel variation is determined by the 5% to 95% model range in ̅ i. We consider the highest performing VIC model simulation for each sub-basin from Chapter 3. Further, we evaluate the climate change impacts on several water balance terms simulated by VIC, i.e., evapotranspiration (ET), surface runoff, baseflow, and streamflow over the LNRB. Climate change impact assessment is carried out at annual and seasonal (winter: DJF, spring: MAM, summer: JJA, and autumn: SON) scales for the 2030s and 2050s 102 against historical (1990s) simulations. Moreover, to evaluate the projected changes in surface and sub-surface soil temperatures, we analyse the seasonal and annual mean of key outputs such as soil temperatures at three different depths (Table 4.3). Nonparametric statistical methods are employed to investigate temporal changes in simulated soil temperatures. These methods are generally considered to be more robust compared to parametric ones and are less affected by the presence of outliers or issues of non-normality (Burford et al. 2009; Déry and Wood 2005a; Gan and Kwong 1992). The Mann-Kendall test (Kendall, 1975; Mann, 1945) is used to examine the statistical significance of trends and Sen’s method (Sen, 1968) is used to obtain estimates of the trend magnitude. The presence of serial correlation can confound the identification of true trends; therefore, the data are examined for significant autocorrelation and then decorrelated (Yue et al., 2002) before applying both Sen’s and Mann-Kendall’s methods. We assume a 5% significance level to assess the significance of temporal trends throughout the study. 4.4. Results 4.4.1. VIC Performance The NSE, KGE, and r values show generally good agreement between monthly simulated and observed streamflow, exhibiting NSE and KGE scores >0.6 for most subwatersheds (Figures 4.3 and 4.4). The calibrated streamflow, including the timing of runoff peaks and low-flow volumes, match closely those in observations for all sub-watersheds, except the Footprint, Gunisao, and Sapochi (Figure 4.5). Model performance is optimal at the Burntwood, Kettle, Limestone, Odei, Taylor, and Weir, with acceptable NSE, KGE, and r values in both calibration and evaluation (Figures 4.3 and 4.4). The LNRB sub-basins are well-calibrated considering the high-performance metrics, whereas relatively larger sub- 103 watersheds with high PBIAS, especially the Grass and Gunisao, are the only sub-basins with low NSE values. Additionally, poor model performance in the southern LNRB may reflect the irregularity in input forcings where interpolation procedures commonly do not take available stations into account due to discontinuous records and data availability (Lilhare et al., 2019). Moreover, discrepancies in simulated flows may be associated with the lake and wetland algorithm represented by the VIC model version adopted for this study. However, there are eight out of ten LNRB sub-basins that have comparatively fewer such features, exhibit satisfactory results (NSE > 0.5) and modelling capacity of the VIC model (Figures 4.3 and 4.4) (Lilhare et al., 2020). Table 4.3 The VIC model calibrated (1990s) values of the D1, D2, and D3 soil layer thicknesses for all subbasins of the LNRB. Sub-basin D1 (m) D2 (m) D3 (m) Burntwood 0.30 0.31 0.43 Footprint 0.30 0.41 0.69 Grass 0.30 0.63 1.49 Gunisao 0.30 1.05 1.33 Kettle 0.30 0.31 0.35 Limestone 0.30 0.32 0.36 Odei 0.30 0.31 0.36 Sapochi 0.30 0.30 0.39 Taylor 0.30 0.30 0.40 Weir 0.30 0.31 0.40 Mean 0.30 0.43 0.62 104 Figure 4.3 Spatial distribution of monthly calibration (1981–1985, 1995–1999) performance metrics, NSE, KGE, r (p < 0.05) and PBIAS, for ten selected sub-watersheds within the LNRB based on the HydroGFD-VIC simulations. 105 Figure 4.4 Spatial distribution of monthly validation (1986–1994) performance metrics, NSE, KGE, r (p < 0.05) and PBIAS, for ten selected sub-watersheds within the LNRB based on the HydroGFD-VIC simulations. 106 Figure 4.5 Simulated and observed daily runoff (mm day-1) averaged over water years 1981–2010 for the LNRB’s ten unregulated sub-basins. An external routing model is used to calculate runoff for the HydroGFDVIC simulations. Note that y-axis scales vary between panels. 4.4.2. CMIP5 Precipitation and Air Temperature Climatology The bias-corrected and spatially disaggregated CMIP5 GCMs show satisfactory baseline precipitation and air temperature mean climatology by its pattern (Figure 4.6). The MME mean precipitation shows a projected increase in winter (~18%), spring (~26%), and autumn (~6%) over the northern and north-eastern LNRB for the 2050s while summer 107 precipitation is projected to decrease over the same period, particularly in the southern LNRB under RCP 8.5 (Figure 4.7). The spatial patterns of precipitation change are reasonably similar under RCP 4.5 relative to RCP 8.5, with lower magnitudes overall (figure not shown). The area-averaged mean daily climatology shows precipitation between November and April, approaching a maximum of 1.5, 3.0, and 2.5 mm day-1 in winter, summer, and autumn, respectively, for the 1990s (Figure 4.6). The maximum intra-seasonal variability for CMIP5 MME precipitation takes place in summer and autumn when daily rates range from nearly 0 to 1.5 mm day-1. Daily precipitation rates are projected to change less in winter and late spring than in summer and autumn, when rates are projected to increase, particularly in June and July. The range of intermodel variation for peak precipitation during summer varies between 1.5 and 3.5 mm day-1. Figure 4.6 Area-averaged water year cycles of CMIP5 Multi-Model Ensemble (MME) mean daily precipitation (a, b) and air temperature (c, d) over the LNRB for the 1990s (black), 2030s (blue), and 2050s (red) (RCP 4.5, left and RCP 8.5, right). Shading corresponds to intermodel ranges. 108 The MME mean air temperature is projected to increase with higher magnitudes under RCP 8.5 than RCP 4.5, during all seasons, in the 2050s (Figure 4.6). The mean air temperature increases in all seasons under RCP 8.5 but with lower magnitudes and similar patterns to the RCP 4.5 scenario (Figure 4.8). Absolute increases are more prominent in winter than in spring and summer, over the entire LNRB, for both pathways; however, projected warming is especially notable in the northern LNRB. Over the 1990s, air temperature remains below -10.0°C between November and March on average, then rises above 0.0°C in late spring inducing snowmelt across the LNRB (Figure 4.6). In the 2030s and 2050s, while the annual variation of air temperature remains similar to the baseline, a consistent increase from 2.0 to 4.0°C is projected for the basin average air temperatures. Figure 4.7 Projected changes (2050s-1990s; %) in the spatial distribution of precipitation during (a) winter, (b) spring, (c) summer, and (d) autumn for RCP 8.5. 109 Figure 4.8 Projected changes (2050s-1990s; °C) in the spatial distribution of air temperature during (a) winter, (b) spring, (c) summer, and (d) autumn for RCP 8.5. 4.4.3. Projected Changes in VIC Simulated Variables For the BaySys project, we were interested in analysing the projected changes in water balance variables and soil temperature under a more extreme climate scenario (i.e., RCP 8.5). The spatial patterns of VIC simulated variables such as rainfall, snowfall, SWE, snow cover, and runoff for RCP 4.5 are similar to those for RCP 8.5, but with slightly lower magnitudes in the 2030s and 2050s (not shown); therefore, here we analysed only simulations under the RCP 8.5 scenario. While the total precipitation substantially increases (~15%) in the 2050s, the simulations project future increases of ~16% in mean annual rainfall and nearly 6% decreases in mean annual snowfall (Figure 4.9). The partitioning of MME mean daily precipitation into rainfall and snowfall exhibits substantial increases in summer rainfall towards the 2050s across the LNRB (Figure 4.10: a, b). The increase in rainfall persists, especially in June and July, with >6.0 mm day-1 in the latter half of the 21st century (Figure 4.10a). Daily snowfall shows higher values during the early winter months 110 (November and December). Warming temperatures and reduced snowfall produce considerable changes in snow accumulation and ablation seasons (Figure 4.10: c, d). Figure 4.9 Future changes (2050s-1990s; %) in the spatial distribution of mean annual (a) rainfall, (b) snowfall, and (c) total precipitation for RCP 8.5. Day-to-day SWE accumulation declines while its seasonality shifts over the 2050s, with more noticeable changes relative to the 1990s across the LNRB (Figure 4.10c). The length of the snow accumulation season is nearly 15 to 30 days shorter in the 2050s, relative to the 1990s. The magnitude and seasonality of snowmelt, which is responsible for generating high flows typically in April or May, show earlier snowmelt freshets in the future and increased snowmelt volume (Figure 4.10d). The increase in daily runoff is much more evident in the 2050s due to earlier snowmelt, which contributes to spring and summer flows. The future evolution of runoff seasonality shows that while the dominant snowmeltgenerated peak flow shifts earlier by ~15 days, spring and summer runoffs are most pronounced during the 2050s. Furthermore, the low-flow period occurs historically in winter in all LNRB sub-basins (Figure 4.11). Streamflow starts to increase in spring, peaks in summer, and eventually declines in autumn. While this overall seasonal pattern continues in the future, the timing of seasonal flows shows a shift of approximately 15 to 20 days, and 111 flow occurs earlier for RCP 8.5 with the intermodel spread of nearly five days (Figure 4.11). Flows start increasing in mid-spring, followed by an earlier and steeper recession in midsummer. Daily mean runoff of the Footprint, Kettle, Limestone, Odei, Sapochi, Taylor, and Weir sub-basins exhibits similar features; however, a noticeable difference in the Grass and Gunisao is the change in timing of runoff peak (Figure 4.11). The peak in these two rivers varies between 15 and 20 days (median values) for RCP 8.5 with a range of nearly 10 days. Figure 4.10 VIC partitioned CMIP5 MME (1981-2070) mean total precipitation into (a) daily rainfall and (b) snowfall; CMIP5-VIC-simulated (c) mean ΔSWE (daily SWE rate) and (d) snowmelt, in (c) values greater than 0 represent snow accumulation, while those below 0 indicate snow ablation; VIC-simulated MME (e) mean daily runoff for the LNRB. All values are spatial averages over the basin under RCP 8.5 scenario. Units are mm day-1. Runoff units are an equivalent regional average rainfall rate rather than a discharge rate. 112 Figure 4.11 The VIC-CMIP5 MME mean daily runoff for the LNRB’s sub-basins. Black, blue, and red curves represent the daily climatology for the base period (1990s), 2050s RCP 4.5, and 2050s RCP 8.5, respectively. Shading corresponds to intermodel uncertainties. 4.4.4. Intermodel Uncertainties in the Projected Water Balance Terms In all LNRB sub-basins, ET is projected to increase by 5% and 15% in the 2030s and 2050s, respectively under both scenarios with 5-10% intermodel uncertainty (Figure 4.12: a, b). The large intermodel uncertainties (both in the 2030s and 2050s) in ET estimates, especially in the Grass, Gunisao, and Weir, highlight less robust changes with respect to 113 historic climate. Under the RCP 8.5 scenario, several basins are projected to experience >10% increase in ET with considerable intermodel variation, which can be attributed to substantial warming, while the modest rise in rainfall reveals excess water over these subbasins. Moreover, the ensemble mean surface runoff is projected to increase in all sub-basins in the 2030s and 2050s under both scenarios (Figure 4.12: c, d). Projected surface runoff is higher for the Gunisao and Grass with more substantial intermodel variation than that of other rivers. The ensemble mean surface runoff is projected to increase more than 25% in all sub-basins under the RCP 8.5 scenario by the 2050s. Similar to the ensemble mean projected change in surface runoff, annual streamflow is projected to increase in the majority of subbasins (Figure 4.12: e, f). However, in the 2030s, there are considerable intermodel variations for projected changes in surface runoff and streamflow. For instance, under both scenarios, only three river basins (Gunisao, Grass, and Weir) are projected to experience large increases (50-100%) with substantial intermodel uncertainty, whereas others are projected to experience modest (10-15%) increases in annual streamflow. Results show that the majority of LNRB sub-basins are projected to experience a rise of 20-50% in annual streamflow, with more robust increases likely in the 2050s. 114 Figure 4.12 The MME mean projected changes (%) in (a, b) ET, (c, d) surface runoff, and (e, f) streamflow at annual scale under RCP 4.5 and 8.5 scenarios for the 2030s and 2050s with respect to the base period (1981– 2010). Bars show MME mean change, while error bars show intermodel variation using one standard deviation. 115 4.4.5. Projected Changes in Soil Thermal Regime In this section, we analyse changes in soil temperature profiles simulated by the VIC model at three different depths. The calibrated soil depths D2 and D3 vary for each sub-basin within the range of 0.3-1.5 m, while D1 (0.3 m) remains static across the LNRB (Table 4.3). The VIC model simulated monthly mean soil temperatures at D1, D2, and D3 over a 10 km grid are evaluated against point observations (Sutton et al., 1999) that show satisfactory results (monthly NSE ranges from 0.6 to 0.8 and r > 0.9) (Figure 4.13). However, VIC version 4.2.d exhibits issues with the computation of soil layer temperatures when soil thermal nodes do not reach the bottom of the soil column (Hamman, Nijssen, Bohn, Gergel, & Mao, 2018). Consequently, modelled soil temperatures for D2 and D3 match the seasonality reasonably well but do not capture the peaks and troughs when compared to the observations (Figure 4.13). Furthermore, the VIC simulated average annual soil temperature climatology is consistent with observations at 0.3 to 0.6 m depths that show ~2.0°C during the 1990s over northern Manitoba and in the LNRB (Beltrami, Gosselin, & Mareschal, 2003; Rouse, 1984; Zhang, Chen, Smith, Riseborough, & Cihlar, 2005). The MME mean annual soil temperatures for all three layers show changes of 1.0°C and 2.2°C during the 2030s and 2050s, respectively (reference period 1990s) under RCP 8.5 across the study domain (Table 4.4). The soil temperatures show the greatest changes in summer (1.8-3.4°C) followed by autumn (1.5-2.4°C), spring (0.7-2.7°C) and winter (0.5-1.1°C) during the 2030s and 2050s (Table 4.4). These long-term soil temperature changes (1.0-2.2°C) exhibit a positive relationship with baseflow changes (14.8%-19.1%) across the study domain (Tables 4.4, 4.5 and 4.6). 116 Table 4.4 The MME mean projected seasonal and annual changes, 2030s-1990s and 2050s-1990s, in VIC simulated soil temperatures at D1, D2, and D3 under RCP 8.5 across the study domain. (2030s-1990s) (℃) (2050s-1990s) (℃) D1 D2 D3 D1 D2 D3 Winter 0.5 0.6 0.5 1.2 1.1 0.8 Spring 1.6 1.1 0.7 2.7 2.0 1.2 Summer 1.9 2.0 1.8 3.2 3.4 3.0 Autumn 1.5 1.5 1.5 2.2 2.2 2.4 Annual 1.4 1.3 1.1 2.3 2.2 1.9 Timeline Table 4.5 The MME mean projected annual trend magnitudes (30 years) of VIC simulated D3 soil temperatures (°C 30-year-1) and baseflows (mm 30-year-1) under RCP 8.5 for the 1990s, 2030s, and 2050s. Trend plots are provided in Figure 4.18 (all significant). D3 soil temperature (°C 30-year-1) Baseflow (mm 30-year -1) 1990s 2030s 2050s 1990s 2030s 2050s Burntwood 0.7 0.9 1.2 16.1 15.1 17.1 Footprint 0.7 1.0 1.2 16.6 16.0 18.1 Gunisao 0.5 1.1 1.4 11.9 12.6 8.6 Grass 0.5 0.9 1.1 19.0 23.0 17.8 Kettle 0.6 0.9 1.1 17.2 22.4 14.0 Limestone 0.6 0.9 1.1 15.2 22.2 12.4 Odei 0.6 1.0 1.2 15.4 19.5 16.0 Sapochi 0.6 1.0 1.2 15.3 19.1 15.9 Taylor 0.6 1.0 1.2 15.4 18.9 16.1 Weir 0.6 0.9 1.2 15.0 22.0 12.2 Mean 0.6 1.0 1.2 15.7 19.1 14.8 Sub-basin 117 Figure 4.13 The VIC model simulated (red) vs observed (black) monthly mean soil temperatures (June 1994 to December 1998) at three depths (D1, D2, and D3) over one grid cell (latitude: 55.8°N, longitude: 98.7°W) within the LNRB. Black dashed line shows 0°C. 118 Figure 4.14 The MME mean projected (a) winter (b) summer and (c) annual trends of VIC simulated soil temperatures at D1 (red), D2 (black), and D3 (blue) under the RCP 8.5 scenario over the LNRB sub-basins. The solid (dashed) trend lines show significant (insignificant) trends for the 1990s, 2030s, and 2050s and shading corresponds to intermodel uncertainties. Black dashed lines show 0°C. Trend magnitudes are provided in Table 4.6. 119 Table 4.6 The MME mean projected seasonal and annual trend magnitudes (all significant) (30-years) of VIC simulated soil temperatures at D1, D2, and D3 under RCP 8.5 across the study domain. Trend plots are provided in Figure 4.14. Timeline D1 (℃ 30-year-1) D2 (℃ 30-year-1) D3 (℃ 30-year-1) 1990s 2030s 2050s 1990s 2030s 2050s 1990s 2030s 2050s Winter 0.4 0.9 1.1 0.4 0.8 0.9 0.4 0.5 0.6 Spring 1.0 1.8 1.7 0.9 1.3 1.4 0.5 0.6 0.8 Summer 0.7 1.5 2.2 0.8 1.6 2.2 0.9 1.5 2.1 Autumn 0.5 1.0 1.0 0.6 1.0 1.1 0.5 1.1 1.2 Annual 0.7 1.3 1.5 0.7 1.2 1.4 0.6 1.0 1.2 A significant increase in the MME mean annual and seasonal soil temperatures appears in all soil layers (Figure 4.14), with the MME mean maximum positive changes of 1.3°C and 2.3°C during the 2030s and 2050s, respectively across the domain (Table 4.6). The annual trend analyses reveal D1 as the warmest among all other soil layers with significant positive trends ranging from 1.4 to 1.6°C 30-year-1 during the 2030s and 2050s (Table 4.6 and Figure 4.14). The MME mean soil temperature under RCP 8.5 projects summer warming of ~5℃ (across all sub-basins), but lesser changes are projected for winter, spring, and autumn. In spring and autumn, D2 and D3 show similar changes and overlap during the projected period, whereas D1 is the warmest soil layer in spring and coldest in autumn (figure not shown). Trend estimates of mean annual soil temperature for the historic period suggest significant warming (range from 0.6 to 0.7°C 30-year-1) for all layers and sub-basins (Table 4.6 and Figure 4.14). However, for future climate conditions, the maximum positive trends are observed for the D1 soil temperature in winter and summer where the trend magnitudes lie within the 0.7-1.3°C and 1.8-2.2°C 30-year-1 range, respectively. The increase in summer soil temperatures exceeds that for winter. The Gunisao sub-basin during winters in the 1990s exhibits insignificant decreasing trends (-0.1°C 30-year-1) for D1 and D2. 120 4.4.6. Long Term Projected Changes in the LNRB’s Hydrology and Soil Temperatures To better understand climate change impacts on the overall hydrology of the domain, the MME mean annual projected changes, relative to the 1990s, in precipitation, air temperature, surface runoff, ET, total soil moisture, soil temperature (D3), and baseflow are then intercompared spatially under RCP 8.5 (Figures 4.15, 4.16, and 4.17). Precipitation and air temperature changes are relatively smaller in the 2030s than 2050s, where the LNRB attains positive changes of 4-8% and 2.0-3.5°C, respectively across the LNRB (Figure 4.15). These positive changes increase in the 2050s for both precipitation (10-15%) and air temperature (3.5-4.5°C) and show wetting and warming trends across the domain. The northern LNRB shows more modest increases (10%) whereas the south-eastern and central sections exhibit 20-40% decreases in surface runoff during the 2030s and 2050s (Figure 4.15: e, f). Moreover, other water balance variables such as ET and total moisture of the soil column show projected increase (5-25%) and decrease (5-15%), respectively across most of the LNRB in the 2030s and 2050s (Figure 4.16). Increasing ET reduces the surface and subsurface moisture content; therefore, total soil moisture decreases >15% over the south-eastern LNRB in the 2030s and 2050s. ET increases from 10-15% to 20-30% whereas soil moisture decreases <10% in the 2050s over north-central and southern sections of the LNRB. The soil temperature climatology for the 1990s ranges from 0.5-1.5°C and 2.0-3.0°C over the northeastern and south-western part of the domain, respectively. Projected changes in soil temperature are more modest in the 2030s than in the 2050s, when the LNRB experiences rises of ~0.5-1.0°C across the domain (Figure 4.17). The D3 soil temperatures increase in the 2050s and range from 2.0-3.5°C across the LNRB. An increase in soil temperature reduces 121 the presence of frozen ground across the domain; therefore, these soil warming trends substantially increase infiltration that result in the highest changes in baseflow of 50-60% and >80% over the north-eastern LNRB in the 2030s and 2050s, respectively. The rest of the domain experiences >40% increases in baseflow during both future timelines. Figure 4.15 The MME mean projected changes (reference period 1990s) (a, b) % for precipitation, (c, d) °C for air temperature, and (e, f) % for surface runoff in the 2030s and 2050s under the RCP 8.5 scenario over the LNRB. Figure 4.16 The MME mean projected % changes (reference period 1990s) in (a, b) ET and (c, d) total soil moisture for the 2030s and 2050s under RCP 8.5 scenario over the LNRB. 122 Over the 2050s, the MME mean annual baseflow continues to increase due to the warmer D3 soil temperatures under the RCP 8.5 scenario that may reduce the presence of frozen ground and degrade permafrost (Figure 4.18). A significant positive relationship is observed between the MME mean annual baseflow and D3 soil temperature over all LNRB sub-basins. The MME mean annual D3 soil temperatures are projected to increase between 1.0℃ to 1.2℃ 30-year-1 during the 2030s and 2050s (Table 4.5). Similarly, baseflows among the LNRB’s sub-watersheds are projected to increase between 13.0 mm to 14.8 mm 30-year-1 during the 2030s and 2050s. The MME mean annual D3 soil temperatures exhibit an overall increase of ~2°C across the study domain by 2070 (Figure 4.18). Late winter and spring flows are expected to increase substantially in the future due to warming soil temperatures that result in increased baseflow (Table 4.7). The simulated MME mean seasonal baseflow increases during the 2030s and 2050s (reference period 1990s) in winter and spring (Table 4.7). With the absence of snowmelt in winter, total runoff derives mainly from subsurface flow, which shows 80.9% and 177.2% increases over the LNRB during the 2030s and 2050s, respectively. Baseflow exhibits the most considerable change in spring relative to other seasons during the 2030s and 2050s attaining 213.5% and 400.3%, respectively (Table 4.7). Snowmelt contributes significantly to total runoff during spring; however, projected runoff increases in this season are mainly due to a projected increase in precipitation and soil temperature. Summer shows minimal projected changes in baseflow during the 2030s and 2050s, where reductions appear for the Odei, Sapochi, and Taylor sub-basins, ranging from -3.3% to -1.9% under the RCP 8.5 scenario (Table 4.7). These changes indicate that in summer, water originates from projected increases in precipitation and resulting surface runoff in these sub-basins. Furthermore, the 123 VIC simulated historic and projected precipitation partitioning in water balance variables suggests relatively greater precipitation partitioning to surface runoff and baseflows than ET in the 2030s and 2050s when compared to the 1990s (Table 4.8). The percentage of evaporative fraction (ET/P) for the historic climate reaches 73% but diminishes to 70% by the 2050s. The percentage of surface runoff (SR/P) and baseflow (BF/P) fractions increase ~2% by the 2050s (Table 4.8). Figure 4.17 The MME mean projected annual climatology (1990s) of VIC simulated (a) D3 soil temperature (°C) and (d) baseflow (mm) under RCP 8.5 scenario over the LNRB. The second and third rows show changes (reference period 1990s), (b, c) °C for soil temperature, and (e, f) % for baseflow, in the 2030s and 2050s, respectively. 124 Table 4.7 The MME mean projected seasonal and annual changes, 2030s-1990s and 2050s-1990s, in VIC simulated baseflow under RCP 8.5 over all sub-basins across the study domain. % Change in baseflow (RCP 8.5), reference period: 1990s Sub-basins Winter Spring Summer Autumn Annual 2030s 2050s 2030s 2050s 2030s 2050s 2030s 2050s 2030s 2050s Burntwood 128.5 277.8 185.2 355.6 1.9 13.4 -2.8 17.0 18.9 47.4 Footprint 137.2 300.3 176.0 335.4 1.8 12.3 0.2 18.4 19.4 46.8 Gunisao 55.7 107.7 158.4 254.7 3.7 13.4 -1.3 3.9 15.4 30.7 Grass 41.5 84.9 272.6 511.1 8.0 21.5 13.3 23.9 21.4 42.2 Kettle 37.3 77.4 241.5 458.3 12.2 29.7 12.2 23.9 22.6 45.8 Limestone 42.6 89.7 241.2 471.2 11.4 34.2 14.0 31.6 22.9 50.7 Odei 105.2 240.7 202.0 375.0 -1.9 8.2 11.7 27.2 18.6 41.8 Sapochi 102.6 235.4 200.2 370.8 -3.0 6.9 11.4 26.7 18.3 41.4 Taylor 102.6 233.1 198.5 366.6 -3.3 6.6 11.1 26.1 18.1 40.9 Weir 55.3 125.6 259.0 503.9 5.9 27.2 11.1 29.4 19.9 47.8 Mean 80.9 177.2 213.5 400.3 3.7 17.3 8.1 22.8 19.6 43.5 Figure 4.18 The MME mean projected annual trends of VIC simulated D3 soil temperature (red) and baseflow (black) under RCP 8.5 scenario over the LNRB sub-basins. The solid trend lines show significant trends for the 1990s, 2030s, and 2050s and shading corresponds to intermodel uncertainties. Trend magnitudes are provided in Table 4.5. 125 Table 4.8 The MME mean historic (1990s) and projected (2030s and 2050s) annual precipitation (P) partitioning in % within water balance components (ET, surface runoff (SR), and baseflow (BF)) under RCP 8.5 for all ten sub-basins. 1990s (%) 2030s (%) 2050s (%) Sub-basins ET/P SR/P BF/P ET/P SR/P BF/P ET/P SR/P BF/P Burntwood 75 16 9 73 17 10 72 17 11 Footprint 72 18 10 70 18 11 70 18 12 Gunisao 81 9 11 79 9 12 77 10 13 Grass 71 16 13 69 16 15 68 16 16 Kettle 72 17 11 70 17 13 68 17 15 Limestone 73 17 10 71 17 12 69 18 13 Odei 71 18 11 69 18 13 68 19 14 Sapochi 72 17 11 70 18 12 68 18 13 Taylor 72 18 11 70 18 12 68 18 14 Weir 73 18 10 71 18 12 69 18 13 Basin average 73 16 11 71 17 12 70 17 13 4.5. Discussion 4.5.1. Projected Changes in the LNRB’s Hydrology Projected precipitation changes are more variable, especially between the southern and north-eastern parts of the LNRB, due to an overlap between model grids. Changes in the partitioning of precipitation between rainfall and snowfall significantly impact the seasonal SWE distribution in the LNRB. During late spring, snowmelt events increase overall in the future, likely due to more frequent warming episodes across the LNRB by the 2050s. The changing phases of projected precipitation in spring and winter induce seasonal shifts in 126 runoff during the 2050s (RCP 8.5). Moreover, the early decline of future SWE storage, due to snowpack depletion, in most of the sub-basins induces seasonal shifts in spring and summer runoffs that, in turn, lead to lower runoff in autumn. Furthermore, a shift from snowfall to rainfall leads to a decrease in surface runoff over the south-eastern LNRB, consistent with other previous studies (Berghuijs, Woods, & Hrachowitz, 2014; McCabe, Wolock, & Valentin, 2018). However, spatial plots indicate an increase (>20%) in surface runoff by the 2050s over some parts of the northern LNRB that is concomitant with increased rainfall and snowfall (Figures 4.9: a, b and 4.15f). The uncertainty in projected water balance estimates increases in the 2050s due to the large intermodel variation. The higher intermodel variations arise from disagreement between GCMs in their simulations of precipitation and air temperature over the study domain, thus leading to less robust changes in surface runoff and streamflow. Previous studies argued that there remains considerable uncertainty in projections of seasonal precipitation across the globe, which could largely be attributed to a lack of understanding in physical processes by climate models (e.g., de Elía et al., 2008; Fekete et al., 2004; Gelfan et al., 2017; Turner and Annamalai, 2012). It is expected that RCP differences will be more significant by the end of the 21st century than they are in the 2050s (Islam, Déry, & Werner, 2017; Mishra & Lilhare, 2016; Poitras, Sushama, Seglenieks, Khaliq, & Soulis, 2011). Despite uncertainties, our results fill existing gaps in projected precipitation transitioning, with relatively less snowfall over the LNRB. In cold regions, such changes are strongly coupled with projected changes in precipitation and air temperature along with latitudinal differences in the sub-basins. Future water availability in most of the sub-basins is likely to be driven by changes in snowmelt and precipitation. Part of this change appears to relate to increased precipitation, air, and soil 127 temperatures. Part of the temperature link may involve increased snowmelt in the headwaters of these river systems. Figure 4.19 Schematic representation of the projected hydrology in warmer and wetter climate over the LNRB. Positive and negative signs show indicate the direction of changes, respectively in the water balance terms and thermodynamic states relative to the 1990s. Future warming may lead to increased atmospheric losses; therefore, relatively higher increases in annual ET can be attributed to significant increases in precipitation and air temperature during late summer and autumn (Figure 4.19). Tam et al. (2019) reported increased water deficits during spring to autumn over the mid-latitudes and Nelson River, and the ET is controlled by both water and energy availability that can limit atmospheric demands. Our analyses show that increases or decreases in surface runoff under the projected climate can be attributed to changes in precipitation type (solid and liquid) and play a crucial role in the LNRB’s hydrology. Our findings over the LNRB are consistent with the existing literature, which indicates that mid-latitudes may become drier in a warmer climate (Cook, Smerdon, Seager, & Coats, 2014; Held & Soden, 2006; Jeong, Sushama, & Naveed Khaliq, 128 2014; Swain & Hayhoe, 2015). However, there are considerable uncertainties in observations and climate model performances over northern Manitoba including the LNRB (de Elía et al., 2008; Lilhare et al., 2019; Tam et al., 2019), which translate to the overall uncertainty in climate change impact assessments. 4.5.2. Warming Impacts on Soil Temperatures and Baseflows Our results show that air temperature mainly drives precipitation partitioning in a warmer and wetter climate. Surface and sub-surface temperature warming may lead to the formation of unfrozen zones across the LNRB, thus providing more baseflow in cold seasons. Soil temperature warming also allows for more infiltration in a projected future climate (Figure 4.19). This results in increased ET and baseflow and contributes to the overall increased river discharge. Results show that changes in runoff depend on how changes in individual water balance components vary spatially, and the extent to which ET varies across different areas of the landscape is energy limited or water limited (Koster and Suarez 1999; Christensen et al. 2008; Gao et al. 2011). ET is energy limited in areas where there is ample supply of water with respect to the available energy: in these areas, increases in available energy in a warmer climate yield increases in ET. However, in many areas across the LNRB, ET is limited by the spatial distribution of lakes and wetlands across the LNRB or the amount of vegetation-available water during the growing season, which is controlled by changes in winter and summer precipitation, rain–snow partitioning, and timing of snowmelt. For example, the Grass sub-basin, which has a high wetland concentration, shows ~12% increase in ET and the highest positive changes in baseflow among all other sub-basins. Overall, soil temperature warming yields increases in ET and streamflow across the domain where such 129 cold and warm season interactions play an important role in partitioning of precipitation between water balance variables. Future warming impacts on the overall water cycle are evident in the present study (Figure 4.19). Increased soil temperature and earlier snowmelt can result in increased stream temperatures, especially during summer, and affect aquatic ecosystems (Ficklin, Stewart, & Maurer, 2013). Our results suggest that, as warming continues, air and soil temperature may become a more significant control of ET, surface runoff, and baseflow for arctic and subarctic watersheds. Comparable results of increasing mean annual and winter flows have been addressed in northern Eurasia, under permafrost and non-permafrost regions (Smith, Sheng, & MacDonald, 2007), and over the permafrost-underlain Yukon River basin (Walvoord & Striegl, 2007) and the Northwest Territories in Canada (St. Jacques & Sauchyn, 2009). Increasing low flows in the LNRB suggest a broad-scale recent mobilization of subsurface water into rivers. This increase of baseflow, probably from soil and groundwater inputs into rivers, suggests permafrost thawing (Smith et al., 2007; Walvoord & Striegl, 2007). These processes are adequately supported by physical evidence of permafrost degradation, such as thermokarst processes, warmer springs and streamflow seasonality, in the Hudson Bay drainage basins including the LNRB (Stadnyk et al., 2019). The LNRB has experienced significantly increasing maximum and minimum seasonal and annual temperatures in the historical and future timelines, which would promote permafrost decay (Zhang, Vincent, Hogg, & Niitsoo, 2000). We quantified that the greater soil temperatures resulting from future warming could alter late spring and summer streamflow by allowing for more infiltration, delaying the runoff response to precipitation events across selected watersheds. Our results suggest that future soil warming may yield increased soil infiltration, stream 130 network interconnectivity, subsurface water movement, and shifts in water storage from lakes and wetlands to soil and groundwater and then to increased baseflow (e.g., Grass sub-basin). Regardless of the exact mechanism in each circumpolar subarctic locality, the findings of increased baseflows throughout much of the LNRB presage a growing role of groundwater processes in the high latitude water cycle under projected climate. As has been suggested by other researchers (e.g., Smith et al., 2007), this could mean shifts from “above-ground” to “below-ground” water storage and increased unsaturated zone storage, soil infiltration, and groundwater movement in permafrost and seasonally frozen ground regions. 4.6. Conclusions Climate change impacts may induce considerable hydro-climatic alterations in midlatitude rivers across the globe. Results obtained in this study provide useful information related to these changes in key hydrological terms under projected warming, rainfall, and snowfall changes over the LNRB. In such a basin, where multiple hydroelectric facilities are situated, understanding of water balance alterations due to climate change is essential. Using the CMIP5 projections and VIC hydrological model, this study demonstrates that in a warmer climate, changes in rainfall and snowfall play a crucial role in altering water cycle components (i.e., ET, surface runoff, and streamflow) at the sub-basin scale. Overall, our results suggest that the projected rainfall and snowmelt increase rapidly, which contribute to substantially more runoff. Despite intermodel uncertainty, many of the LNRB sub-basins are projected to shift in the warmer and wetter climate. The percentage of an evaporative fraction (ET/P) for the current climate is 73% over the LNRB, indicating that water loss in this region is dominated by evaporative processes and not runoff, which accounts for only 27% of the annual precipitation. However, the evaporative fraction percentage decreased from 73% in 131 the current climate to 70% in the future, indicating a trend toward a wetter climate and increased runoff despite an optimistic precipitation scenario. Projected seasonal and annual trends in soil temperature are most pronounced by the 2050s across the LNRB. Warmer temperatures in the soil column will allow more infiltration, contributing to more groundwater movement into river channels in autumn, seen as increased early winter discharge. Streamflow and water availability appear to increase during summer in the 2050s, and water resources managers may experience more significant year-to-year variability and uncertainty in flows due to rapid snowmelt, more rainfall, reduced snowpacks, and warming soil temperatures. While results reported in this study are not precise representations of projected runoff, streamflow, and other water balance changes due to several computational limitations, dataset and modelling uncertainties (i.e., land cover, future climate data, and regulation), they nevertheless provide valuable insight to the projected hydrological state of the LNRB. Moreover, future work is essential to understand the balance between water availability and water demands in different seasons in the LNRB. To reduce uncertainty in future projections, more climate and hydrologic model ensemble members will be desired with better representation of water diversions, and parameterization of permafrost and snowrelated processes over the LNRB. 132 CHAPTER 5: CONCLUSIONS Hydrological modelling has been used globally to assist with decision-making and the implementation of adaptation measures in response to climate change. Streamflow is a major component of the surface water balance, and numerous uncertainties (i.e., input data, model structure, and parameters) remain challenging when applying land surface models for both historical and future projection modelling. These elements have not been commonly integrated in previous studies over a sub-arctic domain across a permafrost gradient, such as in northern Manitoba, Canada (i.e., LNRB). Therefore, this dissertation presents new knowledge on various climate data performances as well as input data selection for hydrological modelling (i.e., the VIC model) studies. Moreover, this study performs comprehensive VIC model parameter sensitivity and uncertainty analyses and projected climate change impacts on overall water balances due to future soil temperature warming and increasing baseflows. Outcomes of this research are anticipated to provide an enhanced understanding in adaptation strategies and decision-making for future water management and regulation developments in the study region. 5.1. Major Findings Findings presented in Chapter 2 of this dissertation (Lilhare et al., 2019) show that annual precipitation from the global reanalysis products ERA-I and WFDEI exceeds that from the ANUSPLIN, IDW, NARR, and HydroGFD datasets over the LNRB during 1981– 2010. The NARR dataset shows greater disagreement, exhibiting a warm bias of ~1.0°C for mean annual air temperature relative to all other datasets. The ERA-I, WFDEI, and HydroGFD (IDW, NARR, and WFDEI) exhibit wet (warm) biases over the central and southern LNRB. There is considerable disparity among all datasets in annual precipitation 133 trends with wetting trends (10.0-30.0 mm decade−1) for all sub-watersheds except ANUSPLIN (1.0-10.0 mm decade-1) and IDW (1.0-5.0 mm decade-1). Mean air temperature trends generally agree across most of the datasets; however, the NARR and IDW show stronger warming relative to other datasets. The Ensemble air temperature does not exhibit any significant warming across the study domain; however, it shows warming of 0.1-0.6°C decade−1 across the LNRB. In this chapter, I conclude that the input dataset selection plays a crucial role in hydro-climatic studies, particularly when assessing hydro-climatic trends for the LNRB. In Chapter 3 (Lilhare et al., 2020) of this dissertation, I evaluate uncertainties propagated through different climate datasets, examined in Chapter 2, where outcomes suggest that the Ensemble-VIC simulations closely match observed streamflow (NSE and KGE ≥ 0.5). In contrast, global reanalysis products such as ERA-I and WFDEI yield high flows (0.5-3.0 mm day-1) with an overestimation of 10-60% in seasonal and annual water balance terms. Overall, the Ensemble-VIC shows a better representation of each water balance term (i.e., total runoff, soil moisture, and ET) followed by the HydroGFD-VIC and NARR-VIC simulations across the study domain. VIC parameters exhibit seasonal importance in VARS analysis, where the second soil layer depth (D2) and variable infiltration curve parameters (binf) emerge as the most sensitive parameters at the annual and seasonal scales. Uncertainty propagated through different climate datasets in the VIC water balance estimation is most evident in the central and southern part of the LNRB due to disagreement between input forcing over these regions. The likelihood distributions of VIC parameters reveal parameter ranges for model calibration in similar hydro-climatic environments, which reduces computational time and modelling efforts. Moreover, the uncertainty envelope in streamflow prediction is presented for all sub- 134 basins, in which input forcing data and model structure are held constant, and the entire uncertainty in streamflow simulation may be attributed to VIC parameters. Chapter 4 of this dissertation investigates the impacts of future warming across a permafrost gradient in the study domain. Results show a 10% decrease and 15-20% increase in the projected annual snowfall and precipitation, respectively over the southern LNRB by the 2050s. Snow accumulation is projected to decline across all sub-basins, particularly in the lower latitudes of the LNRB. Projected uncertainties in major water balance terms (ET, surface runoff, and streamflow) are more substantial in the Gunisao and Grass than other subbasins. The LNRB is projected to experience a warmer and wetter climate, with an evaporative fraction (ET/P) that remains 70% by 2050. The percentage of surface runoff (SR/P) and baseflow (BF/P) fractions increase ~2% in the 2050s. Future warming increases soil temperatures >2.5°C by the 2050s and results in 40-50% more baseflow. Further analyses of soil temperature trends at three different soil depths (D1, D2, and D3) show the most pronounced warming of 1.6°C 30-year-1 in the 2050s in the top soil layer (D1), whereas baseflows increase significantly by 14.0 (26.7) % during the 2030s (2050s). 5.2. Recommendations This research provides a detailed intercomparison of available gridded climate datasets but does not include data homogeneity tests. Moreover, all datasets should be tested over different regions and for different time spans as suggested in both this research as well as previous studies (Choi et al., 2009; Eum et al., 2014; Mesinger et al., 2006). Future work should incorporate other recently developed datasets for intercomparison such as the Canadian Precipitation Analysis (CaPA) (Fortin et al., 2018; Mahfouf et al., 2007), Daymet (daily surface weather data on a 1-km grid for North America, version 3) (Thornton et al., 135 2016), and ERA5-Land (Copernicus Climate Change Service, 2017; Hersbach et al., 2018). The NARR precipitation, WFDEI air temperature, and HydroGFD precipitation and air temperature more closely resemble observations and may perform better in similar topography and hydro-climatic conditions. However, it is recommended that combining available gridded datasets and generating an equally weighted hybrid climate product (i.e., Ensemble) can provide a reliable long-term climate dataset for hydro-climatic studies. Further, assimilation of in situ measurements in gridded climate datasets and advanced interpolation methods may significantly improve the performance of climate data, especially in the LNRB. The VIC parameter sensitivity analyses show D2 and binf as the most sensitive parameters across the LNRB. The sensitivity analyses also recommend parameter ranges, for example, b inf (0.3-0.4), D2 (0.3-0.4 m), D3 (0.6-0.7 m), Ds (0.7-0.8), Dsmax (15-18 mm day-1), and Ws (0.2-0.3), to reduce calibration and modelling efforts in similar hydro-climatic regimes. However, apart from these six soil parameters, it is recommended that the VIC modelling community expands the VARS analysis using other parameters such as for ET and snow. The SA analyses, in this project, show that users should select SA criteria in alignment with their final goals of the modelling application (e.g., flood forecasting, drought analysis, or water balance assessment). Regardless of the metric of choice, often a limited number of parameters controls the majority of model response variations. The SA results provide a list of the most influential VIC parameters over sub-arctic regions to generate reliable modelling results. Therefore, these may reduce calibration efforts and computational burdens in similar hydro-climatic conditions. The VIC users should therefore prioritize influential parameters during model calibration for similar physical and climatic environments. 136 Results obtained in Chapter 4 suggested that projected increases in net precipitation and the future warming of shallow soils across the LNRB may induce more baseflow and overall runoff. It is recommended that these results be tested against additional in situ soil temperature measurements and baseflow observations across the LNRB that would provide more confidence in modelling of soil thermal regimes. There are some known improvements in the VIC parameterization for soil thermodynamics, frozen ground, and permafrost processes, such as developments in soil temperature and thermal nodes computation in addition to improved drainage solutions in updated versions of the VIC model (version 5.0.1) (Andresen et al., 2019; Endalamaw et al., 2017; Gao et al., 2018; Gao et al., 2019; Walvoord & Kurylyk, 2016) that highlight the need for further studies. Future work may focus on quantification of active layer depths, thawing permafrost, and associated soil carbon emissions in north-central Manitoba. In the LNRB, where multiple hydroelectric facilities are situated and ongoing future developments are planned, understanding of future water availability and water balance alterations due to projected changes are essential. While results suggest that increases in surface runoff and streamflow are likely, it will be essential to evaluate the seasonal variability of these changes (Vano, Nijssen, & Lettenmaier, 2015). For instance, there is a need to carefully evaluate excess water availability from streamflow during spring and summer that may be beneficial in later seasons (early and late autumn) for industrial water use and public demand in the LNRB. 5.3. Limitations The availability of gridded datasets over the study domain presents some limitations in this study. For example, recently developed datasets such as Daymet and ERA5-Land that cover the study domain and period, provide additional resources to improve the 137 intercomparison of climate datasets and further VIC modelling analyses. However, dataset selection for this research was based on the BaySys project objectives. This study did not intercompare wind speeds against observations and its sensitivity towards the VIC model response in water balance estimation. The regridding of datasets using bilinear interpolation may also contribute to additional uncertainty since this study did not evaluate these results with other interpolation methods. The sensitivity analyses of VIC parameters are limited to six soil parameters and the VARS method, and do not consider other parameters from lakes, wetlands, land-use, and the routing scheme within the watersheds. Therefore, parameter sensitivity to streamflow generation does not represent possible uncertainty envelopes from additional model parameters. Quantification of climate change impacts on the LNRB hydrology relies on limited GCMs, and do not consider additional climate model outputs in the analyses. This, in turn, leads to gaps in covering potential uncertainty ranges. The results in Chapter 4 provide projected changes in hydrology and linkages between soil temperatures and baseflows across a permafrost gradient in northern Manitoba; however, all selected sub-basins span extensive discontinuous, sporadic discontinuous or isolated patches of permafrost and do not cover any continuous permafrost area. Therefore, findings related to projected changes in soil temperatures and baseflows are limited to discontinuous permafrost areas in a sub-arctic domain. Projected changes in soil thermal regime and hydrology rely on the VIC simulated streamflows, which are calibrated against the observed hydrometric records. However, the VIC simulated soil temperatures are evaluated against observations but the VIC model application does not utilize historic soil temperature observations in model calibration. 138 5.4. Implications of This Work This dissertation presents a comprehensive intercomparison of the multiple climate datasets across the LNRB in northern Manitoba and its utility in the VIC hydrological model to address uncertainty propagated through various input datasets in water balance estimation. The historical intercomparison of annual and seasonal trends is helpful to identify variability across the available datasets, providing a base for climate change impact studies in this region. Moreover, this study extends climate data intercomparisons (e.g., Choi et al., 2009; Eum et al., 2014; Hutchinson et al., 2009; Wong, Razavi, Bonsal, Wheater, & Asong, 2017) to the north-central Canadian watersheds using multiple climate datasets, which have practical and regional importance. Hydrological modelling in this study utilizes a comprehensive sensitivity analysis using the recently developed VARS approach that couples the VIC land surface scheme and extends applicability of VARS to different hydrological models (Razavi & Gupta, 2016b). This research also fulfils the International Association of Hydrological Sciences’ (IAHS’) 23 unsolved problems in hydrology initiatives focused on understanding process changes that control changing runoff response (Blöschl et al., 2019). Findings associated with sensitive parameters and their suggested ranges assist the global hydrological and VIC modelling communities in reducing prediction uncertainty, providing more robust calibrations at lower computational costs in similar hydro-climatic environments. The LNRB, where complex land cover such as lakes, wetlands, permafrost, and various physical characteristics dominate the hydrology, is a prime example of a climate-sensitive basin. This domain has major hydroelectric facilities and thus is of significant environmental and socio-economic importance. An improved understanding of projected changes in the 139 LNRB’s hydrology due to future warming, changes in soil thermal regime, flow timings, and duration have broader implications on industrial water use (Manitoba Hydro) including dams, reservoirs, hydropower generation, community water supply, and use. Our results agree with the previous studies on a general picture of increasing discharge under projected climate scenarios for the LNRB's sub-basins. However, previous studies (Bring & Destouni, 2014; Bring, Shiklomanov, & Lammers, 2017) suggest that the Nelson River basin's upper reaches will dry out in the future; therefore, their contribution to the Nelson River total discharge will diminish over time. Hence, the LNRB may offset decreases in main stem flows through increased future precipitation and runoff. In addition to the earlier efforts, however, our study, by utilizing the information on significant changes but disagreement on the sign, supports Manitoba Hydro on climate change impact assessments where the direction of future discharge changes may be more uncertain, but potentially of considerable magnitude. Moreover, this dissertation and the BaySys project assist Manitoba Hydro on impacts and adaptation strategies for climate change in northern ecosystems as it affects system operations and future hydroelectric developments. The research work carried out by the BaySys project’s Team 2 expands the boundaries of previous studies by focusing on a large, pan-Arctic system experiencing changes in climate (Barber et al., 2012; Bring & Destouni, 2014). This Ph.D. dissertation specifically focuses on the regional importance of the complex hydraulic- and ice-affected regime and its direct linkage and impact on freshwater discharge to Hudson Bay. Results produced by this research serve as input forcing to other teams within the BaySys project looking at biogeochemical and sea ice modelling, ecosystem variability and contaminant transport. The benefit of this research to water resource managers, such as Manitoba Hydro, 140 is to better understand the current and future state of water balance and availability, and associated aquatic environment of the Nelson Watershed and marine environment of Hudson Bay. In addition, this research enhances the understanding of VIC model applicability over the LNRB, various uncertainties, including an investigation of how input forcing impacts climate change assessment and inflow forecasting tools. More broadly, collective outcomes from this research and the BaySys project will benefit Canadian sub-arctic regions, including the Hudson Bay complex, to better understand how seasonal shifts in freshwater eventually affect sediment transport and nutrient delivery in Hudson Bay. Soil temperature warming accompanying projected climate change may enhance the infiltration process, mineralization of carbon stored in the permafrost, and organic matter accumulation in the soil. Thus, future climate change across the Hudson Bay complex may enhance overall streamflow, greenhouse gas emission, and more methane production by the end of 2070s. Further, this project provides new insights on future climate changes in northern Manitoba and the LNRB that contribute to quantifying its impacts on fish productivity and transportation in Hudson Bay. 141 BIBLIOGRAPHY Abebe, N. A., Ogden, F. L., & Pradhan, N. R. (2010). Sensitivity and uncertainty analysis of the conceptual HBV rainfall–runoff model: Implications for parameter estimation. Journal of Hydrology, 389(3–4), 301–310. https://doi.org/10.1016/j.jhydrol.2010.06.007 ACIA. (2004). Impact of a Warming Arctic: Arctic Climate Impact Assessment (p. 139). Cambridge, U.K.; New York, N.Y.: Cambridge University Press. Adam, J. C., Hamlet, A. F., & Lettenmaier, D. P. (2009). Implications of global climate change for snowmelt hydrology in the twenty-first century. Hydrological Processes, 23(7), 962–972. https://doi.org/10.1002/hyp.7201 Adler, R. F., Kidd, C., Petty, G., Morissey, M., & Goodman, H. M. (2001). Intercomparison of global precipitation products: The third precipitation intercomparison project (PIP– 3). Bulletin of the American Meteorological Society, 82(7), 1377–1396. https://doi.org/10.1175/1520-0477(2001)082<1377:IOGPPT>2.3.CO;2 Anderson, J., Chung, F., Anderson, M., Brekke, L., Easton, D., Ejeta, M., … Snyder, R. (2008). Progress on incorporating climate change into management of California’s water resources. Climatic Change, 87, 91–108. https://doi.org/10.1007/s10584-0079353-1 Andresen, C. G., Lawrence, D. M., Wilson, C. J., McGuire, A. D., Koven, C., Schaefer, K., … Zhang, W. (2020). Soil moisture and hydrology projections of the permafrost region: A model intercomparison. The Cryosphere, 14(2), 445–459. https://doi.org/10.5194/tc-14-445-2020, 2020 Anisimov, O. A., Fitzharris, B., Hagan, J. O., Jeffries, R., Marchant, H., Nelson, F. E., … Vaughan, D. G. (2001). Polar Regions (Arctic and Antarctic). In: Climate change 2001: Impacts, adaptation, and vulnerability: Contribution of Working Group II to the third assessment report of the Intergovernmental Panel on Climate Change (pp. 801–841). Cambridge, U.K.; New York: Cambridge University Press. Arora, V. K., & Boer, G. J. (2001). Effects of simulated climate change on the hydrology of major river basins. Journal of Geophysical Research: Atmospheres (1984–2012), 106(D4), 3335–3348. https://doi.org/10.1029/2000JD900620 Arsenault, R., & Brissette, F. P. (2014). Continuous streamflow prediction in ungauged basins: The effects of equifinality and parameter set selection on uncertainty in regionalization approaches. Water Resources Research, 50(7), 6135–6153. https://doi.org/10.1002/2013WR014898 Aster, R. C., Borchers, B., & Thurber, C. H. (2013). Parameter estimation and inverse problems (2nd ed.). Academic Press. Aygün, O., Kinnard, C., & Campeau, S. (2019). Impacts of climate change on the hydrology of northern midlatitude cold regions. Progress in Physical Geography: Earth and Environment. (Sage UK: London, England). https://doi.org/10.1177/0309133319878123 142 Barber, D. G., Asplin, M. G., Papakyriakou, T. N., Miller, L., Else, B. G. T., Iacozza, J., … Wang, F. (2012). Consequences of change and variability in sea ice on marine ecosystem and biogeochemical processes during the 2007–2008 Canadian International Polar Year program. Climatic Change, 115(1), 135–159. https://doi.org/10.1007/s10584-012-0482-9 Barnett, T. P., Adam, J. C., & Lettenmaier, D. P. (2005). Potential impacts of a warming climate on water availability in snow-dominated regions. Nature, 438(7066), 303– 309. https://doi.org/10.1038/nature04141 Beltrami, H., Gosselin, C., & Mareschal, J. C. (2003). Ground surface temperatures in Canada: Spatial and temporal variability. Geophysical Research Letters, 30(10), 1499. https://doi.org/10.1029/2003GL017144 Bennett, K. E., Blanco, J. R. U., Jonko, A., Bohn, T. J., Atchley, A. L., Urban, N. M., & Middleton, R. S. (2018). Global sensitivity of simulated water balance indicators under future climate change in the Colorado Basin. Water Resources Research, 54(1), 132–149. https://doi.org/10.1002/2017WR020471 Bennett, K. E., Werner, A. T., & Schnorbus, M. (2012). Uncertainties in hydrologic and climate change impact analyses in headwater basins of British Columbia. Journal of Climate, 25(17), 5711–5730. https://doi.org/10.1175/JCLI-D-11-00417.1 Berg, P., Donnelly, C., & Gustafsson, D. (2018). Near-real-time adjusted reanalysis forcing data for hydrology. Hydrology and Earth System Sciences, 22(2), 989–1000. https://doi.org/10.5194/hess-22-989-2018 Berghuijs, W. R., Woods, R. A., & Hrachowitz, M. (2014). A precipitation shift from snow towards rain leads to a decrease in streamflow. Nature Climate Change, 4(7), 583– 586. https://doi.org/10.1038/nclimate2246 Bermúdez, M., Neal, J. C., Bates, P. D., Coxon, G., Freer, J. E., Cea, L., & Puertas, J. (2017). Quantifying local rainfall dynamics and uncertain boundary conditions into a nested regional-local flood modeling system. Water Resources Research, 53(4), 2770–2785. https://doi.org/10.1002/2016WR019903 Betts, A. K., Ball, J. H., & Viterbo, P. (2003). Evaluation of the ERA-40 surface water budget and surface temperature for the Mackenzie River basin. Journal of Hydrometeorology, 4(6), 1194–1211. https://doi.org/10.1175/15257541(2003)004<1194:EOTESW>2.0.CO;2 Betts, A. K., & Beljaars, A. C. M. (2017). Analysis of near-surface biases in ERA-Interim over the Canadian Prairies. Journal of Advances in Modeling Earth Systems, 9(5), 2158–2173. https://doi.org/10.1002/2017MS001025 Beven, K. (2006). A manifesto for the equifinality thesis. Journal of Hydrology, 320(1), 18– 36. https://doi.org/10.1016/j.jhydrol.2005.07.007 Beven, K., & Binley, A. (1992). The future of distributed models: Model calibration and uncertainty prediction. Hydrological Processes, 6(3), 279–298. https://doi.org/10.1002/hyp.3360060305 Biemans, H., Hutjes, R. W. A., Kabat, P., Strengers, B. J., Gerten, D., & Rost, S. (2009). Effects of precipitation uncertainty on discharge calculations for main river basins. 143 Journal of Hydrometeorology, 10(4), 1011–1025. https://doi.org/10.1175/2008JHM1067.1 Blenkinsop, S., & Fowler, H. J. (2007). Changes in European drought characteristics projected by the PRUDENCE regional climate models. International Journal of Climatology, 27(12), 1595–1610. https://doi.org/10.1002/joc.1538 Blöschl, G., Bierkens, M. F. P., Chambel, A., Cudennec, C., Destouni, G., Fiori, A., … Zhang, Y. (2019). Twenty-three unsolved problems in hydrology (UPH) – a community perspective. Hydrological Sciences Journal, 64(10), 1141–1158. https://doi.org/10.1080/02626667.2019.1620507 Blöschl, G., & Zehe, E. (2005). On hydrological predictability. Hydrological Processes, 19(19), 3923–3929. https://doi.org/10.1002/hyp.6075 Boden, T., Kaiser, D. P., Sepanski, R. J., & Stoss, F. W. (1994). Trends’ 93: A compendium of data on global change. eds. (p. 984). Oak Ridge National Laboratory, Oak Ridge, Tennessee: ORNL/CDIAC-65, Carbon Dioxide Information Analysis Center. Bosilovich, M. G., Chen, J., Robertson, F. R., & Adler, R. F. (2008). Evaluation of global precipitation in reanalyses. Journal of Applied Meteorology and Climatology, 47(9), 2279–2299. https://doi.org/10.1175/2008JAMC1921.1 Boucher, O., & Best, M. (2010). The WATCH forcing data 1958-2001: A meteorological forcing dataset for land surface-and hydrological-models. (No. 22; p. 41). WATCH technical report. Retrieved from WATCH technical report website: http://www.euwatch.org/media/default.aspx/emma/org/10376311/WATCH+Technical+Report+Nu mber+22+The+WATCH+forcing+data+19582001+A+meteorological+forcing+dataset+for+land+surface-+and+hydrologicalmodels.pdf Bowling, L. C., & Lettenmaier, D. P. (2010). Modeling the effects of lakes and wetlands on the water balance of arctic environments. Journal of Hydrometeorology, 11(2), 276– 295. https://doi.org/10.1175/2009JHM1084.1 Bowling, L. C., Lettenmaier, D. P., Nijssen, B., Graham, L. P., Clark, D. B., El Maayar, M., … others. (2003). Simulation of high-latitude hydrological processes in the Torne– Kalix basin: PILPS Phase 2 (e): 1: Experiment description and summary intercomparisons. Global and Planetary Change, 38(1), 1–30. Bring, A., & Destouni, G. (2014). Arctic Climate and Water Change: Model and Observation Relevance for Assessment and Adaptation. Surveys in Geophysics, 35(3), 853–877. https://doi.org/10.1007/s10712-013-9267-6 Bring, A., Shiklomanov, A., & Lammers, R. B. (2017). Pan-Arctic river discharge: Prioritizing monitoring of future climate change hot spots. Earth’s Future, 5(1), 72– 92. https://doi.org/10.1002/2016EF000434 Bromwich, D. H., Wilson, A. B., Bai, L.-S., Moore, G. W. K., & Bauer, P. (2016). A comparison of the regional Arctic System Reanalysis and the global ERA-Interim Reanalysis for the Arctic. Quarterly Journal of the Royal Meteorological Society, 142(695), 644–658. https://doi.org/10.1002/qj.2527 144 Bukovsky, M. S., & Karoly, D. J. (2007). A brief evaluation of precipitation from the North American Regional Reanalysis. Journal of Hydrometeorology, 8(4), 837–846. https://doi.org/10.1175/JHM595.1 Burford, J. E., Déry, S. J., & Holmes, R. D. (2009). Some aspects of the hydroclimatology of the Quesnel River Basin, British Columbia, Canada. Hydrological Processes, 23(10), 1529–1536. https://doi.org/10.1002/hyp.7253 Centre for Land and Biological Resources Research. (1996). Soil Landscapes of Canada, v.2.2, Research Branch, Agriculture and Agri-Food Canada. Ottawa. Retrieved 23 December 2016, from http://sis.agr.gc.ca/cansis/nsdb/slc/v2.2/index.html Chadburn, S. E., Burke, E. J., Cox, P. M., Friedlingstein, P., Hugelius, G., & Westermann, S. (2017). An observation-based constraint on permafrost loss as a function of global warming. Nature Climate Change, 7(5), 340–344. https://doi.org/10.1038/nclimate3262 Chapin, F. I., Jefferies, R., Reynolds, J., Shaver, G., Svoboda, J., & Chu, E. (1992). Arctic Ecosystems in a Changing Climate—An ecophysiological perspective (1st ed.). San Diego, CA: Academic Press Inc. Retrieved from https://doi.org/10.1016/C2009-002634-8 Chen, J., Brissette, F. P., & Leconte, R. (2011). Uncertainty of downscaling method in quantifying the impact of climate change on hydrology. Journal of Hydrology, 401(3), 190–202. https://doi.org/10.1016/j.jhydrol.2011.02.020 Chen, J., Brissette, F. P., Poulin, A., & Leconte, R. (2011). Overall uncertainty study of the hydrological impacts of climate change for a Canadian watershed. Water Resources Research, 47, W12509. https://doi.org/10.1029/2011WR010602 Cherkauer, K. A., Bowling, L. C., & Lettenmaier, D. P. (2003). Variable infiltration capacity cold land process model updates. Global and Planetary Change, 38(1), 151–159. https://doi.org/10.1016/S0921-8181(03)00025-0 Cherkauer, K. A., & Lettenmaier, D. P. (1999). Hydrologic effects of frozen soils in the upper Mississippi River basin. Journal of Geophysical Research: Atmospheres, 104(D16), 19599–19610. https://doi.org/10.1029/1999JD900337 Cherkauer, K. A., & Lettenmaier, D. P. (2003). Simulation of spatial variability in snow and frozen soil. Journal of Geophysical Research: Atmospheres, 108(D22), 8858. https://doi.org/10.1029/2003JD003575 Choi, W., Kim, S. J., Rasmussen, P. F., & Moore, A. R. (2009). Use of the North American Regional Reanalysis for hydrological modelling in Manitoba. Canadian Water Resources Journal, 34(1), 17–36. https://doi.org/10.4296/cwrj3401017 Christensen, N. S., & Lettenmaier, D. P. (2007). A multimodel ensemble approach to assessment of climate change impacts on the hydrology and water resources of the Colorado River Basin. Hydrology and Earth System Sciences, 11(4), 1417–1434. https://doi.org/10.5194/hess-11-1417-2007 Cook, B. I., Smerdon, J. E., Seager, R., & Coats, S. (2014). Global warming and 21st century drying. Climate Dynamics, 43(9), 2607–2627. https://doi.org/10.1007/s00382-0142075-y 145 Copernicus Climate Change Service. (2017). Retrieved 5 February 2020, from ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service Climate Data Store (CDS) website: https://cds.climate.copernicus.eu/cdsapp#!/home Cosby, B. J., Hornberger, G. M., Clapp, R. B., & Ginn, T. (1984). A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils. Water Resources Research, 20(6), 682–690. https://doi.org/10.1029/WR020i006p00682 Costa, M. H., & Foley, J. A. (1998). A comparison of precipitation datasets for the Amazon basin. Geophysical Research Letters, 25(2), 155–158. https://doi.org/10.1029/97GL03502 Cuo, L., Beyene, T. K., Voisin, N., Su, F., Lettenmaier, D. P., Alberti, M., & Richey, J. E. (2011). Effects of mid-twenty-first century climate and land cover change on the hydrology of the Puget Sound basin, Washington. Hydrological Processes, 25(11), 1729–1753. https://doi.org/10.1002/hyp.7932 Cuo, L., Lettenmaier, D. P., Alberti, M., & Richey, J. E. (2009). Effects of a century of land cover and climate change on the hydrology of the Puget Sound basin. Hydrological Processes, 23(6), 907–933. https://doi.org/10.1002/hyp.7228 Dahri, Z. H., Moors, E., Ludwig, F., Ahmad, S., Khan, A., Ali, I., & Kabat, P. (2018). Adjustment of measurement errors to reconcile precipitation distribution in the highaltitude Indus basin. International Journal of Climatology, 38(10), 3842–3860. https://doi.org/10.1002/joc.5539 de Elía, R., Caya, D., Côté, H., Frigon, A., Biner, S., Giguère, M., … Plummer, D. (2008). Evaluation of uncertainties in the CRCM-simulated North American climate. Climate Dynamics, 30(2–3), 113–132. https://doi.org/10.1007/s00382-007-0288-z Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., … Vitart, F. (2011). The ERA-Interim reanalysis: Configuration and performance of the data assimilation system. Quarterly Journal of the Royal Meteorological Society, 137(656), 553–597. https://doi.org/10.1002/qj.828 Demaria, E. M., Nijssen, B., & Wagener, T. (2007). Monte Carlo sensitivity analysis of land surface parameters using the Variable Infiltration Capacity model. Journal of Geophysical Research: Atmospheres, 112, D11113. https://doi.org/10.1029/2006JD007534 Déry, S. J., & Brown, R. D. (2007). Recent Northern Hemisphere snow cover extent trends and implications for the snow-albedo feedback. Geophysical Research Letters, 34(22), L22504. https://doi.org/10.1029/2007GL031474 Déry, S. J., Hernández‐Henríquez, M. A., Burford, J. E., & Wood, E. F. (2009). Observational evidence of an intensifying hydrological cycle in northern Canada. Geophysical Research Letters, 36, L13402. https://doi.org/10.1029/2009GL038852 Déry, S. J., Hernández-Henríquez, M. A., Owens, P. N., Parkes, M. W., & Petticrew, E. L. (2012). A century of hydrological variability and trends in the Fraser River Basin. 146 Environmental Research Letters, 7(2), 024019. https://doi.org/10.1088/17489326/7/2/024019 Déry, S. J., Mlynowski, T. J., Hernández-Henríquez, M. A., & Straneo, F. (2011). Interannual variability and interdecadal trends in Hudson Bay streamflow. Journal of Marine Systems, 88(3), 341–351. https://doi.org/10.1016/j.jmarsys.2010.12.002 Déry, S. J., Stadnyk, T. A., MacDonald, M. K., Koenig, K. A., & Guay, C. (2018). Flow alteration impacts on Hudson Bay river discharge. Hydrological Processes, 32, 3576– 3587. https://doi.org/10.1002/hyp.13285 Déry, S. J., Stieglitz, M., McKenna, E. C., & Wood, E. F. (2005). Characteristics and trends of river discharge into Hudson, James, and Ungava Bays, 1964-2000. Journal of Climate, 18(14), 2540–2557. https://doi.org/10.1175/JCLI3440.1 Déry, S. J., & Wood, E. F. (2005a). Decreasing river discharge in northern Canada. Geophysical Research Letters, 32(10), L10401. https://doi.org/10.1029/2005GL022845 Déry, S. J., & Wood, E. F. (2005b). Observed twentieth century land surface air temperature and precipitation covariability. Geophysical Research Letters, 32, L21414. https://doi.org/10.1029/2005GL024234 Döll, P., Kaspar, F., & Lehner, B. (2003). A global hydrological model for deriving water availability indicators: Model tuning and validation. Journal of Hydrology, 270(1), 105–134. https://doi.org/10.1016/S0022-1694(02)00283-4 Dore, M. H. (2005). Climate change and changes in global precipitation patterns: What do we know? Environment International, 31(8), 1167–1181. https://doi.org/10.1016/j.envint.2005.03.004 Duan, Q., Sorooshian, S., & Gupta, V. (1992). Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resources Research, 28(4), 1015–1031. https://doi.org/10.1029/91WR02985 Eisner, S., Flörke, M., Chamorro, A., Daggupati, P., Donnelly, C., Huang, J., … Krysanova, V. (2017). An ensemble analysis of climate change impacts on streamflow seasonality across 11 large river basins. Climatic Change, 141(3), 401–417. https://doi.org/10.1007/s10584-016-1844-5 Elsner, M. M., Cuo, L., Voisin, N., Deems, J. S., Hamlet, A. F., Vano, J. A., … Lettenmaier, D. P. (2010). Implications of 21st century climate change for the hydrology of Washington State. Climatic Change, 102(1–2), 225–260. https://doi.org/10.1007/s10584-010-9855-0 Endalamaw, A., Bolton, W. R., Young-Robertson, J. M., Morton, D., Hinzman, L., & Nijssen, B. (2017). Towards improved parameterization of a macroscale hydrologic model in a discontinuous permafrost boreal forest ecosystem. Hydrology and Earth System Sciences, 21(9), 4663–4680. https://doi.org/10.5194/hess-21-4663-2017 Environment and Climate Change Canada. (2016). Canadian Climate Normals 1981-2010 Station Data—Climate—Environment Canada. Retrieved 23 December 2016, from http://climate.weather.gc.ca/climate_normals/results_1981_2010_e.html 147 Environment Canada. (1995). The states of Canada’s climate: Monitoring variability and change (No. 95–1; p. 52). Minister of Public Works and Government Services Canada. Essou, G. R., Sabarly, F., Lucas-Picher, P., Brissette, F., & Poulin, A. (2016). Can precipitation and temperature from meteorological reanalyses be used for hydrological modeling? Journal of Hydrometeorology, 17(7), 1929–1950. https://doi.org/10.1175/JHM-D-15-0138.1 Eum, H.-I., Dibike, Y., Prowse, T., & Bonsal, B. (2014). Inter-comparison of high-resolution gridded climate data sets and their implication on hydrological model simulation over the Athabasca Watershed, Canada. Hydrological Processes, 28(14), 4250–4271. https://doi.org/10.1002/hyp.10236 Fatichi, S., Vivoni, E. R., Ogden, F. L., Ivanov, V. Y., Mirus, B., Gochis, D., … Tarboton, D. (2016). An overview of current applications, challenges, and future trends in distributed process-based models in hydrology. Journal of Hydrology, 537, 45–60. https://doi.org/10.1016/j.jhydrol.2016.03.026 Fekete, B. M., Vörösmarty, C. J., Roads, J. O., & Willmott, C. J. (2004). Uncertainties in precipitation and their impacts on runoff estimates. Journal of Climate, 17(2), 294– 304. https://doi.org/10.1175/1520-0442(2004)017<0294:UIPATI>2.0.CO;2 Fekete, B. M., Wisser, D., Kroeze, C., Mayorga, E., Bouwman, L., Wollheim, W. M., & Vörösmarty, C. (2010). Millennium ecosystem assessment scenario drivers (1970– 2050): Climate and hydrological alterations. Global Biogeochemical Cycles, 24(4), GB0A12. https://doi.org/10.1029/2009GB003593 Ficklin, D. L., Stewart, I. T., & Maurer, E. P. (2013). Effects of climate change on stream temperature, dissolved oxygen, and sediment concentration in the Sierra Nevada in California. Water Resources Research, 49(5), 2765–2782. https://doi.org/10.1002/wrcr.20248 Fortin, V., Roy, G., Stadnyk, T., Koenig, K., Gasset, N., & Mahidjiba, A. (2018). Ten years of science based on the Canadian Precipitation Analysis: A CaPA system overview and literature review. Atmosphere-Ocean, 56(3), 178–196. https://doi.org/10.1080/07055900.2018.1474728 Fowler, H. J., & Archer, D. R. (2006). Conflicting signals of climatic change in the Upper Indus Basin. Journal of Climate, 19(17), 4276–4293. https://doi.org/10.1175/JCLI3860.1 Fowler, H. J., Ekström, M., Blenkinsop, S., & Smith, A. P. (2007). Estimating change in extreme European precipitation using a multimodel ensemble. Journal of Geophysical Research: Atmospheres, 112, D18104. https://doi.org/10.1029/2007JD008619 Fowler, H. J., & Kilsby, C. G. (2007). Using regional climate model data to simulate historical and future river flows in northwest England. Climatic Change, 80(3–4), 337–367. https://doi.org/10.1007/s10584-006-9117-3 Frey, H. C., & Patil, S. R. (2002). Identification and review of sensitivity analysis methods. Risk Analysis, 22(3), 553–578. https://doi.org/10.1111/0272-4332.00039 148 Gan, T. Y., & Kwong, Y. T. (1992). Identification of warming trends in northern Alberta and southern Northwest Territories by the non-parametric Kendall’s test. G.W. Kite & K.D. Harvey (Eds), Using Hydrometric Data to Detect and Monitor Climatic Change. Proceedings of NHRI Workshop No. 8, 43–56. National Hydrology Research Institute, Saskatoon, SK, Canada. Gan, Y., Duan, Q., Gong, W., Tong, C., Sun, Y., Chu, W., … Di, Z. (2014). A comprehensive evaluation of various sensitivity analysis methods: A case study with a hydrological model. Environmental Modelling & Software, 51, 269–285. https://doi.org/10.1016/j.envsoft.2013.09.031 Gao, B., Yang, D., Qin, Y., Wang, Y., Li, H., Zhang, Y., & Zhang, T. (2018). Change in frozen soils and its effect on regional hydrology, upper Heihe basin, northeastern Qinghai–Tibetan Plateau. The Cryosphere, 12(2), 657–673. https://doi.org/10.5194/tc-12-657-2018 Gao, H., Tang, Q., Shi, X., Zhu, C., Bohn, T. J., Su, F., … Wood, E. F. (2010). Water budget record from Variable Infiltration Capacity (VIC) model. Algorithm Theoretical Basis Document for Terrestrial Water Cycle Data Records. Retrieved from https://www.researchgate.net/profile/Xiaogang_Shi/publication/268367169_Water_B udget_Record_from_Variable_Infiltration_Capacity_VIC_Model/links/55715dee08ae e701d61cc286.pdf Gao, J., Xie, Z., Wang, A., Liu, S., Zeng, Y., Liu, B., … Xie, J. (2019). A new frozen soil parameterization including frost and thaw fronts in the community land model. Journal of Advances in Modeling Earth Systems, 11(3), 659–679. https://doi.org/10.1029/2018MS001399 Gebremichael, M., Krajewski, W. F., Morrissey, M. L., Huffman, G. J., & Adler, R. F. (2005). A detailed evaluation of GPCP 1 daily rainfall estimates over the Mississippi River Basin. Journal of Applied Meteorology, 44(5), 665–681. https://doi.org/10.1175/JAM2233.1 Gelfan, A., Gustafsson, D., Motovilov, Y., Arheimer, B., Kalugin, A., Krylenko, I., & Lavrenov, A. (2017). Climate change impact on the water regime of two great Arctic rivers: Modeling and uncertainty issues. Climatic Change, 141(3), 499–515. https://doi.org/10.1007/s10584-016-1710-5 Gemmer, M., Becker, S., & Jiang, T. (2004). Observed monthly precipitation trends in China 1951–2002. Theoretical and Applied Climatology, 77(1–2), 39–45. https://doi.org/10.1007/s00704-003-0018-3 Gerten, D., Schaphoff, S., Haberlandt, U., Lucht, W., & Sitch, S. (2004). Terrestrial vegetation and water balance—Hydrological evaluation of a dynamic global vegetation model. Journal of Hydrology, 286(1), 249–270. https://doi.org/10.1016/j.jhydrol.2003.09.029 Gocic, M., & Trajkovic, S. (2013). Analysis of changes in meteorological variables using Mann-Kendall and Sen’s slope estimator statistical tests in Serbia. Global and Planetary Change, 100, 172–182. https://doi.org/10.1016/j.gloplacha.2012.10.014 Goldberg, D. E. (1989). Genetic Algorithms in Search, Optimization, and Machine Learning. Addison Wesley Publishing Company, Reading, MA. 149 Goulden, M. L., Wofsy, S. C., Harden, J. W., Trumbore, S. E., Crill, P. M., Gower, S. T., … Sutton, D. J. (1998). Sensitivity of boreal forest carbon balance to soil thaw. Science, 279(5348), 214–217. https://doi.org/10.1126/science.279.5348.214 Guo, D., & Wang, H. (2016). CMIP5 permafrost degradation projection: A comparison among different regions. Journal of Geophysical Research: Atmospheres, 121(9), 4499–4517. https://doi.org/10.1002/2015JD024108 Gupta, H. V., Kling, H., Yilmaz, K. K., & Martinez, G. F. (2009). Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. Journal of Hydrology, 377(1–2), 80–91. https://doi.org/10.1016/j.jhydrol.2009.08.003 Haddeland, I., Heinke, J., Biemans, H., Eisner, S., Flörke, M., Hanasaki, N., … Wisser, D. (2014). Global water resources affected by human interventions and climate change. Proceedings of the National Academy of Sciences, 111(9), 3251–3256. https://doi.org/10.1073/pnas.1222475110 Haddeland, I., Skaugen, T., & Lettenmaier, D. P. (2006). Anthropogenic impacts on continental surface water fluxes. Geophysical Research Letters, 33(8), L08406. https://doi.org/10.1029/2006GL026047 Hamman, J. J., Nijssen, B., Bohn, T. J., Gergel, D. R., & Mao, Y. (2018). The Variable Infiltration Capacity model version 5 (VIC-5): Infrastructure improvements for new applications and reproducibility. Geoscientific Model Development, 11(8), 3481– 3496. https://doi.org/10.5194/gmd-11-3481-2018 Held, I. M., & Soden, B. J. (2006). Robust responses of the hydrological cycle to global warming. Journal of Climate, 19(21), 5686–5699. https://doi.org/10.1175/JCLI3990.1 Hersbach, H., de Rosnay, P., Bell, B., Schepers, D., Simmons, A., Soci, C., … Zuo, H. (2018). Operational global reanalysis: Progress, future directions and synergies with NWP (No. 27). Retrieved from https://www.ecmwf.int/node/18765 Hidalgo, H. G., Das, T., Dettinger, M. D., Cayan, D. R., Pierce, D. W., Barnett, T. P., … Nozawa, T. (2009). Detection and attribution of streamflow timing changes to climate change in the western United States. Journal of Climate, 22(13), 3838–3855. https://doi.org/10.1175/2009JCLI2470.1 Hill, M. C., & Tiedeman, C. R. (2007). Effective groundwater model calibration: With analysis of data, sensitivities, predictions, and uncertainty. Hoboken, NJ: John Wiley & Sons. Retrieved from https://onlinelibrary.wiley.com/doi/book/10.1002/0470041080 Hinzman, L. D., Bettez, N. D., Bolton, W. R., Chapin, F. S., Dyurgerov, M. B., Fastie, C. L., … Huntington, H. P. (2005). Evidence and implications of recent climate change in northern Alaska and other arctic regions. Climatic Change, 72(3), 251–298. https://doi.org/10.1007/s10584-005-5352-2 Hively, W. D., Gérard-Marchant, P., & Steenhuis, T. S. (2006). Distributed hydrological modeling of total dissolved phosphorus transport in an agricultural landscape, part II: Dissolved phosphorus transport. Hydrology and Earth System Sciences, 10(2), 263– 276. https://doi.org/10.5194/hess-10-263-2006 150 Hong, Y., Adler, R. F., Huffman, G. J., & Pierce, H. (2010). Applications of TRMM-based multi-satellite precipitation estimation for global runoff prediction: Prototyping a global flood modeling system. In Satellite Rainfall Applications for Surface Hydrology (pp. 245–265). Springer. Retrieved from https://link.springer.com/chapter/10.1007/978-90-481-2915-7_15 Hopkinson, R. F., McKenney, D. W., Milewska, E. J., Hutchinson, M. F., Papadopol, P., & Vincent, L. A. (2011). Impact of aligning climatological day on gridding daily maximum–minimum temperature and precipitation over Canada. Journal of Applied Meteorology and Climatology, 50(8), 1654–1665. https://doi.org/10.1175/2011JAMC2684.1 Horton, P., Schaefli, B., Mezghani, A., Hingray, B., & Musy, A. (2006). Assessment of climate-change impacts on alpine discharge regimes with climate model uncertainty. Hydrological Processes, 20(10), 2091–2109. https://doi.org/10.1002/hyp.6197 Hostetler, S. W. (1991). Simulation of lake ice and its effect on the late-Pleistocene evaporation rate of Lake Lahontan. Climate Dynamics, 6(1), 43–48. https://doi.org/10.1007/BF00210581 Hostetler, S. W., & Bartlein, P. J. (1990). Simulation of lake evaporation with application to modeling lake level variations of Harney-Malheur Lake, Oregon. Water Resources Research, 26(10), 2603–2612. https://doi.org/10.1029/WR026i010p02603 Huisman, J. A., Breuer, L., Bormann, H., Bronstert, A., Croke, B. F. W., Frede, H. G., … Kite, G. (2009). Assessing the impact of land use change on hydrology by ensemble modeling (LUCHEM) III: Scenario analysis. Advances in Water Resources, 32(2), 159–170. https://doi.org/10.1016/j.advwatres.2008.06.009 Hundecha, Y., St-Hilaire, A., Ouarda, T., El Adlouni, S., & Gachon, P. (2008). A nonstationary extreme value analysis for the assessment of changes in extreme annual wind speed over the Gulf of St. Lawrence, Canada. Journal of Applied Meteorology and Climatology, 47(11), 2745–2759. https://doi.org/10.1175/2008JAMC1665.1 Hutchinson, M. F., McKenney, D. W., Lawrence, K., Pedlar, J. H., Hopkinson, R. F., Milewska, E., & Papadopol, P. (2009). Development and testing of Canada-wide interpolated spatial models of daily minimum–maximum temperature and precipitation for 1961–2003. Journal of Applied Meteorology and Climatology, 48(4), 725–741. https://doi.org/10.1175/2008JAMC1979.1 Iman, R. L., & Helton, J. C. (1988). An investigation of uncertainty and sensitivity analysis techniques for computer models. Risk Analysis, 8(1), 71–90. https://doi.org/10.1111/j.1539-6924.1988.tb01155.x Immerzeel, W. W., Van Beek, L. P. H., Konz, M., Shrestha, A. B., & Bierkens, M. F. P. (2012). Hydrological response to climate change in a glacierized catchment in the Himalayas. Climatic Change, 110(3–4), 721–736. https://doi.org/10.1007/s10584011-0143-4 Islam, S. U., & Déry, S. J. (2017). Evaluating uncertainties in modelling the snow hydrology of the Fraser River Basin, British Columbia, Canada. Hydrology and Earth System Sciences, 21(3), 1827–1847. https://doi.org/10.5194/hess-21-1827-2017 151 Islam, S. U., Déry, S. J., & Werner, A. T. (2017). Future climate change impacts on snow and water resources of the Fraser River Basin, British Columbia. Journal of Hydrometeorology, 18(2), 473–496. https://doi.org/10.1175/JHM-D-16-0012.1 Jeong, D. I., Sushama, L., & Naveed Khaliq, M. (2014). The role of temperature in drought projections over North America. Climatic Change, 127(2), 289–303. https://doi.org/10.1007/s10584-014-1248-3 Jiang, T., Chen, Y. D., Xu, C., Chen, X., Chen, X., & Singh, V. P. (2007). Comparison of hydrological impacts of climate change simulated by six hydrological models in the Dongjiang Basin, South China. Journal of Hydrology, 336(3), 316–333. https://doi.org/10.1016/j.jhydrol.2007.01.010 Jorgenson, M. T., Racine, C. H., Walters, J. C., & Osterkamp, T. E. (2001). Permafrost degradation and ecological changes associated with a warming climate in central Alaska. Climatic Change, 48(4), 551–579. https://doi.org/10.1023/A:1005667424292 Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., … Woollen, J. (1996). The NCEP/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society, 77(3), 437–471. https://doi.org/10.1175/15200477(1996)077<0437:TNYRP>2.0.CO;2 Kang, D. H., Gao, H., Shi, X., Islam, S. U., & Déry, S. J. (2016). Impacts of a rapidly declining mountain snowpack on streamflow timing in Canada’s Fraser River Basin. Scientific Reports, 6, 19299. https://doi.org/10.1038/srep19299 Kavetski, D., Kuczera, G., & Franks, S. W. (2006). Calibration of conceptual hydrological models revisited: 2. Improving optimisation and analysis. Journal of Hydrology, 320(1), 187–201. https://doi.org/10.1016/j.jhydrol.2005.07.013 Kay, A. L., Davies, H. N., Bell, V. A., & Jones, R. G. (2009). Comparison of uncertainty sources for climate change impacts: Flood frequency in England. Climatic Change, 92(1–2), 41–63. https://doi.org/10.1007/s10584-008-9471-4 Kendall, M. G. (1975). Rank correlation methods (4th ed.) (Vol. 8). London: Charles Griffin. Kettles, I. M., Tarnocai, C., & Bauke, S. D. (1997). Predicted permafrost distribution in Canada under a climate warming scenario. Current Research, Geological Survey of Canada, 1997-E, 109–115. Kim, K.-Y., Kim, J., Boo, K.-O., Shim, S., & Kim, Y. (2019). Intercomparison of precipitation datasets for summer precipitation characteristics over East Asia. Climate Dynamics, 52(5), 3005–3022. https://doi.org/10.1007/s00382-018-4303-3 Kistler, R., Collins, W., Saha, S., White, G., Woollen, J., Kalnay, E., … Kousky, V. (2001). The NCEP–NCAR 50–year reanalysis: Monthly means CD–ROM and documentation. Bulletin of the American Meteorological Society, 82(2), 247–267. https://doi.org/10.1175/1520-0477(2001)082<0247:TNNYRM>2.3.CO;2 Kochendorfer, J., Rasmussen, R., Wolff, M., Baker, B., Hall, M. E., Meyers, T., … Leeper, R. (2017). The quantification and correction of wind-induced precipitation measurement errors. Hydrology and Earth System Sciences, 21(4), 1973–1989. https://doi.org/10.5194/hess-21-1973-2017 152 Koutsoyiannis, D. (2010). HESS Opinions ‘A random walk on water’. Hydrology and Earth System Sciences, 14(3), 585–601. https://doi.org/10.5194/hess-14-585-2010, 2010 Kouwen, N. (1988). WATFLOOD: A micro-computer based flood forecasting system based on real-time weather radar. Canadian Water Resources Journal, 13(1), 62–77. https://doi.org/10.4296/cwrj1301062 Koven, C. D., Riley, W. J., & Stern, A. (2013). Analysis of permafrost thermal dynamics and response to climate change in the CMIP5 Earth System Models. Journal of Climate, 26(6), 1877–1900. https://doi.org/10.1175/JCLI-D-12-00228.1 Krishnamurti, T. N., Kishtawal, C. M., LaRow, T. E., Bachiochi, D. R., Zhang, Z., Williford, C. E., … Surendran, S. (1999). Improved weather and seasonal climate forecasts from multimodel superensemble. Science, 285(5433), 1548–1550. https://doi.org/10.1126/science.285.5433.1548 Kumar, A., Singh, R., Jena, P. P., Chatterjee, C., & Mishra, A. (2015). Identification of the best multi-model combination for simulating river discharge. Journal of Hydrology, 525, 313–325. https://doi.org/10.1016/j.jhydrol.2015.03.060 Langlois, A., Kohn, J., Royer, A., Cliche, P., Brucker, L., Picard, G., … Willemet, J. M. (2009). Simulation of Snow Water Equivalent (SWE) using thermodynamic snow models in Québec, Canada. Journal of Hydrometeorology, 10(6), 1447–1463. https://doi.org/10.1175/2009JHM1154.1 Leeper, R. D., & Kochendorfer, J. (2015). Evaporation from weighing precipitation gauges: Impacts on automated gauge measurements and quality assurance methods. Atmospheric Measurement Techniques, 8(6), 2291–2300. https://doi.org/10.5194/amt8-2291-2015 Li, L., & Xu, C.-Y. (2014). The comparison of sensitivity analysis of hydrological uncertainty estimates by GLUE and Bayesian method under the impact of precipitation errors. Stochastic Environmental Research and Risk Assessment, 28(3), 491–504. https://doi.org/10.1007/s00477-013-0767-1 Liang, X., & Guo, J. (2003). Intercomparison of land-surface parameterization schemes: Sensitivity of surface energy and water fluxes to model parameters. Journal of Hydrology, 279(1), 182–209. https://doi.org/10.1016/S0022-1694(03)00168-9 Liang, X., Lettenmaier, D. P., Wood, E. F., & Burges, S. J. (1994). A simple hydrologically based model of land surface water and energy fluxes for general circulation models. Journal of Geophysical Research: Atmospheres, 99(D7), 14415–14428. https://doi.org/10.1029/94JD00483 Liang, X., Wood, E. F., & Lettenmaier, D. P. (1996). Surface soil moisture parameterization of the VIC-2L model: Evaluation and modification. Global and Planetary Change, 13(1), 195–206. https://doi.org/10.1016/0921-8181(95)00046-1 Liang, X., Wood, E. F., & Lettenmaier, D. P. (1999). Modeling ground heat flux in land surface parameterization schemes. Journal of Geophysical Research: Atmospheres, 104(D8), 9581–9600. https://doi.org/10.1029/98JD02307 Liang, X., Xie, Z., & Huang, M. (2003). A new parameterization for surface and groundwater interactions and its impact on water budgets with the variable infiltration capacity 153 (VIC) land surface model. Journal of Geophysical Research: Atmospheres, 108(D16), 8613. https://doi.org/10.1029/2002JD003090 Lilhare, R., Déry, S. J., Pokorny, S., Stadnyk, T. A., & Koenig, K. A. (2019). Intercomparison of multiple hydro-climatic datasets across the Lower Nelson River Basin, Manitoba, Canada. Atmosphere-Ocean, 57(4), 262–278. https://doi.org/10.1080/07055900.2019.1638226 Lilhare, R., Pokorny, S., Déry, S. J., Stadnyk, T. A., & Koenig, K. A. (2020). Sensitivity analysis and uncertainty assessment in water budgets simulated by the Variable Infiltration Capacity model for Canadian sub-arctic watersheds. Hydrological Processes, 34(9), 2057–2075. https://doi.org/10.1002/hyp.13711 Liu, Y., & Gupta, H. V. (2007). Uncertainty in hydrologic modeling: Toward an integrated data assimilation framework. Water Resources Research, 43, W07401. https://doi.org/10.1029/2006WR005756 Liu, Y., Tian, F., Hu, H., & Sivapalan, M. (2014). Socio-hydrologic perspectives of the coevolution of humans and water in the Tarim River basin, Western China: The Taiji– Tire model. Hydrology and Earth System Sciences, 18(4), 1289–1303. https://doi.org/10.5194/hess-18-1289-2014 Lohmann, D., Nolte-Holube, R., & Raschke, E. (1996). A large-scale horizontal routing model to be coupled to land surface parametrization schemes. Tellus A, 48(5), 708– 721. https://doi.org/10.1034/j.1600-0870.1996.t01-3-00009.x Lohmann, D., Raschke, E., Nijssen, B., & Lettenmaier, D. P. (1998). Regional scale hydrology: I. Formulation of the VIC-2L model coupled to a routing model. Hydrological Sciences Journal, 43(1), 131–141. https://doi.org/10.1080/02626669809492107 MacDonald, M. K., Stadnyk, T. A., Déry, S. J., Braun, M., Gustafsson, D., Isberg, K., & Arheimer, B. (2018). Impacts of 1.5 and 2.0 °C warming on pan-arctic river discharge into the Hudson Bay complex through 2070. Geophysical Research Letters, 45(15), 7561–7570. https://doi.org/10.1029/2018GL079147 Maggioni, V., Sapiano, M. R., Adler, R. F., Tian, Y., & Huffman, G. J. (2014). An error model for uncertainty quantification in high-time-resolution precipitation products. Journal of Hydrometeorology, 15(3), 1274–1292. https://doi.org/10.1175/JHM-D-130112.1 Mahfouf, J.-F., Brasnett, B., & Gagnon, S. (2007). A Canadian precipitation analysis (CaPA) project: Description and preliminary results. Atmosphere-Ocean, 45(1), 1–17. https://doi.org/10.3137/ao.v450101 Manitoba Hydro. (2015). Climate Change Report—Fiscal Year 2014-2015 (p. 48). Manitoba Hydro. Retrieved from Manitoba Hydro website: https://www.hydro.mb.ca/environment/pdf/climate_change_report_2014_15.pdf Mann, H. B. (1945). Nonparametric tests against trend. Econometrica: Journal of the Econometric Society, 13, 245–259. https://www.jstor.org/stable/1907187 154 Mao, D., & Cherkauer, K. A. (2009). Impacts of land-use change on hydrologic responses in the Great Lakes region. Journal of Hydrology, 374(1), 71–82. https://doi.org/10.1016/j.jhydrol.2009.06.016 Maurer, E. P., Wood, A. W., Adam, J. C., Lettenmaier, D. P., & Nijssen, B. (2002). A longterm hydrologically based dataset of land surface fluxes and states for the conterminous United States. Journal of Climate, 15(22), 3237–3251. https://doi.org/10.1175/1520-0442(2002)015<3237:ALTHBD>2.0.CO;2 McCabe, G. J., Wolock, D. M., & Valentin, M. (2018). Warming is Driving Decreases in Snow Fractions While Runoff Efficiency Remains Mostly Unchanged in SnowCovered Areas of the Western United States. Journal of Hydrometeorology, 19(5), 803–814. https://doi.org/10.1175/JHM-D-17-0227.1 McCarthy, J. J., Canziani, O. F., Leary, N. A., Dokken, D. J., & White, K. S. (2001). Climate Change 2001: Impacts, Adaptation and Vulnerability, Contribution of Working Group II to the Third Assessment Report of the Intergovernmental Panel on Climate Change (IPCC). Cambridge University Press. Retrieved from https://library.harvard.edu/collections/ipcc/docs/27_WGIITAR_FINAL.pdf McKenney, D. W., Hutchinson, M. F., Papadopol, P., Lawrence, K., Pedlar, J., Campbell, K., … Owen, T. (2011). Customized spatial climate models for North America. Bulletin of the American Meteorological Society, 92(12), 1611–1622. https://doi.org/10.1175/2011BAMS3132.1 Mearns, L. O., Sain, S., Leung, L. R., Bukovsky, M. S., McGinnis, S., Biner, S., … Sloan, L. (2013). Climate change projections of the North American Regional Climate Change Assessment Program (NARCCAP). Climatic Change, 120(4), 965–975. https://doi.org/10.1007/s10584-013-0831-3 Mekis, E., & Hogg, W. D. (1999). Rehabilitation and analysis of Canadian daily precipitation time series. Atmosphere-Ocean, 37(1), 53–85. https://doi.org/10.1080/07055900.1999.9649621 Mekonnen, G. B., Matula, S., Doležal, F., & Fišák, J. (2015). Adjustment to rainfall measurement undercatch with a tipping-bucket rain gauge using ground-level manual gauges. Meteorology and Atmospheric Physics, 127(3), 241–256. https://doi.org/10.1007/s00703-014-0355-z Mesinger, F., DiMego, G., Kalnay, E., Mitchell, K., Shafran, P. C., Ebisuzaki, W., … Shi, W. (2006). North American regional reanalysis. Bulletin of the American Meteorological Society, 87(3), 343–360. https://doi.org/10.1175/BAMS-87-3-343 Middelkoop, H., Daamen, K., Gellens, D., Grabs, W., Kwadijk, J. C., Lang, H., … Wilke, K. (2001). Impact of climate change on hydrological regimes and water resources management in the Rhine basin. Climatic Change, 49(1–2), 105–128. https://doi.org/10.1023/A:1010784727448 Milewska, E. J., Hopkinson, R. F., & Niitsoo, A. (2005). Evaluation of geo-referenced grids of 1961–1990 Canadian temperature and precipitation normals. Atmosphere-Ocean, 43(1), 49–75. https://doi.org/10.3137/ao.430104 155 Miller, D. A., & White, R. A. (1998). A conterminous United States multilayer soil characteristics dataset for regional climate and hydrology modeling. Earth Interactions, 2(2), 1–26. https://doi.org/10.1175/10873562(1998)002<0001:ACUSMS>2.3.CO;2 Milly, P. C. D., Dunne, K. A., & Vecchia, A. V. (2005). Global pattern of trends in streamflow and water availability in a changing climate. Nature, 438(7066), 347–350. https://doi.org/10.1038/nature04312 Mishra, V., & Lilhare, R. (2016). Hydrologic sensitivity of Indian sub-continental river basins to climate change. Global and Planetary Change, 139, 78–96. https://doi.org/10.1016/j.gloplacha.2016.01.003 Modarres, R., & Sarhadi, A. (2009). Rainfall trends analysis of Iran in the last half of the twentieth century. Journal of Geophysical Research: Atmospheres, 114, D03101. https://doi.org/10.1029/2008JD010707 Mondal, A., Kundu, S., & Mukhopadhyay, A. (2012). Rainfall trend analysis by MannKendall test: A case study of north-eastern part of Cuttack district, Orissa. International Journal of Geology, Earth and Environmental Sciences, 2(1), 70–78. Mpelasoka, F. S., & Chiew, F. H. S. (2009). Influence of rainfall scenario construction methods on runoff projections. Journal of Hydrometeorology, 10(5), 1168–1183. https://doi.org/10.1175/2009JHM1045.1 Myneni, R. B., Ramakrishna, R., Nemani, R., & Running, S. W. (1997). Estimation of global leaf area index and absorbed PAR using radiative transfer models. IEEE Transactions on Geoscience and Remote Sensing, 35(6), 1380–1393. https://doi.org/10.1109/36.649788 Nash, J., & Sutcliffe, J. V. (1970). River flow forecasting through conceptual models part I— A discussion of principles. Journal of Hydrology, 10(3), 282–290. https://doi.org/10.1016/0022-1694(70)90255-6 Natural Resources Canada. (2010, December 31). Permafrost. Retrieved 23 December 2016, from http://geogratis.gc.ca/api/en/nrcan-rncan/ess-sst/dc7107c0-8893-11e0-aa106cf049291510.html Natural Resources Canada. (2014). Regional, national and international climate modeling | Forests | Natural Resources Canada: Introduction. Retrieved 3 November 2016, from http://cfs.nrcan.gc.ca/projects/3?lang=en_CA Natural Resources Canada. (2016). Potential Evapotranspiration—Open Government Portal. Retrieved 9 April 2018, from https://open.canada.ca/data/en/dataset/a650fb95-1104504b-a344-0fe05877b7cf Nelson, F. E. (2003). (Un) frozen in time. Science, 299(5613), 1673–1675. https://doi.org/10.1126/science.1081111 Newbury, R., & Malaher, G. (1973). The destruction of Manitoba’s last great river. Ottawa: Canadian Nature Federation. Retrieved from http://www.arlis.org/docs/vol2/hydropower/APA_DOC_no._887.pdf 156 Newlands, N. K., Davidson, A., Howard, A., & Hill, H. (2011). Validation and intercomparison of three methodologies for interpolating daily precipitation and temperature across Canada. Environmetrics, 22(2), 205–223. https://doi.org/10.1002/env.1044 Nijssen, B., & Lettenmaier, D. P. (2004). Effect of precipitation sampling error on simulated hydrological fluxes and states: Anticipating the Global Precipitation Measurement satellites. Journal of Geophysical Research: Atmospheres, 109, D02103. https://doi.org/10.1029/2003JD003497 Nijssen, B., Lettenmaier, D. P., Liang, X., Wetzel, S. W., & Wood, E. F. (1997). Streamflow simulation for continental-scale river basins. Water Resources Research, 33(4), 711– 724. https://doi.org/10.1029/96WR03517 Nijssen, B., Schnur, R., & Lettenmaier, D. P. (2001). Global retrospective estimation of soil moisture using the variable infiltration capacity land surface model, 1980-93. Journal of Climate, 14(8), 1790–1808. https://doi.org/10.1175/15200442(2001)014<1790:GREOSM>2.0.CO;2 North American Land Change Monitoring System. (2010). 2010 North American Land Cover at 250 m spatial resolution. Produced by Natural Resources Canada/Canadian Center for Remote Sensing (NRCan/CCRS), United States Geological Survey (USGS); Insituto Nacional de Estadística y Geografía (INEGI), Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) and Comisión Nacional Forestal (CONAFOR). Retrieved 14 April 2017, from The USGS Land Cover Institute (LCI) website: http://www.cec.org/tools-and-resources/map-files/land-cover2010 Oechel, W. C., Hastings, S. J., Vourlrtis, G., Jenkins, M., Riechers, G., & Grulke, N. (1993). Recent change of Arctic tundra ecosystems from a net carbon dioxide sink to a source. Nature, 361(6412), 520–523. https://doi.org/10.1038/361520a0 O’Neil, H. C. L., Prowse, T. D., Bonsal, B. R., & Dibike, Y. B. (2017). Spatial and temporal characteristics in streamflow-related hydroclimatic variables over western Canada. Part 1: 1950–2010. Hydrology Research, 48(4), 915–931. https://doi.org/10.2166/nh.2016.057 Patterson, J. C., & Hamblin, P. F. (1988). Thermal simulation of a lake with winter ice cover. Limnology and Oceanography, 33(3), 323–338. https://doi.org/10.4319/lo.1988.33.3.0323 Pavelsky, T. M., & Smith, L. C. (2006). Intercomparison of four global precipitation data sets and their correlation with increased Eurasian river discharge to the Arctic Ocean. Journal of Geophysical Research: Atmospheres, 111, D21112. https://doi.org/10.1029/2006JD007230 Phillips, D. (1990). The climates of Canada (Vol. 19). Canadian Government Publishing Centre, Environment Canada, Ottawa, Ontario. Poitras, V., Sushama, L., Seglenieks, F., Khaliq, M. N., & Soulis, E. (2011). Projected changes to streamflow characteristics over western Canada as simulated by the Canadian RCM. Journal of Hydrometeorology, 12(6), 1395–1413. https://doi.org/10.1175/JHM-D-10-05002.1 157 Poulin, A., Brissette, F., Leconte, R., Arsenault, R., & Malo, J.-S. (2011). Uncertainty of hydrological modelling in climate change impact studies in a Canadian, snowdominated river basin. Journal of Hydrology, 409(3–4), 626–636. https://doi.org/10.1016/j.jhydrol.2011.08.057 Razavi, S., & Gupta, H. V. (2016a). A new framework for comprehensive, robust, and efficient global sensitivity analysis: 1. Theory. Water Resources Research, 52(1), 423–439. https://doi.org/10.1002/2015WR017558 Razavi, S., & Gupta, H. V. (2016b). A new framework for comprehensive, robust, and efficient global sensitivity analysis: 2. Application. Water Resources Research, 52(1), 440–455. https://doi.org/10.1002/2015WR017559 Razavi, S., Sheikholeslami, R., Gupta, H. V., & Haghnegahdar, A. (2019). VARS-TOOL: A toolbox for comprehensive, efficient, and robust sensitivity and uncertainty analysis. Environmental Modelling & Software, 112, 95–107. https://doi.org/10.1016/j.envsoft.2018.10.005 Reed, S., Koren, V., Smith, M., Zhang, Z., Moreda, F., Seo, D.-J., & DMIP Participants. (2004). Overall distributed model intercomparison project results. Journal of Hydrology, 298(1), 27–60. https://doi.org/10.1016/j.jhydrol.2004.03.031 Renard, B., Kavetski, D., Kuczera, G., Thyer, M., & Franks, S. W. (2010). Understanding predictive uncertainty in hydrologic modeling: The challenge of identifying input and structural errors. Water Resources Research, 46(5), W05521. https://doi.org/10.1029/2009WR008328 Rood, S. B., Samuelson, G. M., Weber, J. K., & Wywrot, K. A. (2005). Twentieth-century decline in streamflows from the hydrographic apex of North America. Journal of Hydrology, 306(1), 215–233. https://doi.org/10.1016/j.jhydrol.2004.09.010 Rouse, W. R. (1984). Microclimate of arctic tree line 2. Soil microclimate of tundra and forest. Water Resources Research, 20(1), 67–73. https://doi.org/10.1029/WR020i001p00067 Sabarly, F., Essou, G., Lucas-Picher, P., Poulin, A., & Brissette, F. (2016). Use of four reanalysis datasets to assess the terrestrial branch of the water cycle over Quebec, Canada. Journal of Hydrometeorology, 17(5), 1447–1466. https://doi.org/10.1175/JHM-D-15-0093.1 Sauchyn, D., Vanstone, J., & Perez-Valdivia, C. (2011). Modes and forcing of hydroclimatic variability in the Upper North Saskatchewan River Basin Since 1063. Canadian Water Resources Journal, 36(3), 205–217. https://doi.org/10.4296/cwrj3603889 Schewe, J., Heinke, J., Gerten, D., Haddeland, I., Arnell, N. W., Clark, D. B., … Kabat, P. (2014). Multimodel assessment of water scarcity under climate change. Proceedings of the National Academy of Sciences, 111(9), 3245–3250. https://doi.org/10.1073/pnas.1222460110 Schneider, U., Fuchs, T., Meyer-Christoffer, A., & Rudolf, B. (2008). Global Precipitation Analysis Products of the GPCC [Technical]. Deutscher Wetterdienst, Offenbach, Germany. 158 Seager, R., Neelin, D., Simpson, I., Liu, H., Henderson, N., Shaw, T., … Cook, B. (2014). Dynamical and thermodynamical causes of large-scale changes in the hydrological cycle over North America in response to global warming. Journal of Climate, 27(20), 7921–7948. https://doi.org/10.1175/JCLI-D-14-00153.1 Sen, P. K. (1968). Estimates of the regression coefficient based on Kendall’s tau. Journal of the American Statistical Association, 63(324), 1379–1389. https://doi.org/10.1080/01621459.1968.10480934 Serreze, M. C., Walsh, J. E., Chapin, F. S., Osterkamp, T., Dyurgerov, M., Romanovsky, V., … Morison, J. (2000). Observational evidence of recent change in the northern highlatitude environment. Climatic Change, 46, 159–207. https://doi.org/10.1023/A:1005504031923 Shafii, M., Tolson, B., & Matott, S. L. (2015). Addressing subjective decision-making inherent in GLUE-based multi-criteria rainfall–runoff model calibration. Journal of Hydrology, 523, 693–705. https://doi.org/10.1016/j.jhydrol.2015.01.051 Sharma, A. R., & Déry, S. J. (2016). Elevational dependence of air temperature variability and trends in British Columbia’s Cariboo Mountains, 1950–2010. Atmosphere-Ocean, 54(2), 153–170. https://doi.org/10.1080/07055900.2016.1146571 Sharma, T., Vittal, H., Chhabra, S., Salvi, K., Ghosh, S., & Karmakar, S. (2018). Understanding the cascade of GCM and downscaling uncertainties in hydro-climatic projections over India. International Journal of Climatology, 38(S1), e178–e190. https://doi.org/10.1002/joc.5361 Sheffield, J., & Wood, E. F. (2007). Characteristics of global and regional drought, 1950– 2000: Analysis of soil moisture data from off-line simulation of the terrestrial hydrologic cycle. Journal of Geophysical Research: Atmospheres, 112(D17), D17115. https://doi.org/10.1029/2006JD008288 Shen, Y., Xiong, A., Wang, Y., & Xie, P. (2010). Performance of high-resolution satellite precipitation products over China. Journal of Geophysical Research: Atmospheres, 115, D02114. https://doi.org/10.1029/2009JD012097 Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 1968 23rd ACM National Conference, 517–524. ACM. Shi, X., Wood, A. W., & Lettenmaier, D. P. (2008). How Essential is Hydrologic Model Calibration to Seasonal Streamflow Forecasting? Journal of Hydrometeorology, 9(6), 1350–1363. https://doi.org/10.1175/2008JHM1001.1 Shrestha, R. R., Schnorbus, M. A., Werner, A. T., & Berland, A. J. (2012). Modelling spatial and temporal variability of hydrologic impacts of climate change in the Fraser River basin, British Columbia, Canada. Hydrological Processes, 26(12), 1840–1860. https://doi.org/10.1002/hyp.9283 Simmons, A. J., & Poli, P. (2015). Arctic warming in ERA-Interim and other analyses. Quarterly Journal of the Royal Meteorological Society, 141(689), 1147–1162. https://doi.org/10.1002/qj.2422 Simmons, A. J., Poli, P., Dee, D. P., Berrisford, P., Hersbach, H., Kobayashi, S., & Peubey, C. (2014). Estimating low-frequency variability and trends in atmospheric 159 temperature using ERA-Interim. Quarterly Journal of the Royal Meteorological Society, 140(679), 329–353. https://doi.org/10.1002/qj.2317 Singh, V. P. (2018). Hydrologic modeling: Progress and future directions. Geoscience Letters, 5(1), 15. https://doi.org/10.1186/s40562-018-0113-z Skinner, W. R., & Gullett, D. W. (1993). Trends of daily maximum and minimum temperature in Canada during the past century. Climatological Bulletin, 27(2), 63–77. Skok, G., Žagar, N., Honzak, L., Žabkar, R., Rakovec, J., & Ceglar, A. (2016). Precipitation intercomparison of a set of satellite- and raingauge-derived datasets, ERA Interim reanalysis, and a single WRF regional climate simulation over Europe and the North Atlantic. Theoretical and Applied Climatology, 123(1), 217–232. https://doi.org/10.1007/s00704-014-1350-5 Slater, A. G., & Lawrence, D. M. (2013). Diagnosing present and future permafrost from climate models. Journal of Climate, 26(15), 5608–5623. https://doi.org/10.1175/JCLI-D-12-00341.1 Smith, A., Delavau, C., & Stadnyk, T. (2015). Identification of geographical influences and flow regime characteristics using regional water isotope surveys in the lower Nelson River, Canada. Canadian Water Resources Journal, 40(1), 23–35. https://doi.org/10.1080/07011784.2014.985512 Smith, L. C., Sheng, Y., & MacDonald, G. M. (2007). A first pan-Arctic assessment of the influence of glaciation, permafrost, topography and peatlands on northern hemisphere lake distribution. Permafrost and Periglacial Processes, 18(2), 201–208. https://doi.org/10.1002/ppp.581 Smith, S. L., & Burgess, M. (2004). Sensitivity of permafrost to climate warming in Canada. Ottawa, ON (Canada): Geological Survey of Canada, Ottawa, ON (Canada). Smith, S. L., & Burgess, M. M. (1998). Mapping the response of permafrost in Canada to climate warming. Current Research - Geological Survey of Canada, 163–191. St. Jacques, J.-M., & Sauchyn, D. J. (2009). Increasing winter baseflow and mean annual streamflow from possible permafrost thawing in the Northwest Territories, Canada. Geophysical Research Letters, 36, L01401. https://doi.org/10.1029/2008GL035822 Stadnyk, T. A., Déry, S. J., MacDonald, M. K., & Koenig, K. A. (2019). Theme I physical environment: The freshwater system. In From science to policy in the greater Hudson Bay marine region: An Integrated Regional Impact Study (IRIS) of climate change and modernization (pp. 113–153). Québec City, Québec, Canada: ArcticNet. Retrieved from https://arcticnet.ulaval.ca/docs/IRIS3-HudsonBay-VF.pdf Stedinger, J. R., Vogel, R. M., Lee, S. U., & Batchelder, R. (2008). Appraisal of the generalized likelihood uncertainty estimation (GLUE) method. Water Resources Research, 44(12), W00B06. https://doi.org/10.1029/2008WR006822 Suh, M. S., Oh, S. G., Lee, D. K., Cha, D. H., Choi, S. J., Jin, C. S., & Hong, S. Y. (2012). Development of new ensemble methods based on the performance skills of regional climate models over South Korea. Journal of Climate, 25(20), 7067–7082. https://doi.org/10.1175/JCLI-D-11-00457.1 160 Sushama, L., Laprise, R., Caya, D., Frigon, A., & Slivitzky, M. (2006). Canadian RCM projected climate-change signal and its sensitivity to model errors. International Journal of Climatology, 26(15), 2141–2159. https://doi.org/10.1002/joc.1362 Sutton, D., Goulden, M. L., Bazzaz, A., Daube, B. C., Fan, S. M., Munger, J. W., & Wofsy, S. C. (1999). BOREAS TF-03 NSA-OBS Tower Flux, Meteorological, and Soil Temperature Data. ORNL Distributed Active Archive Center. https://doi.org/10.3334/ORNLDAAC/361 Swain, S., & Hayhoe, K. (2015). CMIP5 projected changes in spring and summer drought and wet conditions over North America. Climate Dynamics, 44(9), 2737–2750. https://doi.org/10.1007/s00382-014-2255-9 Tam, B. Y., Szeto, K., Bonsal, B., Flato, G., Cannon, A. J., & Rong, R. (2019). CMIP5 drought projections in Canada based on the Standardized Precipitation Evapotranspiration Index. Canadian Water Resources Journal, 44(1), 90–107. https://doi.org/10.1080/07011784.2018.1537812 Tanarhte, M., Hadjinicolaou, P., & Lelieveld, J. (2012). Intercomparison of temperature and precipitation data sets based on observations in the Mediterranean and the Middle East. Journal of Geophysical Research: Atmospheres, 117(D12). https://doi.org/10.1029/2011JD017293 Tapiador, F. J., Turk, F. J., Petersen, W., Hou, A. Y., García-Ortega, E., Machado, L. A., … Huffman, G. J. (2012). Global precipitation measurement: Methods, datasets and applications. Atmospheric Research, 104, 70–97. https://doi.org/10.1016/j.atmosres.2011.10.021 Taubenböck, H., Wurm, M., Netzband, M., Zwenzner, H., Roth, A., Rahman, A., & Dech, S. (2011). Flood risks in urbanized areas–multi-sensoral approaches using remotely sensed data for risk assessment. Natural Hazards and Earth System Sciences, 11, 431–444. https://doi.org/10.5194/nhess-11-431-2011 Taylor, K. E., Stouffer, R. J., & Meehl, G. A. (2011). An overview of CMIP5 and the experiment design. Bulletin of the American Meteorological Society, 93(4), 485–498. https://doi.org/10.1175/BAMS-D-11-00094.1 Teutschbein, C., & Seibert, J. (2010). Regional climate models for hydrological impact studies at the catchment scale: A review of recent modeling strategies. Geography Compass, 4(7), 834–860. https://doi.org/10.1111/j.1749-8198.2010.00357.x Teutschbein, C., Wetterhall, F., & Seibert, J. (2011). Evaluation of different downscaling techniques for hydrological climate-change impact studies at the catchment scale. Climate Dynamics, 37(9–10), 2087–2105. https://doi.org/10.1007/s00382-010-0979-8 Thornton, P. E., Thornton, M. M., Mayer, B. W., Wei, Y., Devarakonda, R., Vose, R. S., & Cook, R. B. (2016). Daymet: Daily Surface Weather Data on a 1-km Grid for North America, Version 3. ORNL DAAC, Oak Ridge, Tennessee, USA. https://doi.org/10.3334/ORNLDAAC/1328 Tobin, C., Nicotina, L., Parlange, M. B., Berne, A., & Rinaldo, A. (2011). Improved interpolation of meteorological forcings for hydrologic applications in a Swiss Alpine 161 region. Journal of Hydrology, 401(1), 77–89. https://doi.org/10.1016/j.jhydrol.2011.02.010 Turk, F. J., Arkin, P., Sapiano, M. R. P., & Ebert, E. E. (2008). Evaluating high-resolution precipitation products. Bulletin of the American Meteorological Society, 89(12), 1911–1916. https://doi.org/10.1175/2008BAMS2652.1 Turner, A. G., & Annamalai, H. (2012). Climate change and the South Asian summer monsoon. Nature Climate Change, 2(8), 587–595. https://doi.org/10.1038/nclimate1495 United States Geological Survey. (2013). Shuttle Radar Topography Mission (SRTM) 1 ArcSecond Global | The Long Term Archive. Retrieved 6 April 2017, from https://lta.cr.usgs.gov/SRTM1Arc USACE. (1956). “Snow Hydrology” Summary Report of the Snow Investigations. Portaland, Oregon: North Pacific Division, US Army Corps of Engineers. Van Griensven, A., Meixner, T., Grunwald, S., Bishop, T., Diluzio, M., & Srinivasan, R. (2006). A global sensitivity analysis tool for the parameters of multi-variable catchment models. Journal of Hydrology, 324(1–4), 10–23. https://doi.org/10.1016/j.jhydrol.2005.09.008 van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., … Rose, S. K. (2011). The representative concentration pathways: An overview. Climatic Change, 109(1), 5. https://doi.org/10.1007/s10584-011-0148-z Vano, J. A., Nijssen, B., & Lettenmaier, D. P. (2015). Seasonal hydrologic responses to climate change in the Pacific Northwest. Water Resources Research, 51(4), 1959– 1976. https://doi.org/10.1002/2014WR015909 Velázquez, J. A., Schmid, J., Ricard, S., Muerth, M. J., St-Denis, B. G., Minville, M., … Turcotte, R. (2013). An ensemble approach to assess hydrological models’ contribution to uncertainties in the analysis of climate change impact on water resources. Hydrology and Earth System Sciences, 17(2), 565–578. https://doi.org/10.5194/hess-17-565-2013 Vrugt, J. A., Diks, C. G., Gupta, H. V., Bouten, W., & Verstraten, J. M. (2005). Improved treatment of uncertainty in hydrologic modeling: Combining the strengths of global optimization and data assimilation. Water Resources Research, 41(1), W01017. https://doi.org/10.1029/2004WR003059 Vrugt, J. A., Gupta, H. V., Bouten, W., & Sorooshian, S. (2003). A Shuffled Complex Evolution Metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters. Water Resources Research, 39(8), 1201. https://doi.org/10.1029/2002WR001642 Wagener, T., & Gupta, H. V. (2005). Model identification for hydrological forecasting under uncertainty. Stochastic Environmental Research and Risk Assessment, 19(6), 378– 387. https://doi.org/10.1007/s00477-005-0006-5 Walvoord, M. A., & Kurylyk, B. L. (2016). Hydrologic impacts of thawing permafrost—A review. Vadose Zone Journal, 15(6), 1–20. https://doi.org/10.2136/vzj2016.01.0010 162 Walvoord, M. A., & Striegl, R. G. (2007). Increased groundwater to stream discharge from permafrost thawing in the Yukon River basin: Potential impacts on lateral export of carbon and nitrogen. Geophysical Research Letters, 34, L12402. https://doi.org/10.1029/2007GL030216 Wang, A., Bohn, T. J., Mahanama, S. P., Koster, R. D., & Lettenmaier, D. P. (2009). Multimodel ensemble reconstruction of drought over the continental United States. Journal of Climate, 22(10), 2694–2712. https://doi.org/10.1175/2008JCLI2586.1 Wang, C., Wang, Z., Kong, Y., Zhang, F., Yang, K., & Zhang, T. (2019). Most of the Northern Hemisphere Permafrost Remains under Climate Change. Scientific Reports, 9(1), 3295. https://doi.org/10.1038/s41598-019-39942-4 Water Survey of Canada. (2016). Home—WaterOffice—Environment Canada. Retrieved 4 January 2017, from http://wateroffice.ec.gc.ca/ Weedon, G. P., Balsamo, G., Bellouin, N., Gomes, S., Best, M. J., & Viterbo, P. (2014). The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanalysis data. Water Resources Research, 50(9), 7505–7514. https://doi.org/10.1002/2014WR015638 Westermann, S., Langer, M., Boike, J., Heikenfeld, M., Peter, M., Etzelmüller, B., & Krinner, G. (2016). Simulating the thermal regime and thaw processes of ice-rich permafrost ground with the land-surface model CryoGrid 3. Geoscientific Model Development, 9(2), 523–546. https://doi.org/10.5194/gmd-9-523-2016 Westmacott, J. R., & Burn, D. H. (1997). Climate change effects on the hydrologic regime within the Churchill-Nelson River Basin. Journal of Hydrology, 202(1), 263–279. https://doi.org/10.1016/S0022-1694(97)00073-5 Wilby, R. L., & Harris, I. (2006). A framework for assessing uncertainties in climate change impacts: Low-flow scenarios for the River Thames, UK. Water Resources Research, 42(2), W02419. https://doi.org/10.1029/2005WR004065 Wilks, D. S. (2011). Statistical methods in the atmospheric sciences (3rd ed., Vol. 100). Academic Press, Elsevier. Williams, G. P., & Gold, L. W. (1976). Ground temperatures. CBD-180. Institute for Research in Construction. Originally Published July. Wong, J. S., Razavi, S., Bonsal, B. R., Wheater, H. S., & Asong, Z. E. (2017). Intercomparison of daily precipitation products for large-scale hydro-climatic applications over Canada. Hydrology and Earth System Sciences, 21(4), 2163–2185. https://doi.org/10.5194/hess-21-2163-2017 Woo, M.-K., & Thorne, R. (2006). Snowmelt contribution to discharge from a large mountainous catchment in subarctic Canada. Hydrological Processes, 20(10), 2129– 2139. https://doi.org/10.1002/hyp.6205 Wu, Z., Lu, G., Wen, L., Lin, C. A., Zhang, J., & Yang, Y. (2007). Thirty‐five year (1971– 2005) simulation of daily soil moisture using the variable infiltration capacity model over China. Atmosphere-Ocean, 45(1), 37–45. https://doi.org/10.3137/ao.v450103 163 Xie, Z., & Yuan, F. (2006). A parameter estimation scheme of the land surface model VIC using the MOPEX databases. IAHS Publication, 307, 169–179. Xu, C., Widén, E., & Halldin, S. (2005). Modelling hydrological consequences of climate change—Progress and challenges. Advances in Atmospheric Sciences, 22(6), 789– 797. https://doi.org/10.1007/BF02918679 Yapo, P. O., Gupta, H. V., & Sorooshian, S. (1998). Multi-objective global optimization for hydrologic models. Journal of Hydrology, 204(1–4), 83–97. https://doi.org/10.1016/S0022-1694(97)00107-8 Yue, S., Pilon, P., Phinney, B., & Cavadias, G. (2002). The influence of autocorrelation on the ability to detect trend in hydrological series. Hydrological Processes, 16(9), 1807–1829. https://doi.org/10.1002/hyp.1095 Yulianti, J. S., & Burn, D. H. (1998). Investigating links between climatic warming and low streamflow in the prairies region of Canada. Canadian Water Resources Journal / Revue Canadienne Des Ressources Hydriques, 23(1), 45–60. https://doi.org/10.4296/cwrj2301045 Zhang, J. L., Li, Y. P., Huang, G. H., Wang, C. X., & Cheng, G. H. (2016). Evaluation of uncertainties in input data and parameters of a hydrological model using a Bayesian framework: A case study of a snowmelt–precipitation-driven watershed. Journal of Hydrometeorology, 17(8), 2333–2350. https://doi.org/10.1175/JHM-D-15-0236.1 Zhang, Q., Sun, P., Singh, V. P., & Chen, X. (2012). Spatial-temporal precipitation changes (1956–2000) and their implications for agriculture in China. Global and Planetary Change, 82, 86–95. https://doi.org/10.1016/j.gloplacha.2011.12.001 Zhang, T., Barry, R. G., Knowles, K., Heginbottom, J. A., & Brown, J. (1999). Statistics and characteristics of permafrost and ground‐ice distribution in the Northern Hemisphere. Polar Geography, 23(2), 132–154. https://doi.org/10.1080/10889379909377670 Zhang, X., Vincent, L. A., Hogg, W. D., & Niitsoo, A. (2000). Temperature and precipitation trends in Canada during the 20th century. Atmosphere-Ocean, 38(3), 395–429. https://doi.org/10.1080/07055900.2000.9649654 Zhang, Y., Chen, W., & Riseborough, D. W. (2008a). Disequilibrium response of permafrost thaw to climate warming in Canada over 1850–2100. Geophysical Research Letters, 35, L02502. https://doi.org/10.1029/2007GL032117 Zhang, Y., Chen, W., & Riseborough, D. W. (2008b). Transient projections of permafrost distribution in Canada during the 21st century under scenarios of climate change. Global and Planetary Change, 60(3–4), 443–456. https://doi.org/10.1016/j.gloplacha.2007.05.003 Zhang, Y., Chen, W., Smith, S. L., Riseborough, D. W., & Cihlar, J. (2005). Soil temperature in Canada during the twentieth century: Complex responses to atmospheric climate change. Journal of Geophysical Research: Atmospheres, 110, D03112. https://doi.org/10.1029/2004JD004910 164 APPENDICES Appendix A I obtain gridded forcing datasets (i.e., Australian National University spline interpolation (ANUSPLIN), North American Regional Reanalysis (NARR), European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim), European Union Water and Global Change (WATCH) Forcing Data ERA-Interim (WFDEI), and Hydrological Global Forcing Data (HydroGFD), and an Inverse Distance Weighted (IDW) dataset constructed from 14 Environment and Climate Change Canada (ECCC) meteorological stations across the LNRB) for an intercomparison in our companion paper (Lilhare et al., 2019) and to further utilize these datasets in VIC model simulations over the LNRB. These datasets were selected based on the availability of specific variables (daily precipitation, minimum and maximum air temperature, and wind speed) that are required as input forcing for the VIC model. Moreover, this work forms a contribution to the Hudson Bay System (BaySys) project, where these six datasets are selected for further data analysis, intercomparison, and hydrological modelling using other land surface models, including VIC. Thus data availability (daily values from 1981 to 2010) and spatial resolution are other critical factors considered during the selection process of the datasets used for this study. The daily gridded, gauge-based dataset ANUSPLIN, developed by Natural Resources Canada (NRCan), is available for the entire Canadian landmass at 10 km spatial resolution (Hopkinson et al., 2011; McKenney et al., 2011). Preliminary analysis shows the updated version (1950–2015) of ANUSPLIN exhibits a consistent dry bias in precipitation over the LNRB; therefore, I retained the earlier version (1950–2011) of this dataset (Hopkinson et al., 2011). This product uses a trivariate thin-plate smoothing spline technique and includes daily 165 total precipitation (mm), minimum and maximum air temperature (°C) at 10 km spatial resolution based on 7514 meteorological stations (1950–2011) across the Canadian landmass south of 60°N. The ANUSPLIN dataset is readily available, relatively easy to apply, and accounts for spatially varying elevation relationships; however, it does not include daily wind speeds, simulates only elevation relationships and performs relatively poorly in sharp spatial elevation gradients (Newlands, Davidson, Howard, & Hill, 2011). The NARR reanalysis product is a high-resolution, regional extension of the National Centres for Environmental Prediction (NCEP)/National Centre for Atmospheric Research (NCAR) global reanalysis data (Kalnay et al., 1996; Kistler et al., 2001). It is developed at 32 km spatial and 3-hourly temporal resolution by utilizing a version of the Eta model and its three-dimensional variational data assimilation system (Mesinger et al., 2006) for the North American continent, available from 1979 to present. The NARR dataset is readily available at high resolution and provides accurate estimates of daily wind speeds, which are required for VIC model forcing; however, it exhibits warm biases in near-surface air temperatures over the study domain and lacks Canadian precipitation data assimilation after 2004 (Lilhare et al., 2019). The ERA-Interim dataset is a global reanalysis product from the European Centre for Medium-Range Weather Forecasts (ECMWF). Originally developed at 0.8° spatial resolution, data are also available at different spatial (0.125° to 3° grids) and 3-hourly temporal resolutions for January 1979 through near real-time (Dee et al., 2011). The ERA-I encompasses multiple variables, including precipitation, minimum and maximum air temperatures, and wind speeds as required by the VIC model, at a high spatial and temporal resolution. It has improved low-frequency variability and stratospheric circulation; however, 166 this product represents intense water cycling (precipitation, evaporation) over the oceans and positive biases in air temperature and humidity in the Arctic. The ERA-I is the wettest among all other datasets and overestimates daily precipitation, especially in the southern LNRB, and does not capture low-level inversions over the study domain (Lilhare et al., 2019). The WFDEI reanalysis dataset relies on a method developed by the European Union's Water and Global Change (WATCH) project (http://www.eu-watch.org) and incorporates in situ observations in the reanalysis (Weedon et al., 2014). The WFDEI dataset was derived from ERA-Interim data (Dee et al., 2011) and was improved by an elevation correction for numerous variables. Further, to ensure that the monthly statistics remain similar to the in situ observations of the Global Precipitation Climatology Center (GPCC), an undercatch correction is adopted whereby the daily variability of the reanalysis product is conserved (Weedon et al., 2014). I obtained the WFDEI-GPCC precipitation, maximum and minimum air temperature, and wind speed data at approximately 55 km spatial and daily temporal scale for this study. The recently developed HydroGFD dataset combines different reanalysis datasets to produce near-real-time, updated hydrologic forcing of precipitation and air temperature (Berg et al., 2018). The HydroGFD data use the already established WFDEI method but with updated climatological observations and interim products to produce near-real-time estimates of precipitation, air temperature, and wind speed at 3- and 6-hourly temporal and 0.50° spatial resolutions (Berg et al., 2018). The HydroGFD dataset is currently limited to precipitation, air temperature, and wind speeds, whereas WFDEI produces several additional variables (Weedon et al., 2014). 167 The IDW dataset of daily precipitation and mean air temperature is derived from 14 ECCC meteorological stations using spatial interpolation by applying the IDW interpolation method in a grid format at 0.10° (∼10 km) horizontal resolution for the LNRB. The IDW dataset is readily available, provides higher temporal resolution (i.e., hourly) data, is relatively easy to apply, and accounts for distance relationships. However, it does not include daily wind speeds due to incomplete records for our study period (1981–2010), accounts for distance effects only (i.e., neglects elevation effects), and limits heterogeneity of the domain. Moreover, meteorological stations are often located in open areas (e.g., airports), and thus may not be representative of forested environments or open bodies of water. Spatially regridded datasets (IDW, ANUSPLIN, NARR, ERA-Interim, WFDEI, and HydroGFD) at daily temporal and 10 km spatial resolutions are assembled to produce an ensemble mean dataset (i.e., referred to as the Ensemble) from 1981 to 2010. As the VIC model requires daily records of precipitation, minimum and maximum air temperatures, and wind speeds, I therefore prepared the Ensemble dataset using the equally weighted average method (Krishnamurti et al., 1999). Daily precipitation, maximum and minimum air temperatures are derived from the equally weighted average of all six gridded products. The daily wind speed ensemble is calculated from four reanalysis products (NARR, ERA-I, WFDEI and HydroGFD) as the other two datasets (IDW and ANUSPLIN) do not have such records. I assign equal weights to each dataset and so ignore the prior knowledge of their modelling capacity. This is one of the most common methods that assumes an equally weighted ensemble produces the best estimates of contemporary and future climate conditions since each model/dataset is equally likely to represent the truth (Krishnamurti et al., 1999; Suh et al., 2012). The concept of a multi-product ensemble has been used widely 168 over global and regional domains to examine uncertainty in forcings and provide input to land surface models under historical and projected future climate conditions (Fowler et al., 2007; Fowler & Kilsby, 2007; Mishra & Lilhare, 2016; A. Wang et al., 2009). 169 Appendix B Area averaged multidata-VIC simulated annual water balance (mm) precipitation (PCP, a), total runoff (TR, b), evapotranspiration (ET, c) and soil moisture (SM, d), represented by different columns, for the LNRB’s ten unregulated sub-basins based on IDW-VIC, ANUSPLIN-VIC, NARR-VIC, ERA-I-VIC, WFDEI-VIC and HydroGFD-VIC simulations (1981–2010). Red bars show multidata-VIC simulations mean, black error bars show inter-VIC simulations variation using standard deviations, while black dots represent the area averaged water balance estimations from Ensemble-VIC. 170 Appendix C Spatial differences of seasonal evapotranspiration (ET) (mm) for the LNRB’s ten unregulated sub-basins based on Ensemble-VIC (ENSEM) minus (1 st row) IDW-VIC, (2nd row) ANUSPLIN-VIC, (3rd row) NARR-VIC, (4th row) ERA-I-VIC, (5th row) WFDEI-VIC, and (6 th row) HydroGFD-VIC simulations, water years 1981–2010, for winter (DJF), spring (MAM), summer (JJA) and autumn (SON). 171 Appendix D Same as Appendix C but for the soil moisture (SM). 172 Appendix E Ratio of factor sensitivity (%) of IVARS 50 at a seasonal scale (1981–2010) for each parameter over all subwatersheds of the LNRB for the three model performance metrics (a) KGE, (b) NSE, and (c) PBIAS. Ratio of factor sensitivity is estimated by normalizing IVARS 50 values in each case, so they add up to 100% for all parameters.