UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analysis and forecasting of extreme temperature and precipitation across the complex terrain of British… Ivo Odon, Pedro 2019

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2019_may_ivoodon_pedro.pdf [ 23.5MB ]
Metadata
JSON: 24-1.0378419.json
JSON-LD: 24-1.0378419-ld.json
RDF/XML (Pretty): 24-1.0378419-rdf.xml
RDF/JSON: 24-1.0378419-rdf.json
Turtle: 24-1.0378419-turtle.txt
N-Triples: 24-1.0378419-rdf-ntriples.txt
Original Record: 24-1.0378419-source.json
Full Text
24-1.0378419-fulltext.txt
Citation
24-1.0378419.ris

Full Text

Analysis and Forecasting of ExtremeTemperature and PrecipitationAcross the Complex Terrain of British ColumbiabyPedro Ivo OdonA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate Studies(Atmospheric Sciences)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)May 2019c© Pedro Ivo Odon 2019The following individuals certify that they have read, and recommendto the Faculty of Graduate and Postdoctoral Studies for acceptance, thedissertation entitled:Analysis and Forecasting of Extreme Temperature and PrecipitationAcross the Complex Terrain of British Columbiasubmitted by Pedro Ivo Odon in partial fulfillment of the requirements forthe degree of Doctor of Philosophyin Atmospheric SciencesExamining Committee:Roland Stull, Atmospheric SciencesSupervisorGregory West, Atmospheric SciencesSupervisory Committee MemberIan McKendry, GeographySupervisory Committee MemberPhillip Austin, Atmospheric SciencesUniversity ExaminerBrett Eaton, GeographyUniversity ExamineriiAbstractThe ability to forecast extreme temperature and precipitation events notonly helps satisfy the public’s desire to better prepare for such events, butcan also provide valuable information about the future risks of such eventsto emergency managers, regional planners, and policy-makers at all levelsof government. Additionally, understanding extreme weather events in thecontext of climate change can help provide confidence and insight into riskassessment.This dissertation advances extreme weather forecasting over the com-plex topography of British Columbia (BC) while accounting for changes inintensity and frequency of extreme events due to nonstationarity.First the problem of finding a dataset to provide climatological distri-butions is addressed. Weather station data coverage, quality, and temporalcompleteness across BC degrade outside of population centres, and as onegoes back in time. This data paucity motivates the search for the best re-analysis to serve as a climatological reference dataset. The 2-m temperatureand daily accumulated precipitation of the reanalyses are compared with ob-servations from meteorological stations distributed over the complex terrainof British Columbia. The observations are separated into climate regionsby Principal Component Analysis (PCA) and K-means clustering, and newverification metrics are introduced to evaluate the best performing reanaly-sis. Upon thorough evaluation, the Japanese 55-year Reanalysis (JRA-55)was found to be best. Its biases are largely explained by the inability of thecoarse-resolution reanalysis to represent terrain characteristics.The second component of this works combines, downscales and bias cor-rects the best performing reanalysis using the high-spatial-resolution Parameter-Elevation Regressions on Independent Slopes Model (PRISM) dataset andusing surface weather station observations. This results in a high-resolution,long-term gridded dataset that is spatially and temporally complete, yieldinga very-high-resolution surface analysis (VHRSA). Biases and mean absoluteerrors are substantially improved over the JRA-55. The VHRSA not onlyrenders a solution to the paucity of observational data across BC, it alsohas the potential to be a valuable dataset for research and operational useiiiAbstractin meteorology, climate studies, and hydrology, to name a few.Next, this dataset is used to create a high-resolution, bias-correctedensemble forecast using the North American Ensemble Forecast System(NAEFS). The post-processed NAEFS is more skillful than the raw NAEFSforecast out to a forecast lead time of 10 days for both 2-m temperature,and daily accumulated precipitation, according to the verification metricsevaluated.Statistical temporal stationarity of extreme values of precipitation andtemperature are assessed for the 60-year VHRSA period. Statistical nonsta-tionarity is tested to determine if important temporal changes are requiredto characterize present-day extreme levels. It is determined that nonstation-ary distributions should be used to represent annual minima values of dailyminimum 2-m temperature during summer months and late winter.Finally, an extreme, or situational awareness index is presented: theParametric Extreme Index (PEI). It can be used to alert forecasters andother end users of future extreme temperature and precipitation events. Theindex is shown to be more skillful than the Standardized Anomaly (SA) asthe rarity of the event increases, with a higher number of hits and a lowernumber of misses. The PEI also takes into account the nonstationarity signalof minimum 2-m temperature.ivPrefaceThis dissertation is composed of five chapters which resulted in the publica-tion of three journal articles. These articles provide the contents of Chapters1, 2, and 3. The details of each article are described in more detail below.Chapter 4 is an article being written as of the time of this writing. Theconclusion is unique to this dissertation.Chapter 1A portion of Chapter 1 has been published in the Canadian Meteorologicaland Oceanographic Society (CMOS) Bulletin:• Odon, P., West, G., and Stull, R. (2017). Vancouver Fall and Winter2016/17: How Bad Was It? CMOS Bulletin SCMO, 45(4):9-12.The objective of the paper is to motivate studying the extreme weatherthat British Columbia (BC) faces throughout fall and winter, and the im-pacts of such events, in particular in Southwest BC where most of the pop-ulation resides. P. Odon carried out the research, analyzed the results, andwrote the paper. Professor Stull provided the computational and financialresources for the study. Professor Stull and Dr. West provided backgroundknowledge on meteorology and insights into the impacts on BC Hydro. Theyalso helped editing the paper.Chapter 2Chapter 2 has been published in Journal of Applied Meteorology and Cli-matology:• Odon, P., West, G., and Stull, R. (2018). Evaluation of reanalysesover British Columbia. Part I: Daily and Extreme 2-m temperature.Journal of Applied Meteorology and Climatology, 57(9):2091-2112.vChapter 3The paper evaluates how different the reanalysis datasets represent dailyand extreme 2-m temperature across the complex terrain of BC. P. Odoncarried out the research, analyzed the results, and wrote the manuscript.Professor Stull provided the computational and financial resources for thestudy and together with Dr. West, edited the contents of the manuscript.Additionally, Dr West provided input into the direction of the research.Chapter 3This is the second part of a two-part study. Chapter 2 is Part I. Chapter 3has been published in Journal of Applied Meteorology and Climatology:• Odon, P., West, G., and Stull, R. (2019). Performance of Reanalysesacross British Columbia. Part II: Evaluation of Daily and ExtremePrecipitation, (in press).The objective of this paper is to evaluate how the different reanaly-sis datasets represent daily and extreme precipitation across BC. P. Odoncarried out the research, analyzed the results, and wrote the manuscript.Professor Stull provided the computational and financial resources for thestudy and together with Dr West, edited the contents of the manuscript.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Understanding Extreme Weather Events and its Impacts —A Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.1 Fall . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.2 Winter . . . . . . . . . . . . . . . . . . . . . . . . . . 91.1.3 Anatomy of an Arctic Outbreak . . . . . . . . . . . . 141.1.4 Final Considerations . . . . . . . . . . . . . . . . . . 141.2 Extreme Data Analysis . . . . . . . . . . . . . . . . . . . . . 161.3 A Historical Record . . . . . . . . . . . . . . . . . . . . . . . 171.4 A Gridded Climatological Dataset . . . . . . . . . . . . . . . 181.5 Forecasting Extreme Weather Events . . . . . . . . . . . . . 191.6 Dissertation Organization . . . . . . . . . . . . . . . . . . . . 192 Performance of Reanalyses across British Columbia. Part I:Evaluation of Daily and Extreme 2-m Temperature . . . . 212.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.2 Data and Methodology . . . . . . . . . . . . . . . . . . . . . 222.2.1 Weather station data . . . . . . . . . . . . . . . . . . 252.2.2 CFSR . . . . . . . . . . . . . . . . . . . . . . . . . . . 27viiTable of Contents2.2.3 ERA-Interim . . . . . . . . . . . . . . . . . . . . . . . 272.2.4 JRA-55 . . . . . . . . . . . . . . . . . . . . . . . . . . 282.2.5 MERRA-2 . . . . . . . . . . . . . . . . . . . . . . . . 282.3 Climate zones . . . . . . . . . . . . . . . . . . . . . . . . . . 292.4 Verification Metrics . . . . . . . . . . . . . . . . . . . . . . . 322.4.1 Daily T2M . . . . . . . . . . . . . . . . . . . . . . . . 322.4.2 Extreme T2M . . . . . . . . . . . . . . . . . . . . . . 332.4.3 Analysis of Variance (ANOVA) . . . . . . . . . . . . 352.5 Nonstationarity . . . . . . . . . . . . . . . . . . . . . . . . . 352.5.1 Daily T2M . . . . . . . . . . . . . . . . . . . . . . . . 362.5.2 Extreme T2M . . . . . . . . . . . . . . . . . . . . . . 372.6 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 372.6.1 Daily T2M . . . . . . . . . . . . . . . . . . . . . . . . 372.6.2 Extreme T2M . . . . . . . . . . . . . . . . . . . . . . 482.7 Stationarity . . . . . . . . . . . . . . . . . . . . . . . . . . . 572.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623 Performance of Reanalyses across British Columbia. PartII: Evaluation of Daily and Extreme Precipitation . . . . . 643.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.2 Data and Methodology . . . . . . . . . . . . . . . . . . . . . 653.2.1 Weather station data . . . . . . . . . . . . . . . . . . 683.2.2 CFSR . . . . . . . . . . . . . . . . . . . . . . . . . . . 713.2.3 ERA-Interim . . . . . . . . . . . . . . . . . . . . . . . 713.2.4 JRA-55 . . . . . . . . . . . . . . . . . . . . . . . . . . 713.2.5 MERRA-2 . . . . . . . . . . . . . . . . . . . . . . . . 723.2.6 PRISM dataset . . . . . . . . . . . . . . . . . . . . . 723.3 Climate zones . . . . . . . . . . . . . . . . . . . . . . . . . . 733.4 Verification Metrics . . . . . . . . . . . . . . . . . . . . . . . 753.4.1 Daily PCP . . . . . . . . . . . . . . . . . . . . . . . . 753.4.2 Extreme PCP . . . . . . . . . . . . . . . . . . . . . . 763.4.3 Kruskal-Wallis Analysis . . . . . . . . . . . . . . . . . 793.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 803.5.1 Daily PCP . . . . . . . . . . . . . . . . . . . . . . . . 803.5.2 Extreme PCP . . . . . . . . . . . . . . . . . . . . . . 933.6 Stationarity . . . . . . . . . . . . . . . . . . . . . . . . . . . 993.6.1 Daily PCP . . . . . . . . . . . . . . . . . . . . . . . . 993.6.2 Extreme PCP . . . . . . . . . . . . . . . . . . . . . . 1053.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105viiiTable of Contents4 Analysis and Forecast of High-resolution Extreme Weather 1094.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 1094.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1104.2.1 Weather station data . . . . . . . . . . . . . . . . . . 1104.2.2 PRISM . . . . . . . . . . . . . . . . . . . . . . . . . . 1134.2.3 NAEFS . . . . . . . . . . . . . . . . . . . . . . . . . . 1134.3 The VHRSA: Downscaling and bias-correcting the JRA-55 . 1134.3.1 Verification of the VHRSA Daily and Extreme T2M . 1174.3.2 Verification of the VHRSA Daily and Extreme PCP . 1194.4 Bias-corrected NAEFS . . . . . . . . . . . . . . . . . . . . . 1234.4.1 Verification of the Forecast . . . . . . . . . . . . . . . 1294.4.2 Verification of Daily PCP . . . . . . . . . . . . . . . . 1304.4.3 Verification of Daily Maximum and Minimum T2M . 1334.5 Stationarity . . . . . . . . . . . . . . . . . . . . . . . . . . . 1344.5.1 Extreme T2M . . . . . . . . . . . . . . . . . . . . . . 1374.5.2 Extreme PCP . . . . . . . . . . . . . . . . . . . . . . 1394.6 The Parametric Extreme Index . . . . . . . . . . . . . . . . . 1394.6.1 Verification of the PEI . . . . . . . . . . . . . . . . . 1414.6.2 PEI Performance across T2M and PCP . . . . . . . . 1474.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1485 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156AppendicesA Surface weather stations . . . . . . . . . . . . . . . . . . . . . 180B Train and test stations . . . . . . . . . . . . . . . . . . . . . . 184ixList of Tables2.1 Overview of the four reanalysis datasets examined in this study. 242.2 Seasonally averaged systematic error (SE in ◦C) and randomerror (RE in ◦C) of mean values, and Kolmogorov-Smirnovstatistic (KS) of daily maximum T2M by climate zone . . . . 392.3 Seasonally averaged systematic error (SE in ◦C) and randomerror (RE in ◦C) of mean values, and Kolmogorov-Smirnovstatistic (KS) of daily minimum T2M by climate zone . . . . 422.4 Seasonally averaged systematic error of 2- and 30-year returnlevels of daily maximum T2M by climate zone (SE2 and SE30respectively in ◦C) . . . . . . . . . . . . . . . . . . . . . . . . 492.5 Seasonally averaged systematic error of 2- and 30-year returnlevels of daily minimum T2M by climate zone (SE2 and SE30respectively in ◦C) . . . . . . . . . . . . . . . . . . . . . . . . 523.1 Overview of the four reanalysis datasets examined in this study. 673.2 Averaged systematic error of: monthly precipitation total(M), two-sample χ2-statistic (χ2), 31-day-window percentageof wet days (W), and 30-year return levels of 1- (1D30), 3-(3D30), 7- (7D30), and 14-day (14D30) accumulated precipi-tation across wetter climate zones Northwest, Maritime West,Maritime East and Southwest, drier climate zones North,South Central, Central and Southeast, and all climate zonesin BC (all systematic errors but two-sample χ2-statistic in %) 814.1 Averaged systematic error and MAE of monthly precipita-tion total and monthly mean daily maximum and minimumT2M (MSE and MMAE respectively), of 2-year return levels(SE2 and MAE2 respectively), averaged systematic error oftwo-sample χ2-statistic (χ2) and KS statistics, across all teststations (all T2M errors but KS statistic in ◦C; all PCP errorsbut two-sample χ2 -statistic in % ) . . . . . . . . . . . . . . . 1224.2 Contingency table. . . . . . . . . . . . . . . . . . . . . . . . . 129xList of TablesA.1 Surface weather stations description. Station name and ab-breviation, longitude and latitude in degrees, elevation in me-ters, variable as either 2-m temperature (T) or 1-day accumu-lated precipitation (p), network and whether it is a SYNOPstation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180B.1 Train surface weather stations description. Station name andabbreviation, longitude and latitude in degrees, elevation inmeters, variable as either 2-m temperature (T) or 1-day ac-cumulated precipitation (p), network station. . . . . . . . . . 184B.2 Test surface weather stations description. Station name andabbreviation, longitude and latitude in degrees, elevation inmeters, variable as either 2-m temperature (T) or 1-day ac-cumulated precipitation (p), network station. . . . . . . . . . 186xiList of Figures1.1 A fine-resolution computer model wind forecast for December12th, 2017 at 1300 UTC. Vector orientation shows direction,colour shows speed. Terrain is grey-shaded. Geographic lo-cations mentioned in text shown with red labels, highwayslabelled in red. . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2 Fifteen-day accumulated precipitation anomaly centered onNovember 8th, 2017, relative to 1979-2017 period. Data fromthe the ECMWF interim reanalysis (ERA-Interim) (Dee et al.,2011). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.3 Cumulative number of rainy days at Vancouver InternationalAirport. Time series for all years between 1937/38 and 2016/17in grey, 1998/99 in red and 2016/17 in black. . . . . . . . . . 71.4 Consecutive rainy days at Vancouver International Airportfor the past 5 years in grey, with fall/winter 2016/17 in black. 81.5 Fifteen-day averaged minimum 2-m temperature anomaly cen-tered on December 11th, 2016, relative to 1979-2017 period. . 101.6 Daily minimum 2-m temperature time series for Vancouver,smoothed using a 5-day rolling window for readability. Timeseries for all years in grey, 1937/38-2016/17 average in blue,2012/13-2016/17 average in red, 1992/92 in orange and 2016/17in black. Days with snowfall in 2016/17 indicated with greendashed lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.7 (a) Public transit buses paralyzed by road conditions duringthe December 11th, 2016 arctic outbreak (Fritz, 2016). (b)A Google Maps screen capture, where orange and red coloursindicate widespread very slow or stopped traffic during theFebruary 3rd, 2017 arctic outbreak. . . . . . . . . . . . . . . . 13xiiList of Figures1.8 (a) 50.0-kPa geopotential height (km) and 2-m temperature(◦C), (b) sea level pressure (kPa) and 6-hr precipitation (mm)for December 4th, 2016 at 1200 UTC. (c) and (d), as in (a)and (b), but for December 6th, 2016 at 1200 UTC. Datafrom the the ECMWF interim reanalysis (ERA-Interim) (Deeet al., 2011). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1 Location of ECCC (red) and BC Hydro (blue) weather sta-tions and British Columbia population distribution (orange).Upside-down triangles indicate valley stations; upright tri-angles indicate non-valley stations. Dashed lines delineatethe dominant climate zones: North, Central, Maritime West,Southwest and Southeast. . . . . . . . . . . . . . . . . . . . . 262.2 Correlation matrix and five dominant climate zone clusters inBC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.3 Observed and reanalysis mean, systematic error, and randomerror magnitude for daily maximum T2M, averaged over sta-tions in the Central climate zone. The vertical dashed linesindicate the change in seasons and the colored dots representthe seasonal average. . . . . . . . . . . . . . . . . . . . . . . . 402.4 Observed and reanalysis mean, systematic error, and randomerror magnitude for daily minimum T2M, averaged over sta-tions in the North climate zone. The vertical dashed linesindicate the change in seasons and the colored dots representthe seasonal average. . . . . . . . . . . . . . . . . . . . . . . . 432.5 a) Kolmogorov-Smirnov (KS) statistic results for daily max-imum T2M, for all four reanalyses by time of year, averagedacross all locations in the Southeast climate zone. The verti-cal dashed lines indicate the change in seasons and the coloreddots represent the seasonal average. KS statistic values closerto zero are better. b) As an example, KS-statistic results fordaily maximum T2M on March 1 at Cranbrook (YXC) areshown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 452.6 Terrain elevation error in meters for ERA-Interim. Otherreanalyses’ errors are similar. . . . . . . . . . . . . . . . . . . 46xiiiList of Figures2.7 Mean systematic and mean random error of daily maximum(TMAX) and minimum (TMIN) T2M for each of the 57 sta-tions as a function of terrain reanalysis error. The solid linesshow the linear regression fits. The nearly horizontal linesshow the random error; the two lines are nearly on top ofeach other. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472.8 30-year return level and systematic error for maximum T2M,for all four reanalyses for Southeast climate zone. The verticaldashed lines indicate the change in seasons and the coloreddots represent the seasonal average. . . . . . . . . . . . . . . . 502.9 30-year return level and systematic error for minimum T2M,for all four reanalyses for Central climate zone. The verticaldashed lines indicate the change in seasons and the coloreddots represent the seasonal average. . . . . . . . . . . . . . . . 532.10 Mean systematic error of 2- and 30-year return levels (RL)of extreme maximum (TMAX) and minimum (TMIN) T2Mfor each of the 57 stations as a function of terrain reanalysiserror. The solid lines show the linear regression fits. . . . . . 552.11 MAE (◦C) of daily maximum and minimum T2M, and of ex-treme (2- and 30-year return levels) maximum and minimumT2M. MAE is averaged across all 57 stations. Smaller values(closer to the origin) are better. . . . . . . . . . . . . . . . . . 562.12 (a) Mean linear trend of daily minimum T2M; (b) Standarddeviation linear trend of daily minimum T2M; (c) Days whenN(µ(t), σ) is statistically significant. The vertical dashedlines indicate the change in seasons and the horizontal dashedlines delineate from top to bottom the North, Central, Mar-itime, Southwest and Southeast climate zones respectively. . . 582.13 (a) Mean linear trend of daily maximum T2M; (b) Standarddeviation linear trend of daily maximum T2M; (c) Days whenN(µ(t), σ) is statistically significant. The vertical dashedlines indicate the change in seasons and the horizontal dashedlines delineate from top to bottom the North, Central, Mar-itime, Southwest and Southeast climate zones respectively. . . 59xivList of Figures2.14 (a) Location linear trend of extreme maximum T2M; (b) Dayswhen GEV (µ(t), σ, κ) for extreme maximum T2M is statis-tically significant; (c) Location linear trend of extreme mini-mum T2M; (d) Days when GEV (µ(t), σ, κ) for extreme min-imum T2M is statistically significant. The vertical dashedlines indicate the change in seasons and the horizontal dashedlines delineate from top to bottom the North, Central, Mar-itime, Southwest and Southeast climate zones respectively. . . 613.1 Location of ECCC (red) and BC Hydro (blue) weather sta-tions used for precipitation analysis and British Columbiapopulation distribution (orange). Upside-down triangles in-dicate valley stations; upright triangles indicate non-valleystations. Dashed lines delineate the dominant climate zonesNorth, Northwest, Central, South Central, Maritime West,Maritime East, Southwest and Southeast. . . . . . . . . . . . 693.2 Cramer correlation matrix and eight dominant climate zoneclusters in BC. Crosses indicate stations where differences inprecipitation distribution are statistically significant. . . . . 743.3 (a) Diagram illustrating the 31-day centered rolling windowfor July 9th, performed over the 31-year period; (b) GEVmodel fitted over the 31-day, 31-year centered window for1-day precipitation total at Vancouver International Airport(YVR); (c) Quantile plot for the GEV fitted model for thesame day, location, and variable; If the GEV is a reasonablemodel, the points on the quantile plot should lie close to theunit diagonal; (d) Return level, return period plot for thesame day, location, and variable; showing the precipitationthat corresponds to a given return period. . . . . . . . . . . . 783.4 Observed and reanalysis running centered 31-day precipita-tion totals and systematic error averaged over stations in theNorthwest climate zone. The vertical dashed lines indicatethe change in seasons and the colored dots represent the sea-sonal average. . . . . . . . . . . . . . . . . . . . . . . . . . . . 823.5 Observed and reanalysis running centered 31-day precipita-tion totals and systematic error averaged over stations in theCentral climate zone. The vertical dashed lines indicate thechange in seasons and the colored dots represent the seasonalaverage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84xvList of Figures3.6 Percentage of (a) wet days; (b) two-sample χ2-statistic; per-centage of (c) dry days, (d) very light, (e) light (f) moderate,(g) heavy and (h) extreme precipitation events for MaritimeWest climate zone. The vertical dashed lines indicate thechange in seasons and the colored dots represent the seasonalaverage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 863.7 Percentage of (a) ”wet days”; (b) Two-sample χ2-statistic;Percentage of (c) ”dry days”, (d) “very light”, (e) “light”(f) “moderate”, (g) “heavy” and (h) “extreme” precipitationevents for North climate zone. The vertical dashed lines indi-cate the change in seasons and the colored dots represent theseasonal average. . . . . . . . . . . . . . . . . . . . . . . . . . 883.8 a) Winter precipitation totals of (a) PRISM. Winter precip-itation totals bilinearly interpolated to PRISM grid of (b)CFSR, (c) ERA-Interim, (d) JRA-55 and (e) MERRA-2. Sys-tematic error of (f) CFSR, (g) ERA-Interim, (h) JRA-55 and(i) MERRA-2. The dots represent weather station location. . 903.9 JRA-55 mean systematic error of 31-day precipitation totalsfor each of the 66 stations as a function of PRISM meansystematic error. The solid lines show the linear regression fits. 923.10 30-year return level and systematic error of 1-day precipita-tion total, for all four reanalyses for Southwest climate zone.The vertical dashed lines indicate the change in seasons andthe colored dots represent the seasonal average. . . . . . . . . 943.11 30-year return level and systematic error of 1-day precipita-tion total, for all four reanalyses for North climate zone. Thevertical dashed lines indicate the change in seasons and thecolored dots represent the seasonal average. . . . . . . . . . . 963.12 MAE of 31-day precipitation totals, χ2−statistic, ”wet days”,and 2-year and 30-year return levels of 1-, 3-, 7-, and 14-dayprecipitation totals. MAE is averaged across all 66 stations.Values closer to 0 at the origin of the plot are better. . . . . . 98xviList of Figures3.13 Seasonal precipitation total at Vancouver International Air-port (YVR) for July 9th (grey), 10-year left moving aver-age (black), decadal averages of seasonal precipitation total(1981-1990, 1991-2000 and 2001-2010) (blue) and 30-year sea-son precipitation total averaged over the study period 1981-2010 (green). The red line shows the linear regression fit tothe three decadal mean values. Individual 90% confidence in-tervals (CI; green dotted line) and multiple (1 − αWalker) ×100% CI (red dotted line) are drawn; α0 = 0.10. This ex-ample indicates that a stationary distribution is appropriatefor Jul 9th because the three decadal average values all fallwithin both CI. . . . . . . . . . . . . . . . . . . . . . . . . . . 1013.14 (a) Mean linear trend of season accumulated precipitation;(b) Mean linear trend of days with precipitation. The verticaldashed lines indicate the change in seasons and the horizontaldashed lines delineate from top to bottom the North, Central,Northwest, Maritime West, Maritime East, Southwest, SouthCentral and Southeast climate zones respectively. Units arepercent change over the 30-year period. . . . . . . . . . . . . 1033.15 Mean linear trend of (a) “very light” , (b) “light”, (c) “mod-erate” and (d) “extreme” precipitation events. The verticaldashed lines indicate the change in seasons and the horizontaldashed lines delineate from top to bottom the North, Central,Northwest, Maritime West, Maritime East, Southwest, SouthCentral and Southeast climate zones respectively. Units arepercent change over the 30-year period. . . . . . . . . . . . . 1044.1 Location of ECCC (blue) and BC Hydro (yellow) T2M weatherstations and British Columbia population distribution (red). . 1114.2 Location of ECCC (blue) and BC Hydro (yellow) PCP weatherstations and British Columbia population distribution (red). . 112xviiList of Figures4.3 (a) 30-year monthly mean of daily maximum T2M PRISMclimate on August 28th. (b) 2-year return level systematicerror between training stations and reanalysis interpolated totraining stations as a function of monthly mean daily max-imum T2M bias between JRA-55 and PRISM interpolatedto training stations. The solid line shows the linear regres-sion fit. (c) Downscaling difference derived by subtracting(a) from (e). (d) JRA-55 daily maximum T2M on 28 Aug2017. (e) Monthly mean of daily maximum T2M JRA-55climate billinearly interpolated to PRISM grid. (f) VHRSAdaily maximum T2M on 28 Aug 2017. . . . . . . . . . . . . . 1164.4 (a) Observed, JRA-55 and VHRSA monthly mean daily max-imum T2M averaged over test stations. (b) 2-year returnlevels of observed, JRA-55 and VHRSA maximum T2M av-eraged over test stations. (c) Monthly mean bias of JRA-55and VHRSA daily maximum T2M averaged over test stations.(d) As in (c) but for 2-year return levels. (e) As in (c) butfor MAE. (f) As in (e) but for 2-year return levels. (g) KSstatistic for daily maximum T2M. The vertical dashed linesindicate the change in seasons and the coloured dots representthe errors seasonal averages. In (a) and (b), values closer toobservations (black line) are better. In (c)-(g), values closerto zero are better. . . . . . . . . . . . . . . . . . . . . . . . . 1184.5 (a) Observed, JRA-55 and VHRSA monthly mean daily min-imum T2M averaged over test stations. (b) 2-year returnlevels of observed, JRA-55 and VHRSA minimum T2M av-eraged over test stations. (c) Monthly mean bias of JRA-55and VHRSA daily minimum T2M averaged over test stations.(d) Same as in (c) but for 2-year return levels. (e) Same asin (c) but for MAE. (f) Same as in (e) but for 2-year returnlevels. (g) KS statistic for daily minimum T2M. The verticaldashed lines indicate the change in seasons and the coloureddots represent the errors seasonal averages. Values of bias,MAE, and KS Statistic closer to zero are better. . . . . . . . 120xviiiList of Figures4.6 (a) Observed, JRA-55 and VHRSA mean of monthly PCPtotals averaged over test stations. b) 2-year return levels ofobserved, JRA-55 and VHRSA PCP averaged over test sta-tions. (c) Monthly mean bias of bias-corrected JRA-55 PCPaveraged over test stations. (d) Same as in (c) but for 2-yearreturn levels. (e) Same as in (c) but for MAE. (f) Same as in(e) but for 2-year return levels. g) χ2-statistic of PCP. Thevertical dashed lines indicate the change in seasons and thecoloured dots represent the errors seasonal averages. Valuesof bias, MAE, and χ2-statistic closer to zero are better. . . . 1214.7 (a) NAEFS raw Canadian control member forecast of dailyminimum T2M on previous day 2 Jan, 2017; (b) Same as in(a) but for VHRSA; (c) Ensemble mean raw NAEFS forecaston current day 3 Jan, 2017; (d) Ensemble mean bias-correctedNAEFS forecast on current day 3 Jan, 2017. DMB methoddone separately for each NAEFS member and then averagedto generate the bias-corrected ensemble mean. (e) VHRSAdaily minimum T2M on 3 Jan, 2017 presented for comparisonwith (c) and (d). The dots represent weather station location. 1264.8 (a) NAEFS raw Canadian control member forecast of 1-dayaccumulated PCP on previous day 7 Nov, 2016; (b) Sameas in (a) but for VHRSA; (c) Ensemble mean raw NAEFSforecast on current day 8 Nov, 2016; (d) Ensemble mean bias-corrected NAEFS forecast on current day 8 Nov, 2016. DMBmethod done separately for each NAEFS member and thenaveraged to generate the bias-corrected ensemble mean. (e)VHRSA 1-day accumulated PCP on 8 Nov, 2016 presentedfor comparison with (c) and (d). The dots represent weatherstation location. . . . . . . . . . . . . . . . . . . . . . . . . . 1274.9 (a) NAEFS raw Canadian control member forecast of dailymaximum T2M on previous day 27 Aug, 2017; (b) Same asin (a) but for VHRSA; (c) Ensemble mean raw NAEFS fore-cast on current day 28 Aug, 2017; (d) Ensemble mean bias-corrected NAEFS forecast on current day 28 Aug, 2017. DMBmethod done separately for each NAEFS member and thenaveraged to generate the bias-corrected ensemble mean. (e)VHRSA daily maximum T2M on 28 Aug, 2017 presented forcomparison with (c) and (d). The dots represent weatherstation location. . . . . . . . . . . . . . . . . . . . . . . . . . 128xixList of Figures4.10 a) TS for events above the 95th percentile during Fall 2016;(b) Same as in a) but for H; (c) Same as in a) but for F ; d)Systematic error of daily PCP; (e) Same as in a) but for b;(f) Same as in a) but for d ; (g) Same as in a) but for QBS;(h) Same as in a) but for CRPS. Values of TS and H closerto one are better. Values of F and systematic error closer tozero are better. Combined values of b = 0 and d = 1 arebetter. Values of QBS and CRPS closer to zero are better. . 1324.11 a) TS for events above the 90th percentile during Summer2017; (b) Same as in a) but for H; (c) Same as in a) but for F; d) Systematic error of daily maximum T2M. Values of TSand H closer to one are better. Values of F and systematicerror closer to zero are better. . . . . . . . . . . . . . . . . . 1334.12 a) Systematic error of daily minimum T2M during Winter2016/17; (b) Same as in a) but for QBS; (c) Same as in a)but for b ; (d) Same as in a) but for d ; (e) Same as in a) butfor CRPS. Values of systematic error, QBS and CRPS closerto zero are better. Combined values of b = 0 and d = 1 arebetter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1354.13 (a) Trend of extreme minimum T2M in July; (b) Locationswhere trend is statistically significant; (c) Locations whereGEV (µ(t), σ, κ) for extreme minimum T2M is statisticallysignificant. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1384.14 (a) Trend of extreme maximum T2M in August; (b) Locationswhere trend is statistically significant. . . . . . . . . . . . . . 1404.15 (a) The NCDF (purple) and the GCDF (black) show the prob-ability (y-axis) that the extreme value (F−1c (pc); x-axis) willnot be exceeded. (b) The NAEFS-derived probability thatF−1c (pc) will not be exceeded is given by Ff (F−1c (pc)). Thecomplement of the NCDF, given by 1 − Ff (F−1c (pc)), com-putes the NAEFS-derived probability that F−1c (pc) will beexceeded. The PEI is related to the area above the NCDFbounded by pc = 0 and pc = Fc(F−1f (1)) (pc = 0.12 in thefigure). PEI=0.02 indicating low extreme values have a lowprobability of exceedance. . . . . . . . . . . . . . . . . . . . . 1424.16 Day 1 PEI over BC on 3 Jan 2017. . . . . . . . . . . . . . . . 1434.17 Day 1 PEI over BC on 8 Nov 2016. . . . . . . . . . . . . . . . 1444.18 Day 1 PEI over BC on 28 Aug 2017. . . . . . . . . . . . . . . 145xxList of Figures4.19 EDS as function of 1 - BR for Day 1 PEI and SA across allstations during summer of 2017. Rarer events are to the right.Values of EDS closer to 1 are better. . . . . . . . . . . . . . . 147A.1 Elevation of 2-m temperature and 1-day accumulated precip-itation stations. . . . . . . . . . . . . . . . . . . . . . . . . . . 183xxiAcknowledgementsII thank my research supervisor, Dr. Roland Stull, who provided means andfunding for my dissertation, helped writing and editing journal papers, andprovided me with a great environment for my studies. I also wish to thankthe other members of my supervisory committee, Drs. Alex Cannon, GregWest and Ian McKendry.I am grateful for the positive and supportive environment in the office,and the help of many of the members of the Weather Forecast ResearchTeam. Without it I’m not sure my transition into Atmospheric Scienceswould have been possible. I must specifically thank Maggie Campbell andDavid Siuta for their patience, always allowing me to stop by their desksand ask questions about Atmospheric Sciences or coding; Dr. Greg Westfor help writing and editing journal papers, and many conversations whichhelped spark many of the ideas presented in this dissertation.Dr. Henryk Modzelewski, Roland Schigas and Tim Chui provided excel-lent assistance with IT questions, computing resources and coding.BC Hydro, the Canadian Natural Science and Engineering ResearchCouncil (NSERC), and Mitacs are acknowledged for providing funds forthis research.xxiiDedicationTo my wife Andreia. This degree would not be possible without your un-conditional love and support. This degree is yours as much as it is mine.xxiiiChapter 1IntroductionThe goal of this research is to provide accurate forecasts of extreme dailymaximum and minimum 2-m temperature (hereafter T2M) and 1-day ac-cumulated precipitation (hereafter PCP) over the mountainous terrain ofBritish Columbia (BC), Canada. Forecasting extreme events as well as fu-ture trends requires accurately identifying such patterns in space and time.During the October-March storm season, southwest BC is affected by ex-treme rainfall events in low-elevations, typically falling as snow at mid tohigh elevations (Odon et al., 2017). One type of storm process, known as anatmospheric river, makes landfall on the Coast Mountains of BC, producingheavy to extreme precipitation due to orographic lifting (Ralph et al., 2004,2005). During winter, arctic air outbreaks can result in record-breakinglow temperatures. These and other types of extreme cool season weatherevents cause floods, property and infrastructure damage, power outages, andbusiness and public-transportation disruption (Odon et al., 2017). Duringsummer, long droughts (Odon et al.) in combination with record-breakingsummer temperatures across BC result in high electricity demand, livestocksuffering, wildfire fighting costs that exceeded half a billion dollars in 2017alone, and insured property losses that reach hundreds of millions of dollars(Phillips, 2018).Complicating matters further, global climate models project an increasein the frequency and intensity of extreme temperature and precipitationevents (Kharin and Zwiers, 2000; Zwiers and Kharin, 2005; Zwiers andWehner, 2013). This will make the impacts of such events, such as in-creases in winter runoff leading to flooding, and overwhelmed drainage andsewage-system capacity (White et al., 2016; Sun et al., 2018b), more difficultto manage and mitigate.This necessitates increasing our understanding and ability to forecastextreme weather events. Such information can provide earlier and moreaccurate detection of such events, which in turn can help community re-sponders, government stakeholders and even the media to take appropriateaction to mitigate damage and reduce casualties (Herring et al., 2017).Although numerical weather prediction has greatly improved over the11.1. Understanding Extreme Weather Events and its Impacts — A Case Studypast decades, forecasting extreme T2M and PCP events, especially extremePCP beyond two days, remains challenging (Lalaurette, 2003).This dissertation has four goals: (1) Employ an appropriate statisticalmodel that specifically addresses extreme T2M and PCP events. Such amodel provides accurate estimation of future extreme levels only if suffi-cient historical data is available. Extreme data modelling, however, leads toa waste of information when implemented in practice, since only the maxi-mum observed values over a specified period are used. This creates a scarcityof data. Even without the scarcity problem, the simple availability of accu-rate and complete historical data can be an impediment. This leads to thenext two goals: (2) evaluate which reanalysis best represents observationsto serve as a surrogate climatological dataset with respect to daily and ex-treme minimum and maximum T2M, and daily and extreme PCP over BC;and (3) based on that evaluation, improve the best reanalysis if possible tocreate a more accurate gridded climatological dataset to better understandT2M and PCP extreme levels. Finally, (4) combine this gridded datasetwith gridded forecasts to create bias-corrected forecasts that provide moreaccurate information about extreme forecast events.Motivation as to why this research is important and its potential benefitsare provided in the next section, via a case-study overview of the extremecool season of 2016/2017 in southwest BC. Each of the four goals are intro-duced in the four following sections, with details presented in subsequentchapters.1.1 Understanding Extreme Weather Events andits Impacts — A Case StudyThe fall and winter seasons of 2016/17 were noteworthy for the cold, wetweather they brought to regions of British Columbia (BC). Although allof BC experienced some impacts, the South Coast region saw particularlylarge impacts from an exceptionally wet fall that quickly transitioned intoa persistently cold, snowy winter.1.1.1 FallOctober saw significant snow accumulations at mid to high elevations, andrecord rainfall amounts in many low-elevation areas across the South Coast.Temperatures warmed to record levels ahead of an atmospheric river eventon November 8th, 2016. Many daily records were set, with new monthly21.1. Understanding Extreme Weather Events and its Impacts — A Case Studyrecords of 19.4◦C in Vancouver and 22.4◦C in Abbotsford (see Fig. 1.1 forall mentioned locations). The combination of record warmth and rain-on-snow led to large snowmelt contributions to runoff, in addition to the heavyrainfall itself.31.1.UnderstandingExtremeWeatherEventsanditsImpacts—ACaseStudyFigure 1.1: A fine-resolution computer model wind forecast for December 12th, 2017 at 1300 UTC. Vector orien-tation shows direction, colour shows speed. Terrain is grey-shaded. Geographic locations mentioned in text shownwith red labels, highways labelled in red.41.1. Understanding Extreme Weather Events and its Impacts — A Case StudyFigure 1.2 illustrates the fifteen-day precipitation anomalies for the firsthalf of November. In several locations of the Coast it precipitated twicemore than normal.For the province’s primary electric utility, BC Hydro, the situation madefor challenging reservoir operations, balancing dam safety requirements withminimizing downstream flooding. Although most reservoirs across the SouthCoast were full or near full, the most concerning was the Campbell Riversystem. The upper part of the watershed received the normal Novembermonthly precipitation amount compressed into the first eight days of themonth (BC Hydro, 2016). Furthermore, one-week and two-week reservoirinflows set new records, and were estimated to exceed 1-in-100 year volumes.In emergency coordination calls and meetings, BC Hydro worked withnearby towns, regional districts, and the province, deciding to increase dis-charges to a record-tying 600 m3 · s−1, enough to fill an olympic-sized swim-ming pool every four seconds. This was done to mitigate the risk of over-topping the dam in subsequent storms, which would’ve meant passing thefull reservoir inflows (1,100 m3 ·s−1), flooding communities downstream (BCHydro, 2016).At Vancouver International Airport (hereafter, Vancouver) there wereonly three days without rain in October and only two in November. Infact, the fall and winter 2016/17 period had the second highest frequencyof rain since 1937/38, with 121 days of rain during the 182-day (6-month)period (Fig. 1.3). Only the 1998/99 period was higher, with 131 days ofrain. Furthermore, this 2016/17 period featured 22 consecutive days of rain(October 12th - November 2nd; Fig. 1.4). This was the longest stretch inthe past six years, and the ninth longest since 1937/38.51.1. Understanding Extreme Weather Events and its Impacts — A Case StudyFigure 1.2: Fifteen-day accumulated precipitation anomaly centered onNovember 8th, 2017, relative to 1979-2017 period. Data from the theECMWF interim reanalysis (ERA-Interim) (Dee et al., 2011).61.1.UnderstandingExtremeWeatherEventsanditsImpacts—ACaseStudyFigure 1.3: Cumulative number of rainy days at Vancouver International Airport. Time series for all years between1937/38 and 2016/17 in grey, 1998/99 in red and 2016/17 in black.71.1.UnderstandingExtremeWeatherEventsanditsImpacts—ACaseStudyFigure 1.4: Consecutive rainy days at Vancouver International Airport for the past 5 years in grey, with fall/winter2016/17 in black.81.1. Understanding Extreme Weather Events and its Impacts — A Case Study1.1.2 WinterIn December the pattern abruptly became much cold and drier. Five arcticair outbreaks brought anomalous cold temperatures, and snow to BritishColumbia in December through February.Figure 1.5 shows the 15-day average minimum 2-m temperature anomalyacross British Columbia centered on second arctic outbreak on December11th, 2016. It shows the difference between what British Columbians wouldexpect for that time of year and what actually happened. Temperatureswere much colder than average, with much of BC more than 10◦C belownormal.In Vancouver, daily minimum temperatures were below freezing for mostof December, January and the first part of February. This period was colderthan winters in recent memory (5-year average, red line in Fig. 1.6) andthe long term average (80-year average, blue line in Fig. 1.6). There were54 days of below freezing daily minimum temperature during the entire 6-month period, more than any such period since 1992/93 (orange line in Fig.1.6).The five arctic outbreaks can be identified in the smoothed minimumtemperature time series by the large departures from the 80-year average,occurring on approximately December 5th and 11th, 2016, January 1st and10th, 2017 and February 2nd, 2017; most followed by snow days (greendashed line in Fig. 1.6).These weather events directly impacted the general population as wellas government and industry. Garbage was left uncollected for weeks inneighbourhoods of Vancouver, Burnaby, Surrey and Delta (Correia, 2016;McElroy, 2016) due to persistent snow and ice cover in alleys and lanes.Across Metro Vancouver, many side streets and sidewalks were left un-cleared and unsalted. Residents had difficulties simply getting around thecity and questioned the Vancouver Mayor’s commitment to the issues (Laanela,2016). Nearby cities also experienced salt shortages, rationing their sup-plies. Stores and wholesalers were also having trouble meeting the demand(of Vancouver, 2017).By the end of December, the city of Vancouver had spent $2.5 million onsnow and ice reduction, triple the amount used in the previous two winterperiods combined (of Vancouver, 2017). Even after the unusual cold winter,Vancouver still needed to address the nearly 15,000 potholes (almost doublethan normal) affecting drivers (Vancouver, 2017).91.1. Understanding Extreme Weather Events and its Impacts — A Case StudyFigure 1.5: Fifteen-day averaged minimum 2-m temperature anomaly cen-tered on December 11th, 2016, relative to 1979-2017 period.101.1.UnderstandingExtremeWeatherEventsanditsImpacts—ACaseStudyFigure 1.6: Daily minimum 2-m temperature time series for Vancouver, smoothed using a 5-day rolling windowfor readability. Time series for all years in grey, 1937/38-2016/17 average in blue, 2012/13-2016/17 average in red,1992/92 in orange and 2016/17 in black. Days with snowfall in 2016/17 indicated with green dashed lines.111.1. Understanding Extreme Weather Events and its Impacts — A Case StudyLong waits, cancellations and delays also left commuters questioningVancouver’s public transportation system readiness for such cold winterweather (McElroy, 2017a) (Fig. 1.7a). Additionally, strong arctic outflowwinds over the Salish Sea during the outbreaks led B.C. Ferries to cancelsailings, impacting ferry commuters (Brend, 2016).During the arctic outbreaks, approximately $5 million was spent to op-erate an ice-clearing cable collar system on the Port Mann Bridge —theprovince’s primary east-west corridor for both commercial and commutertraffic. By contrast, just $300,000 was spent to operate it in 2015/16. Laneclosure due to crews clearing the bridge and Highway 1 led to major trafficproblems during the last arctic outbreak on February 3rd, 2017 (Saltman,2017) (Fig. 1.7b).During the third arctic outbreak on January 3rd, 2017, BC Hydro seta new record for power consumption, breaking the old record set on Nov29th, 2006 (Paetkau, 2017). Finally, in the last outbreak in early February,a storm cycle brought record-breaking snow and freezing rain to the FraserValley. Abbotsford observed 57.8 cm of snow and there were reports of up to80 cm in Chilliwack (MacMahon, 2017). Freezing rain accumulated 2-4 cm(Meuse, 2017). 361,000 BC Hydro customers lost power (Luymes, 2017) andVancouver was cut off from the rest of BC with highways 1, 3, 5 and 99 allclosed (McElroy, 2017b). The large-scale features providing the atmosphericcontext of these intense arctic outbreaks will be described next.121.1.UnderstandingExtremeWeatherEventsanditsImpacts—ACaseStudyFigure 1.7: (a) Public transit buses paralyzed by road conditions during the December 11th, 2016 arctic outbreak(Fritz, 2016). (b) A Google Maps screen capture, where orange and red colours indicate widespread very slow orstopped traffic during the February 3rd, 2017 arctic outbreak.131.1. Understanding Extreme Weather Events and its Impacts — A Case Study1.1.3 Anatomy of an Arctic OutbreakThe five arctic outbreak events were each comprised of a cold continentalairmass associated with a strong arctic high pressure system in Alaska andthe Yukon, combined with a warmer, moist air mass associated with a Pacificlow pressure system. The lows played an important role in producing periodsof snow across the South Coast.In far northern latitudes, during winter the sun angle is low or below thehorizon. This means incoming solar radiation is very limited, while outgoingradiation from the earth’s surface continues unabated. This creates a nega-tive surface energy budget and the air cools, building very cold airmasses atthe surface. Such airmasses are associated with strong high pressure such asthe ones seen over Alaska and the Yukon during the first outbreak betweenDecember 4th, 2016 and December 6th, 2016 (blue shading in Fig. 1.8a,and red sea level pressure contours in Fig. 1.8b).The cold, dense, stable airmass flows through valleys, fjords, and straitson its way to Vancouver, partially blocked by the higher mountainous terrainof British Columbia (Figs. 1.8c and d). Where the valleys, fjords andstraits widen, the cold air spreads and thins, accelerating into arctic outflowwinds (Jackson, 1996). An example is shown in a fine-resolution computermodel forecast for the second arctic outbreak on December 12th (Fig. 1.1).These outflow winds have impacts of their own, like the cancelled BC Ferriessailings in the Strait of Georgia mentioned in the previous section.Often an arctic outbreak event will draw to a close with the approach ofa Pacific low pressure system ushering in moist, mild air. In the December5th outbreak, an upper-level low over the Gulf of Alaska moved southwarddown the coast (Figs. 1.8a and c). With low-level arctic air in place, theupper-level trough and associated surface low pressure brought 5-15 cm ofsnow across the South Coast, with 5.4 cm measured in Vancouver (Figs.1.8c and d). While this is relatively little snow compared to snowfalls inother parts of Canada, it caused major disruptions.1.1.4 Final ConsiderationsThe fall season was abnormally wet for multiple months, especially in termsof frequency of rain. This culminated in an early November storm cyclethat caused headaches for emergency management personnel, and featuredgreater than 1-in-100-year cumulative flows for parts of Vancouver Island. Aseries of five arctic outbreaks led to a well-below-normal winter, and an ab-normally large number of days with below-freezing minimum temperatures.141.1. Understanding Extreme Weather Events and its Impacts — A Case StudyFigure 1.8: (a) 50.0-kPa geopotential height (km) and 2-m temperature (◦C),(b) sea level pressure (kPa) and 6-hr precipitation (mm) for December 4th,2016 at 1200 UTC. (c) and (d), as in (a) and (b), but for December 6th, 2016at 1200 UTC. Data from the the ECMWF interim reanalysis (ERA-Interim)(Dee et al., 2011).151.2. Extreme Data AnalysisThe result was a 6-month period of weather that, by the metrics discussedherein, was the worst in recent memory, and among the worst in history forthe South Coast of BC. Chapter 4 gives more details about how extremethese seasons were, and how well they were forecast.1.2 Extreme Data AnalysisOver the last few decades, techniques that can provide forecast users withadvanced warning of extreme weather events have been developed. To avoiddata waste, these techniques work by comparing model forecasts of vari-ables to historical means and climate distributions (Lalaurette, 2003; Dutraet al., 2013). For instance, the extreme forecast index (EFI) (Lalaurette,2003; Zsoter, 2006) developed at the European Centre for Medium-RangeWeather Forecasts (ECMWF) indicates how far the forecast distribution foran event deviates from the climate distribution. Another example is the sit-uational awareness table developed by National Oceanic and AtmosphericAdministration and the National Weather Service (NOAA/NWS) (Grummand Hart, 2001; Graham et al., 2013). It provides standardized anomalies,percentiles, and return periods of forecasts by comparing them to a reanal-ysis or to a so-called “model climate”.Such techniques provide forecasters and users with more information onwhich to base their decisions by making extreme events easier to identify andplacing them in historical context. However, there is room for improvementbecause they do not utilize extreme data analysis methods.Extreme levels can be estimated by fitting an appropriate distributionto a sample of observed extreme values recorded over a specified period.According to the extremal types theorem, in the limit as the sample size ofobserved extreme values increases, the distribution stabilizes and approachesone of the three extreme value distributions (EVD), regardless of the distri-bution from which the original extreme values were drawn from (Gnedenko,1943; De Haan, 1976). The three EVD distributions can be simplified to oneGeneralized Extreme Value (GEV) distribution function (Jenkinson, 1955).This theorem is analogous to the central limit theorem, which states that inthe limit as the sample size of mean values becomes large, the shape of thesampling distribution converges to a Gaussian distribution. Both theoremssuggest extreme and mean values behave differently and as such should beanalyzed separately.The above explanation might seem unnecessary at first, but to make thepoint clear, it implies an appropriate statistical method of extremes needs161.3. A Historical Recordto be employed, rather than defining an extreme value as being a certaindistance away from mean values in order to avoid data waste. Second, theconvergence behaviour of the EVD is quite important. It implies that a longhistorical record is needed for extreme data analysis to be correctly employedand justify data waste. Therefore, to lay the ground work for an improvedstatistical model of extreme weather events, first the best historical datasetneeds to be identified to accurately place extreme T2M and PCP into ahistorical context. The challenge is a lack of long-term complete historicalsurface weather station records at many locations across the province (Karlet al., 1993; Odon et al., 2018).A new methodology to address extreme weather events is introduced inChapter 2 and used thereafter throughout the entire dissertation.1.3 A Historical RecordWeather station coverage across BC is limited outside of the southwest BCpopulation centres (Karl et al., 1993; Odon et al., 2018), and often the qualityand completeness of station records degrades as well. This data paucitymotivates the search for the best model climate from a gridded reanalysisdataset. This would provide a long, continuous record with complete spatialcoverage over BC.A reanalysis is a 3-D gridded dataset created by combining rerun modelforecasts with observed weather data via data assimilation. The result is along-term gridded dataset that is spatially complete and physically consis-tent. They contain a large variety of atmospheric variables including manythat are not directly observed (Dee et al., 2011), potentially rendering afeasible solution to the paucity of observational data over BC.Some of the most well-known reanalysis datasets are the Climate Fore-cast System Reanalysis (CFSR) from the National Centers for EnvironmentPrediction (NCEP), the ECMWF interim reanalysis (ERA-Interim), theJapanese 55-year Reanalysis (JRA-55) from the Japanese MeteorologicalAgency (JMA), and the Modern Era Retrospective-Analysis for Researchand Applications (MERRA-2) from the National Aeronautics and Space Ad-ministration (NASA). They represent the current generation of reanalyseswith improvements over previous generation reanalyses R1 (NCEP/NCARReanalysis I) and R2 (NCEP/DOE Reanalysis II), ERA-15 ( ECMWF 15-year Re-Analysis) and ERA-40 ( ECMWF 40-year Re-Analysis), JRA-25(Japanese 25-year Reanalysis), and MERRA (Modern Era Retrospective-Analysis for Research and Applications) (Dee et al., 2011; Saha et al., 2010).171.4. A Gridded Climatological DatasetEven though the observational data assimilated by the four reanalysesare largely identical, their assimilation methods are not, resulting in signifi-cant differences (Mooney et al., 2011; Saha et al., 2010; Koster et al., 2016).Such differences can have an impact on how well extreme weather events arerepresented.New techniques to assess the performance of the new generation reanal-yses — CFSR, ERA-Interim, JRA-55 and MERRA-2 — with respect todaily and extreme maximum and minimum T2M, and daily and extremePCP over mountainous BC are described in Chapters 2 and 3 respectively.1.4 A Gridded Climatological DatasetReanalyses are often downscaled to a higher resolution (Cosgrove, 2003;Juang and Kanamitsu, 1994; Rasmussen et al., 2011; Stefanova et al., 2012;Abatzoglou, 2013). High-resolution gridded temperature and precipitationdatasets are widely used because of the benefits of spatial completeness andfine-scale features due to topography. Such details are important for mod-elling applications in fields like hydrology, ecology, and agriculture (Thorn-ton et al., 1997; Mote et al., 2005; Abatzoglou, 2013; Stoklosa et al., 2015).Across BC there are two types of state-of-the-art datasets suited forlong-term analysis: the very high-resolution Parameter-Elevation Regres-sions on Independent Slopes Model climatology (PRISM; Daly et al. (1994,1997, 2002)), which provides access to a 30-arc-second (∼ 800 m) griddedmonthly-mean and annual-mean maximum and minimum temperature, andprecipitation values for the 1981-2010 period (PCIC et al., 2014); and theCanadian homogenized stations dataset, which accounts for non-climaticchanges in the data such as changes in station siting, instrumentation, timeof observation and procedure (Vincent, 1998; Vincent and Gullett, 1999;Vincent et al., 2002; Mekis, 2005; Mekis and Brown, 2010; Mekis and Vin-cent, 2011). However, PRISM lacks temporal resolution and the stationslack spatial completeness.The temporal resolution of the reanalysis, the spatial resolution of PRISM,and the homogeneity and ground truth of the stations are therefore combinedto create a new, very-high-resolution surface analyses (hereafter referred toas the VHRSA). The methodology is described in Chapter 4.181.5. Forecasting Extreme Weather Events1.5 Forecasting Extreme Weather EventsMountain ranges can dramatically impact precipitation distribution due toorographic lifting and blocking, among other processes. Similarly, tempera-ture exhibits large spatial variability in complex terrain. This is in part dueto typical changes in temperature with height (i.e., the lapse rate), but alsodue to mountain-specific processes like cold air pooling, heating and cool-ing of near-slope air, diabatic processes involved in cross-barrier airmasstransformation (Smith et al., 2003), and blocking of temperature advection,among others. All these processes make the spatial distribution of precip-itation and temperature quite complicated in the complex terrain of BC(Deng et al., 2005; Astsatryan et al., 2015). Lack of model resolution to re-solve terrain features and processes make accurate forecasts in mountainousregions particularly challenging (Junker et al., 1992; Kunz and Kottmeier,2006; Smith et al., 2010; Haren et al., 2015).Although numerical weather prediction models are far less frequentlyverified over mountainous regions, several studies investigated whether high-resolution models performed better than their lower-resolution counterpartsacross complex terrain. ?, ? and ? suggest that at higher resolutions(< 4 km) there are still major problems associated with precipitation andtemperature forecasts. The model microphysis parameterizations lead tooverprediction along the steep windward slopes and underprediction in thelee of major barriers. In contrast, Schirmer and Jamieson (2015), Weusthoffet al. (2010), Garvert et al. (2005) and Ikeda et al. (2010) concluded thathigher-resolution models were equal or better at simulating orographic in-fluences compared to the low-resolution models. However, all studies agreethat higher resolutions models better resolve terrain features.The forecast methodology described in Chapter 4 both bias corrects anddownscales the North American Ensemble Forecast System (NAEFS) usingthe VHRSA dataset.1.6 Dissertation OrganizationThe outline of the dissertation is as follows. Chapters 2 and 3 evaluatereanalyses to identify the one that yields the best long-term and completeT2M and PCP historical dataset over BC, respectively. A new methodologyis also introduced in both chapters to address how extremes in T2M andPCP vary across space and time (around the calendar year). Chapter 4introduces a new statistical model that is used to forecast extreme T2M and191.6. Dissertation OrganizationPCP, that accounts for the issues of model resolution and non-stationarity(changes due to climate change, urbanization, and other effects). The resultsare summarized in the conclusion.20Chapter 2Performance of Reanalysesacross British Columbia.Part I: Evaluation of Dailyand Extreme 2-mTemperature2.1 IntroductionPerformance of the new generation reanalyses — CFSR, ERA-Interim, JRA-55 and MERRA-2 — is assessed with respect to daily and extreme maximumand minimum 2-m temperature (hereafter T2M) over mountainous BC. Thepurpose of this Chapter is to identify which reanalysis best represents ob-servations as a surrogate climatological dataset. The end goal is for thisreanalysis to serve as a tool for the creation of a new extreme forecast indexin Chapter 4.Previous studies have shown differences in T2M between the ERA-Interimand CFSR, and previous generation reanalyses ERA-40 and R1 over the Ti-betan Plateau in central Asia (Bao and Zhang, 2013). Betts et al. (2009)showed that the seasonal cycle of mean T2M is higher than observationsover the Mississippi River basin for both ERA-40 and ERA-Interim. Berget al. (2003) found that ERA-40 has substantial biases in T2M and dewpointtemperatures over land in North America. Wang et al. (2011) evaluated themean T2M performance over the entire globe from a mix of new and pre-vious generation reanalyses (ERA-40, R1, R2 and CFSR). More recently,Lindsay et al. (2014) evaluated the performance of mean T2M over the Arc-tic for CFSR, MERRA-1, ERA-Interim, and JRA-25. All studies agree thatnewer reanalyses outperform previous generation reanalyses.Both daily and extreme maximum and minimum T2M are analyzed,(defined in the Data and Methodology section), and compared with weather212.2. Data and Methodologystations in BC for the period 1980-2010. T2M is examined because of itsimportance and impact on ecosystems (Parmesan et al., 2000), agriculture(Rosenzweig et al., 2001), tourism (Patz et al., 2005; White et al., 2016),urbanization and health (Curriero et al., 2001; White et al., 2016; Odonet al., 2017), and due to observational data availability (Jones et al., 1985;Peterson and Vose, 1997; Jones et al., 1999).Second, trends in both daily and extreme T2M during the study pe-riod are examined in order to determine if significant statistical changesoccurred over time. This is to assess whether or not a stationary climato-logical distribution is appropriate to represent present-day daily and extremedistributions.The outline of the Chapter is as follows. In section 2.2, a brief de-scription of the different reanalyses and of the weather-station observationsis given. In sections 2.3- 2.5 a methodology for dividing BC into climatezones, the various metrics used for evaluating daily and extreme reanalysisT2M, and a method for assessing statistical nonstationarity are described.In section 2.6, the daily and extreme T2M from the reanalyses are evaluated.In section 2.7, the trends of both daily and extreme T2M are examined, andthe corresponding changes in return levels and return periods of extremesare discussed. The results are summarized in section 2.8.2.2 Data and MethodologyDaily maximum/minimum T2M from 57 weather stations from 1 Jan 1980 to31 Dec 2010 are used in this study for evaluation of the CFSR, ERA-Interim,JRA-55 and MERRA-2 reanalyses. A description of the weather stationdataset is given below. The authors chose the 1980-2010 as the overlappingperiod for comparison because the MERRA-2 began in 1980 (Gelaro, 2015);and because in 2011 the CFSR was extended using NCEP’s Climate ForecastSystem Version 2 (CFSv2) operational model, and differences between themodels used to produce CFSR and the operational CFSv2 may affect dataevaluation across the 2010/2011 boundary period (Saha et al., 2014).A description of the different reanalyses and of the weather stationdataset is given below. A summary of the reanalyses atmospheric modelsand configurations is presented in Table 2.1.T2M from the weather stations used in this study are not assimilated bythe CFSR Wang et al. (2011), or the MERRA-2 (Bosilovich et al., 2015).The SYNOP stations (surface synoptic observations) are directly assimi-lated by the ERA-Interim (Dee et al., 2011) and the JRA-55 (Kobayashi222.2. Data and Methodologyet al., 2015) (see Appendix A for assimilated stations). Therefore evaluat-ing against these observations provides a dependent measure of accuracy. Adiscussion of how the reanalyses perform across non-assimilated stations isgiven in section 2.6.232.2.DataandMethodologyTable 2.1: Overview of the four reanalysis datasets examined in this study.Institution Reanalysis Model AssimilationmethodPeriod Downloadgrid(lat× lon)TimeintervalReferenceNCEP/NCAR CFSR CFST382/L64(globalhorizontalresolution∼ 38 km)3D-Var GSI 1979-2011(current asCFSv2)0.5◦ × 0.5◦(∼ 50 km)0000,0600,1200and 1800UTCSaha et al.(2010)ECMWF ERA-InterimIFST255/L60(globalhorizontalresolution∼ 79 km)4D-Var 1979-current 0.5◦ × 0.5◦(∼ 50 km)0000,0600,1200and 1800UTCDee et al.(2011)JMA JRA-55 JMAT319/L60(globalhorizontalresolution∼ 55 km)4D-Var 1958-2012(current asJCDAS)0.5616◦ ×0.5616◦(∼ 55 km)0000,0600,1200and 1800UTCEbita et al.(2009)NASA MERRA-2GEOS-5.12.4AGCM(lat × lon)0.5◦ ×0.625◦/L723D-Var GSI 1980-current 0.5◦ ×0.625◦0030,0130,......2330UTCGelaro(2015);Gelaroet al.(2017)242.2. Data and Methodology2.2.1 Weather station dataThe stations were initially selected because of their proximity to populationcentres and so that they would be geographically dispersed around BC. Theyrepresent conditions for a mixture of coastal, intermountain and continentalclimates. Ideally, variations in climatological time series should be causedonly by changes in weather and climate. However, decades-long time seriescan be affected by inhomogeneities such as missing data, changes in instru-mentation, station relocation, changes in the local environment surroundingthe weather station such as urbanization or land use cover, and changes intime of observation, to name a few (Peterson et al., 1998; Mekis and Hogg,1999; Jones et al., 1985). Some changes cause sudden discontinuities whileother changes, such as changes in the local environment around the station,cause gradual trends in the data.Of the 74 stations initially selected for the study, 8 were eliminated dueto either station relocation or discontinuation during the 1980 to 2010 studyperiod. Additionally, eight stations with more than 4% missing data wereexcluded, and one station was eliminated due to poor data quality. Ofthe remaining 57 stations, 48 are from Environment and Climate ChangeCanada (ECCC) and nine are from BC Hydro, the primary electric utilitycompany for BC.Figure 2.1 shows the locations of all 57 stations overlaid with the popu-lation distribution across BC, and with the province’s political boundaries.Red and blue stations represent ECCC and BC Hydro stations, respectively,and their corresponding three-letter abbreviations that are referenced inthe paper. Fifty-two of the stations are located in valleys (indicated byupside-down triangles), and five are in non-valley locations (upright trian-gles). Of the 4.8 million people living in BC, about half live in the Vancouvermetropolitan area in southwest BC.252.2.DataandMethodologyFigure 2.1: Location of ECCC (red) and BC Hydro (blue) weather stations and British Columbia populationdistribution (orange). Upside-down triangles indicate valley stations; upright triangles indicate non-valley stations.Dashed lines delineate the dominant climate zones: North, Central, Maritime West, Southwest and Southeast.262.2. Data and MethodologyTemporal data homogeneity for ECCC stations was assessed by iden-tifying non-climatic shifts in the annual and monthly means of the dailymaximum and minimum T2M using techniques based on regression models(Vincent, 1998; Wang et al., 2007). The non-climatic shifts were mainlydue to station relocation, or changes in observing practices and automation(Vincent and Gullett, 1999). Monthly and daily maximum and minimumtemperatures were then adjusted for the shifts identified in the temperatureannual series (Vincent et al., 2002).BC Hydro station data were manually quality controlled based on rangelimits, spatiotemporal consistency, present weather conditions, and are usedto verify reanalysis performance at non-valley stations. Initial stationar-ity analysis indicated potentially spurious trends in the data, so they areexcluded from the analysis of nonstationarity in daily and extreme T2M.Furthermore, ECCC stations with more than 1% missing data are also ex-cluded, leaving 26 ECCC stations for the nonstationarity analysis.2.2.2 CFSRIn 2010, NCEP introduced the CFSR. Previous NCEP reanalyses have beenamong the most used NCEP products in history. Many known errors inthe assimilation of observational data and execution of previous reanalyseswere corrected in the CFSR, resulting in a superior product in most respects(Saha et al., 2010).CFSR uses NCEP’s global coupled atmosphere-ocean model. It consistsof a spectral triangular atmospheric grid (Saha et al., 2006) at a horizontalresolution of T382 (∼ 38 km) and a hybrid sigma-pressure system with 64vertical levels extending from the surface to approximately 0.26 hPa.CFSR was the first NCEP global reanalysis to directly assimilate satelliteradiances, and to use three-dimensional variational data assimilation (3D-Var) in a Gridpoint Statistical Interpolation (GSI) scheme rather than aSpectral Statistical Interpolation (SSI) scheme.The T2M variable in the CFSR is derived primarily from satellite radi-ances and radiosonde information. Observations of T2M from land stationsare not assimilated (Saha et al., 2010).2.2.3 ERA-InterimThe ECMWF introduced the ERA-Interim in 2011, in part to replace ERA-40 (Dee et al., 2011). The ERA-Interim configuration uses a spectral T255272.2. Data and Methodologyhorizontal resolution (∼ 79 km), and a hybrid sigma-pressure system with60 vertical levels with the top of the atmosphere located at 0.1 hPa.Some of the improvements over ERA-15 and ERA-40 include a newhumidity analysis; new model physics where the prognostic equations aresolved using a semi-Lagrangian scheme, a variational bias correction of satel-lite radiances (Dee et al., 2011); and improvements on various technical as-pects of reanalysis such as data selection, quality control, bias correction,and performance monitoring; each of which can have a major impact onthe quality of the reanalysis. Additionally, the use of 4D-Var for the at-mospheric analysis in ERA-Interim is a major step forward. The improved4D-Var data assimilation scheme allows the cost function to be minimizedover a 12-h assimilation time interval rather than a single time, as is thecase for ERA-40 and CFSR.The T2M variable is directly assimilated by the ERA-Interim from SYNOPstations (Dee et al., 2011).2.2.4 JRA-55In 2011 the JMA produced the Japanese 55-year Reanalysis (JRA-55) (Ebitaet al., 2009). The JRA-55 extends 55 years from 1958 to 2012, and willbe continued in real time as the JMA Climate Data Assimilation System(JCDAS). JRA-55 uses a spectral model integrated at a TL319 (∼ 55 km)horizontal resolution with 60 vertical levels up to 0.1 hPa in hybrid sigma-pressure coordinates.It employs a 4D-Var scheme which seeks the initial condition that bestfits the forecast to the observations within a 6-h assimilation interval ratherthan a single time. The reanalysis also contains a new radiation scheme,variational bias correction for satellite radiances, an update on dynamicaland physical processes such as the prognostic equations being solved in asemi-Lagrangian form rather than Eulerian (Ebita et al., 2011; Takeuchiet al., 2013). These upgrades significantly reduce model biases versus theJMA’s previous generation reanalysis, enhance the dynamical consistencyof analysis fields, and advance the handling of satellite radiances.Similarly to ERA-Interim, T2M is directly assimilated from SYNOPstations (Ebita et al., 2011).2.2.5 MERRA-2The new Modern Era Retrospective-Analysis for Research and Applications(MERRA-2) was produced by NASA in 2015. The grid used for MERRA-2282.3. Climate zonesis 0.5◦ latitude 0.625◦ longitude (∼ 55 km) with 72 vertical levels in hybridsigma-pressure coordinates, from the surface to 0.01 hPa.It uses an upgraded 3D-Var assimilation scheme based on the GSI witha 6-h update cycle. Some other improvements over NASA’s previous gen-eration reanalysis include an updated physics model, aerosol assimilation,corrections in precipitation for land surface and imbalances in water andenergy cycles, to name a few (Rienecker et al., 2011; Gelaro, 2015; Gelaroet al., 2017).Similarly to CFSR, the T2M is derived primarily from satellite radi-ances, and radiosonde, aircraft and wind-profiler information. Observationsof T2M from land stations are not assimilated (Gelaro et al., 2017).2.3 Climate zonesThe location of BC immediately east of the Pacific Ocean, and its com-plex topography, produce distinctive climate zones that vary with elevation,location relative to ocean and mountains, and latitude (Peel et al., 2007).Climate zones vary from the wettest in Canada to the hottest and driest,with the formidable barrier of the Coast Mountains serving as the dividing-line between the two.Due to temperature ranges and variability that are unique to each cli-mate zone (Moore et al., 2008), reanalysis performance is evaluated sepa-rately for each zone.To determine the different climate zones, Principal Component Analysis(PCA) is conducted on the daily maximum and minimum T2M correlationmatrix for each station. The corresponding station correlation is calculatedfrom the daily temperature values for the entire study period. The highcorrelation values (Fig. 2.2) are owed to annual seasonal effects, which insection 2.5 becomes an important factor as the authors want to evaluate thedegree to which the reanalyses capture the annual cycle across the differentclimate zones.The Principal Components are ordered by the size of their respectiveeigenvalues — their rank corresponds to their relative importance in de-scribing temperature variations. 93% of the variability of both maximumand minimum temperature can be explained by the first four components,which are retained, leading to a significant reduction of data while retainingmost of the variance.A K-means clustering analysis is then performed on the four componentsto find a natural grouping of the data, where each component belongs to292.3. Climate zonesFigure 2.2: Correlation matrix and five dominant climate zone clusters inBC.302.3. Climate zonesthe cluster with the nearest centroid. One to 10 clusters were tested and afive-cluster solution was chosen (Fig. 2.2). After five clusters, the additionaldecrease in the within-cluster dissimilarity is not substantial, and the num-ber of stations per cluster becomes too low. The five climate zones (encircledby dashed lines in Fig. 2.1) are subsequently matched to those identified byChilton (1981) and Moore et al. (2008):1. Islands and Coast Mountains (hereafter Maritime):Heavily influenced by the Pacific Ocean, this zone includes island sta-tions on the immediate coast, windward of the Coast Mountains. Ithas mild winters and cool summers. Fall, winter, and spring featurefrequent landfalling low pressure systems, while summers are fairlydry. High elevations accumulate very deep snowpacks in the winterand are extensively glaciated.2. Coast Mountains (hereafter Southwest):The climate of the Coast Mountains is similar to the Maritime cli-mate, with a slightly less pronounced maritime influence. This moreprotected region has slightly more sunshine and greater temperaturevariability relative to the Maritime climate zone.3. Interior Plateau (hereafter Central):On the lee side of the Coast Mountains, the Interior Plateau is abroad, elevated region broken occasionally by narrow valleys. In thisdrier and more continental climate, seasonal and diurnal differences intemperature are much greater than at the coast. Summers tend to behot and dry; winters cooler and less moist. The Okanagan Valley isthe southernmost, hottest, and driest part of the province and in allof Canada.4. Columbia Mountains and Southern Rockies (hereafter Southeast):Farther east, westerly winds again ascend the Columbia and RockyMountains. These ranges also restrict easterly flow of cold, continen-tal arctic air into the region. High elevations are wet and cool. Valleybottoms are semi-arid with hot summers, and frequent, cold temper-atures below temperature inversions in the cool season.5. The Interior Plains, Northern and Central Plateaus and Mountains(hereafter North):312.4. Verification MetricsIn the Northern and Central Plateaus and Mountains, winters arecolder and drier than in the South, due to frequent influxes of conti-nental arctic air, and fewer Pacific storms. Summers are short, warm,and wetter than southern BC due to the northward migration of thestorm track.To the east of the northern Rocky Mountains lies an extension of theGreat Plains. It sees limited precipitation due to the many topographicbarriers between it and the Pacific Ocean, and very frequent influxesof continental arctic air due to proximity and a lack of topographicbarriers to the north and east. This area experiences a continentalclimate with long, cold, dry winters and short, warm summers.2.4 Verification MetricsThe statistical behaviors of daily and extreme T2M are compared betweenobserved weather station data, and their corresponding location in the re-analyses. In order to evaluate the agreement between observations at sta-tion locations and reanalyses at grid points, the methods Nearest Neighbor,Inverse Distance Weighting (IDW), Bilinear and Bicubic interpolation areinvestigated for comparison with the observed values. No method is supe-rior across all locations and IDW is adopted. A more detailed description ofeach interpolation method can be found in Mooney et al. (2011) and Stahlet al. (2006).2.4.1 Daily T2MTo obtain a smooth climatology of daily maximum and minimum T2M av-eraged over 31 years for each calendar day, a 31-day centered rolling windowwas applied to both reanalysis and observed temperatures. The rolling av-erage is the unweighted mean of the 961 sample values centered around acalendar day. This window length was chosen because it smooths out noiserelated to extreme weather events while still capturing monthly variations indaily T2M. For the reanalyses, the daily maximum (minimum) T2M is de-fined as the highest (lowest) value of the six-hourly T2M outputs for CFSR,ERA-Interim and JRA-55; and as the highest (lowest) value of the hourlyT2M outputs for MERRA-2 in a calendar day. The Canadian meteorologi-cal convention is followed, where it defines a calendar day to be from 0601UTC of the current day to 0600 UTC of the following day (MeteorologicalService of Canada, 2015) .322.4. Verification MetricsThe bias (or systematic error), the random error, and the two-sampleKolmogorov-Smirnov (KS) statistic are then computed to estimate how ac-curately each reanalysis captures T2M.Bias is the average difference between the reanalysis and the observa-tion. The random error, measured by centered root mean squared error(CRMSE), expresses how concentrated the errors are from the mean biason a given calendar day. The two-sample KS statistic is used to determinethe largest absolute difference between the reanalysis empirical cumulativedistribution function (ECDF) and the observation ECDF. It detects thegeneral difference between the two distributions.2.4.2 Extreme T2MIn addition to societal impacts, extreme temperatures are a concern forpower utilities primarily because of their effect on load. That is, duringperiods of extreme heat, electrical load (demand) is high due to air condi-tioning. The same goes for extreme cold and heating. Other impacts includetransmission line thermal ratings — wherein power transmission can be lim-ited during periods of hot weather, limiting a utility’s ability to supply powerfor consumption or market trading. It is important for utilities to identifyand anticipate such events, so that they can plan optimal system operationsin light of abnormal load considerations.A return level is a value of T2M that is expected to occur, on averageonce every return period. For extreme T2M, the 2- and 30-yr return lev-els are examined, which represent significant and extreme departures fromnormal, respectively (2- and 30-year return levels are also known as 2- and30-year recurrence intervals, or 0.5 and 0.03 Annual Exceedance Probabili-ties (AEP)).To do this, a 31-day centered rolling window is used over the 31-yearperiod to obtain 31 values of annual maximum (minimum) T2M for eachcalendar day, one value for each of the 31 years. In other words, the 31 an-nual maximum (minimum) T2M in a given 31-day centered rolling windoware the highest (lowest) values in each of the 31 years, which then com-prise 31 values of extreme T2M for each calendar day. As with daily T2M,the window length struck a balance between having a large enough samplesize while ensuring that each day within the window would have a similarclimatological distribution.Return levels can be estimated by fitting an appropriate distribution to asample of extreme values. According to the theorem of extremal types, in thelimit as the sample size of extreme values becomes large, the distribution sta-332.4. Verification Metricsbilizes and approaches one of the three extreme value distributions (EVD),regardless of the distribution from which the original T2M was drawn (Gne-denko, 1943; De Haan, 1976). This theorem for extreme values is analogousto the central limit theorem for mean values, which states that in the limitas the sample size of mean values becomes large, the sampling distributionconverges to a Gaussian distribution. Both theorems suggest extreme andmean values behave differently and as such should be analyzed differently.The three EVD distributions can be simplified to one Generalized ExtremeValue (GEV) distribution function (Jenkinson, 1955).A GEV dresses these 31 sample values for each calendar day by themethod of L moments. The method of L moments is chosen because itproduces more efficient estimates of return levels for small sample sizes thanthe method of maximum likelihood (Hosking et al., 1985; Hosking, 1990).The small sample size of 31 extremes values is a concern because: 1) the threeEVD distributions converge at different rates (Davis, 1982), and 2) the dailyT2M from which the extreme values are sampled are normally distributedand serially correlated, resulting in a slower convergence rate to either of thethree distributions (Leadbetter et al., 1983). Thus, a goodness-of-fit test isconducted to evaluate the closeness of the fitted GEV distributions to thedata.A Lilliefors test compares the largest absolute difference between thefitted GEV cumulative distribution function (CDF) and the observationECDF. The null hypothesis is that the observed data is drawn from a GEVdistribution. A sufficiently large critical difference between the fitted GEVCDF and the observation ECDF results in the null hypothesis being rejected.A parametric bootstrap procedure determines this critical difference.Namely, 100 samples of size 31 are generated from the fitted GEV distribu-tion for each calendar day at each station. Then, 100 critical differences arederived from the comparison of each generated sample ECDF and the fittedGEV CDF. The 90th percentile of the resulting collections of critical differ-ences is used as the critical difference for the rejection of the null hypothesisthat the sample originates from a GEV distribution, which corresponds toaccepting or rejecting the null hypothesis at the α0 = 0.10 significance level.At this 10% significance level, 10% of the tests are expected to exceeda critical difference, rejecting the null hypothesis. However, there are 57weather stations and 365 calendar days totaling N0 = 20805 independentLilliefors tests. To reduce the probability of incorrectly rejecting one ormore of the N0 true null hypotheses, Walker’s criterion with αWalker =1 − (1 − α0)1/N0 = 5.06 × 10−6 is regarded as significant. Namely, the nullhypothesis is rejected when the critical difference between the fitted GEV342.5. NonstationarityCDF and the observations ECDF exceeds the (1 − αWalker)th percentile ofthe resulting collections of critical differences between the fitted GEV CDFand the generated ECDF of the samples (Wilks, 2016).Less than 4% (6%) of the locations and calendar days are rejected duringthe 1980-2010 study period, suggesting maximum (minimum) extreme T2Mbehaviour can be described by a GEV distribution.The 2- and 30-year return levels are then estimated from the GEV dis-tribution for each calendar day. The bias of the 2- and 30-year return levelsare calculated to estimate how well each reanalysis captures extreme T2M.The 2- and 30-yr return levels of maximum and minimum extreme T2Mare chosen because less than 1% of their 90% respective confidence intervalsoverlap, indicating the difference between the two return levels is statisticallysignificant. The sampling uncertainty of the estimates of the 2- and 30-yr return levels are determined by a parametric bootstrap procedure. 100samples of size 31 are generated from the fitted GEV distribution for eachcalendar day at each station. Then, the 2- and 30-year return levels areestimated from each generated sample. The 5th and 95th percentiles of theresulting collection of 2- and 30-year return levels is used as lower and upperbounds of the 90% confidence intervals for the true 2- and 30-year returnlevels.2.4.3 Analysis of Variance (ANOVA)The mean systematic error of daily maximum and minimum T2M, and of 2-and 30-year return levels are calculated for each location from all calendarday systematic errors. Comparisons between mean systematic errors of thereanalyses (CFSR, ERA-Interim, JRA-55 and MERRA-2) and of the ob-served levels of T2M (daily maximum, daily minimum, extreme maximumand extreme minimum) are made using a two-way ANOVA. Tukey’s honestsignificant difference (HSD) is applied following statistical significance in theANOVA to identify differences in pairwise comparisons of reanalyses meansystematic errors.2.5 NonstationarityStatistical nonstationarity is tested to determine if there are important tem-poral changes in daily and extreme T2M. Specifically, is nonstationarity sub-stantial enough to require a more complex characterization of extreme levelsthat takes into account temporal changes, or is a comparatively simpler sta-tionary model accurate enough to represent daily and extreme T2M?352.5. Nonstationarity2.5.1 Daily T2MFor daily T2M, a 91-day centered rolling window is used to obtain 2821 val-ues of maximum and minimum T2M for each calendar day over the 31 years.The rolling window is extended to 91 days to increase the sample size whilestill accounting for seasonal variations. A Gaussian distribution dresses the2821 sample values, where the parameters of the Gaussian distribution areestimated by the method of maximum likelihood.An advantage of maximum likelihood over L moments for parameter es-timation is its adaptability to changes in model. The maximum likelihoodestimation (MLE) of nested models leads to a simple likelihood ratio test(LRT) procedure of any one nonstationary model where the mean and stan-dard deviation vary with time, against a stationary model where the meanand standard deviation are constant. Because the method of L momentsare defined for only identically distributed random samples, the use of time-dependent parameters in a Gaussian distribution prevent us from using themethod of L moments.First, a Shapiro-Wilks test tests the null hypothesis that the observeddata is drawn from a Gaussian distribution at the αFDR = 0.05 level of sig-nificance. The False Discovery Rate (FDR) procedure is similar to Walker’scriteria, but less strict. The collection of p-values from the N0 = 20805Shapiro-Wilks tests are arranged in ascending order. The null hypothesesare rejected if their respective p-values are lower than a threshold level pFDRthat depends on the distribution of ordered p-values (Wilks, 2016).For the locations and calendar days where the Shapiro-Wilks test istested during the 1980-2010 study period, more than 90% of the resultssuggest both maximum and minimum daily T2M can be described by aGaussian distribution.Second, the null hypothesis that all 31 years of daily maximum/minimumT2M are drawn from the same Gaussian distribution are tested against thealternate hypothesis that each of the 31 years of the daily temperature data isdrawn from a different Gaussian distribution due to time dependencies in theparameters. This alternate hypothesis is tested against the null hypothesisusing a LRT.To perform the LRT, it is necessary to fit 31 Gaussian distributionsseparately to the each year of data, and compare these 31 distributions withthe single Gaussian distribution fit using the full data set for each calendarday at each station. In principle, time dependence can be assumed on bothparameters. However, if the changes in either of the parameters is negligible,it is advantageous to keep them constant.362.6. Results and discussionThe Gaussian distribution parameters are assumed to vary with time asµ(t) = µ0 + µ1(t− t0), and σ(t) = σ0 + σ1(t− t0). The slope coefficients µ1and σ1 represent the annual rate of change in mean and standard deviationof daily T2M respectively. The null hypotheses µ1 = 0 and σ1 = 0 are testedagainst the alternate hypotheses that µ1 6= 0 and σ1 6= 0.With models M0 = N(µ, σ) and M1 = N(µ(t), σ(t)), the LRT comparesthe difference between the maximized log-likelihoods under models M0 andM1 respectively. A sufficiently large difference indicates that the nonstation-ary model explains substantially more of the variation in the data than thestationary model. Small differences suggest that the increase in model sizedoes not bring worthwhile improvements in the capacity of the model to ex-plain the data. This critical difference is determined by the χ2 distribution.Model M0 is rejected at the αFDR = 0.05 level of significance.2.5.2 Extreme T2MFor extreme maximum and minimum T2M, the rolling 31-day centeredrolling window is maintained for the 31 year period. A GEV dresses these31 sample values by the method of maximum likelihood. A nonstationaryGEV distribution is compared, where only the location parameter is allowedto exhibit trend, with a stationary GEV distribution with constant location,scale and shape parameters.The GEV distribution location parameter is assumed to vary with timeas µ(t) = µ0 +µ1(t− t0), where the slope coefficient µ1 represent the annualrate of change in location values of extreme T2M.With models M1 = GEV (µ(t), σ, κ) and M0 = GEV (µ, σ, κ), the alter-nate hypothesis is tested against the null hypothesis that the extremes oftemperature are drawn from the same GEV distribution using a LRT at theαFDR = 0.05 level of significance. Namely, the alternate hypothesis is thata nonstationary model explains substantially more of the variation in thedata, and consequently changes in return levels should be accounted for.2.6 Results and discussion2.6.1 Daily T2MThe Central climate zone is representative of the overall performance of thedifferent reanalyses for daily maximum T2M across BC (Figure 2.3). Theseasonal cycle of temperature is well captured by all reanalyses (Fig. 2.3). Allfour reanalyses show a negative (cold) systematic error (Fig. 2.3; Table 2.2)372.6. Results and discussionthroughout the entire calendar year. ERA-Interim and JRA-55 outperformCFSR and, to a lesser degree, MERRA-2. Additionally, MERRA-2 has aless consistency in performance, and usually has a higher systematic errorduring spring. MERRA-2 and CFSR also tend to have a higher (worse)random error during spring (Fig. 2.3). The Maritime climate zone has thelowest systematic across all reanalyses all year long, and the lowest season-to-season temperature variation (Figure not shown; Table 2.2). The Northclimate zone has the highest random error across all reanalyses during wintertime, and the highest season-to-season temperature variation. This might beexpected to some extent since higher observed variability typically makes forpoorer model forecast random error, and the reanalyses are based on modelforecasts.382.6.ResultsanddiscussionTable 2.2: Seasonally averaged systematic error (SE in ◦C) and random error (RE in ◦C) of mean values, andKolmogorov-Smirnov statistic (KS) of daily maximum T2M by climate zoneWinter Spring Summer FallSE RE KS SE RE KS SE RE KS SE RE KSNorthCFSR -2.05 1.44 0.29 -5.59 1.12 0.53 -4.22 0.91 0.40 -3.86 0.86 0.43ERA-Interim -1.88 1.33 0.23 -3.47 0.69 0.37 -3.71 0.54 0.38 -3.21 0.79 0.36JRA-55 -1.67 1.56 0.28 -3.80 0.78 0.39 -3.20 0.53 0.34 -3.01 0.77 0.35MERRA-2 -3.32 1.52 0.30 -6.11 1.12 0.57 -4.24 0.75 0.42 -3.76 0.85 0.40CentralCFSR -2.15 1.11 0.30 -4.80 0.99 0.50 -4.36 0.81 0.40 -3.68 0.73 0.40ERA-Interim -1.55 0.97 0.18 -2.67 0.56 0.30 -3.15 0.58 0.62 -2.67 0.64 0.30JRA-55 -2.20 0.96 0.25 -2.75 0.52 0.30 -2.98 0.54 0.29 -2.88 0.64 0.31MERRA-2 -3.61 1.12 0.37 -4.98 1.01 0.57 -2.84 0.64 0.42 -2.89 0.80 0.40MaritimeCFSR -0.81 0.76 0.25 -1.40 0.79 0.39 -1.35 0.76 0.41 -1.02 0.65 0.35ERA-Interim -0.05 0.72 0.25 -0.05 0.74 0.32 -0.29 0.81 0.30 -0.27 0.67 0.31JRA-55 -0.37 0.80 0.29 -0.38 0.88 0.33 -1.03 0.98 0.36 -0.62 0.79 0.34MERRA-2 -0.56 0.73 0.24 -0.66 0.77 0.32 -0.45 0.81 0.30 -0.40 0.66 0.28SouthwestCFSR -2.84 0.75 0.39 -3.23 0.90 0.39 -2.57 0.97 0.28 -2.74 0.79 0.34ERA-Interim -2.57 0.72 0.35 -2.28 0.67 0.31 -2.16 0.74 0.27 -2.54 0.68 0.32JRA-55 -2.44 0.73 0.33 -2.13 0.76 0.29 -2.49 0.77 0.28 -2.53 0.77 0.32MERRA-2 -3.15 0.78 0.44 -3.15 0.98 0.39 -1.94 0.75 0.24 -1.98 0.76 0.26SoutheastCFSR -4.87 0.98 0.47 -7.83 1.19 0.67 -7.02 1.01 0.52 -5.91 0.86 0.51ERA-Interim -4.19 1.03 0.39 5.99 0.72 0.53 -6.28 0.71 0.47 -5.25 0.83 0.46JRA-55 -4.63 0.94 0.43 -5.52 0.73 0.51 -5.54 0.67 0.43 -4.95 0.84 0.44MERRA-2 -5.61 1.06 0.14 -7.51 1.10 0.65 -5.64 0.76 0.46 -4.58 1.04 0.42392.6.ResultsanddiscussionFigure 2.3: Observed and reanalysis mean, systematic error, and random error magnitude for daily maximumT2M, averaged over stations in the Central climate zone. The vertical dashed lines indicate the change in seasonsand the colored dots represent the seasonal average.402.6. Results and discussionDaily minimum T2M and their errors are illustrated in Figure 2.4 for theNorth climate zone, which is representative of other zones. Similar to dailymaximum T2M, all four reanalyses capture season-to-season temperaturevariations regardless of the temperature differences between the differentclimate zones, but generally exhibit a cold bias (Figure not shown; Table2.3). In contrast, the systematic error is much smaller for daily minimumthan for daily maximum T2M across all zones. Again, ERA-Interim andJRA-55 outperform CFSR and MERRA-2 for systematic and random errors,both of which are higher during winter months. For the North, systematicerror is positive (warm) during winter months for the ERA-Interim andJRA-55, and in the Maritime it is positive all year across all reanalyses(Figure not shown; Table 2.3).412.6.ResultsanddiscussionTable 2.3: Seasonally averaged systematic error (SE in ◦C) and random error (RE in ◦C) of mean values, andKolmogorov-Smirnov statistic (KS) of daily minimum T2M by climate zoneWinter Spring Summer FallSE RE KS SE RE KS SE RE KS SE RE KSNorthCFSR 0.48 1.78 0.20 -1.42 1.17 0.24 -1.51 0.66 0.28 -0.62 1.01 0.25ERA-Interim 1.13 1.37 0.19 -0.32 0.88 0.17 -0.48 0.53 0.17 -0.58 0.86 0.21JRA-55 1.93 1.71 0.25 -0.51 0.82 0.17 -1.02 0.53 0.26 -0.67 0.86 0.26MERRA-2 -3.03 1.56 0.23 -3.89 1.01 0.22 -2.24 0.66 0.37 -3.35 1.09 0.35CentralCFSR 0.84 1.56 0.14 -0.77 0.99 0.17 -1.60 0.68 0.40 -0.18 0.90 0.15ERA-Interim 0.98 1.24 0.14 -0.01 0.75 0.18 -0.17 0.61 0.62 -0.09 0.83 0.15JRA-55 1.16 1.34 0.16 -0.24 0.74 0.16 -0.64 0.63 0.29 -0.37 0.83 0.15MERRA-2 -3.54 1.49 0.24 -2.86 0.89 0.33 -1.59 0.71 0.29 -2.39 1.08 0.27MaritimeCFSR 2.40 0.87 0.35 2.10 0.66 0.37 1.11 0.63 0.29 2.08 0.80 0.37ERA-Interim 2.71 0.88 0.33 2.74 0.68 0.40 1.96 0.31 0.38 2.52 0.83 0.37JRA-55 2.79 0.93 0.39 2.89 0.70 0.47 1.72 0.67 0.37 2.67 0.83 0.43MERRA-2 1.57 0.86 0.29 1.58 0.63 0.31 0.90 0.63 0.26 1.58 0.80 0.32SouthwestCFSR -0.65 0.79 0.19 -0.94 0.68 0.26 -1.79 0.69 0.27 -0.82 0.70 0.22ERA-Interim -1.08 0.84 0.21 -0.40 0.59 0.22 -0.19 0.59 0.20 -0.66 0.65 0.21JRA-55 -0.14 0.74 0.22 0.23 0.59 0.26 -0.23 0.62 0.19 -0.05 0.69 0.26MERRA-2 -3.08 0.91 0.22 -1.94 0.69 0.34 -1.00 0.61 0.29 -1.58 0.75 0.28SoutheastCFSR -3.43 1.41 0.24 -3.89 1.16 0.41 -4.22 0.80 0.54 -3.03 1.08 0.33ERA-Interim -2.94 1.19 0.24 -2.62 0.83 0.32 -2.41 0.74 0.37 -2.57 1.05 0.29JRA-55 -2.22 1.10 0.25 -2.52 0.81 0.33 -2.32 0.81 0.35 -2.18 0.97 0.28MERRA-2 -7.83 1.30 0.55 -7.83 1.13 0.60 -3.84 0.79 0.51 -4.69 1.36 0.46422.6.ResultsanddiscussionFigure 2.4: Observed and reanalysis mean, systematic error, and random error magnitude for daily minimumT2M, averaged over stations in the North climate zone. The vertical dashed lines indicate the change in seasonsand the colored dots represent the seasonal average.432.6. Results and discussionFinally, a KS statistic is used to determine how well the reanalyses cap-tures the ECDF of observed daily maximum T2M across the Southeastclimate zone (Fig. 2.5). Values closer to zero are better, which would indi-cate no difference between the observation ECDF and that of the reanaly-sis. ERA-Interim and JRA-55 generally outperform CFSR and MERRA-2.MERRA-2 once again exhibits high variability in performance throughoutthe seasons, followed closely by CFSR (Table 2.2). MERRA-2 performsbest in the Maritime climate zone. For daily minimum T2M the dispar-ity between reanalyses is more pronounced, with Era-Interim and JRA-55consistently outperforming CFSR and MERRA-2 (Table 2.3). Furthermore,daily minimum T2M is consistently better captured by the reanalyses thanmaximum T2M.A reality of the relatively sparse surface station network in the com-plex terrain of BC is that most stations are located in valleys, and higher-elevation stations are lacking (Fig. 2.6). This is especially true for stationswith quality, long records. With this caveat in mind, the effects of re-analysis terrain error (reanalysis terrain elevation minus real-world stationelevation) are examined. Due to the relatively coarse resolution of the re-analyses, they tend to reduce terrain amplitude, exhibiting negative terrainerror over mountain ridges, and positive errors over valleys. Reanalysis dailyT2M systematic error is strongly related to reanalysis terrain error, but notclimate zone (Fig. 2.7). Namely, larger positive terrain errors are associatedwith larger negative temperature biases. That is, typically valleys exhibit acold bias in reanalyses, and a warm bias over mountain ridges.The regression line is steeper for CFSR and MERRA-2, indicating thatERA-Interim and JRA-55 suffer from terrain-related biases to a lesser ex-tent. Additionally, minimum T2M have less terrain-related systematic errorthan maximum T2M. The magnitudes of the random errors are smaller anduncorrelated with terrain error across all reanalyses.442.6. Results and discussionFigure 2.5: a) Kolmogorov-Smirnov (KS) statistic results for daily maximumT2M, for all four reanalyses by time of year, averaged across all locations inthe Southeast climate zone. The vertical dashed lines indicate the changein seasons and the colored dots represent the seasonal average. KS statisticvalues closer to zero are better. b) As an example, KS-statistic results fordaily maximum T2M on March 1 at Cranbrook (YXC) are shown.452.6.ResultsanddiscussionFigure 2.6: Terrain elevation error in meters for ERA-Interim. Other reanalyses’ errors are similar.462.6.ResultsanddiscussionFigure 2.7: Mean systematic and mean random error of daily maximum (TMAX) and minimum (TMIN) T2Mfor each of the 57 stations as a function of terrain reanalysis error. The solid lines show the linear regression fits.The nearly horizontal lines show the random error; the two lines are nearly on top of each other.472.6. Results and discussion2.6.2 Extreme T2MThe Southeast climate zone has the highest maximum T2M 30-year returnlevels, and the performance of the reanalyses here is representative of theoverall performance across BC (Fig. 2.8; Table 2.4). All four reanalysesshow negative (cold) systematic error (Fig. 2.8) throughout the entire calen-dar year. ERA-Interim and JRA-55 have a consistent systematic error andoutperform MERRA-2 and CFSR. MERRA-2 once again has a larger vari-ation in performance. CFSR also generally performs worse but has a moreconsistent bias. The Maritime climate zone has the lowest 30-year returnlevels, with values below 30◦C during summer (Figure not shown).482.6.ResultsanddiscussionTable 2.4: Seasonally averaged systematic error of 2- and 30-year return levels of daily maximum T2M by climatezone (SE2 and SE30 respectively in ◦C)Winter Spring Summer FallSE2 SE30 SE2 SE30 SE2 SE30 SE2 SE30NorthCFSR -4.48 -7.66 -6.76 -6.28 -3.66 -3.30 -5.31 -5.51ERA-Interim -2.89 -4.38 -3.94 -3.68 -3.65 -3.60 -3.99 -4.46JRA-55 -3.43 -4.83 -4.10 -3.69 -3.03 -2.72 -3.60 -3.98MERRA-2 -4.90 -8.35 -8.25 -8.66 -5.05 -4.68 -4.93 -5.72CentralCFSR -5.21 -6.20 -5.82 -5.08 -4.01 -3.84 -4.49 -5.11ERA-Interim -2.76 -3.64 -2.97 -2.61 -2.94 -3.22 -2.90 -3.66JRA-55 -3.68 -4.65 -3.29 -3.27 -2.93 -3.23 -3.37 -4.00MERRA-2 -6.16 -8.20 -6.60 -6.40 -3.29 -3.15 -3.58 -3.78MaritimeCFSR -2.30 4.72 -2.21 -3.72 -1.85 -4.06 -2.27 -4.01ERA-Interim -1.39 -3.83 -1.18 -1.95 -1.38 -2.66 -1.74 -3.16JRA-55 -1.90 -4.42 -1.71 -2.53 -2.19 -4.00 -2.32 -4.02MERRA-2 -1.95 -3.95 -1.69 -3.05 -1.28 -3.01 -1.52 -3.07SouthwestCFSR -3.78 -4.69 -3.59 -3.79 -2.69 -3.52 -3.05 -3.71ERA-Interim -3.31 -4.12 -2.62 -2.34 -2.12 -2.21 -2.75 -2.98JRA-55 -3.24 -4.55 -2.55 -2.73 -2.92 -3.53 -3.35 -3.84MERRA-2 -4.60 -4.84 -4.48 -4.09 -2.33 -2.17 -1.90 -1.79SoutheastCFSR -6.61 -8.73 -8.35 -6.84 -6.81 -6.81 -6.09 -5.76ERA-Interim -4.60 -5.36 -5.72 -5.05 -6.08 -5.98 -4.99 -4.74JRA-55 -5.49 -6.40 -5.80 -5.48 -5.65 -5.81 -5.35 -5.38MERRA-2 -7.18 -8.98 -8.57 -7.17 -6.22 -5.90 -4.43 -3.72492.6.ResultsanddiscussionFigure 2.8: 30-year return level and systematic error for maximum T2M, for all four reanalyses for Southeastclimate zone. The vertical dashed lines indicate the change in seasons and the colored dots represent the seasonalaverage.502.6. Results and discussionFor minimum T2M 30-year return levels, results are mixed. The JRA-55and ERA-Interim show generally positive systematic errors (larger duringthe winter) in the North and Central zones, while MERRA-2 and CFSRshow mostly negative systematic errors (Fig. 2.9; Table 2.5). This indicatesa warm bias in extreme minimum T2M for the JRA-55 and ERA-Interim,and a cold bias for the MERRA-2 and CFSR. In contrast, the Southeastand Southwest climate zones generally show negative systematic error inextreme minimum T2M in all four reanalyses, indicating they are typicallytoo cold (Figure not shown; Table 2.5). The Maritime climate zone has thelowest season-to-season extreme T2M variations (Figure not shown).512.6.ResultsanddiscussionTable 2.5: Seasonally averaged systematic error of 2- and 30-year return levels of daily minimum T2M by climatezone (SE2 and SE30 respectively in ◦C)Winter Spring Summer FallSE2 SE30 SE2 SE30 SE2 SE30 SE2 SE30NorthCFSR -1.02 -3.88 -1.75 -4.42 -1.14 -1.80 -0.16 -2.69ERA-Interim 2.20 1.73 0.95 1.76 0.11 0.04 -0.04 0.11JRA-55 3.93 4.95 1.05 1.49 -0.59 -0.63 0.11 0.35MERRA-2 -3.72 -4.74 -3.25 -3.97 -1.11 -1.26 -3.41 -4.50CentralCFSR 0.24 -3.04 -1.19 -2.76 -1.34 -1.66 0.57 -0.13ERA-Interim 2.56 3.74 0.73 1.30 0.73 1.17 0.63 1.37JRA-55 4.16 5.99 0.98 2.18 0.17 0.35 0.67 1.23MERRA-2 -3.02 -2.83 -2.06 -2.54 0.02 0.31 -1.51 -1.24MaritimeCFSR 3.09 1.90 3.44 3.12 2.38 3.07 3.55 3.41ERA-Interim 3.14 2.31 3.83 3.94 3.33 3.90 3.60 3.83JRA-55 3.83 2.74 4.51 4.51 3.43 3.99 4.46 4.17MERRA-2 1.91 0.87 3.15 3.08 2.64 3.39 3.19 2.99SouthwestCFSR -1.25 -2.86 -1.32 -2.71 -2.32 -2.58 -0.82 -1.89ERA-Interim -2.55 -3.84 -0.98 -2.36 -0.08 0.12 -1.09 -2.42JRA-55 0.15 -0.67 0.74 0.64 0.31 0.69 0.73 -0.03MERRA-2 -5.32 -6.64 -2.10 -3.30 -0.38 0.01 -1.66 -3.12SoutheastCFSR -6.91 -8.54 -5.11 -8.83 -4.37 -4.85 -3.41 -6.00ERA-Interim -3.69 -2.83 -2.77 -3.28 -2.35 -2.16 -2.89 -3.53JRA-55 -1.16 -0.27 -2.09 -2.32 -2.00 -1.82 -1.88 -2.40MERRA-2 -8.91 -10.19 -5.78 -7.64 -3.18 -2.84 -4.49 -5.81522.6.ResultsanddiscussionFigure 2.9: 30-year return level and systematic error for minimum T2M, for all four reanalyses for Central climatezone. The vertical dashed lines indicate the change in seasons and the colored dots represent the seasonal average.532.6. Results and discussionERA-Interim and JRA-55 extreme minimum T2M systematic errors aresmaller than CFSR and MERRA-2, which have larger errors during wintermonths. MERRA-2 once again exhibits the largest variation in performancethroughout the calendar year and across different climate zones.Similar to daily T2M, the systematic errors of 30- and 2-yr return levelT2M are related to reanalysis terrain elevation error (Fig. 2.10). Namely, thelarger the positive reanalysis terrain elevation error relative to the station’sreal world elevation, the larger the negative systematic error of the reanal-ysis temperature. Furthermore, the regression lines are vertically shifted inCFSR and MERRA-2, indicating larger errors across all elevations despitesimilar horizontal grid spacing (which should give them similar abilities toresolve terrain).An ANOVA indicates differences in mean systematic error of daily andextreme T2M between the reanalyses at the 1% significance level. After mul-tiple comparisons by Tukey’s HSD, daily and extreme minimum T2M areshown to be significantly better captured than daily and extreme maximumT2M. ERA-Interim and JRA-55 consistently significantly outperform CFSRand MERRA-2 for both daily and extreme T2M. If the four reanalyses areseparated into two groups — the more accurate JRA-55 and ERA-Interim,and the less accurate CFSR and MERRA-2 — the difference in mean sys-tematic error between groups is statistically significant. Within groups, theJRA-55 (CFSR) mean systematic error is smaller than the ERA-Interim(MERRA-2) mean systematic error, but this difference is not statisticallysignificant.Finally, Figure 2.11 shows the mean absolute error (MAE) of the reanal-yses daily and extreme (2- and 30-year return levels) T2M averaged overthe entire study period and all stations. The smaller the enclosed area fora given reanalysis, the better its performance. Daily minimum and maxi-mum T2M have smaller MAE than extreme minimum and maximum T2M.Daily and extreme minimum T2M are better captured by the reanalysesthan daily and extreme maximum T2M. Additionally, the errors of 2-yearand 30-year return levels are of similar magnitude. Results are similar fornon-assimilated stations. ERA-Interim and JRA-55 outperform CFSR andto a lesser degree MERRA-2 (not shown).542.6.ResultsanddiscussionFigure 2.10: Mean systematic error of 2- and 30-year return levels (RL) of extreme maximum (TMAX) andminimum (TMIN) T2M for each of the 57 stations as a function of terrain reanalysis error. The solid lines showthe linear regression fits.552.6.ResultsanddiscussionFigure 2.11: MAE (◦C) of daily maximum and minimum T2M, and of extreme (2- and 30-year return levels)maximum and minimum T2M. MAE is averaged across all 57 stations. Smaller values (closer to the origin) arebetter.562.7. Stationarity2.7 StationarityIn selecting a nonstationary model for the Gaussian distribution, a decisionmust be made about whether both the mean and standard deviation, oronly one of the parameters should be treated as varying in time.Variations through time of mean and standard deviation values for eachcalendar day at each station are modelled as a linear trend for daily maxi-mum and minimum T2M. The null hypothesis that there is no trend, whichindicates the parameters do not vary with time, is tested against the alter-nate hypotheses that the parameters are time-dependent at the αFDR = 0.05significance level.Figure 2.12a indicates a warming trend of mean values of daily minimumT2M across BC during summer, fall and winter, and a weak cooling trendduring spring across North and Central climate zones. Furthermore, Fig-ure 2.12b suggests a weak decreasing trend in standard deviation of dailyminimum T2M for fall and winter season. However, none of trends arestatistically significant.Under models M0 = N(µ, σ) and M1 = N(µ(t), σ), however, the criticaldifference for comparing these two models is statistically significant at theαFDR = 0.05 significance level, mostly near urban areas in the Southwestduring summer (Fig . 2.12c). This implies that a linear trend componentexplains a substantial amount of the variation in the data, and is likely to bea genuine effect in the warming trend process rather than a chance featurein the observed data.Figure 2.13a suggests a warming trend of mean values of daily maximumT2M across the Southeast and Southwest climate zones during summer, falland winter. Figure 2.13b shows a decreasing trend in standard deviation ofdaily maximum T2M during fall and winter across BC, and increasing trendin the spring.None of the trends in the daily maximum T2M parameters are statis-tically significant at the αFDR = 0.05 significance level, and a stationaryGaussian distribution with model M0 = N(µ, σ) is accurate enough to rep-resent daily maximum T2M during the study period with few exceptions atthe Southeast climate zone (Fig. 2.13c).Analogously, Fig. 2.14a shows strong warming of extreme maximum T2Mduring summer, and strong warming of extreme minimum T2M during win-ter (Fig. 2.14c). The LRT deviance statistic suggests the evidence sup-porting such a trend is weak, implying no significant improvement over thestationary model at the αFDR = 0.05 significance level (Fig. 2.14b,d). Dur-ing the 31 years of the study period, changes in return levels of extreme572.7. StationarityFigure 2.12: (a) Mean linear trend of daily minimum T2M; (b) Standarddeviation linear trend of daily minimum T2M; (c) Days when N(µ(t), σ) isstatistically significant. The vertical dashed lines indicate the change in sea-sons and the horizontal dashed lines delineate from top to bottom the North,Central, Maritime, Southwest and Southeast climate zones respectively.582.7. StationarityFigure 2.13: (a) Mean linear trend of daily maximum T2M; (b) Standarddeviation linear trend of daily maximum T2M; (c) Days when N(µ(t), σ) isstatistically significant. The vertical dashed lines indicate the change in sea-sons and the horizontal dashed lines delineate from top to bottom the North,Central, Maritime, Southwest and Southeast climate zones respectively.592.7. StationarityT2M are generally too weak for nonstationarity to be required in order tocharacterize extreme levels.602.7.StationarityFigure 2.14: (a) Location linear trend of extreme maximum T2M; (b) Days when GEV (µ(t), σ, κ) for extrememaximum T2M is statistically significant; (c) Location linear trend of extreme minimum T2M; (d) Days whenGEV (µ(t), σ, κ) for extreme minimum T2M is statistically significant. The vertical dashed lines indicate thechange in seasons and the horizontal dashed lines delineate from top to bottom the North, Central, Maritime,Southwest and Southeast climate zones respectively.612.8. Conclusion2.8 ConclusionThe skill of reanalysis daily and extreme T2M were investigated acrossBritish Columbia (BC) during the 1980-2010 study period. To identify thebest reanalyses among CFSR, ERA-Interim, JRA-55 and MERRA-2, thesystematic error, random error and KS statistic were used to compare dailymaximum and minimum T2M. To evaluate extreme maximum and minimumT2M, systematic error of 2- and 30-year return levels were compared.ERA-Interim and JRA-55 had the lowest systematic and random error,and more consistent lower errors throughout the year and across all climatezones in BC. One possible explanation besides T2M being assimilated bythese reanalyses is differences in the data assimilation methods. In regionswhere data is sparse, a reanalysis solution is more dependent on the modelstructure than on observations (Lindsay et al., 2014). ERA-Interim andJRA-55 use 4D-Var while CFSR and MERRA-2 use 3D-Var. 4D-Var betterinterprets physical information implicit in the meteorological model equa-tions where observations are sparse, and extracts and transfers informationinto data-void areas in a more consistent way (The´paut, 2006; Whitakeret al., 2009; Dee et al., 2011). Both issues are relevant for BC as weatherstation coverage is sparse outside Southwestern BC, and the Pacific Oceanto the west is a well-known data void (Hacker et al., 2003; Spagnol, 2005).Furthermore, the previous generation JRA-25, which has a 3D-Var assimi-lation method, shows a larger variation in performance in T2M throughoutthe calendar year (Lindsay et al., 2014) when compared to the results ofJRA-55 shown here.Minimum temperatures are consistently better captured by the reanaly-ses than maximum temperatures, whether it is daily or extreme T2M. Thisconsistent difference in skill between maximum and minimum temperaturessuggest they should be treated separately, rather than assessing hot or coldextreme weather events by looking at different tails of a common distributionof mean T2M (Graham et al., 2013; Hart and Grumm, 2001).Furthermore, reanalysis daily and extreme maximum T2M generallyhave a cold bias, as do daily and extreme minimum T2M. Previous studiesshow a cold bias in T2M in CFSR and ERA-Interim at different pressurelevels (Bao and Zhang, 2013), and in MERRA-1, JRA-55 and CFSR overAntarctica (Jones et al., 2016).ERA-Interim extends as far back as 1979, whereas JRA-55 extends backto 1958, giving the latter a longer record and a larger sample size. A longerhistorical record is desirable, as a natural difficulty of extreme weather statis-tics is the limited amount of extreme data. Additionally, the standard errors622.8. Conclusionof the estimated parameters and of the return levels are expected to decreaseas the sample size increases, whether they were obtained by MLE or methodof L moments (Hosking et al., 1985; Hosking, 1990; Cai and Hames, 2010).There is a noticeable warming trend in the Southwest and Southeastclimate zones of BC during summer months for daily maximum T2M; andduring summer months for extreme maximum T2M; during summer, falland winter months for daily minimum T2M; and during winter for extrememinimum T2M. Vincent and Mekis (2006) showed that the number of coldevents has significantly decreased while the number of warm events hassignificantly increased over BC. This warming trend could in part be causedby urbanization surrounding the weather stations (Gunn, 2010; Jones et al.,1985). The population of BC has grown faster than the Canadian nationalrate, with Southwest BC leading the growth rate (Demography Division,2016). Despite some evidence of warming and cooling trends, most werenot statistically significant. Thus, for this 31-year study period, temporallyvarying distributions are not needed to represent T2M over BC. It is possiblethat, were a longer period of record used, more widespread significant trendsmight be found.While random error and KS statistic cannot be bias corrected or easilyimproved otherwise, this study found that systematic errors in daily andextreme T2M are largely explained by the reanalysis terrain elevation error,and thus could be largely corrected/eliminated. Systematic T2M errors are alarger component of the total error than random errors, thus a bias-correctedreanalysis T2M would be a substantially improved dataset, particularly overcomplex terrain. The methodology presented here should be able to be usedin other parts of the world.Such a correction will be discussed in Chapter 4. First, Chapter 3 willevaluate daily and extreme accumulated precipitation produced by thesereanalyses over BC.63Chapter 3Performance of Reanalysesacross British Columbia.Part II: Evaluation of Dailyand Extreme Precipitation3.1 IntroductionIn Chapter 2 performance of latest-generation reanalyses with respect todaily and extreme maximum and minimum 2-m temperature (T2M) overmountainous BC is discussed. In this chapter performance with respectto daily and extreme precipitation is assessed (hereafter daily and extremePCP; defined in section 3.2), and trends in both daily and extreme PCPduring the study period are examined to determine if significant statisticalchanges occurred over the timespan of the dataset.An added difficulty is that verification of precipitation extremes is morechallenging than that of extreme temperatures (Bhend and Whetton, 2013;van Oldenborgh et al., 2013). Because of the small spatial and time scalesof precipitation, in general, numerical models do not simulate precipitationas well as they do temperature (Kendon et al., 2014; Ravishankar et al.,2016), and have limited ability to faithfully represent extreme precipitationevents (Zhu et al., 2014). For instance, differences in the reanalyses pa-rameterizations of convection and other physical processes can impact howwell extreme precipitation events are represented (Dee et al., 2011; Lindsayet al., 2014).Previous studies have evaluated reanalysis performance with respect toprecipitation and moisture fields. Flux tower observations over the NorthernHemisphere of temperature, wind speed, precipitation, downward shortwaveradiation, net surface radiation, and latent and sensible heat fluxes were usedto evaluate the performance of CFSR, ERA-Interim, ERA-40, and MERRA,where ERA-40 was found to have the lowest precipitation bias and ERA-643.2. Data and MethodologyInterim best captured precipitation variability (Decker et al., 2012). Hodgeset al. (2011) explored how well JRA-25, ERA-Interim, MERRA, and CFSRidentify extratropical cyclones over the Northern and Southern Hemispheres,and found that the latest-generation reanalyses better represent cyclones,especially in the Southern Hemisphere. Berg et al. (2003) found that ERA-40 has positive biases in precipitation over land in North America, andRuiz-Barradas and Nigam (2005) found that ERA-40 also has positive biasesin evapotranspiration during the warm season over the U.S. Great Plains.Finally, Bosilovich et al. (2015) found that MERRA-2 has a more consistentglobal precipitation average than MERRA, and a lower global precipitationbias than JRA-55 and CFSR when compared to the Global PrecipitationClimatology Project (GPCP; Adler et al. (2003)).Precipitation is examined here because of its various financial, societal,and environmental impacts such as hydroelectric power generation (Odonet al., 2017), flooding and water management (White et al., 2016; Odonet al., 2017; Sun et al., 2018b), agriculture (Rosenzweig et al., 2001; Sunet al., 2018b), tourism (Patz et al., 2005; White et al., 2016), health (Currieroet al., 2001; Patz et al., 2005), and flora and fauna (Parmesan et al., 2000).Furthermore, several studies have noted increases in the frequency andintensity of extreme precipitation events in various parts of the world (Mannand Emanuel, 2006; Krishnamurthy et al., 2009; Donat et al., 2013; Westraet al., 2013; Ravishankar et al., 2016). An increase in extreme precipita-tion may lead to other impacts such as increase in winter runoff, which inturn may lead to flooding, and challenges in drainage and sewage systemscapacities (White et al., 2016; Sun et al., 2018b).In section 3.2, a brief description of the different reanalyses and of theweather station observations is given. In sections 3.3 and 3.4, we describethe methodology for dividing BC into climate zones, and the various metricsused for evaluating daily and extreme reanalysis PCP. In section 3.5, dailyand extreme PCP from the reanalyses are evaluated. In section 3.6, themethods for assessing statistical nonstationarity are introduced, and trendsof both daily and extreme PCP are examined. Results are summarized inthe conclusion.3.2 Data and MethodologyDaily accumulated precipitation from 66 weather stations from 1 Jan 1980to 31 Dec 2010 are used in this study to evaluate the CFSR, ERA-Interim,JRA-55 and MERRA-2 reanalyses. The 1980-2010 study period is chosen653.2. Data and Methodologybecause it is the longest overlap between the four reanalyses. MERRA-2began in 1980 (Gelaro, 2015; Gelaro et al., 2017) and CFSR ended in 31 Dec2010. From January 1, 2011 forward, the CFSR was extended using NCEPClimate Forecast System Version 2 (CFSv2) operational model. Differencesbetween the model used to produce the CFSR and the operational CFSv2may affect data evaluation past the extension date (Saha et al., 2014).A description of the different reanalyses and of the weather stationdataset is given below. A summary of the reanalyses atmospheric mod-els and configurations are presented in Table 3.1. A broader description andcomparison of the latest and previous generation reanalyses can be found inChapter 2.663.2.DataandMethodologyTable 3.1: Overview of the four reanalysis datasets examined in this study.Institution Reanalysis Model AssimilationmethodPeriod Downloadgrid(lat× lon)TimeintervalReferenceNCEP/NCAR CFSR CFST382/L64(globalhorizontalresolution∼ 38 km)3D-Var GSI 1979-2011(current asCFSv2)0.5◦ × 0.5◦(∼ 50 km)6-haccumu-lationat 0000,0600,1200and 1800UTCSaha et al.(2010)ECMWF ERA-InterimIFST255/L60(globalhorizontalresolution∼ 79 km)4D-Var 1979-current 0.5◦ × 0.5◦(∼ 50 km)6-h and12-haccumu-lationat 0000and 1200UTCDee et al.(2011)JMA JRA-55 JMAT319/L60(globalhorizontalresolution∼ 55 km)4D-Var 1958-2012(current asJCDAS)0.5616◦ ×0.5616◦(∼ 55 km)3-haccumu-lationat 0000,0300, ...,and 2100UTCEbita et al.(2009)NASA MERRA-2GEOS-5.12.4AGCM(lat × lon)0.5◦ ×0.625◦/L723D-Var GSI 1980-current 0.5◦ ×0.625◦1-haccumu-lationat 0030,0130, ...,and 2330UTCGelaro(2015);Gelaroet al.(2017)673.2. Data and MethodologyPrecipitation from the weather stations used in this study are not as-similated by the ERA-Interim (Dee et al., 2011), the MERRA-2 (Bosilovichet al., 2015) or the JRA-55 (Kobayashi et al., 2015); but are indirectlyassimilated by the CFSR (see subsection 2b; Wang et al. (2011)). There-fore evaluating against these observations provides a reasonably independentmeasure of accuracy.3.2.1 Weather station dataDue to the higher spatial variability of precipitation compared to temper-ature, more precipitation stations are included in this study. Of the 111geographically-dispersed precipitation stations initially selected, 45 stationswith more than 4% missing data were excluded. Of the remaining 66 sta-tions, 57 are from Environment and Climate Change Canada (ECCC) and9 are from BC Hydro (Table A.1).Figure 3.1 shows the locations of all 66 stations overlaid with populationdistribution across BC. Fifty-seven stations are located in valleys (indicatedby upside-down triangles), and nine are in non-valley locations (uprighttriangles).683.2.DataandMethodologyFigure 3.1: Location of ECCC (red) and BC Hydro (blue) weather stations used for precipitation analysis andBritish Columbia population distribution (orange). Upside-down triangles indicate valley stations; upright trian-gles indicate non-valley stations. Dashed lines delineate the dominant climate zones North, Northwest, Central,South Central, Maritime West, Maritime East, Southwest and Southeast.693.2. Data and MethodologyLong-term time series often contain variations caused by changes in theenvironment surrounding the gauges, instrumentation, observing proceduresincluding the time of observation, station location, or discontinuation of thestation. As a result, variations unrelated to changes in weather and climatecan be introduced into the time series. Different adjustment techniquesfor precipitation have been developed to address impacts on climate datahomogenization [e.g., Jones et al. (1985); Peterson et al. (1998); Mekis andHogg (1999); Vincent et al. (2002)].The methods to adjust daily rainfall and snowfall for ECCC stations aredescribed in Mekis and Vincent (2011). For each rain gauge type, correc-tions were implemented to account for undercatch due to wind, evaporation,and gauge-specific wetting losses. A complete description of gauges can befound in Metcalfe et al. (1997) and Devine and Mekis (2008). For snow-fall, density corrections based on coincident ruler and Nipher measurementswere applied to all snow ruler measurements (Mekis and Brown, 2010). Adetailed description of trace (non-measurable precipitation amount) relatedissues and adjustments are given in Mekis (2005) and Mekis and Vincent(2011).Daily total precipitation was calculated by adding a station’s daily rainand snow observations together. In case of station relocation, a new iden-tification number is often given to the new location and observations fromthe two stations are combined to create a longer time series. Adjustmentsare applied to join the two datasets, based on standardized ratios betweenthe sites and neighbouring sites, or overlapping observation periods (Vincentet al., 2009).Data homogeneity for BC Hydro stations was assessed by BC Hydrousing Double-Mass Curves (DMC) (Searcy and Hardison, 1960). The the-ory of the double-mass curve is based on the graph of the precipitation ata station against precipitation of surrounding reference stations during thesame period. A break in the slope of the double-mass curve means thata change in the constant of proportionality between the station and sur-rounding reference stations has occurred. The data before the date that thechange occurred is modified to match the historic relationship between thestation and its reference stations.Reliable data is required in order to detect trends in daily and extremePCP. When dealing with trends of daily data, it is important that the datasetis nearly complete during the studied period. Furthermore, when analyzingdecade-long trends, it is important that years with many missing data, ifthey occur, are relatively few and not clustered during a certain time inter-val, as this period might have had an anomalous climate (Moberg and Jones,703.2. Data and Methodology2005; Vicente-Serrano et al., 2010). Finally, the reliability of frequency ofextreme precipitation events is closely related to the sample size used duringthe study period (Hosking et al., 1985; Hosking, 1990; Porth et al., 2001; Caiand Hames, 2010). Stations with more than 1% missing data are excludedfrom our nonstationarity analysis, leaving twenty-nine ECCC and six BCHydro stations.3.2.2 CFSRPrecipitation is generated by the model during the direct assimilation oftemperature and humidity information from satellite radiances (Saha et al.,2010). Then, the model-generated precipitation is corrected with the CPC(the National Oceanic and Atmospheric Administration (NOAA) ClimatePrediction Center) Merged Analysis of Precipitation (CMAP) product (Xieand Arkin, 1997) — which defines 5-day mean precipitation on a 2.5◦× 2.5◦latitude-longitude grid over the globe by merging information derived fromgauge and satellite observations — and the CPC Unified Gauge-Based Anal-ysis of Global Daily Precipitation (CPCU) product — which interpolatesquality-controlled rain gauge reports collected from the Global Telecommu-nication System (GTS) and many other national and international archives(Saha et al., 2010) on a 0.5◦ × 0.5◦ latitude-longitude grid over the globe.Finally, an algorithm accounts for orographic enhancements in precipitation(Xie et al., 2007).3.2.3 ERA-InterimPrecipitation is generated by the model during the variational analysis ofupper-air atmospheric fields such as temperature, wind, humidity and ozone,in combination with direct assimilation of 2-m temperature, 2-m relativehumidity and 10-m winds from land stations, and upper-air temperatures,wind, and specific humidity from radiosonde data (Dee et al., 2011).3.2.4 JRA-55Precipitation is model-generated during direct assimilation of upper-air tem-peratures and humidity information from satellite radiances, and direct as-similation of surface pressure, 2-m temperature, 2-m relative humidity, 10-mwinds from land stations, and upper-air temperatures, winds, and specifichumidity from radiosonde data (Ebita et al., 2011).713.2. Data and Methodology3.2.5 MERRA-2There are two kinds of precipitation fields in the MERRA-2 system. The pre-cipitation generated by the model during the assimilation procedure (Bloomet al., 1996; Reichle et al., 2017) (PRECTOT is the variable name in the dataproduct file), and the corrected precipitation that is seen by the MERRA-2land surface and that modulates aerosol wet deposition over land and ocean(PRECTOTCORR is the variable name in the data product file). As men-tioned, for a consistent independent evaluation of all reanalyses performance,the former is used in this study.The precipitation is generated by the model during the direct assim-ilation of temperature and humidity information from satellite radiances(Bloom et al., 1996; Koster et al., 2016).3.2.6 PRISM datasetPrecipitation is not evenly distributed over weather station areas in com-plex terrain (Taesombat and Sriwongsitanon, 2009). In order to estimateareal precipitation, it is preferable to have as many weather stations as pos-sible. However, spatial and temporal coverage is a limiting factor (Karlet al., 1993; Odon et al., 2018), as is accuracy and reliability of precipitationrecords (Metcalfe et al., 1997; Serreze et al., 2005). Additionally, in orderto evaluate the agreement between observations and reanalyses, it is impor-tant to realize that each grid point in a reanalysis represents an averagecentered on the geographical coordinates of the grid point. By contrast, anobservation represents a point measurement within a reanalysis grid-box,which may or may not be representative of the grid box average. Further-more, grid resolution and location in each reanalysis dataset is different,and the location of a station can vary from the centre to the edge of thegrid box. Interpolation techniques may produce inaccurate results becauseof the effects of topographical variation and the limited number of avail-able rainfall stations (Taesombat and Sriwongsitanon, 2009). In order toidentify regional and terrain biases and to improve the accuracy of arealrainfall estimation, the reanalyses are billinearly interpolated to the samehigh-resolution grid as the Parameter-Elevation Regressions on IndependentSlopes Model ( PRISM; Daly et al. (1994, 1997, 2002)) climatology for gridcomparison.PRISM climatology, produced by the Pacific Climate Impact Consortium(PCIC), the Pacific Institute for Climate Solutions (PICS) and the BC Gov-ernment provides access to a 30-arc-second (∼ 800 m) gridded precipitation723.3. Climate zonesclimatology for the 1981-2010 climate normal period, for land-surface areasof BC (PCIC et al., 2014). PRISM has been tested and verified throughoutthe United States and has been applied in numerous countries across theglobe including western Canada previously for the 1961-1990 period. In thisstudy, the 1981-2010 climate period is used.3.3 Climate zonesPrincipal component analysis (PCA) is employed to separate stations intogroups with similar precipitation climates. Because precipitation has bothfrequency and magnitude, contingency tables are employed. Daily PCP(mm) is divided into seven non-overlapping bins for the entire study period:[0, 0.25], (0.25, 1.0], (1.0, 2.5], (2.5, 5.0], (5.0, 10.0], (10.0, 20.0] and (20.0,∞).The number of events falling into each bin is recorded in a contingency ta-ble for the entire study period . A χ2-test compares the pairwise differencebetween the distributions of the 66 stations. A sufficiently large differencebetween the distribution of the daily PCP amounts over the seven bins in-dicate the stations have different precipitation characteristics and thereforebelong to different climate zones. Small differences suggest that the sta-tions have similar climates. This critical difference is determined by theχ2-distribution. Due to the large number of stations pair comparisons (2141pairs), the null hypothesis that the stations belong to the same climate zoneis rejected at the αWalker = 1−(1−α0)1/N0 = 2.31×10−5 level of significance,where α0 = 0.05 and N0 = 66 (Fig. 3.2).Cramer’s correlation is computed for each pair χ2-statistic, and PCA isconducted on the Cramer’s correlation matrix of daily PCP. The first 41components explain 90% of the variability in the data, and are retained.The number of components is larger here than in Chapter 2 of this studybecause of the higher variability in precipitation. For temperature, wherethe Pearson correlation between the stations is high (see section 2.3 formore details), it is possible to capture most of their variance using a smallernumber of principal components. For precipitation, where the Cramer cor-relation between the stations is lower (Fig. 3.2), more principal componentsare needed to capture most of their variance.733.3.ClimatezonesFigure 3.2: Cramer correlation matrix and eight dominant climate zone clusters in BC. Crosses indicate stationswhere differences in precipitation distribution are statistically significant.743.4. Verification MetricsA K-means clustering analysis is then performed on the componentsto arrange the data into groups. One to 10 clusters were tested and aneight-cluster solution was chosen (Fig. 3.2). This analysis yields 8 climatezones (North, Central, Northwest, Maritime West, Maritime East, South-west, South Central, and Southeast) that roughly matched to those identifiedin Chapter 2.The Maritime climate zone from Chapter 2 has been divided into Mar-itime East and Maritime West, while the Southeast climate zone has beendivided into Southeast and South Central. These additional climate zoneshighlight the differences between windward and leeward regions, which seeenhanced upslope precipitation and rain shadows, respectively.3.4 Verification MetricsThe statistical behaviours of daily and extreme PCP are compared betweenobserved weather station data, and their corresponding location in the re-analyses and PRISM. As in Chapter 2, four horizontal interpolation meth-ods were trialed for reanalysis output. Also as in Chapter 2, the methodsNearest Neighbour, Inverse Distance Weighting (IDW), Bilinear and Bicubicinterpolation perform very similarly; IDW is used. Furthermore, IDW fromland-only grid points (omitting sea grid points) was also tested for coastalstations in the Maritime West, Maritime East, Southwest and Northwest cli-mate zones ( a grid point is classified as a land point based on each reanalysisland-sea mask). Some small variations between the two IDW methods werefound, but overall the land-only IDW is similar in result to the IDW.3.4.1 Daily PCPThe Canadian meteorological convention in defining a calendar day to befrom 0601 UTC of that day to 0600 UTC of the following day is followed.For CFSR, JRA-55 and MERRA-2, precipitation accumulation intervals aresummed over this window to get daily PCP (Table 3.1). For ERA-Interim,precipitation is accumulated in the forecast sense, i.e., reset to zero at 0000and 1200 UTC. Six-hour accumulated precipitation for the 6-h periods pre-ceding 0000, 0600, 1200 and 1800 UTC are obtained in the following manner:for 0600 UTC the 0000-UTC 6-h accumulated precipitation is used; for 1200UTC the 0000-UTC 6-h accumulated precipitation is subtracted from the0000-UTC 12-h accumulated precipitation; for 1800 UTC the 1200-UTC 6-haccumulated precipitation is used; and for 0000 UTC the 1200-UTC 6-h pre-cipitation is subtracted from the 1200-UTC 12-h accumulated precipitation.753.4. Verification MetricsA 31-day centered rolling accumulation window is used to obtain smoothmonthly mean precipitation for each calendar day. These 31-day accumu-lated values for each calendar day are then averaged over the 31-year evalu-ation period (1980-2010). This is done for each station and reanalysis data.The percentage bias (or systematic error) is then computed to estimate howaccurately each reanalysis captures monthly precipitation (31-day precipi-tation totals).For the same 31-day centered rolling window, days with daily PCP below1.0 mm are classified (and hereafter referred to) as ”dry” days, and daysequal or above 1.0 mm are classified (and hereafter referred to) as ”wet” days.In climate studies, this delineation is typically made at trace amounts or 1.0mm (Vincent and Mekis, 2006; Werner and Cannon, 2016). In this study 1.0mm was chosen because coarse resolution models tend to overforecast thefrequency and spatial extent of light precipitation events (e.g., Zhu and Luo(2015), and shown later in this study). On every 31-day window, systematicerror is computed from the number of wet days at a station location inthe reanalysis with respect to the actual number of wet days observed at astation.Finally, wet days are divided into the five non-overlapping intervals:[1.0 mm, 50th), [50th, 75th), [75th, 90th), [90th, 95th) and [95th, 100th]. Thepercentiles are calculated from the entire wet-day climatological distribution,centered on each calendar day, using each station’s observed data.This allows for an evaluation of a reanalysis’ ability to correctly capturethe frequency and distribution of precipitation intensities. If the reanal-ysis distribution is very close to that of the observation, the number ofexpected occurrences in each bin will be very close. This difference betweenthe number of “light” ([1.0 mm, 50th)), “light” ([50th, 75th)), “moderate”([75th, 90th)), “heavy” ([90th, 95th)) and “extreme” ([95th, 100th]) precipi-tation events in each dataset is given by the two-sample χ2-statistic.3.4.2 Extreme PCPThe definition of an extreme precipitation event varies widely. One possi-bility is to define it as an event in which precipitation over some specifiedperiod exceeds some threshold, either at a point measured by a single raingauge, or spatially averaged over some region. The choice of threshold alsovaries. Some studies use fixed absolute thresholds while others use a fixedpercentile based on the distribution specific to a given location, so that it isspecific to the location climatology. In this study, the 2- and 30-yr returnlevels for every station for each calendar day in the 1980-2010 study period763.4. Verification Metricsare estimated.Furthermore, end users of precipitation forecasts, such as hydrologists,are concerned with both peak flows and total volume of flows, especiallywhen they deviate from climatology. Correspondingly, they require accu-rate estimates of precipitation intensity and accumulation over a range oftime scales. For the flashy reservoirs of the BC South Coast, the timescalesof their total volume concerns typically range from 1 to 14 days. Hydropowerfacilities and their associated operating procedures are designed assumingestimated minimum and maximum volumes that could be possibly expectedat various timescales. Extreme precipitation accumulations, whether accu-mulated over 1 or 14 days, will cause extreme total flow volumes over thosetime scales, pushing or possibly exceeding the limits of a hydropower facil-ity. The effects of heavy or extreme precipitation events can be compoundedif components of a hydropower facility (e.g., spill gates) are out of servicefor planned or unplanned maintenance. Given that flow volumes are im-portant at a range of time scales, 2- and 30-year return levels of 1-, 3-, 7-,and 14-day accumulated precipitation are examined (2- and 30-year returnlevels are also known as 2- and 30-year recurrence intervals, or 0.5 and 0.03Annual Exceedance Probabilities (AEP)).To do this, the maximum 1-, 3-, 7-, and 14-day accumulated precipi-tation within a 31-day centered rolling window, for each calendar day arecalculated (Fig. 3.3a). This is done for each year from 1980-2010 inclusive,yielding 31 annual maximum values of 1-, 3-, 7-, and 14-day accumulatedprecipitation for each calendar day. A 31-day window was chosen so that allvalues within the window are from the same time of year, and would havesimilar climatological precipitation distributions.A Generalized Extreme Value distribution (GEV) fitted by the methodof L moments dresses these 31 sample values of 1-, 3-, 7-, and 14-day accu-mulated precipitation for each calendar day (Fig. 3.3b-d). A GEV is chosenbecause of the interest in the statistical behaviour of the 31 annual maximumvalues of 1-, 3-, 7-, and 14-day accumulated precipitation at each calendarday. Estimates of the 2- and 30-year return levels of annual maximum arethen obtained from the fitted GEV.The Lilliefors (Wilks, 2011) test compares the largest difference, in ab-solute value, between the theoretical GEV cumulative distribution function(CDF) and the observed empirical cumulative distribution function (ECDF).The null hypothesis is that the observed data is drawn from the distributionbeing tested (i.e., the observation ECDF and GEV CDF are indistinguish-able), and a sufficiently large critical difference results in the null hypothesisbeing rejected.773.4. Verification MetricsFigure 3.3: (a) Diagram illustrating the 31-day centered rolling window forJuly 9th, performed over the 31-year period; (b) GEV model fitted over the31-day, 31-year centered window for 1-day precipitation total at VancouverInternational Airport (YVR); (c) Quantile plot for the GEV fitted model forthe same day, location, and variable; If the GEV is a reasonable model, thepoints on the quantile plot should lie close to the unit diagonal; (d) Returnlevel, return period plot for the same day, location, and variable; showingthe precipitation that corresponds to a given return period.783.4. Verification MetricsA parametric bootstrap procedure determines the critical value that re-sults in the rejection of the null hypothesis. 100 samples of size 31 are gener-ated from the fitted GEV distribution for each calendar day at each station,and a critical value is derived from each of the 100 generated samples. Thatis, the 90th percentile of the resulting collections of critical values is thenused as the critical value for the rejection of the null hypothesis — that thesample originates from the GEV distribution at the 10% significance level(α0 = 0.10). However, there are 66 weather stations and 365 calendar daystotaling N0 = 24090 independent Lilliefors tests. Due to very large numberof tests, the αWalker = 1− (1−α0)1/N0 = 4.37× 10−6 is instead regarded assignificant, and the (1− αWalker)th percentile of the resulting collections ofcritical values is then used as the critical difference (Wilks, 2016).Less than 1% of the locations and calendar days where the null hypoth-esis is tested for all different accumulations totals were rejected during the1980-2010 study period, suggesting the annual extremes in fact can be de-scribed by a GEV distribution.The systematic error of the 2- and 30-year return levels were then calcu-lated to estimate how each reanalysis captures annual extremes of 1-, 3-, 7-,and 14-day accumulated precipitation. These two return levels are chosenbecause less than 1% of their 90% confidence intervals overlap, indicatingthe difference between the two return levels is statistically significant, andbecause the 30-year return level is the most extreme verifiable value giventhe length of the data.To calculate their 90% confidence intervals, 100 samples of size 31 aregenerated from the fitted GEV distribution for each calendar day at eachstation, and the 2- and 30-year return levels are estimated from each gener-ated sample. Then, the 5th and 95th percentiles of the resulting collectionof 2- and 30-year return levels are used as the lower and upper bounds ofthe 90% confidence intervals for the true 2- and 30-year return levels.3.4.3 Kruskal-Wallis AnalysisThe mean systematic error of monthly precipitation total, two-sample χ2-statistic, 31-day-window percentage of wet days, and of 2- and 30-year returnlevels of 1-, 3-, 7-, and 14-day accumulated precipitation are calculated foreach station from all calendar days systematic errors. Comparisons betweenthese eleven mean systematic errors of each reanalyses (CFSR, ERA-Interim,JRA-55 and MERRA-2) are performed using eleven independent Kruskal-Wallis nonparametric tests. Eleven independent Kruskal-Wallis tests areused because the mean systematic errors are skewed, and due to the differ-793.5. Results and discussionent magnitude and variability of each type of mean systematic error. Finally,Nemenyi’s test (Hollander et al., 2013) is applied following statistical signif-icance at the αWalker = 9.53× 10−3 level (α0 = 0.10) in the Kruskal-Wallisresults to identify significant performance differences in pairwise compar-isons between the reanalyses mean systematic errors.3.5 Results and discussion3.5.1 Daily PCPFirst, performance of daily PCP across the climate zones in the reanalysesare investigated. Reanalysis performance in the Northwest climate zone(Fig. 3.4) is representative of performance across the wetter climate zones(Northwest, Maritime West, Maritime East and Southwest; latter three notshown). The seasonal cycle and magnitudes of 31-day precipitation totalsare fairly well captured. All four reanalyses exhibit similar seasonal cyclesand differ mostly in magnitude of annual bias (Table 3.2). They showa positive (wet) bias all year long for the Northwest and Maritime Eastclimate zones. In the Maritime West (the wettest climate zone in BC and inall of Canada) and Southwest zones, JRA-55, ERA-Interim and CFSR havea negative (dry) bias throughout the year (not shown).803.5.ResultsanddiscussionTable 3.2: Averaged systematic error of: monthly precipitation total (M), two-sample χ2-statistic (χ2), 31-day-window percentage of wet days (W), and 30-year return levels of 1- (1D30), 3- (3D30), 7- (7D30), and 14-day(14D30) accumulated precipitation across wetter climate zones Northwest, Maritime West, Maritime East andSouthwest, drier climate zones North, South Central, Central and Southeast, and all climate zones in BC (allsystematic errors but two-sample χ2-statistic in %)M χ2 W 1D30 3D30 7D30 14D30WetCFSR 8.57 26.16 6.44 -31.35 -21.24 -13.74 -8.71ERA-Interim 14.98 26.44 5.57 -11.83 -3.74 0.80 2.57JRA-55 11.50 18.92 4.29 -15.79 -6.51 -2.32 -0.03MERRA-2 23.22 18.10 9.04 11.30 -4.80 0.79 3.46DryCFSR 98.56 15.90 20.41 -2.15 12.21 26.01 38.01ERA-Interim 55.48 16.23 14.74 -11.99 -1.05 7.01 15.67JRA-55 59.30 13.93 12.63 -6.53 5.59 13.06 19.98MERRA-2 69.45 12.44 17.65 2.89 10.61 17.04 24.46AllCFSR 48.11 21.65 12.58 -18.52 -6.54 3.73 11.82ERA-Interim 32.77 21.96 9.60 -11.90 -2.56 3.53 8.33JRA-55 32.50 16.73 7.95 -11.72 -1.19 4.43 8.76MERRA-2 43.53 15.61 12.83 -5.05 1.98 7.93 12.69813.5.ResultsanddiscussionFigure 3.4: Observed and reanalysis running centered 31-day precipitation totals and systematic error averagedover stations in the Northwest climate zone. The vertical dashed lines indicate the change in seasons and thecolored dots represent the seasonal average.823.5. Results and discussionOverall, JRA-55 and ERA-Interim outperform CFSR and MERRA-2.The latter two exhibit higher bias and higher variability in bias throughoutthe wetter climate zones and seasons.In contrast, in the South Central zone, with the driest locations in allCanada, all reanalyses have a large wet bias all year long (not shown). Inthe remaining drier climate zones, North, Southeast and Central, reanalyseshave a large wet bias for most of the year and the largest wet bias duringspring (Fig. 3.5 for Central climate zone). Again, JRA-55 and ERA-Interimoutperform CFSR, and to a lesser extent MERRA-2 (Table 3.2). The generalwet bias in most zones could be a result of the low resolution of the reanalysesthat tends to spread out precipitation into drier portions of grid cells, failingto capture the locally drier climates of lower (valley) elevations, where moststations are located.833.5.ResultsanddiscussionFigure 3.5: Observed and reanalysis running centered 31-day precipitation totals and systematic error averagedover stations in the Central climate zone. The vertical dashed lines indicate the change in seasons and the coloreddots represent the seasonal average.843.5. Results and discussionThe percentages of “wet” and “dry” days, and the two-sample χ2-statisticare illustrated in Figure 3.6 (a, b, and c) for the Maritime West climatezone, which is representative of the wetter climate zones (Maritime West,Southwest, Northwest and Maritime East; latter three not shown). Thereanalyses capture the annual cycle of precipitation frequency better thanthey do the precipitation amounts (Fig. 3.6a, c; Table 3.2). The two-sampleχ2-statistic (Fig. 3.6b) is used to determine how well the reanalyses cap-ture the histogram of observed precipitation events across all bins. Lowervalues are better, and 0 would indicate no difference between the observa-tion histogram and that of the reanalysis. The MERRA-2 is the best dueto its consistent low values of χ2-statistic across the wetter climate zones,followed closely by JRA-55 (Table 3.2). Finally, looking at the precipita-tion percentile bins, reanalyses overestimate the number of “very light” and“light” precipitation events (Fig. 3.6d, e), and underestimate the numberof “moderate”, “heavy” and “extreme” precipitation events for the Mar-itime West and Southwest zones (Fig. 3.6f-h). The opposite is true for theMaritime East and Northwest (not shown).853.5.ResultsanddiscussionFigure 3.6: Percentage of (a) wet days; (b) two-sample χ2-statistic; percentage of (c) dry days, (d) very light, (e)light (f) moderate, (g) heavy and (h) extreme precipitation events for Maritime West climate zone. The verticaldashed lines indicate the change in seasons and the colored dots represent the seasonal average.863.5. Results and discussionReanalysis performance in the North region (Fig. 3.7) is representative ofthe drier climate zones (North, South Central, Central and Southeast; latterthree not shown). All reanalyses substantially overestimate the percentageof “wet” days, and the opposite for “dry” days (Fig. 3.7a, c; Table 3.2).MERRA-2 and JRA-55 outperform CFSR and ERA-Interim for two-sampleχ2-statistic (Fig. 3.7b; table 3.2). Reanalyses overestimate the occurrence of“very light” precipitation events (Fig. 3.7d-h), somewhat captures “light”,“moderate”, and “heavy” events (Fig. 3.7d-g), and underestimates “ex-treme” precipitation events (Fig. 3.7h) for the North Southeast climatezones. For the Central and South Central all precipitation type events arewell captured.Overall, across all climate zones and reanalyses, “very light” and “light”precipitation events are overestimated. This is expected as lower-resolutionmodels tend to overforecast such events; and extreme precipitation eventsare underestimated. Lower-resolution models and reanalyses typically arenot able to resolve extreme precipitation maxima, especially in the wettestzones Maritime West and Southwest.873.5.ResultsanddiscussionFigure 3.7: Percentage of (a) ”wet days”; (b) Two-sample χ2-statistic; Percentage of (c) ”dry days”, (d) “verylight”, (e) “light” (f) “moderate”, (g) “heavy” and (h) “extreme” precipitation events for North climate zone. Thevertical dashed lines indicate the change in seasons and the colored dots represent the seasonal average.883.5. Results and discussionLong-term, homogenized stations in mountainous BC are mostly locatedin valleys (Fig. 3.1). To identify regional and upper-elevation biases, eachreanalysis grid is bilinearly interpolated to the PRISM grid for comparisonto PRISM’s 30-year mean seasonal values. Similar to station-to-station com-parison, JRA-55 and ERA-Interim outperform CFSR and MERRA-2 (Fig.3.8). All four reanalyses show a dry bias during winter along the windwardand upper elevations of the Islands, Coast Mountains and Rocky Mountainswhere the Maritime West, Northwest and Southeast climate zones are lo-cated; and a wet bias along the Interior Plateau and leeward, lower elevationregions of Vancouver Island and Lower Mainland where the Maritime East,Southwest and Central climate zones are located. The North climate zonehas the lowest systematic error across the better reanalyses. Fall results aresimilar to Winter results (not shown). Summer has the lowest systematicerrors across the better reanalyses (JRA-55 and ERA-Interim). MERRA-2and CFSR exhibit a wet bias across the entire province (not shown). Finally,Spring has the largest systematic errors across all four reanalyses with wetbiases in the Maritime East, Southwest, Central and North climate zones,and dry biases across the Maritime West, Northwest and Southeast climatezones. CFSR followed by MERRA-2 shows a larger wet bias across most ofthe province (not shown).893.5.ResultsanddiscussionFigure 3.8: a) Winter precipitation totals of (a) PRISM. Winter precipitation totals bilinearly interpolated toPRISM grid of (b) CFSR, (c) ERA-Interim, (d) JRA-55 and (e) MERRA-2. Systematic error of (f) CFSR, (g)ERA-Interim, (h) JRA-55 and (i) MERRA-2. The dots represent weather station location.903.5. Results and discussionAn examination of the magnitudes of the biases of each metric at eachstation shows that reanalysis precipitation systematic error relative to weatherstation observations is highly correlated to reanalysis precipitation system-atic error relative to PRISM (PRISM and all four reanalyses are bilinearlyinterpolated to station locations for station-to-station comparison). PRISMwas developed to create a climatological precipitation on a regularly spacegrid that addresses spatial scales and patterns of orographic precipitation(Daly et al., 1994). This is an indication that the reanalysis biases in 31-day precipitation totals (shown for JRA-55 in Fig. 3.9) can be largelyexplained by topographic and synoptic parameters such as terrain steep-ness, exposure, elevation, location of barriers, and wind speed and direction,that are incorporated into PRISM. Strong correlations are also obtained forCFSR (0.78 ≤ R ≤ 0.94), ERA-Interim (0.79 ≤ R ≤ 0.93) and MERRA-2(0.77 ≤ R ≤ 0.91; not shown). Hence, the PRISM accurately representsstation observations, and can be used for bias correction and downscaling ofreanalyses.913.5.ResultsanddiscussionFigure 3.9: JRA-55 mean systematic error of 31-day precipitation totals for each of the 66 stations as a functionof PRISM mean systematic error. The solid lines show the linear regression fits.923.5. Results and discussion3.5.2 Extreme PCPThe Southwest climate zone has the highest population density, and reservoirsizes in the region are small relative to the magnitude of heavy and extremeprecipitation events. These two factors make it among the most sensitiveto extreme precipitation events. It is also one of the wettest regions —it has one of the highest 30-year return levels of 1-, 3-, 7-, and 14-dayaccumulated precipitation, second only to the Maritime West zone. Similarto the results of daily PCP for the wetter climate zones, the Southwest (Fig.3.10) and Maritime West climate zones show that all reanalyses are typicallytoo dry for extreme precipitation events, and too wet for the Northwest andMaritime East climate zones.933.5.ResultsanddiscussionFigure 3.10: 30-year return level and systematic error of 1-day precipitation total, for all four reanalyses forSouthwest climate zone. The vertical dashed lines indicate the change in seasons and the colored dots representthe seasonal average.943.5. Results and discussionMERRA-2 performs best in this regard, followed by ERA-Interim andJRA-55, and then the CFSR. The percent magnitude of the 1-day dry biasesin 30-year return levels are similar to those accumulated over 3, 7, and 14days (not shown; Table 3.2).For the wetter climates zones, the highest values of extreme precipitationoccur during storm season (October to February for southwest BC; Fig.3.10).These tend to be associated with non-convective, synoptic systems.For the drier climate zones the highest values of extreme precipitationoccur during summer. Biases are smaller and typically associated with thun-derstorm convection. In these zones, all reanalyses generally exhibit a drybias (e.g., North, Fig. 3.11), with the exceptions that MERRA-2 has a wetbias during the summer peak, and JRA-55 a near zero bias. Furthermore, allreanalyses have smaller biases when compared to the wetter climate zones(cf. Figs. 3.10, 3.11), indicating that 30-year return levels of of 1-, 3-,7-, and 14-day precipitation totals are fairly well captured all year long forall accumulation periods (not shown) in drier zones. This is notably differ-ent from the inability of the reanalyses to capture the annual cycle of dailyprecipitation in drier zones (Fig. 3.5). This is somewhat surprising sinceone might expect relatively coarse-resolution reanalyses to capture monthlyaccumulated precipitation more accurately than extreme convective events.953.5.ResultsanddiscussionFigure 3.11: 30-year return level and systematic error of 1-day precipitation total, for all four reanalyses for Northclimate zone. The vertical dashed lines indicate the change in seasons and the colored dots represent the seasonalaverage.963.5. Results and discussionA Kruskal-Wallis analysis indicates significant differences in mean sys-tematic error of daily and extreme PCP between the four reanalysis datasetsat the αWalker = 9.53×10−3 level. After multiple comparisons by Nemenyi’stest, for daily PCP the ERA-Interim, JRA-55 and MERRA-2 significantlyoutperform CFSR. For extreme PCP, the MERRA-2 and JRA-55 reanalysessignificantly outperform ERA-Interim and CFSR.All of these results are summarized in Fig. 3.12 where mean absoluteerror (MAE) of the reanalyses daily and extreme PCP are averaged overthe entire study period and all stations. The closer the value is to 0 fora given reanalysis, the better its performance. MERRA-2 and JRA-55 arethe better reanalyses, outperforming CFSR for all metrics, and to a lesserextent, ERA-Interim for daily PCP. This averaging also hides the greatervariability in bias of the poorer performing reanalyses, which is harder tocorrect for. For extreme PCP, the difference between MERRA-2 and JRA-55, and ERA-Interim is more noticeable with the former two outperformingthe latter. Additionally, the errors of 30-year return levels are smaller thanthose of 2-year return levels. Although MERRA-2 outperforms JRA-55 inmost metrics on average, the differences are not significant.973.5.ResultsanddiscussionFigure 3.12: MAE of 31-day precipitation totals, χ2−statistic, ”wet days”, and 2-year and 30-year return levels of1-, 3-, 7-, and 14-day precipitation totals. MAE is averaged across all 66 stations. Values closer to 0 at the originof the plot are better.983.6. Stationarity3.6 Stationarity3.6.1 Daily PCPStationarity in daily and extreme PCP distributions is assessed to deter-mine if temporal changes are significant. If there are significant temporalchanges, that would mean a traditional, stationary distribution (based onthe 1981-2010 period) would not be appropriate to represent the presentday expected precipitation distribution. First, variations through time —for each calendar day at each station — for season precipitation totals, num-ber of ”wet” days, “very light”, “light”, “moderate”, “heavy” and “extreme”precipitation events are modeled using linear regression to identify patterns.Due to large year-to-year background variability in precipitation, thestudy period is divided into three decades (1981-1990, 1991-2000 and 2001-2010). A 91-day centered rolling window is used to obtain seasonal precipi-tation totals for each calendar day; one value per decade per calendar day,yielding 3 values for each calendar day. The seasonal precipitation total fora given calendar day is assumed to vary with time as µ(t) = µ0 +µ1(t− t0),where the slope coefficient µ1 represents the average change from one decadeto the next (Fig. 3.13).Over the same 91-day centered rolling window, the number of ”wet” daysis assumed to vary with time as µ(t) = µ0 + µ1(t− t0), where µ1 representsthe 10-year rate of change of the number of wet days in the 91-day centeredrolling window. Namely, the number of ”wet” days is recorded per decadeper calendar day, yielding 3 values for each calendar day. The delineationof ”wet” days is lowered from 1.0 mm to trace because stationarity is beingassessed on weather station data only (excessive occurrences of very lightmodel precipitation are not an issue).Finally, still over the same 91-day centered rolling window ”wet” daysare recorded into five non-overlapping intervals: [Trace, 50th), [50th, 75th),[75th, 90th), [90th, 95th) and [95th, 100th]th for “very light”, “light”, “mod-erate”, “heavy” and “extreme” precipitation events, respectively. The num-ber of events in each bin is assumed to vary with time as µ(t) = µ0+µ1(t−t0),where the slope coefficient µ1 represents the average change in the numberof events from one decade to the next.Second, confidence intervals (CI) (Sun et al., 2018a) are implementedto assess significant changes of 10-year means from one decade to the nextof season precipitation totals, number of ”wet” days, “very light”, “light”,“moderate”, “heavy” and “extreme” precipitation events.A sufficiently large change in means from one decade to the next indicates993.6. Stationaritythere is a trend, and therefore nonstationarity is required to characterizedaily PCP. Small changes suggest a simpler, stationary model is accurateenough to represent precipitation. This critical difference is determinedby the (1 − αWalker) × 100% CI due to multiple testing, where αWalker =4.37× 10−6 and α0 = 0.10. Namely, variations in the mean over successivedecades are large enough to be considered statistically significant if they falloutside the lower and upper bounds of the (1−αWalker)% CI of the resultingdifference between the means of the 10-year periods and the 30-year period(Fig. 3.13). The individual 90% CI is adopted instead of the 95% CI toreduce the probability of making a type II error — that is, to reduce theprobability of failing to see that there is a trend.1003.6.StationarityFigure 3.13: Seasonal precipitation total at Vancouver International Airport (YVR) for July 9th (grey), 10-yearleft moving average (black), decadal averages of seasonal precipitation total (1981-1990, 1991-2000 and 2001-2010)(blue) and 30-year season precipitation total averaged over the study period 1981-2010 (green). The red line showsthe linear regression fit to the three decadal mean values. Individual 90% confidence intervals (CI; green dottedline) and multiple (1− αWalker)× 100% CI (red dotted line) are drawn; α0 = 0.10. This example indicates that astationary distribution is appropriate for Jul 9th because the three decadal average values all fall within both CI.1013.6. StationarityFigure 3.14a indicates a noticeable drying trend of seasonal precipi-tation totals across southern BC during spring and summer (late summerin South Central zone), and a weaker wet trend during late summer, falland early winter across most of BC. The weather stations are organized byclimate zones with North on top, followed by Central, Northwest, and south-ern regions Maritime West, Maritime East, Southwest, South Central andSoutheast at the bottom. Similarly, precipitation frequency (Fig . 3.14b)suggests an increase in the number of dry days during spring and summeracross most of BC (particularly for the South Central zone in summer),and a weaker increase in the frequency of ”wet” days during fall and win-ter across southern BC. None of the trends are significant at the multiple(1−αWalker)× 100% CI. At the individual 90% CI, the trends are generallynot statistically significant, with exceptions in the Northwest and MaritimeWest climate zone during summer.Furthermore, Figure 3.15a, b and c, indicates a weak increase in thenumber of “very light”, “light” and “moderate” precipitation events acrossMaritime West, East and Southwest BC during fall (and winter for “verylight” events), and a weak decrease in such events during spring and summer(substantially stronger for “very light” events in the South Central zone).Again they are not significant — none of the trends fall outside either themultiple (1 − αWalker) × 100% CI or the individual 90% CI. “heavy” (Fig.3.15d) and in particular “extreme” (not shown) precipitation events do notshow any clear pattern.1023.6. StationarityFigure 3.14: (a) Mean linear trend of season accumulated precipitation;(b) Mean linear trend of days with precipitation. The vertical dashed linesindicate the change in seasons and the horizontal dashed lines delineate fromtop to bottom the North, Central, Northwest, Maritime West, MaritimeEast, Southwest, South Central and Southeast climate zones respectively.Units are percent change over the 30-year period.1033.6.StationarityFigure 3.15: Mean linear trend of (a) “very light” , (b) “light”, (c) “moderate” and (d) “extreme” precipitationevents. The vertical dashed lines indicate the change in seasons and the horizontal dashed lines delineate fromtop to bottom the North, Central, Northwest, Maritime West, Maritime East, Southwest, South Central andSoutheast climate zones respectively. Units are percent change over the 30-year period.1043.7. ConclusionDespite some clear trends, none of the changes of 10-year means from onedecade to the next are statistically significant at the multiple (1−αWalker)×100% CI, and only a few isolated cases are statistically significant at theindividual 90% CI. During the 1981-2010 study period, a stationary 30-yearmean is accurate enough to represent mean values of season precipitationtotals, number of ”wet” days, “very light”, “light”, “heavy”, “heavy” and“extreme” precipitation events.3.6.2 Extreme PCPFor extreme PCP, the rolling 91-day centered rolling window is maintainedfor 1-, 3-, 7-, and 14-day precipitation totals. A GEV dresses the 31 an-nual maximum values for each calendar day by the method of maximumlikelihood. A nonstationary GEV distribution is compared, where only thelocation parameter is allowed to exhibit trend, with a stationary GEV dis-tribution with constant location, scale and shape parameters.The GEV distribution location parameter is assumed to vary with timeas µ(t) = µ0 + µ1(t − t0), where the slope coefficient µ1 represent the 10-year rate of change in location values of extreme 1-, 3-, 7-, and 14-dayprecipitation totals.None of the variations in the average of annual maxima over successivedecades fall outside either the (1 − αWalker) × 100% CI, or the 90% CI.Furthermore, no clear trend in extreme PCP is discernible (not shown).During the 31 years of the study period, a stationary GEV distributionwith model M0 = GEV (µ, σ, κ) is accurate enough to represent 1-, 3-, 7-,and 14-day precipitation totals.3.7 ConclusionReanalysis performance for daily and extreme precipitation (PCP) is evalu-ated across British Columbia (BC) during the 1980-2010 study period. Tocompare daily PCP among CFSR, ERA-Interim, JRA-55 and MERRA-2,the systematic error of 31-day precipitation total, wet days, and two-sampleχ2−statistic are calculated. To identify performance of extreme PCP, thesystematic error of 2- and 30-year return levels of 1-, 3-, 7- and 14- dayaccumulated precipitation are compared.In reanalyses, precipitation is generally better represented over areaswell-covered with accurate, complete, and coherent observations of all nu-merical forecast variables, which are used to correct the reanalyses. In rel-atively data-sparse areas such as BC (Hacker et al., 2003; Spagnol, 2005),1053.7. Conclusionreanalysis precipitation relies mostly on the underlying model output ratherthan observations (Dee et al., 2011; Lindsay et al., 2014). Namely, despitemodel estimation of precipitation being based on temperature and humid-ity information derived from the assimilated observations, approximationsused in the models representation of moist processes over data-sparse areasstrongly affect the quality and consistency of the hydrological cycle (Deeet al., 2011).The model generated-precipitation in CFSR traditionally has a wet bias(Saha et al., 2010). In this study, CFSR wasn’t the best reanalysis in anyof the metrics. Results suggest a wet bias across all climate zones in BCmatching previous results which indicate a wet bias over the western Pacificand in mid-high latitudes (Saha et al., 2010; Wang et al., 2011).JRA-55 has the most consistent, and among the smallest, systematicerror throughout the year and across the different climate zones in BC.Previous studies concluded that the quality of the JRA-55 improved signifi-cantly when compared with that of JRA-25 (Kobayashi et al., 2015; Haradaet al., 2016), especially in the Pacific Ocean north of 30◦ N (Harada, 2018).ERA-Interim exhibits the largest variation in performance throughoutthe calendar year and across the different climate zones, with lower system-atic errors across the drier climates zones than the wetter climate zones; cap-turing the wettest months in the dry climate zones, and missing the correctamount of precipitation during storm season in the wet climate zones. How-ever, ERA-Interim performed fairly well for daily PCP across BC. Uppalaet al. (2005) explains the various difficulties encountered in ERA-40 with theassimilation of humidity information, which led to a generally poor repre-sentation of the global transport of moisture in the atmosphere. Accordingto our results, those problems seem to have been corrected on ERA-Interim.Additionally, Uppala et al. (2005) concludes that there was an improvementover previous generation reanalysis ERA-40 to ERA-Interim in precipitationover higher latitudes.A previous MERRA study (Bosilovich et al., 2015) suggests that thesparse coverage of precipitation gauges in high latitudes may lead to sig-nificant biases. Studies have documented the difficulty of conserving at-mospheric dry mass while guaranteeing that the net source of water fromprecipitation and surface evaporation equals the change in the total atmo-spheric water (Trenberth and Smith, 2005; Bosilovich et al., 2008; Berrisfordet al., 2011). Reconsideration of these issues were taken into account duringthe development of MERRA-2. In this study, MERRA-2 performed well,particularly for extreme PCP.In summary, JRA-55 and MERRA-2 better capture precipitation distri-1063.7. Conclusionbution across BC all year, and have the lowest systematic errors across thewet climate zones during storm season. This makes them the better choicesfor a gridded climatological dataset of daily precipitation over BC. For dailyPCP, MERRA-2 and JRA-55 are the better reanalyses followed closely byERA-Interim. For extreme PCP, MERRA-2 and JRA-55 are the better re-analyses, with the lowest systematic errors throughout the year and acrossdifferent climate zones.According to Chapter 2 and this study, ERA-Interim performs betterfor daily and extreme 2-m temperature than it does for daily and extremePCP, even though the two fields should influence one another. A possibleexplanation is that many reanalyses do not directly assimilate 2-m air tem-perature observations, whereas ERA-Interim does. In contrast, PCP is amodel-produced field, influenced indirectly by surface and upper-air tem-perature and humidity observations (Dee et al., 2011; Lader et al., 2016).There is a noticeable drying trend in precipitation total during springand summer months across southern BC, and a wet trend during early fallfor northern and southwestern BC. These patterns also manifest themselvesin dry- and wet-day frequencies. The strongest signal is drying in the SouthCentral zone in summer. These findings add more information to previousstudies. Vincent and Mekis (2006) showed that the the number of dayswith precipitation per year also have significantly increased from 1950-2003across BC, and Zhang et al. (2000) showed a distinct drying pattern in thesouthern regions of BC during summer and spring during the second half ofthe twentieth century. Finally, our analysis shows that spring and summerhave been getting drier for much of BC, which is in line with future climateprojections (Haughian et al., 2012).The number of light and moderate precipitation events has generally in-creased during fall and winter months, and decreased mostly during springand summer across BC. Despite the clear evidence of a dry trend for springand summer months, and a wet trend during fall and winter months, thetrends are not as discernible for precipitation intensity. Finally, there is nodiscernible pattern in changes in return levels of extreme PCP, or frequencyof heavy and extreme precipitation events. A different study also showedno consistent trends in the number of precipitation extremes during the lastcentury (Zhang et al., 2001). By contrast, Groisman et al. (2005) showed anincrease in heavy and very heavy precipitation events south of 55◦N from1910-2001 across BC. The lack of consistency between periods and method-ology for computing the trends have made it difficult to compare resultsacross different studies. One possibility is that changes in precipitationwere occurring too slowly to be discerned in the 31-year study period, given1073.7. Conclusionthe the considerable year-to-year variability in precipitation. Therefore, itis possible that with a longer record discernible trends may be found. Forthis 31-year study period, apart from some isolated cases, no statisticallysignificant trends are found. Thus, a stationary distribution is sufficient torepresent daily and extreme PCP over BC.Chapter 2 shows that, across mountainous BC, ERA-Interim and JRA-55 are the most consistent and accurate reanalyses for daily and extremetemperature. This chapter shows that JRA-55 and MERRA-2 are the mostconsistent reanalyses for daily and extreme precipitation. More consistentbiases are favoured, as they are more easily removed by bias correction. Itis important to note that, as expected, the results for daily and extremetemperatures are more conclusive than the results for daily and extremeprecipitation since models do not simulate precipitation as well as they dotemperature (Kendon et al., 2014; Ravishankar et al., 2016), and have diffi-culties to represent extreme precipitation (Zhu et al., 2014). Due to highervariability in precipitation across BC, there is a large variation in perfor-mance —even for the better reanalyses JRA-55 and MERRA-2 – across thedifferent climate zones and seasons.Furthermore, the longer JRA-55 record is advantageous in that the stan-dard errors (of the estimated parameters used in extremes modelling and theresulting return levels) are expected to decrease as the sample size increases(Hosking et al., 1985; Hosking, 1990; Cai and Hames, 2010).Based on these findings, and the temperature findings in Chapter 2,the JRA-55 is recommended as the most accurate reanalysis over BC. Thischapter concludes that for daily PCP, the JRA-55 systematic error relativeto weather station observations is highly correlated to JRA-55 systematicerror relative to PRISM. It suggests that the biases can be explained by to-pographic and synoptic parameters— parameters that were implemented inthe development of PRISM. In the next chapter the bias corrections based onerror dependencies found in Chapter 2 and this study will be implemented inthe JRA-55 to create an even more accurate gridded climatological dataset.This will then be used, in conjunction with a probabilistic forecast datasetto create an extreme temperature and precipitation forecast index.108Chapter 4Analysis and Forecast ofHigh-resolution ExtremeWeather4.1 IntroductionThis Chapter has four objectives. First, the high-resolution PRISM clima-tology, the JRA-55 and the homogenized weather station dataset are com-bined to downscale and bias correct the JRA-55 and create a 30-arc-second(∼ 800 m) very-high-resolution surface analysis (VHRSA) of daily maximumand minimum 2-m temperature, and 1-day accumulated precipitation.VHRSA temperature and precipitation datasets have the potential ben-efit of rendering a feasible solution to the paucity of observational dataacross BC due to its inherent spatial and temporal completeness. Addition-ally, such a dataset has the ability to resolve fine-scale topographic featuresthat are important for a wide variety meteorological, climatological, andhydrological studies.The accuracy of the VHRSA is compared to the JRA-55 — the bestperforming reanalyses across BC. Performance is assessed with respect todaily maximum and minimum 2-m temperature, and 1-day accumulatedprecipitation over mountainous BC (hereafter daily and extreme T2M asdefined in section 2.3; daily and extreme PCP as defined in section 3.2).Then, the VHRSA is used to statistically downscale and bias correctNorth American Ensemble Forecast System (NAEFS) forecasts, creating avery-high-resolution probabilistic forecast across BC.Finally, trends in extreme T2M and PCP are examined to determinewhether or not a stationary climatological distribution is appropriate to rep-resent present-day extreme distributions. The stationarity results are usedto inform the creation of a new extreme weather forecast index. The objec-tive of the index to improve upon existing tools (that heighten operationalawareness of potentially extreme events) by more accurately calculating the1094.2. Dataextremity of the event, and producing fewer false alarms.In section 4.2, a brief description of the different datasets and of theweather station observations are given. In section 4.3, the methodology tostatistically downscale and bias correct the JRA-55 reanalysis, and the vari-ous metrics used to evaluate daily and extreme T2M and PCP are presented.In section 4.4, the methodology to statistically downscale and bias correctthe NAEFS forecast is presented. The metrics used to evaluate daily T2Mand PCP with respect to the raw NAEFS forecast are also presented. Insection 4.5, the methods for assessing statistical nonstationarity are intro-duced, and trends of extreme T2M and PCP are examined. Section 4.6introduces the Parametric Extreme Index (PEI) and compares it to Stan-dardized Anomalies (SA). Results are summarized in the conclusion.4.2 DataObservational data were obtained from 62 surface weather stations for dailymaximum and minimum T2M, and from 69 stations for daily accumulatedPCP, for the period 1 Jan 1958 to 31 Dec 2017. These data are used inconjunction with the PRISM dataset to bias-correct and downscale the JRA-55, which in turn is used to bias-correct and downscale the North AmericanEnsemble Forecast System (NAEFS) forecast.A detailed description of the T2M and PCP weather station datasetsare given in chapter 2 and 3, respectively. Information on the JRA-55 at-mospheric models and configurations are presented in tables 2.1 and 3.1. Amore detailed description of the JRA-55 can be found in Chapter 2 or inEbita et al. (2009, 2011) and Takeuchi et al. (2013).4.2.1 Weather station dataOf the 72 and 118 geographically-dispersed T2M and PCP stations initiallyselected for this study respectively, 10 T2M and 49 PCP stations with morethan 10% missing data were excluded. The upper bound threshold wasraised to 10% in this Chapter to strike a balance between having enoughstations with a complete temporal record, while still having representativedata in each climate zone in the province. A broad description of climatezones can be found in Chilton (1981) and Moore et al. (2008). A more de-tailed description and formal statistical derivation of T2M and PCP climatezones can be found in Chapters 2 and 3, respectively.Figures 4.1 and 4.2 shows the locations of all Environment and ClimateChange Canada (ECCC) and BC Hydro stations and their corresponding1104.2. DataFigure 4.1: Location of ECCC (blue) and BC Hydro (yellow) T2M weatherstations and British Columbia population distribution (red).three-letter abbreviations, overlaid with the population distribution acrossBC (these maps are shown again because they have a different selection ofstations from those shown in Chapters 2 and 3). Station elevation variesfrom sea-level to alpine, with most stations located in valleys with differentorientations, slopes and elevations, and some stations located on mountainslopes with different slope angles (see Appendix A for more details on sta-tions).For nonstationary analysis, stations with more than 1% missing dataare excluded, leaving eighteen ECCC T2M stations and eight ECCC PCPstations. BC Hydro stations were excluded due to apparently spurious tem-poral trends. A detailed discussion of station homogeneity can be found inChapter 2 for T2M and Chapter 3 for PCP.1114.2. DataFigure 4.2: Location of ECCC (blue) and BC Hydro (yellow) PCP weatherstations and British Columbia population distribution (red).1124.3. The VHRSA: Downscaling and bias-correcting the JRA-554.2.2 PRISMIn this study, gridded daily T2M and PCP climatology for the 1981-2010climate normal period are used to downscale and improve the accuracy ofthe JRA-55 by identifying high-resolution regional and terrain biases in theJRA-55 daily T2M and PCP. A more detailed description of the PRISM canbe found in Chapter 3 or in Daly et al. (1994, 1997, 2002).4.2.3 NAEFSThe North American Ensemble Forecast System (NAEFS) is the amalga-mation of the Meteorological Service of Canada (MSC) Global EnsemblePrediction System (GEPS) (Charron et al., 2010) and the United StatesNational Centers for Environmental Prediction (NCEP) Global EnsembleForecast System (GEFS) (Toth and Kalnay, 1997). The GEPS and GEFSare each comprised of 20 perturbed and 1 control members. Both ensem-ble systems employ ensemble Kalman filtering to generate perturbed initialconditions, and perturbed physical parameterizations (Wang et al., 2013;Hou et al., 2015; Wei et al., 2008; Zhou et al., 2017). When combined, theNAEFS provides a 42-member forecasts out to 16 days that is of higherquality than either ensemble alone (Zhu et al., 2014). The available archiveof NAEFS data was downloaded at a horizontal grid spacing of 1.0◦ × 1.0◦every 6 hours (0000, 0600, 1200 and 1800 UTC) out to 14 days.4.3 The VHRSA: Downscaling andbias-correcting the JRA-55The statistical behaviour of extreme T2M and PCP are compared betweenobserved weather station data, and the JRA-55 and PRISM interpolatedto the stations locations. As in Chapters 2 and 3, the Inverse DistanceWeighting (IDW) is used to interpolate reanalysis and PRISM output. Amore detailed description of each interpolation method can be found inMooney et al. (2011) and Stahl et al. (2006).For JRA-55, precipitation accumulation intervals are summed over the0601UTC-0600UTC window to get daily PCP (Table 3.1). The daily max-imum (minimum) T2M is defined as the highest (lowest) value of the six-hourly T2M outputs in the same window (Table 2.1).The highest (lowest) value of daily maximum (minimum) T2M withina 31-day centered rolling window is selected for each calendar day, for eachyear in the 60-year dataset. Thus, each calendar day has 60 values of annual1134.3. The VHRSA: Downscaling and bias-correcting the JRA-55maximum (minimum) T2M and PCP. This is done for both the station dataand JRA-55 interpolated to the stations.A Generalized Extreme Value distribution is fitted over the 60 samplevalues for each calendar day by the method of L moments. As in previouschapter, estimates of return levels derived by the method of L moments aremore reliable than the method of maximum likelihood (Hosking et al., 1985;Hosking, 1990).A Lilliefors test is conducted to evaluate the goodness-of-fit test of thefitted GEV distribution to the observed data. The null hypothesis that theobserved data is drawn from a GEV distribution is rejected in favour of thealternate hypothesis if a large difference between the fitted GEV cumula-tive distribution function (CDF) and the observation empirical cumulativedistribution function (ECDF) is detected.This large critical difference is determined by a parametric bootstrapprocedure. Namely, 100 samples of size 60 are generated from the fittedGEV distribution for each calendar day at each station. Next, 100 criti-cal differences are computed from the comparison of each generated sampleECDF and the fitted GEV CDF. Accepting or rejecting the null hypothesisat the 10% significance level is equivalent to comparing the 90th percentileof the resulting collections of critical differences with the actual differencebetween the fitted GEV CDF and the observation ECDF. At this 10% sig-nificance level, 10% of the tests are expected to exceed a critical difference,rejecting the null hypothesis. However, to reduce the probability of incor-rectly rejecting one or more of the true null hypotheses due to multipletesting, the global 10% significance level is regarded as significant. Namely,the null hypothesis is rejected when the critical difference between the fittedGEV CDF and the observations ECDF exceeds the global 90th percentile ofthe resulting collections of critical differences between the fitted GEV CDFand the generated ECDF of the samples (Wilks, 2016).At the global 10% significance level, less than 8% of the stations andcalendar days are rejected during the 1958-2017 study period. It indicatesextreme T2M and PCP values can be described by a GEV distribution.The 2-year return level is then estimated from the GEV distribution foreach calendar day and station.Next, the monthly mean PRISM values are assumed to be valid on the15th of each month. The twelve monthly values are linearly interpolated intime to obtain 365 calendar day values of monthly means (example shownfor maximum T2M, Fig. 4.3a).For the JRA-55, the same 31-day centered rolling window is used toaccumulate PCP values and obtain monthly precipitation for each calendar1144.3. The VHRSA: Downscaling and bias-correcting the JRA-55day. These 31-day accumulated values are then averaged over the 30-yearPRISM climate period (1981-2010). For T2M, the 31-day window is usedto obtain monthly mean values of daily maximum and minimum T2M foreach calendar day, which are also averaged over the 30-year PRISM climateperiod.Then, for each calendar day and variable, a linear regression summarizesthe relationship between the bias of the JRA-55 2-year return level valuerelative to the station 2-year return level value, and the bias of the monthlymean JRA-55 value relative to the monthly mean PRISM value (additivebias for T2M, multiplicative bias for PCP). This yields 365 regression modelsand 730 parameters — two parameters per calendar day (the slope and y-intercept; Fig. 4.3b). The 2-year return level is chosen because it is themedian of a GEV distribution and therefore a robust representation of thecentre of the data. Two-thirds of the stations are randomly selected to trainthe linear regression model and the remaining 1/3 of the stations are usedto test it (see Appendix B for a list of train and test stations). This crossvalidation is needed to test how well the model will generalize to grid pointsthroughout BC.Finally, for each calendar day, the JRA-55 is bilinearly interpolated (Fig.4.3e) to the PRISM grid (Fig. 4.3a) and the monthly bias between the JRA-55 and the PRISM is calculated (additive bias for T2M, multiplicative biasfor PCP; Fig. 4.3c).The VHRSA is generated by starting with the JRA-55 field (Fig. 4.3d),and bias correcting it applying the linear regression equation (Fig. ??b) tothe monthly mean bias (Fig. 4.3c). The final result is the VHRSA field(Fig. 4.3f). That is, the VHRSA value is computed by subtracting (multi-plying for PCP) the downscaling difference for T2M (ratio for PCP) fromthe JRA-55 value.The downscaling differences for T2M in Figure 4.3c are typically lessthan 0◦C in valleys and greater than 0◦C in mountainous regions, leadingto a dataset with larger vertical temperature differences. Similarly, for PCPthe downscaling ratio is typically greater than 1 in valleys and less than 1 inmountainous regions, leading to a dataset that is drier in valleys and wetterin ridges and upper elevation regions (see Figure 3.9 and subsection 3.5.1for more details).1154.3.TheVHRSA:Downscalingandbias-correctingtheJRA-55Figure 4.3: (a) 30-year monthly mean of daily maximum T2M PRISM climate on August 28th. (b) 2-year returnlevel systematic error between training stations and reanalysis interpolated to training stations as a function ofmonthly mean daily maximum T2M bias between JRA-55 and PRISM interpolated to training stations. The solidline shows the linear regression fit. (c) Downscaling difference derived by subtracting (a) from (e). (d) JRA-55daily maximum T2M on 28 Aug 2017. (e) Monthly mean of daily maximum T2M JRA-55 climate billinearlyinterpolated to PRISM grid. (f) VHRSA daily maximum T2M on 28 Aug 2017.1164.3. The VHRSA: Downscaling and bias-correcting the JRA-55To evaluate whether the VHRSA is statistically better as a griddedclimatological dataset than the JRA-55, the monthly mean bias (or sys-tematic error) and mean absolute error (MAE; or random error) are com-puted to estimate how accurately each dataset captures T2M and PCP. Thebias measures the average difference between the datasets and the observa-tion. The MAE measures the mean of the absolute differences between thedatasets and observations. The two-sample Kolmogorov-Smirnov (KS) andχ2-statistics are computed to estimate how accurately each dataset capturesthe distribution of T2M and PCP respectively. The two-sample KS statis-tic determines the largest absolute difference between the gridded datasetvalues of the T2M ECDF, and the observed ECDF values. The χ2-statisticmeasures the difference between the number of “very light” ([1.0 mm, 50th)),“moderate” ([50th, 75th)), “moderate” ([75th, 90th)), “heavy” ([90th, 95th))and “extreme” ([95th, 100th]) PCP events in each dataset.Finally, the bias and MAE of the 2-year return levels are calculated toestimate how well each gridded analysis captures extreme T2M and PCP.For more details on these metrics, refer to sections 2.4 and 3.4.Given that the test stations belong to different climate zones, the medianis used to average the monthly mean and 2-year return levels biases, MAE,KS and χ2-statistics across all test stations on each calendar day (Figures4.4, 4.5 and 4.6).4.3.1 Verification of the VHRSA Daily and Extreme T2MFor daily and extreme maximum T2M, the seasonal cycle of temperature isvery well captured and maintained over the VHRSA (Fig. 4.4a,b). Addi-tionally, the VHRSA clearly outperforms the JRA-55 throughout the entirecalendar year for both bias and MAE of daily and extreme T2M, and KS-statistic (Fig. 4.4c-g). The cold bias in the JRA-55 has mostly been removedfor both daily and extreme T2M (Fig. 4.4c,d), and the VHRSA better cap-tures the observations ECDF (Fig. 4.4g). The coloured points represent theseasonal error of each metric. They indicate that the errors in the VHRSAare consistently substantially smaller than those of the JRA-55 across allmetrics evaluated.The results for daily and extreme minimum T2M are nearly as good asthose for daily maximum T2M. Both the JRA-55 and the VHRSA capturethe seasonal cycle of daily and extreme minimum T2M well (Fig. 4.5a,b).The JRA-55 exhibits high variability in bias throughout the seasons for bothmean and 2-year return levels of daily minimum T2M (Fig. 4.5c,d). Thecold bias from mid-spring to mid-fall, and the warm bias during winter have1174.3. The VHRSA: Downscaling and bias-correcting the JRA-55Figure 4.4: (a) Observed, JRA-55 and VHRSA monthly mean daily max-imum T2M averaged over test stations. (b) 2-year return levels of ob-served, JRA-55 and VHRSA maximum T2M averaged over test stations.(c) Monthly mean bias of JRA-55 and VHRSA daily maximum T2M aver-aged over test stations. (d) As in (c) but for 2-year return levels. (e) As in(c) but for MAE. (f) As in (e) but for 2-year return levels. (g) KS statisticfor daily maximum T2M. The vertical dashed lines indicate the change inseasons and the coloured dots represent the errors seasonal averages. In(a) and (b), values closer to observations (black line) are better. In (c)-(g),values closer to zero are better. 1184.3. The VHRSA: Downscaling and bias-correcting the JRA-55been reduced in the VHRSA. The VHRSA consistently outperforms theJRA-55 for daily MAE, averaging about 2◦ C (Fig. 4.5e,f), similar to thatof the VHRSA maximum T2M. Since the VHRSA is corrected relative to2-year return values, the MAE is even lower here, around 1◦ C annually —again similar to that of the VHRSA maximum T2M. Finally, the VHRSAbetter captures the observed ECDF throughout the entire calendar year(Fig. 4.5g).4.3.2 Verification of the VHRSA Daily and Extreme PCPThe results for PCP, while not quite as good as those for T2M, still generallyshow substantial improvements after statistical downscaling and bias correc-tion. The seasonal cycle is well captured and the wet bias is reduced overthe VHRSA (Fig. 4.6a,b). The VHRSA consistently outperforms the JRA-55 throughout the seasons for monthly bias and MAE (Fig 4.6c,e). Duringthe most important part of the year, storm season, the VHRSA better cap-tures the observations distribution of precipitation (Fig 4.6g). Although thebias in extreme PCP is not improved (Fig. 4.6d), the VHRSA consistentlyoutperforms the JRA-55 for extreme PCP MAE (Fig. 4.6f). One possibleexplanation is that the biases across BC are cancelling each other resultingin a smaller bias but not a smaller MAE.The monthly mean systematic errors and MAE of daily maximum andminimum T2M and PCP, and the 2-year return level systematic errors andMAE are calculated from all calendar day systematic errors (Table 4.1).Similarly, the mean KS statistics and χ2-statistic are computed for dailymaximum and minimum T2M, and PCP respectively. Comparisons betweenthe fifteen mean systematic errors of the JRA-55 and the VHRSA are madeusing fifteen independent Mann-Whitney nonparametric tests (Hollanderet al., 2013). Fifteen independent Mann-Whitney tests are used becausethe MAE is inherently skewed as it is left-bounded at 0, and due to thedifferences in magnitude and variability of each type of error.A sufficiently large difference between the fifteen mean systematic errorsof the JRA-55 and the VHRSA indicate that daily and extreme T2M andPCP are more accurate in the VHRSA than in the JRA-55. Small differ-ences in the mean systematic errors suggest that the datasets have similaraccuracy. Due to the fact that there are fifteen error pair comparisons be-tween the JRA-55 and the VHRSA, the null hypothesis that the datasetshave similar accuracy is rejected at the αWalker = 1− (1−α0)1/N0 = 0.0034level of significance, where α0 = 0.05 and N0 = 15.The Mann-Whitney tests indicate that all errors of daily and extreme1194.3. The VHRSA: Downscaling and bias-correcting the JRA-55Figure 4.5: (a) Observed, JRA-55 and VHRSA monthly mean daily min-imum T2M averaged over test stations. (b) 2-year return levels of ob-served, JRA-55 and VHRSA minimum T2M averaged over test stations.(c) Monthly mean bias of JRA-55 and VHRSA daily minimum T2M aver-aged over test stations. (d) Same as in (c) but for 2-year return levels. (e)Same as in (c) but for MAE. (f) Same as in (e) but for 2-year return levels.(g) KS statistic for daily minimum T2M. The vertical dashed lines indicatethe change in seasons and the coloured dots represent the errors seasonalaverages. Values of bias, MAE, and KS Statistic closer to zero are better.1204.3. The VHRSA: Downscaling and bias-correcting the JRA-55Figure 4.6: (a) Observed, JRA-55 and VHRSA mean of monthly PCP totalsaveraged over test stations. b) 2-year return levels of observed, JRA-55 andVHRSA PCP averaged over test stations. (c) Monthly mean bias of bias-corrected JRA-55 PCP averaged over test stations. (d) Same as in (c) butfor 2-year return levels. (e) Same as in (c) but for MAE. (f) Same as in(e) but for 2-year return levels. g) χ2-statistic of PCP. The vertical dashedlines indicate the change in seasons and the coloured dots represent theerrors seasonal averages. Values of bias, MAE, and χ2-statistic closer tozero are better.1214.3. The VHRSA: Downscaling and bias-correcting the JRA-55Table 4.1: Averaged systematic error and MAE of monthly precipitationtotal and monthly mean daily maximum and minimum T2M (MSE andMMAE respectively), of 2-year return levels (SE2 and MAE2 respectively),averaged systematic error of two-sample χ2-statistic (χ2) and KS statistics,across all test stations (all T2M errors but KS statistic in ◦C; all PCP errorsbut two-sample χ2 -statistic in % )MSE MMAE SE2 MAE2 KS χ2T2M MAXVHRSA 0.32 1.87 0.12 0.76 0.11JRA-55 -3.91 4.21 -4.17 4.17 0.39T2M MINVHRSA -0.21 1.83 0.53 1.02 0.12JRA-55 -0.57 2.59 0.18 2.62 0.28PCPVHRSA 12.26 15.84 -9.94 12.26 16.62JRA-55 30.59 34.63 5.72 23.79 20.58T2M and PCP of the VHRSA are smaller than those of the JRA-55 atthe global 5% significance level, with the exception of extreme PCP bias asshown in Figure 4.6d).In general, the results of maximum T2M are better than those of mini-mum T2M, which in turn are better than the results of PCP. One possibleexplanation is the coefficient of determination R2. The R2 is a standardmeasure of the goodness-of-fit of a regression. It can be interpreted as theproportion of the variation of the bias of the JRA-55 2-year return level valuerelative to the station 2-year return level value that is described or accountedfor by the regression (Fig. 4.3b). Higher R2 values are obtained for maxi-mum T2M (0.62 ≤ R2 ≤ 0.91) than for minimum T2M, (0.69 ≤ R2 ≤ 0.85)which in turn are higher than the R2 values of PCP (0.58 ≤ R2 ≤ 0.81)throughout the year.Finally, looking at the different climate zones across all variables, resultsare very similar, there is always a clear improvement in the VHRSA (notshown). For instance, the large cold systematic errors of daily and extrememaximum T2M across the Southeast climate zone found in Chapter 2 (Ta-bles 2.2 and 2.3) have been removed, and the large wet biases in systematicerrors of daily PCP across the drier climate zones found in Chapter 3 (North,1224.4. Bias-corrected NAEFSSouth Central, Central and Southeast; Table 3.2) have been reduced.4.4 Bias-corrected NAEFSA rolling window of the 30 most recent values (as would be used in real-timeforecasting) is used to bias correct and downscale each member’s daily maxi-mum and minimum T2M, and 1-day accumulated PCP forecast towards theVHRSA. Given that the daily maximum and minimum T2M were not avail-able across in the NAEFS archive used, the daily maximum (minimum) T2Mis defined as the highest (lowest) value of the six-hourly T2M outputs (forthe calendar day 0601-0600 UTC). The 1-day accumulated PCP is summedover the four six-hourly outputs for each day.The NAEFS forecast Ft is bias corrected using an additive (multiplica-tive) degree-of-mass-balance bias-correction factor DMBt for T2M (PCP)(Grubiˇsic´ et al., 2005; McCollor and Stull, 2008a; Bourdin et al., 2014).First, the raw NAEFS forecast (∼ 100-km horizontal grid spacing) is bilin-early interpolated to the VHRSA grid (∼ 800-m grid spacing). Then, thebias-correction factor is calculated and updated daily, and applied to theraw interpolated NAEFS forecast to generate a bias-corrected, downscaledNAEFS forecast Fˆt. The factor DMBt is a combination of the previous-dayDMBt−1 and the difference (ratio for PCP) of the previous day forecast-observation pair (Ft−1−Ot−1 for T2M; Ft−1/Ot−1 for PCP), weighted by atime parameter τ = 30:DMBt =τ − 1τDMBt−1 +1τ(Ft−1 −Ot−1). (4.1)The factor functions as a 30-day exponentially decaying rolling window.Namely, the influence of the forecast-observation pairs decreases with ane-folding time of τ = 30 days from most recent to least recent. The valueof τ = 30 is chosen because McCollor and Stull (2008a) have shown ithas the optimal minimum bias between the bias-corrected forecast and theobservation for T2M and PCP forecasts over BC. The bias correctionFˆt = Ft −DMBt (4.2)is done separately for each NAEFS member and then averaged to generatethe bias-corrected ensemble mean (Fˆt = Ft/DMBt for PCP; Figures 4.7,4.8 and 4.9).The degree-of-mass-balance method is chosen for three reasons. (1) Dueto the 30-arc-second (∼ 800 m) grid resolution of the VHRSA, an algorithm1234.4. Bias-corrected NAEFSis needed to improve the forecast in an operational setting (with associatedtime constraints) without exceeding computing resources. (2) This biascorrection has been well studied, used and proven to work well. (3) Itefficiently updates daily and retains a memory without creating an ever-growing training dataset.Typically, operational raw ensemble forecasts are found to exhibit a lowensemble spread (Buizza, 1997; Hamill, 2001; Hamill et al., 2008; McCol-lor and Stull, 2008b), which in turn leads to overconfidence in probabilityassessment. Namely, ensemble forecasts are underdispersive and produceuncalibrated probabilistic forecasts (that is, the empirically derived proba-bilities are not accurate). A Nonhomogeneous Gaussian regression (NGR) isused to calibrate the forecast Fˆt at the stations for T2M and PCP (Gneitinget al., 2005; Hagedorn et al., 2008). The NGR assumes the forecast errorsform of a Gaussian distribution dressed about the ensemble mean:Nt(F¯t, s2ens). (4.3)Here, F¯t is the bias-corrected ensemble mean calculated from 4.1 and4.2, where each member is equally weighted. The variance of the residuals— which are assumed to be Gaussian distributed— is computed as a linearfunction of the ensemble variance s2ens. The regression parameters are chosento minimize the continuous ranked probability score (CRPS).Previous studies have reported good results calibrating surface temper-ature error distributions (which are approximately Gaussian) using NGR(Hamill et al., 2008; Hagedorn et al., 2008). Thorarinsdottir and Gneiting(2010) extend NGR to handle non-Gaussian distributed variables such aswind speeds. Once again, this method has been chosen to strike a balancebetween computing demands and well reported results.The bias-corrected, calibrated NAEFS is evaluated for daily PCP andminimum T2M during fall (defined here as October, November and De-cember) and winter (defined here as December, January and February) of2016/17, which were abnormally wet and cold seasons (e.g., Figs. 4.7 and4.8); and daily maximum T2M during summer of 2017 (defined here as June,July and August) , which was abnormally hot (e.g., Fig. 4.9). For more de-tails on the severity and impacts refer to subsection 1.1. The month priorto each season is used to spin up the bias correction and calibration of theraw forecast. Again, a 30-day training period is chosen to strike a balancebetween computing demands and reported results. According to McCollorand Stull (2008a), the DMB requires a 40 day training period for precipi-tation and a shorter 14 day training period for temperature to accomplish1244.4. Bias-corrected NAEFSbias reduction. Hagedorn et al. (2008) and Hamill et al. (2008) suggest NGRcan calibrate 2-m temperature and lighter precipitation events with a 30-day training dataset. A total of 3588, 4140, 4416 observations are used toevaluate forecast performance during fall, winter and summer respectively.1254.4.Bias-correctedNAEFSFigure 4.7: (a) NAEFS raw Canadian control member forecast of daily minimum T2M on previous day 2 Jan,2017; (b) Same as in (a) but for VHRSA; (c) Ensemble mean raw NAEFS forecast on current day 3 Jan, 2017; (d)Ensemble mean bias-corrected NAEFS forecast on current day 3 Jan, 2017. DMB method done separately for eachNAEFS member and then averaged to generate the bias-corrected ensemble mean. (e) VHRSA daily minimumT2M on 3 Jan, 2017 presented for comparison with (c) and (d). The dots represent weather station location.1264.4.Bias-correctedNAEFSFigure 4.8: (a) NAEFS raw Canadian control member forecast of 1-day accumulated PCP on previous day 7 Nov,2016; (b) Same as in (a) but for VHRSA; (c) Ensemble mean raw NAEFS forecast on current day 8 Nov, 2016;(d) Ensemble mean bias-corrected NAEFS forecast on current day 8 Nov, 2016. DMB method done separatelyfor each NAEFS member and then averaged to generate the bias-corrected ensemble mean. (e) VHRSA 1-dayaccumulated PCP on 8 Nov, 2016 presented for comparison with (c) and (d). The dots represent weather stationlocation.1274.4.Bias-correctedNAEFSFigure 4.9: (a) NAEFS raw Canadian control member forecast of daily maximum T2M on previous day 27 Aug,2017; (b) Same as in (a) but for VHRSA; (c) Ensemble mean raw NAEFS forecast on current day 28 Aug, 2017;(d) Ensemble mean bias-corrected NAEFS forecast on current day 28 Aug, 2017. DMB method done separatelyfor each NAEFS member and then averaged to generate the bias-corrected ensemble mean. (e) VHRSA dailymaximum T2M on 28 Aug, 2017 presented for comparison with (c) and (d). The dots represent weather stationlocation.1284.4. Bias-corrected NAEFS4.4.1 Verification of the ForecastMost forecast verification metrics have weaknesses (Murphy and Winkler,1987). In order to obtain an informed picture of the skill of the forecast,different forecast attributes need to be analyzed using different metrics.The bias (or systematic error) simply measures the difference betweenthe average forecast and average observation (additive for T2M, multiplica-tive for PCP), and therefore expresses the bias of the forecasts. A positivebias indicates that forecasted value was larger than the observed, which iscalled overforecasting. Conversely, a negative bias indicates that the eventwas underforecast, namely, the forecast value was less than the observed.Metrics based on a 2 × 2 contingency table are also used for T2M andPCP verification (Table 4.2). Stephenson (2000) suggests the hit (H) andfalse alarm (F) rates should be used in combination with the threat score(TS). This is helpful for extreme T2M and PCP values because the TS givesa more reasonable comparison between H and F due to the large number ofcorrect rejections ignored.Forecast ObservedYes NoYes Hit (a) False alarm (b)No Miss (c) Correct rejection (d)Table 4.2: Contingency table.H is the ratio of correct forecasts to the number of times the event hasbeen observed:H =aa+ b. (4.4)F is the ratio of false alarms to the total number of nonoccurrences ofthe event:F =bb+ d. (4.5)TS is the number of hits divided by the total number of occasions onwhich that event was either forecast and/or observed:TS =aa+ b+ c. (4.6)1294.4. Bias-corrected NAEFSThe worst possible forecast has high systematic error, H=0, F=1 andTS=0. Conversely, the best possible forecast has systematic error equalsto zero, H = 1, F = 0 and TS = 1. In order to evaluate higher impactT2M and PCP events, the 90th percentile of daily maximum T2M and dailyPCP, and the 10th percentile of daily minimum T2M are calculated within a31-day centered rolling window for each calendar day. A hit is issued whenboth the observed and the forecast ensemble mean values are above/belowthat calendar day 90th/10th percentile on any given day.For probabilist forecast verification, the quantile Brier Score (QBS) isused to measure the reliability of the ensemble 90th percentile of daily T2Mand daily PCP (Bentzien and Friederichs, 2014). The QBS is connected tothe reliability diagram. A reliable or well-calibrated forecast with QBS =0 has points falling closer to the 1:1 perfect-reliability line. That is, if theobservations fall at or below the 90th percentile forecast 90% of the time,then the probabilistic forecast is reliable for that percentile.For overall reliability of any T2M and PCP quantile, Candille et al.(2007) suggest the bias score (b) is combined with the dispersion score (d).The scores b and d are connected to the rank or PIT histogram (Anderson,1996; Talagrand et al., 1999). A perfectly reliable forecast with b = 0 anddispersion equal to 1 (d = 1) has a flat rank histogram. A large negative(positive) value of b indicates a negative (positive) bias, where a PIT his-togram would slope down towards the right (left). A value of d greater(smaller) than 1 characterizes an underdispersed (overdispersed) ensemblewith U-shaped (bell-shaped) PIT histogram.Finally, the CRPS is analyzed to evaluate forecast sharpness. Forecaststhat are frequently much different from climatology are sharp forecasts. Thebest possible forecast has CRPS=0.A bootstrap procedure calculates the 95% confidence interval for eachmetric. A metric is estimated from 100 generated samples of the ensembleforecast-observation pairs (Candille et al., 2007). Then, the 2.5th and 97.5thpercentiles of the resulting collection of metrics are used as the lower andupper bounds of the 95% confidence intervals for the true metric value.The 95% confidence interval measures the uncertainty around the metricand determines significant statistical differences between the raw and thepost-processed forecasts.4.4.2 Verification of Daily PCPFor daily PCP, the post-processed NAEFS exhibits a statistically signifi-cantly better (higher) TS and H than the raw NAEFS out to a forecast lead1304.4. Bias-corrected NAEFStime of 7 days (Fig. 4.10a,b). However, the post-processed NAEFS alsoproduces more F than the raw NAEFS (Fig. 4.10c). This is likely due togeneral inability of the coarse-resolution NAEFS to produce more extremevalues. As expected, as lead time increases, TS and H worsen (decrease)and F improves (decreases). One possible explanation is that the number ofmisses increase as the systematic error increases with increasing lead time(Fig. 4.10d). Further, ensemble member forecast solutions become morerandom in nature as lead time increases and the ensemble mean typicallytrends towards climatology, decreasing the chances of forecasting (H or F)a 90th event.The systematic error of the post-processed NAEFS is also statisticallysignificantly smaller than those of the raw NAEFS as far out as day 9 (Fig4.10d), indicating the efficacy of the bias correction. As expected, the 95%confidence interval widens with lead time across all metrics. It suggests thepost-processed NAEFS has statistically significantly less bias than the rawNAEFS as far out as day 7 across most metrics. Beyond this lead time,although the post-processed NAEFS still outperforms the raw NAEFS, theuncertainty around the metric is too high to determine whether statisticallysignificant differences exist.The post-processed NAEFS exhibits a much lower QBS, b, d and CRPSvalues (Fig. 4.10e-h) than the raw NAEFS. This indicates a statisticallysignificantly more reliable, calibrated and sharp forecast. The combined band d scores (Fig. 4.10e,f) match the results of the QBS (Fig. 4.10g). Thepost-processed NAEFS is still underdispersed. Namely, too many observa-tions fall in the low and high percentiles. Finally, the CRPS indicates theCDF of the post-processed NAEFS is sharper than the raw NAEFS out today 4, suggesting forecast are frequently much different than climatology.Although the downscaled NAEFS is statistically significantly more reli-able, calibrated and sharp that the raw NAEFS, there is room for improve-ment. QBS, b and d values remain high, suggesting a different calibrationmethod may lead to a more skillful forecast.A Gamma-distributed calibration method was also attempted for PCP(Scheuerer and Hamill, 2015; Baran and Lerch, 2016). Although the versa-tility in shape of the gamma distribution makes it an attractive candidatefor representing precipitation error distributions, the results are unsatisfy-ing in an operational setting due to erratic results. A gamma distributionis more difficult to work with than a Gaussian distribution, and obtainingparameter estimates from data is not as straightforward Wilks (2011). Moreresearch is needed in this area.1314.4. Bias-corrected NAEFSFigure 4.10: a) TS for events above the 95th percentile during Fall 2016; (b)Same as in a) but for H; (c) Same as in a) but for F ; d) Systematic errorof daily PCP; (e) Same as in a) but for b; (f) Same as in a) but for d ; (g)Same as in a) but for QBS; (h) Same as in a) but for CRPS. Values of TSand H closer to one are better. Values of F and systematic error closer tozero are better. Combined values of b = 0 and d = 1 are better. Values ofQBS and CRPS closer to zero are better.1324.4. Bias-corrected NAEFSFigure 4.11: a) TS for events above the 90th percentile during Summer 2017;(b) Same as in a) but for H; (c) Same as in a) but for F ; d) Systematicerror of daily maximum T2M. Values of TS and H closer to one are better.Values of F and systematic error closer to zero are better.4.4.3 Verification of Daily Maximum and Minimum T2MFor daily maximum T2M, the post-processed NAEFS exhibits a statisticallysignificantly better (higher) TS and H than the raw NAEFS out to a forecastlead time of 7 days (Fig. 4.11a,b), at which point the uncertainty aroundthe metric is too high to determine significant differences in skill. Similar todaily PCP, the post-processed NAEFS also produces more false alarms (F)than the raw NAEFS (Fig. 4.11c).The systematic error of the post-processed NAEFS is also statisticallysignificantly smaller than those of the raw NAEFS as far out as day 10 (Fig.4.11d).Finally, the post-processed NAEFS exhibits a much better QBS, b, d,and CRPS than the raw NAEFS for daily maximum T2M (not shown),suggesting a more reliable, calibrated and sharp forecast.The results for daily minimum T2M are not nearly as good as those fordaily maximum T2M and PCP for TS, F and H. This is largely because the1334.5. Stationarityraw NAEFS daily minimum T2M have much less systematic error than theraw NAEFS daily maximum T2M (higher H, TS and F values; not shown),so it is not as easy for the bias correction to improve upon them. TheVHRSA post-processed NAEFS exhibits lower TS, H and F than the rawNAEFS. However, the systematic error of daily minimum T2M is statisti-cally significantly better (lower) than those of the raw NAEFS as far as day10 (Fig 4.12a).The values of QBS, b, d, and CRPS of daily minimum T2M are of similarmagnitude to those of daily maximum T2M (Fig. 4.11b-e). This indicatesa statistically significantly more reliable, calibrated and sharp forecast. Thecombined b and d scores (Fig. 4.11c,d) match the results of the QBS (Fig.4.11b). The post-processed NAEFS is slightly underdispersed. Namely, afew observations fall in the low and high percentiles. Finally, the CRPSindicates the CDF of the post-processed NAEFS is more accurate than theraw NAEFS out to day 7. As expected, as the lead time increases, the QBSand the CRPS increase as the ensemble mean and distribution becomes lessreliable and sharp respectively. b increases in concert with systematic bias,indicating that the mean and distribution follow similar trends. The prob-abilistic spread goes from being too narrow to well-calibrated, as indicatedby d — the raw and calibrated spreads increase quickly with lead time.The downscaled NAEFS is statistically significantly more reliable, cali-brated and sharp that the raw NAEFS for both T2M and PCP. However,QBS, b and d values for T2M suggest the NGR calibration method providesbetter results for T2M which is approximately Gaussian distributed.4.5 StationarityMany studies using either homogenized station datasets or gridded GeneralCirculation Models (GCM) have shown that climate has undergone changeson a multidecadal time-scale (Mekis and Vincent, 2011; Kharin and Zwiers,2000; Zwiers and Kharin, 2005; Zhang et al., 2001; Groisman et al., 2005;Odon et al., 2018).Analysis of multidecadal trends of extremes are an important aspect ofclimate research. Changes in the magnitude and frequency of extremes haveenvironmental and socioeconomical consequences. It is therefore of greatinterest to evaluate changes in extreme T2M and PCP in the VHRSA sinceclimate trends at the regional scale are not easy to detect in a station-basedstudy, as it was done in Chapters 2 and 3, due to large natural variabilityin regional climate (Dettinger et al., 1998; Odon et al.).1344.5. StationarityFigure 4.12: a) Systematic error of daily minimum T2M during Winter2016/17; (b) Same as in a) but for QBS; (c) Same as in a) but for b ; (d)Same as in a) but for d ; (e) Same as in a) but for CRPS. Values of systematicerror, QBS and CRPS closer to zero are better. Combined values of b = 0and d = 1 are better.1354.5. StationarityGenerally, T2M is Gaussian distributed while PCP is not. Accordingly,parametric and nonparametric methods were presented in Chapters 2 and3, respectively. In this Chapter, nonparametric methods are presented asthey are better suited to compare T2M and PCP on equal footing.In order to detect spatial climate signals, a nonparametric robust regres-sion is used in the VHRSA to determine the trends in extreme T2M andPCP. This method has the following advantages when compared to the lin-ear regression used in Chapter 2: 1) the regression is robust against outliersand 2) the data doesn’t need to be Gaussian distributed.The trend is assumed to vary with time as µ(t) = µ0+µ1(t−t0), where theslope coefficient µ1 represents the annual rate of change in extreme T2M andPCP. The magnitude of µ1 is calculated by Theil-Sen single median (Theil,1950) method which computes the slopes of all possible combinations of pairsof the 60 sample values. After calculating these 1475 slopes, the median istaken as the magnitude of µ1.Such an analysis is computationally expensive at every grid point in theVHRSA. Hence, the trend is performed only on the 15th of each monthwhich is assumed to represent the trend of each month.The trend statistical significance is determined using the nonparametricMann-Kendall test (Mann, 1945; Wilks, 2011) which has been applied inhydroclimate trend studies (Kumar et al., 2013).The Mann-Kendall test isrobust against outliers, independent of the T2M or PCP distribution, andprovides a more powerful analysis for non-Gaussian distributed data suchas PCP (Yue et al., 2002; Onoz and Bayazit, 2003).The null hypothesis µ1 = 0 (no trend) is tested against the alternatehypothesis µ1 6= 0 that there is a trend at the αFDR = 0.10 level of signifi-cance. Namely, a sufficiently large change in extremes throughout the studyperiod indicates there is a trend, and therefore nonstationarity is requiredto characterize extreme T2M and PCP. Because the Mann-Kendall test re-quires serially independent data, the effective sample size is used rather thanthe actual sample size.Finally, the rolling 31-day centered rolling window is maintained anda GEV dresses these 60 annual values for the 15th of each month by themethod of maximum likelihood. A nonstationary GEV distribution is com-pared, where only the location parameter is allowed to exhibit trend, witha stationary GEV distribution with constant location, scale and shape pa-rameters. The GEV distribution location parameter is assumed to vary withtime as µ(t) = µ0 + µ1(t− t0), where the slope coefficient µ1 represents theannual rate of change in extreme T2M and PCP.With models GEV (µ(t), σ, κ) and GEV (µ, σ, κ), the alternate hypothe-1364.5. Stationaritysis is tested against the null hypothesis that extreme values of T2M and PCPare drawn from the same GEV distribution using a LRT at the αFDR = 0.10level of significance. Namely, the alternate hypothesis suggests that a non-stationary model explains more of the variation in the time series and thatconsequently changes in return levels should be accounted for. For moredetails on this methodology, refer to subsections 2.5.2 and 3.6.2.4.5.1 Extreme T2MFigures 4.13a,b indicate a statistically significant warming trend of extremevalues of minimum T2M across most of BC during July at the αFDR =0.10 significance level, with exceptions in the Central climate zone. Similarresults are obtained during January, May, June and August (not shown).A non-statistically-significant warming trend of extreme minimum T2M ispresent throughout the rest of the months across all of BC (not shown).Under models GEV (µ, σ, κ) and GEV (µ(t), σ, κ), the LRT suggests theevidence supporting a trend in July is strong, implying a nonstationarymodel brings significant improvements over a stationary model at the αFDR =0.10 significance level (Fig. 4.13c). Namely, during the 60 years of the studyperiod, changes in return levels of extreme minimum T2M significant, andnonstationarity is required to characterize extreme levels. Again, similarresults are obtained during January, May, June and August (not shown).1374.5.StationarityFigure 4.13: (a) Trend of extreme minimum T2M in July; (b) Locations where trend is statistically significant;(c) Locations where GEV (µ(t), σ, κ) for extreme minimum T2M is statistically significant.1384.6. The Parametric Extreme IndexAnalogously, figures 4.14a,b indicate a statistically significant warmingtrend of extreme values of maximum T2M across the South Central climatezone throughout the year (shown for August). None of the LRT are statis-tically significant (not shown). A stationary GEV distribution with modelGEV (µ, σ, κ) is accurate enough to represent extreme maximum T2M.4.5.2 Extreme PCPNone of the trends in extreme PCP are statistically significant. Furthermore,no clear trend in extreme PCP is discernible across the months (not shown).During the 60 years of the study period, a stationary GEV distribution withmodel GEV (µ, σ, κ) is accurate enough to represent extreme PCP.4.6 The Parametric Extreme IndexA new Parametric Extreme Index (PEI) is herein developed and presented,with the goal of alerting forecasters of extreme events. It does this by pro-viding a concise, single index value derived from comparing the downscaledNAEFS forecast cumulative distribution function (hereafter NCDF) to theGEV cumulative distribution function (derived from the VHRSA; hereafterGCDF). Namely, on any given day, the NCDF is compared to the GCDFfor a given location and calendar day. The 31-day rolling window approachis useful for operational forecasting because it allows for assessment of ex-treme weather on days within a similar climate. For instance, comparingthe true annual maximum T2M (which likely occurs in July) to a winter-time NCDF is of no use since it has no chance of occurring during the wintermonths. Finally, the GCDF is a distribution of climatological extremes only,which in turn means that the PEI will alert only if the NCDF forecasts anextreme event, rather than just anomalous events, which occur frequently.The PEI ranges from 0 (no forecasted chance of extreme weather) to 1 (highforecasted probability of very extreme weather).Figure 4.15 illustrates the NCDF and GCDF for the PEI calculation on28 August 2017 at one grid point. The PEI value is low because there is alow probability of exceeding low extreme values. For reference, the 2-yearreturn level, which represents the median of the GCDF, is far to the right.The PEI is calculated according to the formulaPEI =12∫ Fc(F−1f (1))01− Ff (F−1c (pc))√1− pc dpc, (4.7)1394.6. The Parametric Extreme IndexFigure 4.14: (a) Trend of extreme maximum T2M in August; (b) Locationswhere trend is statistically significant.1404.6. The Parametric Extreme Indexwhere pc is the probability that an extreme value (F−1c (pc)) will not be ex-ceeded in the extreme climate record, the inverse CDF F−1f and F−1c aresimply the quantile functions of the forecast and of the extreme climate dis-tributions respectively, and 1−Ff (F−1c (pc)) denotes the complement of theNCDF. Namely, if Ff (F−1c (pc)) represents the NAEFS-derived probabilitythat the extreme value F−1c (pc) will not be exceeded, 1− Ff (F−1c (pc)) rep-resents the probability that the value will be exceeded. pc = 0 is the lowerintegral limit. It represents the lowest extreme value in the GCDF (F−1c (0)),or the 0th percentile. pc = Fc(F−1f (1)) is the upper integral limit. It repre-sents the probability that the highest forecast value issued F−1f (1) will notbe exceeded in the extreme climate distribution. In essence, the integral isbounded by the lowest extreme value in the GCDF and the highest forecastissued in the NCDF. The calculations are made so that more weight is givento the right tail of the distributions (where rarer, more extreme values are);that is, as pc −→ 1,√1− pc −→ 0 and 1√1−pc −→∞.If the NCDF lies below the lowest extreme value in the GCDF (F−1c (0)),the PEI computes its lowest value: PEI = 0. As the magnitude of thehighest forecast value increases (F−1f (1) −→∞) and the less frequent, moreextreme that value is in the GCDF (Fc(F−1f (1)) −→ 1), the PEI approachesits highest value: PEI = +1.Figure 4.16 shows the PEI during an arctic outbreak on 3 January 2017(VHRSA for this date shown in Fig. 4.7e). Since the event happened duringJanuary, the GCDF is given by GEV (µ(t), σ, κ). Namely, a nonstationarydistribution is used to characterize extreme levels. According to the resultsof section 4.5, current extreme levels are warmer than a stationary distribu-tion based on the 60-year record would indicate. If a stationary distributionhad been used, the PEI would yield a lower index (alert) value, as the ex-treme minimum T2M levels would have been colder and therefore harderto reach (not shown). Figures 4.17 and 4.18 illustrate the PEI during anatmospheric river event on 8 November 2016 (Fig. 4.8e), and a heat wave on28 August 2017 (Fig. 4.9e). Both days had large socioeconomical impacts(Odon et al., 2017) (see Chapter 1 for more details).4.6.1 Verification of the PEIExisting situational/extreme awareness tools, such as the Standardized Anomaly(SA) and EFI, do not handle extreme values in a proper manner, somethingthat this dissertation has sought to improve upon. Previous methods detectanomalous weather events by comparing them with historical means and/or1414.6. The Parametric Extreme IndexFigure 4.15: (a) The NCDF (purple) and the GCDF (black) show the prob-ability (y-axis) that the extreme value (F−1c (pc); x-axis) will not be ex-ceeded. (b) The NAEFS-derived probability that F−1c (pc) will not be ex-ceeded is given by Ff (F−1c (pc)). The complement of the NCDF, given by1−Ff (F−1c (pc)), computes the NAEFS-derived probability that F−1c (pc) willbe exceeded. The PEI is related to the area above the NCDF bounded bypc = 0 and pc = Fc(F−1f (1)) (pc = 0.12 in the figure). PEI=0.02 indicatinglow extreme values have a low probability of exceedance. 1424.6. The Parametric Extreme IndexFigure 4.16: Day 1 PEI over BC on 3 Jan 2017.1434.6. The Parametric Extreme IndexFigure 4.17: Day 1 PEI over BC on 8 Nov 2016.1444.6. The Parametric Extreme IndexFigure 4.18: Day 1 PEI over BC on 28 Aug 2017.1454.6. The Parametric Extreme Indexclimate distributions (Lalaurette, 2003; Dutra et al., 2013) rather than ex-treme value distributions as defined above. Little is found in the literatureon the skill of either of these existing indices. Generally, SAs are associatedwith relatively high false alarm rates for medium- and long-range forecasts(Alcott, 2014). In contrast, another study suggests short-range forecastsSAs perform well (Qian et al., 2016). The EFI shows good skill for forecastsbetween the 10th and 90th percentiles (?).It is worth mentioning that the climatological mean and standard devia-tion of the SA are neither robust nor resistant. If the variable distributionsare asymmetric or non-Gaussian distributed (like precipitation or integratedprecipitable water), the mean and the standard deviation will misrepresentthe centre and spread of the data respectively.Caution is advised when verifying high-impact events. Studies haveshown that the skill of the verification metric of choice is proportional tothe frequency of extreme weather events (Schaefer, 1990). This creates themisleading impression that rare events cannot be skillfully forecast (Ferro,2007). For instance, the hit rate converges to zero with increasing rarity ofthe event (Stephenson et al., 2008).As most of the scores vanish with increasing rarity of the event (Ghelliand Primo, 2009; Stephenson et al., 2008), the Extreme Dependency Score(EDS) is proposed as a non-vanishing alternative to verify rare events andcompare the PEI with the SA during the summer of 2017. The bias-correcteduncalibrated NAEFS is used so that both indices are compared on equalfooting since the SA only uses the ensemble mean and cannot fully benefitfrom a calibrated forecast. The EDS is given by:2ln(a+cn)ln(an) − 1, (4.8)where n = a + b + c + d (see Table 4.2 for contingency table). The EDS ispresented with the base rate (BR). BR represents the probability that anevent occurs:BR =a+ cn(4.9)In Figure 4.19 the EDS is plotted as a function of 1 - BR. Thus, largevalues of 1 - BR represent small BR which are interpreted as rare events. Theprobability 1 - BR varies from 0.4 to 0.95 to capture very rare events. Eachpoint in the graph represents BR and EDS values calculated in a contingencytable generated with observations varying from above the 95th percentile tobelow the 99th percentile. That is, as the rarity of the event increases and1464.6. The Parametric Extreme IndexFigure 4.19: EDS as function of 1 - BR for Day 1 PEI and SA across allstations during summer of 2017. Rarer events are to the right. Values ofEDS closer to 1 are better.the observations fall above the 95th percentile to below the 99th percentile,a hit is issued if the SA falls above +1 (or below -1 for minimum T2M);for the PEI, a hit is issued if the PEI falls above 0. SA values above +1(below -1) are considered anomalous because they are at least one standarddeviation away from the mean.The best possible forecast has EDS=1. An EDS=1 happens if the fore-cast is either perfect or if it tends to overforecast the extreme event byavoiding any miss (c=0). In essence, as the number of hits (a) increase, andthe number of misses (c) decrease, the EDS approaches 1, the best possibleforecast. Conversely, the worst possible forecast has negative EDS values.That is, as the number of hits decrease, and the number of misses increase,the EDS computes negative values.4.6.2 PEI Performance across T2M and PCPIn Figure 4.19, the PEI outperforms the SA for all BR values for Summer2017. The shape of the curve determines the degree of dependency of theEDS on the BR. It implies that the number of misses (c) decreases fasterthan the number of hits (a) since the EDS approaches 1 as BR approachesless frequent events (Ghelli and Primo, 2009).Similar results are obtained during fall and winter (not shown). Based1474.7. Conclusionon these findings, the PEI clearly outperforms the SA for alerting users toextreme weather events, because the higher EDS values indicate a highernumber of hits and a lower number of misses across a range of extremeevents (ranging in rarity, or, in other words, return intervals).4.7 ConclusionA new very-high-resolution surface analyses (VHRSA) of daily maximumand minimum 2-m temperature (T2M), and 1-day accumulated precipita-tion (PCP) dataset that spans the period 1958-2017 was created by down-scaling and bias correcting the best performing reanalysis, the JRA-55. TheVHRSA was subsequently evaluated against the JRA-55 across the complexterrain of British Columbia (BC). In the evaluation systematic error, MAEand KS statistics were used to compare daily maximum and minimum T2M.To evaluate extreme maximum and minimum T2M, systematic error of 2-year return levels were compared. To compare daily PCP, the systematicerror of 31-day precipitation total, and two-sample χ2−statistic were calcu-lated. For extreme PCP, the systematic error of 2-year return levels of 1-dayaccumulated precipitation were compared.The VHRSA consistently exhibited better scores across all metrics through-out the year and across BC for both daily and extreme T2M. The cold biasof daily maximum T2M in the JRA-55 was mostly removed in the VHRSA.Although the difference in performance in daily and extreme minimum T2Mis not as apparent as it is in daily and extreme maximum T2M, the VHRSAstill outperforms the JRA-55 throughout the year and across all metrics eval-uated. The cold (warm) bias during fall and winter (spring and summer) inthe JRA-55 has been reduced in VHRSA.The results for PCP, while not quite as good as those for T2M, stillgenerally show substantial improvements in the VHRSA. For monthly andextreme PCP, all metrics indicate the VHRSA consistently outperformsthe JRA-55 throughout the seasons. In general, the downscaling and bias-correcting leads to a dataset with higher temperatures in valleys and lowertemperatures in mountainous regions. Similarly, for PCP the VHRSA isdrier in valleys and wetter in ridges and upper elevation regions. This isan important improvement since different studies have shown that NARRunder-predicts precipitation and temperature in mountainous regions acrossBC (Trubilowicz et al., 2016; Hunter et al., 2019).Thus, the VHRSA will serve as a valuable new spatially and temporallycomplete, high-resolution, 60-year climatological dataset for BC. Although1484.7. Conclusionin the present study the VHRSA was developed for BC, the same method-ology could be used to improve upon reanalysis datasets wherever obser-vations and PRISM or PRISM-like datasets are available — namely, theUnited States. The VHRSA has a wide range of potential applications inmeteorology, climatology, and hydrology. This study developed the VHRSAwith the intention of using it for some such applications.The VHRSA was first used to bias-correct and downscale the NAEFSensemble forecast. In order to conduct a thorough evaluation of the forecast,the systematic error, hit rate, false alarm rate, threat score, quantile BrierScore, bias and dispersion scores, and CRPS were used to compare the post-processed NAEFS to the raw NAEFS forecast.The post-processed NAEFS is generally statistically significantly moreskilful than the raw NAEFS forecast for forecast lead times out to 10 daysfor both T2M and PCP according to the metrics evaluated. Some previousstudies have concluded that higher-resolution models were better or at leastequal in performance to lower-resolution models across mountainous regions(Schirmer and Jamieson, 2015; Weusthoff et al., 2010; Garvert et al., 2005;Ikeda et al., 2010). A skilful very-high-resolution forecast is useful across thecomplex terrain of BC as many studies have shown that mountain rangesplay an important role in the regional and synoptic scale weather features(Deng et al., 2005; Astsatryan et al., 2015), and influence the distributionand intensity of precipitation and temperature (Junker et al., 1992; Kunzand Kottmeier, 2006; Smith et al., 2010; Haren et al., 2015).The Nonhomogeneous Gaussian regression (NGR) probabilistic calibra-tion method resulted in better results for temperature than precipitation.Previous studies have also shown that NGR results in better calibrationfor temperature (Hagedorn et al., 2008) than precipitation (Hamill et al.,2008). Baran and Lerch (2016) and Hamill et al. (2008) suggest differentcalibration techniques for non-Gaussian distributed fields such as precipi-tation and wind speeds. However, finding new robust methods to improvethe forecast calibration using such distributions in an operational settingwithout exacerbating computing demands is an ongoing issue. Nipen andStull (2011), Baran and Lerch (2016) and Hamill et al. (2008) present dif-ferent calibration techniques in datasets of the order of 103 as it was donehere calibrating point forecasts for station locations. Nipen and Stull (2011)evaluated different calibration methods on 1225 grid points across NorthAmerica between 2001 and 2004 for 2-m temperature, mean sea level pres-sure, wind speed, precipitable water and relative humidity. Baran and Lerch(2016) verified wind speed forecast across Germany between 2010 and 2011on 83220 grid points, and Hamill et al. (2008) verified precipitation skill at1494.7. Conclusiona 32-km horizontal resolution grid across North America from 1982 to 2001.However, conducting such a calibration for the number of grid points in thedownscaled NAEFS, which are on the order of 106, is problematic. In Siutaet al. (2017), which also used a dataset three orders of magnitude smaller,a regression method was reported to be computationally cost effective, andgave improved post-calibrated results for wind, which is a non-Gaussian dis-tributed variable. They argue this likely worked because the wind-forecasterror distribution was quasi-Gaussian.The VHRSA is then used to estimate the magnitude and the statisticalsignificance of trends extreme T2M and PCP. There is a noticeable andstatistically significant warming trend in the extreme maximum T2M acrossthe South Central climate zone throughout the year, and in extreme valuesof minimum T2M across most of BC during Summer and late January. Theresults for extreme T2M support the findings of Chapter 2, that there is anoticeable and statistically significant warming trend in the Southwest andSoutheast climate zones of BC during summer months for extreme maximumT2M; during winter months for extreme minimum T2M across BC. Vincentet al. (2012) performed a trend study during a similar study period (1950-2010) for annual maximum and minimum T2M and also reported the resultsfor minimum T2M were more significant than the results of maximum T2M.Finally, the new Parametric Extreme Index (PEI) outperformed an exist-ing situational/extremes awareness index, the Standardized Anomaly (SA).Its substantially superior EDS value across events with a range of rareness(extremity), indicated that the new index had more hits and less misseswhen detecting such events. The SA tool is already widely used and likedby operational forecasters. Since the PEI is substantially more useful, it hasgreat potential to see even more widespread use than the SA, for alertingforecasters to the potential for extreme weather.Although the PEI was developed using the VHRSA over BC, a similarmethodology could be employed to create such an index in other parts ofthe world. The formulation of the PEI alone would likely lead to improvedextremes alerting. Having more accurate and reliable probabilistic forecasts,however, will further improve PEI performance.Future work should examine the performance of alternative calibra-tion methods to further improve the PEI. One possibility is to upscale theVHRSA and probabilistic forecast somewhat, to perhaps 3 km, to reducecomputational demands. ? suggests decreasing grid spacing below 4 kmprovides more detail and structure but has only a limited impact on ac-curacy. This may make calibration of the entire gridded forecast over BCoperationally feasible without losing accuracy. Trends in mean values of1504.7. ConclusionT2M and PCP using the VHRSA should be evaluated as well.151Chapter 5ConclusionThe behaviour of extreme and daily values are inherently different. Addi-tionally, the distribution of daily and extreme maximum and minimum 2-mtemperature (T2M), and 1-day accumulated precipitation (PCP) are alsodifferent. This difference in behaviour and distribution motivated the de-velopment of new methodologies to statistically analyze and forecast them.The end goal is to better alert end users to forecasted extremes.First, in order to identify forecasted values as extreme (i.e., being on thetails of the climatological distribution), the best possible source for climato-logical data needed to be identified. The paucity of surface weather stationdata outside of southwest BC population centres motivated the search for thebest gridded reanalysis dataset to replace observations as a surrogate clima-tological dataset. Performance of the latest-generation reanalyses, the Cli-mate Forecast System Reanalysis (CFSR), the ECMWF interim reanalysis(ERA-Interim), the Japanese Meteorological Agency (JMA), and the Mod-ern Era Retrospective-Analysis for Research and Applications (MERRA-2),were rigorously evaluated with respect to daily and extreme maximum andminimum T2M, and daily and extreme PCP over the complex terrain of BC.New methodologies were developed for statistical evaluation of the re-analyses. These included evaluations based on a daily rolling window, break-ing results down by climate zone (which were also determined via a statis-tically rigorous method), and applying metrics that have rarely if ever beenused in previous meteorological studies. In addition to bringing improvedmethodologies into the field, an effort was made to account for the inherentassumptions in these methods, assumptions which are sometimes ignored.For example, statistical significance had to be determined accounting formultiple testing on a rolling window. Additional efforts were made to findmethods and metrics that would be appropriate specifically for extreme val-ues, such as the Generalized Extreme Value distribution (GEV).The result of the reanalysis evaluation showed the JRA-55 best repre-sents the climatological means, distributions, and distributions of extremeT2M and PCP over BC. While the ERA-Interim showed similar consistentand accurate results for daily and extreme T2M, and the MERRA-2 showed152Chapter 5. Conclusionbetter results for daily and extreme PCP, the JRA-55 was clearly the supe-rior overall reanalysis with consistent and accurate results across all metricsevaluated.In the reanalysis evaluation, it became clear that the coarse-resolutionJRA-55 had biases related to its inability to resolve the complex terrainof BC. Thus, a methodology was developed to downscale and bias correctT2M and PCP from the JRA-55. The temporal resolution of the JRA-55(6-hourly), the very high spatial resolution of the PRISM dataset (∼ 800m), and the homogeneity and ground truth from surface weather stations,were combined to create a spatially and temporally complete very-high-resolution surface analysis (VHRSA). Again, this new VHRSA was rigor-ously evaluated, and was generally found to be significantly less biased, andmore accurate than the JRA-55. Thus, a new, significantly improved, very-high-resolution 60-year daily surface analysis dataset was created for BC.This should prove to be an immensely valuable dataset for research andoperational use by the meteorological, climatological, and hydrological com-munities. For the purposes of this study, it also rendered a feasible solutionto the paucity of observational data across BC.Probabilistic forecasts of extremes requires a probabilistic forecast dataset.The North American Forecast System (NAEFS) is one such widely used en-semble forecast dataset. However, it is relatively low resolution, like thereanalysis datasets. Downscaling and bias correction were needed here aswell. The VHRSA was used to statistically downscale and bias correct theNAEFS. Computational restraints prohibited use of the VHRSA to proba-bilistically calibrate the forecasts, however, calibration was done for pointforecasts at weather station locations. Calibration is needed so that theforecasted probabilities are more accurate (that is, an event with a fore-casted 30% probability of occurrence is observed on average 30 out of 100times). While the bias correction was largely effective at significantly reduc-ing bias, the calibration delivered mixed results. It reported better resultscalibrating surface temperature — which is approximately Gaussian — thanprecipitation, a non-Gaussian distributed variable. Nonetheless, overall thepost-processed NAEFS showed significant improvement over the raw, inmost cases out to a forecast lead time of 10 days.The VHRSA was also used to analyze trends in extreme T2M and PCPduring the 1958-2017 period. A significant warming trend in extreme val-ues of minimum T2M across most of BC during Summer and January werefound — that is, extreme daily minimum temperatures are getting less coldat some times of year. Additionally, there was a significant increasing trendin extreme maximum T2M values across the South Central climate zone153Chapter 5. Conclusionthroughout the year — extremes of daily maximum temperatures are get-ting hotter. While there were some indications of warming trends in othermonths, they were not statistically significant. For extreme PCP, no in-creasing or decreasing trends were apparent.These post-processed probabilistic forecast distributions were then usedto create a new extreme forecast index (sometimes referred to as a situ-ational awareness tool). This so-called Parametric Extreme Index (PEI)differentiates itself and improves upon existing such indexes in that:1. It accounts for nonstationarity when appropriate. This accounts fortime-changing climate distributions, which were herein shown to besignificant for some variables, regions and times of year. It yields amore appropriate gauge of how extreme a value is based on the currentclimate.2. It is based on forecast distributions that are bias corrected and down-scaled onto a very-high-resolution grid. Further, point forecast valuesof the index are probabilistically calibrated. This improves the forecastaccuracy and bias of the index.3. It uses a more appropriate extreme values climate distribution, ratherthan a mean climate distribution. A mean climate distribution is inap-propriate for non-Gaussian distributed variables such as precipitation.The mean and the standard deviation will misrepresent the centreand spread of the data respectively leading to meaningless Standard-ized Anomalies values. Additionally, it leads to more false alarms asforecasts far away from the climate mean occur frequently.Indeed, the PEI was shown to substantially outperform the well knownStandardized Anomalies index across a range of extreme events ranging inrarity, with a higher number of hits and a lower number of misses. Both aredesirable to alert forecast users of future extreme weather events.Collectively, this dissertation developed appropriate statistical methodsto analyze and forecast extreme weather events. These new techniques leadto the creation of a new dataset, and the ability to better forecast extremeweather events. This methodology can be used to improve reanalyses wher-ever observations and PRISM or PRISM-like datasets are available, and beused for modelling in fields like hydrology, ecology, and agriculture besidesmeteorology and climatology. A tool to better forecast extreme weatherevents provides earlier and more accurate detection of such events, which inturn can help community responders, emergency managers, regional plan-154Chapter 5. Conclusionners, government and the media to take appropriate action to mitigate dam-age and reduce casualties.Future work will examine the performance of alternative calibratingmethods to post-process the very-high-resolution NAEFS forecast in a trueoperational setting without exacerbating computing demands. The PEI willbe evaluated across different seasons and forecast lead times. Finally, non-linear trends will be tested for precipitation to evaluate whether a morecomplex nonstationary model is more accurate to represent extreme precip-itation than the ones studied here. Trends in mean values of temperatureand precipitation using the VHRSA will also be evaluated.155BibliographyJ. T. Abatzoglou. Development of gridded surface meteorological data forecological applications and modelling. International Journal of Climatol-ogy, 33(1):121–131, 2013. ISSN 08998418. doi: 10.1002/joc.3413.R. F. Adler, G. J. Huffman, A. Chang, R. Ferraro, P.-P. Xie, J. Janowiak,B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind,P. Arkin, #. And, and E. Nelkin. The Version-2 Global Precipita-tion Climatology Project (GPCP) Monthly Precipitation Analysis (1979-Present). J. Hydrometeor, 4:1147–1167, 2003. URL http://precip.gsfc.nasa.gov.T. Alcott. How I Learned to Stop Worrying and Love the Global Ensembles,2014.J. Anderson. A method for producing and evaluating probabilistic forecastsfrom ensemble model integrations. Journal of Climate, 9:1518–1530, 1996.H. Astsatryan, A. Shakhnazaryan, V. Sahakyan, Y. Shoukourian,V. Kotroni, Z. Petrosyan, R. Abrahamyan, and H. Melkonyan. WRF-ARW model for prediction of high temperatures in south and south eastregions of Armenia. In Proceedings - 11th IEEE International Conferenceon eScience, eScience 2015, pages 207–213, 2015. ISBN 9781467393256.doi: 10.1109/eScience.2015.82.X. Bao and F. Zhang. Evaluation of NCEP-CFSR, NCEP-NCAR, ERA-Interim, and ERA-40 reanalysis datasets against independent soundingobservations over the Tibetan Plateau. Journal of Climate, 26:206–214,2013. ISSN 08948755. doi: 10.1175/JCLI-D-12-00056.1.S. Baran and S. Lerch. Mixture EMOS model for calibrating ensembleforecasts of wind speed. Environmetrics, 27(2):116–130, 2016. ISSN1099095X. doi: 10.1002/env.2380.BC Hydro. Record inflows, but limited downstream flood-ing, thanks to BC Hydro operations and some luck, 2016.156BibliographyURL https://www.bchydro.com/news/conservation/2016/comox-campbell-river-watersheds-break-records.html.S. Bentzien and P. Friederichs. Decomposition and graphical portrayal ofthe quantile score. Quarterly Journal of the Royal Meteorological Society,140(683):1924–1934, 2014. ISSN 1477870X. doi: 10.1002/qj.2284.A. A. Berg, J. Famiglieti, J. Walker, and P. Houser. Impact of bias correctionto reanalysis products on simulations of North American soil moisture andhydrological fluxes. Journal of Geophysical Research, 108(4490), 2003.ISSN 0148-0227. doi: 10.1029/2002JD003334.P. Berrisford, D. Dee, K. Fielding, M. Fuentes, P. Kallberg, S. Kobayashi,and S. Uppala. The ERA-Interim Archive Version 2.0. ERA report series1, 1(1):1–27, 2011.A. K. Betts, M. Ko¨hler, and Y. Zhang. Comparison of river basin hydromete-orology in ERA-Interim and ERA-40 reanalyses with observations. Jour-nal of Geophysical Research Atmospheres, 114(2), 2009. ISSN 01480227.doi: 10.1029/2008JD010761.J. Bhend and P. Whetton. Consistency of simulated and observed re-gional changes in temperature, sea level pressure and precipitation. Cli-matic Change, 118(3-4):799–810, jun 2013. ISSN 0165-0009. doi: 10.1007/s10584-012-0691-2. URL http://link.springer.com/10.1007/s10584-012-0691-2.S. C. Bloom, L. L. Takacs, A. M. da Silva, D. Ledvina, S. C. Bloom,L. L. Takacs, A. M. da Silva, and D. Ledvina. Data Assimilation UsingIncremental Analysis Updates. Monthly Weather Review, 124(6):1256–1271, jun 1996. ISSN 0027-0644. doi: 10.1175/1520-0493(1996)124〈1256:DAUIAU〉2.0.CO;2.M. G. Bosilovich, J. Chen, F. R. Robertson, and R. F. Adler. Evalu-ation of global precipitation in reanalyses. Journal of Applied Mete-orology and Climatology, 47(9):2279–2299, 2008. ISSN 15588424. doi:10.1175/2008JAMC1921.1.M. G. Bosilovich, S. Akella, L. Coy, R. Cullather, C. Draper, R. Gelaro,R. Kovach, Q. Liu, A. Molod, P. Norris, K. Wargan, W. Chao, R. Re-ichle, L. Takacs, Y. Vikhliaev, S. Bloom, A. Collow, S. Firth, G. Labow,157BibliographyG. Partyka, S. Pawson, O. Reale, S. D. Schubert, and M. Suarez. MERRA-2: Initial Evaluation of the Climate. Technical Report Series on GlobalModeling and Data Assimilation, 43(September):1–139, 2015.D. R. Bourdin, T. N. Nipen, and R. B. Stull. Reliable probabilistic forecastsfrom an ensemble reservoir inflow forecasting system. Water ResourcesResearch, 50:3108–3130, 2014. doi: 10.1002/2014WR015462.Received.Y. Brend. BC Ferries cancels sailing due to ’severe weather’ as winterchill hits — CBC News, 2016. URL https://www.cbc.ca/news/canada/british-columbia/weather-ferry-arctic-outbreak-1.3880910.R. Buizza. Potential Forecast Skill of Ensemble Prediction andSpread and Skill Distributions of the ECMWF Ensemble Pre-diction System. Monthly Weather Review, 125(1):99–119,1997. ISSN 0027-0644. doi: 10.1175/1520-0493(1997)125〈0099:PFSOEP〉2.0.CO;2. URL http://journals.ametsoc.org/doi/abs/10.1175/1520-0493{%}281997{%}29125{%}3C0099{%}3APFSOEP{%}3E2.0.CO{%}3B2.Y. Cai and D. Hames. Minimum sample size determination for generalizedextreme value distribution. Communications in Statistics: Simulation andComputation, 40(1):87–98, 2010. ISSN 15324141. doi: 10.1080/03610918.2010.530368.G. Candille, C. Coˆte´, P. L. Houtekamer, and G. Pellerin. Verification of anEnsemble Prediction System against Observations. Monthly Weather Re-view, 135(7):2688–2699, 2007. ISSN 0027-0644. doi: 10.1175/MWR3414.1.URL http://journals.ametsoc.org/doi/abs/10.1175/MWR3414.1.M. Charron, G. Pellerin, L. Spacek, P. L. Houtekamer, N. Gagnon, H. L.Mitchell, and L. Michelin. Toward Random Sampling of Model Errorin the Canadian Ensemble Prediction System. Monthly Weather Review,138(5):1877–1901, 2010. ISSN 0027-0644. doi: 10.1175/2009MWR3187.1.URL http://journals.ametsoc.org/doi/abs/10.1175/2009MWR3187.1.R. R. Chilton. A Summary of Climatic Regimes of British Columbia. B.C.Min. Environ., Victoria, B.C., 1981.C. Correia. Winter conditions testing patience of Lower Mainland resi-dents — CBC News, 2016. URL https://www.cbc.ca/news/canada/british-columbia/lower-mainland-snow-issues-1.3913674.158BibliographyB. A. Cosgrove. Real-time and retrospective forcing in the North Ameri-can Land Data Assimilation System (NLDAS) project. Journal of Geo-physical Research, 108(D22):8842, 2003. ISSN 0148-0227. doi: 10.1029/2002JD003118.F. C. Curriero, J. a. Patz, J. B. Rose, and S. Lele. The Association BetweenExtreme Precipitation and Waterborne Disease Outbreaks in the UnitedStates , 1948 – 1994. American Journal of Public Health, 91(8):1194–1199,2001. ISSN 0090-0036. doi: 10.2105/AJPH.91.8.1194.C. Daly, R. Neilson, and D. Phillips. A Statistical-Topographic Model forMapping Climatological Precipitation over Mountainous Terrain. Journalof Applied Meteorology, 33:140–158, 1994.C. Daly, G. Taylor, and W. Gibson. The Prism approach to mapping precipi-tation and temperature. In 10th AMS Conference on Applied Climatology,number 1, pages 10–12, Reno, NV, 1997.C. Daly, G. H. Taylor, W. P. Gibson, T. W. Parzybok, G. L. Johnson, andP. A. Pasteris. Climate division normals derived from topographically-sensitive climate grids. In 13th AMS Conf. on Applied Climatology, pages177–180, Portland, OR, 2002.R. A. Davis. The Rate of Convergence in Distribution of the Maxima.Statistica Neerlandica, 36(1):31–35, 1982.L. De Haan. Sample Extremes: An Elementary Introduction. StatisticaNeerlandica, 30(4):161–172, 1976. doi: 10.1111/j.1467-9574.1976.tb00275.x.M. Decker, M. A. Brunke, Z. Wang, K. Sakaguchi, X. Zeng, and M. G.Bosilovich. Evaluation of the Reanalysis Products from GSFC, NCEP,and ECMWF Using Flux Tower Observations. Journal of Climate, 25:1916–1944, 2012.D. P. Dee, S. M. Uppala, A. J. Simmons, P. Berrisford, P. Poli, S. Kobayashi,U. Andrae, M. A. Balmaseda, G. Balsamo, P. Bauer, P. Bechtold, A. C. M.Beljaars, L. van de Berg, J. Bidlot, N. Bormann, C. Delsol, R. Dragani,M. Fuentes, A. J. Geer, L. Haimberger, S. B. Healy, H. Hersbach, E. V.Ho´lm, L. Isaksen, P. K˚allberg, M. Ko¨hler, M. Matricardi, A. P. McNally,B. M. Monge-Sanz, J.-J. Morcrette, B.-K. Park, C. Peubey, P. de Rosnay,C. Tavolato, J.-N. The´paut, and F. Vitart. The ERA-Interim reanalysis:configuration and performance of the data assimilation system. Quarterly159BibliographyJournal of the Royal Meteorological Society, 137(656):553–597, apr 2011.ISSN 00359009. doi: 10.1002/qj.828.Demography Division. Canadian Demographics at a Glance. StatisticsCanada, 91(003):1– 81, 2016. ISSN 1916-1832. doi: 91-003-XWE.X. Deng, R. Stull, and B. Columbia. A Mesoscale Analy-sis Method for Surface Potential Temperature in Mountainousand Coastal Terrain. Monthly weather review, 133(2):389–408,2005. ISSN 0027-0644. doi: 10.1175/MWR-2859.1. URL http://journals.ametsoc.org/doi/pdf/10.1175/MWR-2859.1{%}0Ahttp://journals.ametsoc.org/doi/abs/10.1175/MWR-2859.1.M. D. Dettinger, D. R. Cayan, H. F. Diaz, and D. M. Meko. North-South precipitation patterns in western North America on interannual-to-decadal timescales. Journal of Climate, 11(12):3095–3111, 1998. ISSN08948755. doi: 10.1175/1520-0442(1998)011〈3095:NSPPIW〉2.0.CO;2.K. a. Devine and E´. Mekis. Field accuracy of Canadian rain measurements.Atmosphere - Ocean, 46(2):213–227, 2008. ISSN 07055900. doi: 10.3137/ao.460202.M. G. Donat, L. V. Alexander, H. Yang, I. Durre, R. Vose, R. J. H.Dunn, K. M. Willett, E. Aguilar, M. Brunet, J. Caesar, B. Hewitson,C. Jack, A. M. G. Klein Tank, A. C. Kruger, J. Marengo, T. C. Peterson,M. Renom, C. Oria Rojas, M. Rusticucci, J. Salinger, A. S. Elrayah, S. S.Sekele, A. K. Srivastava, B. Trewin, C. Villarroel, L. A. Vincent, P. Zhai,X. Zhang, and S. Kitching. Updated analyses of temperature and precip-itation extreme indices since the beginning of the twentieth century: TheHadEX2 dataset. Journal of Geophysical Research: Atmospheres, 118(5):2098–2118, mar 2013. ISSN 2169897X. doi: 10.1002/jgrd.50150.E. Dutra, M. Diamantakis, I. Tsonevsky, E. Zsoter, F. Wetterhall, T. Stock-dale, D. Richardson, and F. Pappenberger. The extreme forecast indexat the seasonal scale. Atmospheric Science Letters, 14(4):256–262, 2013.ISSN 1530261X. doi: 10.1002/asl2.448.A. Ebita, S. Kobayashi, Y. Ota, and M. Moriya. The new Japanese reanal-ysis project: JRA55. In Proc. 15th GCOS/WCRP Atmospheric Observa-tion Panel for Climate, volume Global Cli, page 12.2, Geneva, Switzer-land, 2009. Climate Prediction Division, JMA.160BibliographyA. Ebita, S. Kobayashi, Y. Ota, M. Moriya, R. Kumabe, O. Kazu-toshi, Y. Harada, S. Yasui, K. Miyaoka, K. Takahashi, H. Kamahori,C. Kobayashi, H. Endo, M. Soma, Y. Oikawa, and T. Ishimizu. TheJapanese 55-year Reanalysis “JRA-55”: An Interim Report. SOLA: Sci-entific Online Letters of the Atmosphere, 7:149–152, 2011.C. a. T. Ferro. A Probability Model for Verifying Deterministic Forecasts ofExtreme Events. Weather and Forecasting, 22(5):1089–1100, 2007. ISSN0882-8156. doi: 10.1175/WAF1036.1.W. Y. Fritz. Winter is no reason to mind your mannersin @WinnieYeo . Twitter, (Dec 20, 2016):Retrieved fromhttps://twitter.com/winnieyeo?lang=en., 2016.M. F. Garvert, B. a. Colle, and C. F. Mass. The 13–14 December2001 IMPROVE-2 Event. Part I: Synoptic and Mesoscale Evolutionand Comparison with a Mesoscale Model Simulation. Journal of theAtmospheric Sciences, 62(10):3474–3492, 2005. ISSN 0022-4928. doi:10.1175/JAS3549.1. URL http://journals.ametsoc.org/doi/abs/10.1175/JAS3549.1.R. Gelaro. MERRA-2 and future reanalysis at GMAO. In NOAA ClimateReanalysis Task Force Technical Workshop, page 15pp., College Park,MD, 2015.R. Gelaro, W. McCarty, M. J. Sua´rez, R. Todling, A. Molod, L. Takacs,C. a. Randles, A. Darmenov, M. G. Bosilovich, R. Reichle, K. Wargan,L. Coy, R. Cullather, C. Draper, S. Akella, V. Buchard, A. Conaty, A. M.da Silva, W. Gu, G. K. Kim, R. Koster, R. Lucchesi, D. Merkova, J. E.Nielsen, G. Partyka, S. Pawson, W. Putman, M. Rienecker, S. D. Schu-bert, M. Sienkiewicz, and B. Zhao. The modern-era retrospective analysisfor research and applications, version 2 (MERRA-2). Journal of Climate,30(14):5419–5454, 2017. ISSN 08948755. doi: 10.1175/JCLI-D-16-0758.1.A. Ghelli and C. Primo. On the use of the extreme dependency score toinvestigate the performance of an NWP model for rare events. Meteoro-logical Applications, 16:537–544, 2009.B. Gnedenko. On the Limiting Distribution of a Supercritical BranchingProcess in a Random Environment. Annals of Mathematics, 44:423–453,1943.161BibliographyT. Gneiting, A. E. Raftery, A. H. Westveld, and T. Goldman. Cali-brated Probabilistic Forecasting Using Ensemble Model Output Statis-tics and Minimum CRPS Estimation. Monthly Weather Review, 133(5):1098–1118, 2005. ISSN 0027-0644. doi: 10.1175/MWR2904.1. URLhttp://journals.ametsoc.org/doi/abs/10.1175/MWR2904.1.R. Graham, T. Alcott, N. Hosenfeld, and R. Grumm. Anticipating a rareevent utilizing forecast anomalies and a situational awareness display: TheWestern U.S. Storms of 18-23 January 2010. Bulletin of the AmericanMeteorological Society, 94(12):1827–1836, 2013. ISSN 00030007. doi: 10.1175/BAMS-D-11-00181.1.P. Groisman, R. Knight, D. Easterling, T. R. Karl, G. C. Hegerl, and V. N.Razuvaev. Trends in precipitation intensity in the climate record. Journalof Climate, 18:1326–1350, 2005. ISSN 08948755. doi: 10.1175/JCLI3339.1.V. Grubiˇsic´, R. K. Vellore, and A. W. Huggins. Quantitative PrecipitationForecasting of Wintertime Storms in the Sierra Nevada: Sensitivity tothe Microphysical Parameterization and Horizontal Resolution. MonthlyWeather Review, 133(10):2834–2859, 2005. ISSN 0027-0644. doi: 10.1175/MWR3004.1. URL http://journals.ametsoc.org/doi/abs/10.1175/MWR3004.1.R. Grumm and R. Hart. Standardized anomalies applied to significant coldseason weather events: Preliminary findings. Weather & Forecasting, 16:736–754, 2001. ISSN 0882-8156. doi: 10.1175/1520-0434(2001)016〈0736:SAATSC〉2.0.CO;2.A. M. Gunn. Climate and Weather Extremes. Greeenwood, 2010. ISBN9780313355684. doi: 10.1016/B978-0-12-386917-3.00010-5.J. P. Hacker, E. Scott, K. And, and R. B. Stull. Ensemble Experiments onNumerical Weather Prediction Error and Uncertainty for a North PacificForecast Failure. Weather and Forecasting, 18:12–31, 2003.R. Hagedorn, T. M. Hamill, and J. S. Whitaker. Probabilistic ForecastCalibration Using ECMWF and GFS Ensemble Reforecasts. Part I: Two-Meter Temperatures. Monthly Weather Review, 136(7):2608–2619, 2008.ISSN 0027-0644. doi: 10.1175/2007MWR2410.1. URL http://journals.ametsoc.org/doi/abs/10.1175/2007MWR2410.1.T. M. Hamill. NOTES AND CORRESPONDENCE Interpretation of RankHistograms for Verifying Ensemble Forecasts. Mon. Wea. Rev., 129:162Bibliography550–560, 2001. URL https://journals.ametsoc.org/doi/pdf/10.1175/1520-0493{%}282001{%}29129{%}3C0550{%}3AIORHFV{%}3E2.0.CO{%}3B2.T. M. Hamill, R. Hagedorn, and J. S. Whitaker. Probabilistic ForecastCalibration Using ECMWF and GFS Ensemble Reforecasts. Part II: Pre-cipitation. Monthly Weather Review, 136(7):2620–2632, 2008. ISSN 0027-0644. doi: 10.1175/2007MWR2411.1. URL http://journals.ametsoc.org/doi/abs/10.1175/2007MWR2411.1.M. Harada. Introduction to reanalysis and JRA-55. In TCC Training Sem-inar on Seasonal Forecast, Tokyo Climate Center - Japan MeteorologicalAgency, number Tokyo Climate Center - Japan Meteorological Agency,pages 17–26, Tokyo, Japan, 2018.Y. Harada, H. Kamahori, C. Kobayashi, H. Endo, S. Kobayashi, Y. Ota,H. Onoda, K. Onogi, K. Miyaoka, and K. Takahashi. The JRA-55 Reanal-ysis: Representation of Atmospheric Circulation and Climate Variability.Journal of the Meteorological Society of Japan. Ser. II, 94(3):269–302,2016. ISSN 0026-1165. doi: 10.2151/jmsj.2016-015.R. V. Haren, R. J. Haarsma, G. J. V. Oldenborgh, and W. Hazeleger. Res-olution dependence of European precipitation in a state-of-the-art atmo-spheric general circulation model. Journal of Climate, 28(13):5134–5149,2015. ISSN 08948755. doi: 10.1175/JCLI-D-14-00279.1.R. E. Hart and R. H. Grumm. Using Normalized Climatological Anomaliesto Rank Synoptic-Scale Events Objectively. Monthly Weather Review,129:2426–2442, 2001. ISSN 0027-0644. doi: 10.1175/1520-0493(2001)129〈2426:UNCATR〉2.0.CO;2.S. R. Haughian, P. J. Burton, S. W. Taylor, and C. L. Curry. Expected ef-fects of climate change on forest disturbance regimes in British Columbia.BC Journal of Ecosystems and Management, 13(1):16–39, 2012. ISSN15945685.S. C. Herring, N. Christidis, A. Hoell, J. P. Kossin, C. J. S. Iii, and P. a.Stott. Explaining Extreme Events of 2016 From a Climate Perspective.Bulletin of the American Meteorological Society, 98(12):S1–S157, 2017.ISSN 0003-0007. doi: 10.1175/BAMS-D-17-0118.1. URL http://www.ametsoc.net/eee/2016/2016{_}bams{_}eee{_}high{_}res.pdf.163BibliographyK. I. Hodges, R. W. Lee, and L. Bengtsson. A comparison of extratropical cy-clones in recent reanalyses ERA-Interim, NASA MERRA, NCEP CFSR,and JRA-25. Journal of Climate, 24:4888–4906, 2011. ISSN 08948755.doi: 10.1175/2011JCLI4097.1.M. Hollander, D. A. Wolfe, and E. Chicken. Nonparametric Statistical Meth-ods. John Wiley and Sons, Hoboken, New Jersey, 3rd edition, 2013. ISBN978-0-471-19045-5.J. Hosking. L-Moments : Analysis and Estimation of Distributions UsingLinear Combinations of Order Statistics. Journal of the Royal StatisticalSociety - Series B, 52(1):105–124, 1990.J. Hosking, J. Wallis, and E. Woo. Estimation of the generalized ex-treme value distribution by the method of probability weighted mo-ments. Technometrics, 27(3):251–261, 1985. ISSN 0040-1706. doi:dx.doi.org/10.1080/00401706.1985.10488049.D. Hou, Y. Zhu, X. Zhou, R. Wobus, J. Peng, Y. Luo, and B. Cui. The2015 Upgrade of NCEP’s Global Ensemble Forecast System (GEFS).In 27th Conf. on Weather Analysis and Forecasting/23rd Conf. on Nu-merical Weather Prediction, pages Amer. Meteor. Soc., 2A.6., Chicago,IL, jun 2015. AMS. URL https://ams.confex.com/ams/27WAF23NWP/webprogram/Paper273406.html.C. Hunter, R. Moore, and I. McKendry. Evaluation of the north americanregional reanalysis (narr) precipitation fields in a topographically complexdomain. Hydrological Sciences Journal, 0(ja):null, 2019.K. Ikeda, R. Rasmussen, C. Liu, D. Gochis, D. Yates, F. Chen, M. Tewari,M. Barlage, J. Dudhia, K. Miller, K. Arsenault, V. Grubiˇsic´, G. Thomp-son, and E. Guttman. Simulation of seasonal snowfall over Colorado.Atmospheric Research, 97(4):462–477, sep 2010. ISSN 0169-8095. doi: 10.1016/J.ATMOSRES.2010.04.010. URL https://www.sciencedirect.com/science/article/pii/S0169809510001043.P. L. Jackson. Surface winds during an intense outbreak of arctic air inSouthwestern British Columbia. Atmosphere-Ocean, 34(2):285–311, jun1996. ISSN 0705-5900. doi: 10.1080/07055900.1996.9649566.A. F. Jenkinson. The frequency distribution of the annual maximum (orminimum) values of meteorological elements. Quarterly Journal of the164BibliographyRoyal Meteorological Society, 81(348):158–171, apr 1955. ISSN 00359009.doi: 10.1002/qj.49708134804.P. Jones, S. Raper, B. Santer, B. Cherry, C. Goodess, R. Bradley, H. Diaz,P. Kelly, and T. Wigley. A grid point surface air temperature data set forthe Northern Hemisphere. Technical Report DOE/EV/10098-2, UnitedStates Department of Energy, Washington, D.C., 1985.P. D. Jones, M. New, D. Parker, S. Martin, and I. G. Rigor. Surface air tem-perature and its changes over the past 150 years. Reviews of Geophysics,37(2):173–199, 1999.R. W. Jones, I. Renfrew, A. Orr, B. G. M. Webber, D. M. Holland, andM. A. Lazzara. Evaluation of four global reanalysis products using in situobservations in the Amundsen Sea Embayment, Antarctica. J. Geophys.Res. Atmos., 121, 2016. ISSN 2169897X. doi: 10.1002/2016JD025476.H.-m. H. Juang and M. Kanamitsu. The NMC Nested Regional SpectralModel. Monthly Weather Review, 122(January):3–26, 1994.N. W. Junker, J. E. Hoke, B. E. Sullivan, K. F. Brill, and F. J. Hughes.Seasonal geographic variations in quantita- tive precipitation predictionby NMC’s nested-grid model and medium-range forecast model. Wea.Forecasting, 7:410–429, 1992.T. R. Karl, P. Y. Groisman, R. W. Knight, and R. R. Heim. Recent vari-ations of snow cover and snowfall in North America and their relationto precipitation and temperature variations. Journal of Climate, 6(7):1327–1344, 1993. ISSN 08948755. doi: 10.1175/1520-0442(1993)006〈1327:RVOSCA〉2.0.CO;2.E. J. Kendon, N. M. Roberts, H. J. Fowler, M. J. Roberts, S. C. Chan, andC. A. Senior. Heavier summer downpours with climate change revealed byweather forecast resolution model. Nature Climate Change, 4(7):570–576,jul 2014. ISSN 1758-678X. doi: 10.1038/nclimate2258.V. V. Kharin and F. W. Zwiers. Changes in the extremes in an ensembleof transient climate simulations with a coupled atmosphere-ocean GCM.Journal of Climate, 13(21):3760–3788, 2000. ISSN 08948755. doi: 10.1175/1520-0442(2000)013〈3760:CITEIA〉2.0.CO;2.S. Kobayashi, Y. Ota, Y. Harada, and A. Ebita. The JRA-55 Reanalysis:General Specifications and Basic Characteristics. Journal of the Meteoro-logical Society of Japan, 93(1):5–48, 2015.165BibliographyR. D. Koster, W. McCarty, L. Coy, R. Gelaro, A. Huang, D. Merkova, E. B.Smith, M. Sienkiewicz, and K. Wargan. MERRA-2 Input Observations:Summary and Assessment. Technical Report NASA/TM–2016-104606,National Aeronautics and Space Administration - NASA, Greenbelt, 2016.C. K. B. Krishnamurthy, U. Lall, H.-H. Kwon, C. K. B. Krishnamurthy,U. Lall, and H.-H. Kwon. Changing Frequency and Intensity of RainfallExtremes over India from 1951 to 2003. Journal of Climate, 22(18):4737–4746, sep 2009. ISSN 0894-8755. doi: 10.1175/2009JCLI2896.1.S. Kumar, V. Merwade, J. L. Kinter, and D. Niyogi. Evaluation of tem-perature and precipitation trends and long-term persistence in CMIP5twentieth-century climate simulations. Journal of Climate, 26(12):4168–4185, 2013. ISSN 08948755. doi: 10.1175/JCLI-D-12-00259.1.M. Kunz and C. Kottmeier. Orographic Enhancement of Precipitation overLow Mountain Ranges. Part II: Simulations of Heavy Precipitation Eventsover Southwest Germany. Journal Of Applied Meteorology And Climatol-ogy, 45:1041–1055, 2006. ISSN 1558-8424. doi: 10.1175/JAM2390.1.M. Laanela. Vancouver mayor absent as councillor calls for snow removalinquiry — CBC News, 2016. URL https://www.cbc.ca/news/canada/british-columbia/vancouver-snow-removal-1.3922282.R. Lader, U. S. Bhatt, J. E. Walsh, T. S. Rupp, and P. a. Bieniek. Two-metertemperature and precipitation from atmospheric reanalysis evaluated forAlaska. Journal of Applied Meteorology and Climatology, 55(4):901–922,2016. ISSN 15588432. doi: 10.1175/JAMC-D-15-0162.1.F. Lalaurette. Early detection of abnormal weather conditions using aprobabilistic extreme forecast index. Quarterly Journal of the RoyalMeteorological Society, 129(594):3037–3057, 2003. ISSN 00359009. doi:10.1256/qj.02.152.M. Leadbetter, G. Lindgren, and H. Rootzen. Extremes and Related Prop-erties of Random Sequences and Processes. Springer Verlag, 1983. ISBN9780387775005. doi: 10.1007/978-0-387-98135-2.R. Lindsay, M. Wensnahan, a. Schweiger, and J. Zhang. Evaluation of sevendifferent atmospheric reanalysis products in the arctic. Journal of Climate,27(7):2588–2606, 2014. ISSN 08948755. doi: 10.1175/JCLI-D-13-00014.166BibliographyG. Luymes. B.C. weather: Winter storm blows over in time for FamilyDay — Vancouver Sun, 2017. URL https://vancouversun.com/news/local-news/winter-storm-blows-over-in-time-for-family-day.M. MacMahon. Fraser Valley hit hardest by snowfall, more on the way -NEWS 1130, 2017. URL https://www.citynews1130.com/2017/02/06/fraser-valley-hit-hardest-snowfall-way/.H. B. Mann. Non-parametric tests against trend. Econometrica, 13(3):245–259, 1945. ISSN 0168-6054. doi: 10.1016/j.annrmp.2004.07.001.M. E. Mann and K. A. Emanuel. Atlantic hurricane trends linked to climatechange. Eos, Transactions American Geophysical Union, 87(24):233, jun2006. ISSN 0096-3941. doi: 10.1029/2006EO240001.D. McCollor and R. Stull. Hydrometeorological Accuracy Enhancementvia Postprocessing of Numerical Weather Forecasts in Complex Ter-rain. Weather and Forecasting, 23:131–144, 2008a. ISSN 0882-8156.doi: 10.1175/2008WAF2222177.1. URL http://journals.ametsoc.org/doi/abs/10.1175/2008WAF2222177.1.D. McCollor and R. Stull. Hydrometeorological Short-Range Ensemble Fore-casts in Complex Terrain. Part I: Meteorological Evaluation. Weatherand Forecasting, 23(4):533–556, 2008b. ISSN 0882-8156. doi: 10.1175/2008WAF2007063.1. URL http://journals.ametsoc.org/doi/abs/10.1175/2008WAF2007063.1.J. McElroy. Extended cold snap putting a strain on ser-vices throughout Metro Vancouver — CBC News, 2016.URL https://www.cbc.ca/news/canada/british-columbia/extended-cold-snap-putting-a-strain-on-services-throughout-metro-vancouver-1.3897474.J. McElroy. Vancouver in its longest cold snap in over 30 years — CBC News,2017a. URL https://www.cbc.ca/news/canada/british-columbia/when-vancouver-had-winter-1.3918910.J. McElroy. Weather closes every highway linking LowerMainland to rest of B.C. — CBC News, 2017b. URLhttps://www.cbc.ca/news/canada/british-columbia/weather-closes-every-highway-linking-lower-mainland-to-rest-of-b-c-1.3975497.167BibliographyE´. Mekis. J3.7 Adjustments for trace measurements in Canada. In 15thConference on Applied Climatology, Savannah, Georgia, USA, 20-24 June2005, 2005.E´. Mekis and R. Brown. Derivation of an adjustment factor map for theestimation of the water equivalent of snowfall from ruler measurements inCanada. Atmosphere - Ocean, 48(4):284–293, 2010. ISSN 07055900. doi:10.3137/AO1104.2010.E. Mekis and W. D. Hogg. Rehabilitation and analysis of Canadian dailyprecipitation time series. Atmosphere-Ocean, 37(1):53–85, 1999. ISSN0705-5900. doi: 10.1080/07055900.1999.9649621.E´. Mekis and L. a. Vincent. An overview of the second generation adjusteddaily precipitation dataset for trend analysis in Canada. Atmosphere -Ocean, 49(2):163–177, 2011. ISSN 07055900. doi: 10.1080/07055900.2011.583910.J. R. Metcalfe, B. Routledge, and K. Devine. Rainfall measurement inCanada: Changing observational methods and archive adjustment pro-cedures. Journal of Climate, 10(1):92–101, 1997. ISSN 08948755. doi:10.1175/1520-0442(1997)010〈0092:RMICCO〉2.0.CO;2.Meteorological Service of Canada. MANOBS - manual of surface weatherobservations. 2015. ISBN 9781100254456.M. Meuse. Fraser Valley weather the worst veteran CBC cameraman hasever seen — CBC News, 2017. URL https://www.cbc.ca/news/canada/british-columbia/freezing-rain-worst-ever-1.3974244.A. Moberg and P. D. Jones. Trends in indices for extremes in daily tem-perature and precipitation in central and western Europe, 1901-99. Inter-national Journal of Climatology, 25(9):1149–1171, 2005. ISSN 08998418.doi: 10.1002/joc.1163.P. a. Mooney, F. J. Mulligan, and R. Fealy. Comparison of ERA-40, ERA-Interim and NCEP/NCAR reanalysis data with observed surface air tem-peratures over Ireland. International Journal of Climatology, 31(4):545–557, 2011. ISSN 08998418. doi: 10.1002/joc.2098.R. D. Moore, D. L. Spittlehouse, P. H. Whitfield, and K. Stahl. Chapter3 - Weather and Climate. Number November. Ministry of Forests andRange, Research Branch, Victoria, B.C. and FORREX Forest ResearchExtension Partnership, Kamloops, BC, 2008.168BibliographyP. W. Mote, A. F. Hamlet, M. P. Clark, and D. P. Lettenmaier. Decliningmountain snowpack in western north America. Bulletin of the AmericanMeteorological Society, 86(1):39–49, 2005. ISSN 00030007. doi: 10.1175/BAMS-86-1-39.A. H. Murphy and R. L. Winkler. A General Framework for Forecast Ver-ification. Monthly Weather Review, 115(7):1330–1338, jul 1987. ISSN0027-0644. doi: 10.1175/1520-0493(1987)115〈1330:AGFFFV〉2.0.CO;2. URL http://dx.doi.org/10.1175/1520-0493(1987)115{%}3C1330:AGFFFV{%}3E2.0.CO2.T. Nipen and R. Stull. Calibrating probabilistic forecasts from an NWPensemble. Tellus, Series A: Dynamic Meteorology and Oceanography, 63(5):858–875, 2011. ISSN 02806495. doi: 10.1111/j.1600-0870.2011.00535.x.P. Odon, G. West, and R. Stull. Reanlyses Performance across BritishColumbia . Part II : Evaluation of Extreme Precipitation. Journal ofApplied Meteorology and Climatology, (in press).P. Odon, G. West, and R. Stull. Vancouver Fall and Winter 2016/17: HowBad Was It? CMOS Bulletin SCMO, 45(4):9–12, 2017.P. Odon, G. West, and R. Stull. Evaluation of reanalyses over BritishColumbia. Part I: Daily and Extreme 2-m temperature. Journal of AppliedMeteorology and Climatology, 57(9):2091–2112, 2018. ISSN 15588432. doi:10.1175/JAMC-D-18-0058.1.C. of Vancouver. City crews continue work to re-move remaining snow and ice — City of Vancou-ver, 2017. URL https://vancouver.ca/news-calendar/city-crews-continue-work-to-remove-remaining-snow-and-ice-jan-2017.aspx.B. Onoz and M. Bayazit. The power of statistical tests for trend detection.Turkish Journal of Engineering and Environmental Sciences, 27(4):247–251, 2003. ISSN 13000160. doi: 10.1016/j.ica.2004.10.030.J. Paetkau. Record breaking power usage as cold snap continues— CBC News, 2017. URL https://www.cbc.ca/news/canada/british-columbia/bc-hydro-cold-winter-power-ted-olynyk-1.3921852.169BibliographyC. Parmesan, T. L. Root, and M. R. Willig. Impacts of extreme weather andclimate on terrestiral biota. Bulletin of the American Meteorological Soci-ety, 81(3):443–450, 2000. ISSN 0003-0007. doi: 10.1175/1520-0477(2000)081〈0443.J. A. Patz, D. Campbell-Lendrum, T. Holloway, and J. A. Foley. Impactof regional climate change on human health. Nature, 438(7066):310–317,2005. ISSN 14764687. doi: 10.1038/nature04188.PCIC, C. Pacific Climate Impacts, V. University of, P. C. Group, and O. S.University. PCIC and Prism Climate Group, 2014: High Resolution Cli-matology, 2014. URL https://pacificclimate.org/.M. C. Peel, B. L. Finlayson, and T. a. McMahon. Updated world mapof the Ko¨ppen-Geiger climate classification. Hydrology and Earth Sys-tem Sciences, 11(5):1633–1644, 2007. ISSN 16077938. doi: 10.5194/hess-11-1633-2007.T. C. Peterson and R. S. Vose. An Overview of the Global HistoricalClimatology Network Temperature Database. Bulletin of the Ameri-can Meteorological Society, 78(12):2837–2849, 1997. ISSN 00030007. doi:10.1175/1520-0477(1997)078〈2837:AOOTGH〉2.0.CO;2.T. C. Peterson, D. R. Easterling, T. R. Karl, P. Groisman, N. Nicholls,N. Plummer, S. Torok, I. Auer, R. Boehm, D. Gullett, L. Vincent,R. Heino, H. Tuomenvirta, O. Mestre, T. Szentimrey, J. Salinger, E. J.Forland, I. Hanssen-Bauer, H. Alexandersson, P. Jones, and D. Parker.Homogeneity Adjustsments of in situ Atmospheric Climate Data: AReview. International Journal of Climatology, 18(13):1493–1517, 1998.ISSN 08998418. doi: 10.1002/(SICI)1097-0088(19981115)18:13〈1493::AID-JOC329〉3.0.CO;2-T.D. Phillips. Canada’s Top Ten Weather Stories 2017. CMOS Bulletin SCMO,46(1):13–17, 2018.L. S. Porth, D. C. Boes, R. a. Davis, C. a. Troendle, and R. M. King.Development of a technique to determine adequate sample size using sub-sampling and return interval estimation. Journal of Hydrology, 251(1-2):110–116, 2001. ISSN 00221694. doi: 10.1016/S0022-1694(01)00442-5.W. Qian, N. Jiang, and J. Du. Anomaly-Based Weather Analysis versusTraditional Total-Field-Based Weather Analysis for Depicting RegionalHeavy Rain Events. Weather and Forecasting, 31(1):71–93, 2016. ISSN170Bibliography0882-8156. doi: 10.1175/WAF-D-15-0074.1. URL http://journals.ametsoc.org/doi/abs/10.1175/WAF-D-15-0074.1.F. M. Ralph, P. J. Neiman, and G. a. Wick. Satellite and CALJET AircraftObservations of Atmospheric Rivers over the Eastern North Pacific Oceanduring the Winter of 1997/98. Monthly Weather Review, 132(7):1721–1745, 2004. ISSN 0027-0644. doi: 10.1175/1520-0493(2004)132〈1721:SACAOO〉2.0.CO;2.F. M. Ralph, P. J. Neiman, and R. Rotunno. Dropsonde Observations inLow-Level Jets over the Northeastern Pacific Ocean from CALJET-1998and PACJET-2001: Mean Vertical-Profile and Atmospheric-River Char-acteristics. Monthly Weather Review, 133(4):889–910, 2005. ISSN 0027-0644. doi: 10.1175/MWR2896.1.R. Rasmussen, C. Liu, K. Ikeda, D. Gochis, D. Yates, F. Chen, M. Tewari,M. Barlage, J. Dudhia, W. Yu, K. Miller, K. Arsenault, V. Grubiˇsic´,G. Thompson, and E. Gutmann. High-resolution coupled climate runoffsimulations of seasonal snowfall over Colorado: A process study of currentand warmer climate. Journal of Climate, 24(12):3015–3048, 2011. ISSN08948755. doi: 10.1175/2010JCLI3985.1.A. Ravishankar, G. Meehl, and L. F. Bosart. Attribution of ExtremeWeather Events in the Context of Climate Change. National Academyof Sciences, page 200, 2016. doi: 10.17226/21852.R. H. Reichle, Q. Liu, R. D. Koster, C. S. Draper, S. P. Mahanama, and G. S.Partyka. Land Surface Precipitation in MERRA-2. Journal of Climate,30(5):1643–1664, 2017. ISSN 08948755. doi: 10.1175/JCLI-D-16-0570.1.M. M. Rienecker, M. J. Suarez, R. Gelaro, R. Todling, and J. Bacmeister.MERRA: NASA’s Modern-Era Retrospective Analysis for Research andApplications. Journal of Climate, 24:3624–3648, 2011.C. Rosenzweig, a. Iglesius, X. B. Yang, P. R. Epstein, and E. Chivian.Climate change and extreme weather events: Implications for food pro-duction, plant diseases, and pests. Global Change and Human Health, 2:90–104, 2001.A. Ruiz-Barradas and S. Nigam. Warm season rainfall variability over theU.S. Great Plains in observations, NCEP and ERA-40 reanalyses, andNCAR and NASA atmospheric model simulations. Journal of Climate,18(11):1808–1830, 2005. ISSN 08948755. doi: 10.1175/JCLI3343.1.171BibliographyS. Saha, S. Nadiga, C. Thiaw, J. Wang, W. Wang, Q. Zhang, H. M.V. D. Dool, H.-L. Pan, S. Moorthi, D. Behringer, D. Stokes, M. Pena,S. Lord, G. White, W. Ebisuzaki, P. Peng, and P. Xie. The NCEP Cli-mate Forecast System. Journal of Climate, 19:3483–3517, 2006. doi:10.1175/2010Bams3001.1.S. Saha, S. Moorthi, H.-L. Pan, X. Wu, J. Wang, S. Nadiga, P. Tripp,R. Kistler, J. Woollen, D. Behringer, H. Liu, D. Stokes, R. Grumbine,G. Gayno, J. Wang, Y.-T. Hou, H.-Y. Chuang, H.-M. H. Juang, J. Sela,M. Iredell, R. Treadon, D. Kleist, P. Van Delst, D. Keyser, J. Derber,M. Ek, J. Meng, H. Wei, R. Yang, S. Lord, H. Van Den Dool, A. Ku-mar, W. Wang, C. Long, M. Chelliah, Y. Xue, B. Huang, J.-K. Schemm,W. Ebisuzaki, R. Lin, P. Xie, M. Chen, S. Zhou, W. Higgins, C.-Z. Zou,Q. Liu, Y. Chen, Y. Han, L. Cucurull, R. W. Reynolds, G. Rutledge, andM. Goldberg. The NCEP Climate Forecast System Reanalysis. Bulletinof the American Meteorological Society, 91(8):1015–1057, aug 2010. ISSN0003-0007. doi: 10.1175/2010BAMS3001.1.S. Saha, S. Moorthi, X. Wu, J. Wang, S. Nadiga, P. Tripp, D. Behringer,Y. T. Hou, H. Y. Chuang, M. Iredell, M. Ek, J. Meng, R. Yang, M. P.Mendez, H. Van Den Dool, Q. Zhang, W. Wang, M. Chen, and E. Becker.The NCEP climate forecast system version 2. Journal of Climate, 27(6):2185–2208, 2014. ISSN 08948755. doi: 10.1175/JCLI-D-12-00823.1.J. Saltman. BY THE NUMBERS: Winter’s impactcostly on Alex Fraser and Port Mann bridges, 2017.URL https://vancouversun.com/news/local-news/by-the-numbers-winters-costly-impact-on-alex-fraser-and-port-mann-bridges.J. Schaefer. The Critical Success as an Indicator of Warning Skill. Weatherand Forecasting, 5:570–575, 1990.M. Scheuerer and T. M. Hamill. Statistical Postprocessing of Ensemble Pre-cipitation Forecasts by Fitting Censored, Shifted Gamma Distributions*.Monthly Weather Review, 143(11):4578–4596, 2015. ISSN 0027-0644. doi:10.1175/MWR-D-15-0061.1. URL http://journals.ametsoc.org/doi/10.1175/MWR-D-15-0061.1.M. Schirmer and B. Jamieson. Verification of analysed and forecasted winterprecipitation in complex terrain. Cryosphere, 9(2):587–601, 2015. ISSN19940424. doi: 10.5194/tc-9-587-2015.172BibliographyJ. K. Searcy and C. H. Hardison. Double-Mass Curves. Manual of Hydrology:Part 1. General Surface-Water Techniques, page 41, 1960.M. C. Serreze, A. P. Barrett, and F. Lo. Northern High-Latitude Precip-itation as Depicted by Atmospheric Reanalyses and Satellite Retrievals.Monthly Weather Review, 133(12):3407–3430, 2005. ISSN 0027-0644. doi:10.1175/MWR3047.1. URL https://journals.ametsoc.org/doi/pdf/10.1175/MWR3047.1.D. Siuta, G. West, R. Stull, and T. Nipen. Calibrated Probabilistic Hub-Height Wind Forecasts in Complex Terrain. Weather and Forecasting,32(2):555–577, 2017. ISSN 0882-8156. doi: 10.1175/WAF-D-16-0137.1.URL http://journals.ametsoc.org/doi/10.1175/WAF-D-16-0137.1.B. L. Smith, S. E. Yuter, P. J. Neiman, and D. E. Kingsmill. Water VaporFluxes and Orographic Precipitation over Northern California Associatedwith a Landfalling Atmospheric River. Monthly Weather Review, 138(1):74–100, 2010. ISSN 0027-0644. doi: 10.1175/2009MWR2939.1. URLhttp://journals.ametsoc.org/doi/abs/10.1175/2009MWR2939.1.R. B. Smith, Q. Jiang, M. G. Fearon, P. Tabary, M. Dorninger, J. D. Doyle,and R. Benoit. Orographic precipitation and air mass transformation: AnAlpine example. Quarterly Journal of the Royal Meteorological Society,129(588 PART B):433–454, 2003. ISSN 00359009. doi: 10.1256/qj.01.212.G. J. Spagnol. Rocketsonde buoy system : observing system simulationexperiments. Retrospective Theses and Dissertations, 1919-2007, 2005.doi: 10.14288/1.0052550.K. Stahl, R. D. Moore, J. a. Floyer, M. G. Asplin, and I. G. McKendry. Com-parison of approaches for spatial interpolation of daily air temperature ina large region with complex topography and highly variable station den-sity. Agricultural and Forest Meteorology, 139(3-4):224–236, 2006. ISSN01681923. doi: 10.1016/j.agrformet.2006.07.004.L. Stefanova, V. Misra, S. Chan, M. Griffin, J. J. O’Brien, and T. J.Smith III. A proxy for high-resolution regional reanalysis for the South-east United States: assessment of precipitation variability in dynami-cally downscaled reanalyses. Climate Dynamics, 38(11-12):2449–2466,jun 2012. ISSN 0930-7575. doi: 10.1007/s00382-011-1230-y. URLhttp://link.springer.com/10.1007/s00382-011-1230-y.173BibliographyD. B. Stephenson. Use of the “Odds Ratio” for Diagnosing Fore-cast Skill. Weather and Forecasting, 15(2):221–232, 2000. ISSN0882-8156. doi: 10.1175/1520-0434(2000)015〈0221:UOTORF〉2.0.CO;2. URL http://journals.ametsoc.org/doi/abs/10.1175/1520-0434{%}282000{%}29015{%}3C0221{%}3AUOTORF{%}3E2.0.CO{%}3B2.D. B. Stephenson, B. Casati, C. Ferro, and C. Wilson. The extreme depen-dency score: a non-vanishing measure for forecast of rare events. Meteo-rological Applications, 15:41–50, 2008. ISSN 13504827. doi: 10.1002/met.J. Stoklosa, C. Daly, S. D. Foster, M. B. Ashcroft, and D. I. Warton. Aclimate of uncertainty: Accounting for error in climate variables for speciesdistribution models. Methods in Ecology and Evolution, 6(4):412–423,2015. ISSN 2041210X. doi: 10.1111/2041-210X.12217.F. Sun, M. L. Roderick, and G. D. Farquhar. Rainfall statistics, stationarity,and climate change. Proceedings of the National Academy of Sciences,(23):201705349, 2018a. ISSN 0027-8424. doi: 10.1073/pnas.1705349115.URL http://www.pnas.org/lookup/doi/10.1073/pnas.1705349115.Q. Sun, C. Miao, Q. Sun, and C. Miao. Extreme Rainfall (R20mm,RX5day) in Yangtze–Huai, China, in June–July 2016: The Role ofENSO and Anthropogenic Climate Change. Bulletin of the AmericanMeteorological Society, 99(1):S102–S106, jan 2018b. ISSN 0003-0007.doi: 10.1175/BAMS-D-17-0091.1. URL http://journals.ametsoc.org/doi/10.1175/BAMS-D-17-0091.1.W. Taesombat and N. Sriwongsitanon. Areal rainfall estimation using spa-tial interpolation techniques. ScienceAsia, 35(3):268–275, 2009. ISSN15131874. doi: 10.2306/scienceasia1513-1874.2009.35.268.Y. Takeuchi, H. Eito, T. Fujieda, and T. Fujita. Outline of the operationalnumerical weather prediction of the Japan Meteorological Agency. Tech-nical Report March, Japan Meteorological Agency, 2013.O. Talagrand, B. Vautard, and B. Strauss. Evaluation of probabilistic pre-diction systems. In Proc. Workshop on Predictability, pages ECMWF,1–25, Reading, UK, 1999.H. Theil. A rank-invariant method of linear and polynomial regression anal-ysis. Proc. Ned. Akad. Wet, 53(2):1397–1412, 1950. ISSN 2008-8140(Print). doi: 10.1007/978-94-011-2546-8.174BibliographyJ.-N. The´paut. Assimilating only surface pressure observations in 3D- and4D-Var. In Proceedings of the ECMWF workshop on atmospheric reanal-ysis, number June, pages 1–44, Reading, UK, 2006.T. L. Thorarinsdottir and T. Gneiting. Probabilistic Forecasts of WindSpeed : Ensemble Model Output Statistics using Heteroskedastic Cen-sored Regression. J.R. Statist. Soc. A, 173:371–388, 2010.P. E. Thornton, S. W. Running, and M. A. White. Generating surfaces ofdaily meteorological variables in complex terrain. Journal of Hydrology,190:214–251, 1997.Z. Toth and E. Kalnay. Ensemble Forecasting at NCEP and theBreeding Method. Monthly Weather Review, 125(12):3297–3319,1997. ISSN 0027-0644. doi: 10.1175/1520-0493(1997)125〈3297:EFANAT〉2.0.CO;2. URL http://journals.ametsoc.org/doi/abs/10.1175/1520-0493{%}281997{%}29125{%}3C3297{%}3AEFANAT{%}3E2.0.CO{%}3B2.K. E. Trenberth and L. Smith. The mass of the atmosphere: A constraint onglobal analyses. Journal of Climate, 18(6):864–875, 2005. ISSN 08948755.doi: 10.1175/JCLI-3299.1.J. W. Trubilowicz, J. M. Shea, G. Jost, and R. D. Moore. Suitability ofNorth American Regional Reanalysis ( NARR ) output for hydrologicmodelling and analysis in mountainous terrain. 2347(March):2332–2347,2016. doi: 10.1002/hyp.10795.S. M. Uppala, P. W. K˚allberg, A. J. Simmons, U. Andrae, V. da Costa Bech-told, M. Fiorino, J. K. Gibson, J. Haseler, a. Hernandez, G. a. Kelly, X. Li,K. Onogi, S. Saarinen, N. Sokka, R. P. Allan, E. Andersson, K. Arpe,M. a. Balmaseda, a. C. Beljaars, L. van de Berg, J. Bidlot, N. Bormann,S. Caires, F. Chevallier, a. Dethof, M. Dragosavac, M. Fisher, M. Fuentes,S. Hagemann, E. Ho´lm, B. J. Hoskins, L. Isaksen, P. a.E.M. Janssen,R. Jenne, a. P. McNally, J. F. Mahfouf, J. J. Morcrette, N. a. Rayner,R. W. Saunders, P. Simon, a. Sterl, K. E. Trenberth, a. Untch, D. Vasil-jevic, P. Viterbo, and J. Woollen. The ERA-40 re-analysis. QuarterlyJournal of the Royal Meteorological Society, 131(612):2961–3012, 2005.ISSN 00359009. doi: 10.1256/qj.04.176.G. J. van Oldenborgh, F. J. Doblas Reyes, S. S. Drijfhout, and E. Hawkins.Reliability of regional climate model trends. Environmental Research Let-ters, 8(1):014055, mar 2013. ISSN 1748-9326. doi: 10.1088/1748-9326/8/175Bibliography1/014055. URL http://stacks.iop.org/1748-9326/8/i=1/a=014055?key=crossref.2aa3fd77fc365abd77e9d8773bb66738.C. Vancouver. Wicked winter leaves behind 15,000 potholes in Vancou-ver — CTV Vancouver News, 2017. URL https://bc.ctvnews.ca/wicked-winter-leaves-behind-15-000-potholes-in-vancouver-1.3348547?autoPlay=true.S. M. Vicente-Serrano, S. Beguer´ıa, J. I. Lo´pez-Moreno, M. a. Garc´ıa-Vera,and P. Stepanek. A complete daily precipitation database for northeastSpain: Reconstruction, quality control, and homogeneity. InternationalJournal of Climatology, 30(8):1146–1163, 2010. ISSN 08998418. doi: 10.1002/joc.1850.L. a. Vincent. A technique for the identification of inhomogeneities in Cana-dian temperature series. Journal of Climate, 11(5):1094–1104, 1998. ISSN08948755. doi: 10.1175/1520-0442(1998)011〈1094:ATFTIO〉2.0.CO;2.L. a. Vincent and D. W. Gullett. Canadian historical and homogeneoustemperature datasets for climate change analyses. International Journalof Climatology, 19(12):1375–1388, 1999. ISSN 08998418. doi: 10.1002/(SICI)1097-0088(199910)19:12〈1375::AID-JOC427〉3.0.CO;2-0.L. a. Vincent and E´. Mekis. Changes in Daily and Extreme Temper-ature and Precipitation Indices for Canada over the Twentieth Cen-tury. Atmosphere-Ocean, 44(2):177–193, 2006. ISSN 0705-5900. doi:10.3137/ao.440205.L. a. Vincent, X. Zhang, B. R. Bonsal, and W. D. Hogg. Homogenization ofdaily temperatures over Canada. Journal of Climate, 15(11):1322–1334,2002. ISSN 08948755. doi: 10.1175/1520-0442(2002)015〈1322:HODTOC〉2.0.CO;2.L. a. Vincent, E. J. Milewska, R. Hopkinson, and L. Malone. Bias inminimum temperature introduced by a redefinition of the climatologi-cal day at the canadian synoptic stations. Journal of Applied Meteo-rology and Climatology, 48(10):2160–2168, 2009. ISSN 15588424. doi:10.1175/2009JAMC2191.1.L. a. Vincent, X. L. Wang, E. J. Milewska, H. Wan, F. Yang, and V. Swail. Asecond generation of homogenized Canadian monthly surface air temper-ature for climate trend analysis. Journal of Geophysical Research Atmo-spheres, 117(17):1–13, 2012. ISSN 01480227. doi: 10.1029/2012JD017859.176BibliographyW. Wang, P. Xie, S. H. Yoo, Y. Xue, A. Kumar, and X. Wu. An assessmentof the surface climate in the NCEP climate forecast system reanalysis.Climate Dynamics, 37(7-8):1601–1620, 2011. ISSN 09307575. doi: 10.1007/s00382-010-0935-7.X. Wang, D. Parrish, D. Kleist, and J. Whitaker. GSI 3DVar-Based En-semble–Variational Hybrid Data Assimilation for NCEP Global Fore-cast System: Single-Resolution Experiments. Monthly Weather Re-view, 141(11):4098–4117, 2013. ISSN 0027-0644. doi: 10.1175/MWR-D-12-00141.1. URL http://journals.ametsoc.org/doi/abs/10.1175/MWR-D-12-00141.1.X. L. Wang, Q. H. Wen, and Y. Wu. Penalized maximal t test for detectingundocumented mean change in climate data series. Journal of AppliedMeteorology and Climatology, 46(6):916–931, 2007. ISSN 15588424. doi:10.1175/JAM2504.1.M. Wei, Z. Toth, R. Wobus, and Y. Zhu. Initial perturbations based on theensemble transform (ET) technique in the NCEP global operational fore-cast system. Tellus, Series A: Dynamic Meteorology and Oceanography, 60A(1):62–79, 2008. ISSN 02806495. doi: 10.1111/j.1600-0870.2007.00273.x.A. T. Werner and A. J. Cannon. Hydrologic extremes - An intercomparisonof multiple gridded statistical downscaling methods. Hydrology and EarthSystem Sciences, 20(4):1483–1508, 2016. ISSN 16077938. doi: 10.5194/hess-20-1483-2016.S. Westra, L. V. Alexander, F. W. Zwiers, S. Westra, L. V. Alexander,and F. W. Zwiers. Global Increasing Trends in Annual Maximum DailyPrecipitation. Journal of Climate, 26(11):3904–3918, jun 2013. ISSN0894-8755. doi: 10.1175/JCLI-D-12-00502.1. URL http://journals.ametsoc.org/doi/abs/10.1175/JCLI-D-12-00502.1.T. Weusthoff, F. Ament, M. Arpagaus, and M. W. Rotach. Assessing theBenefits of Convection-Permitting Models by Neighborhood Verification:Examples from MAP D-PHASE. Monthly Weather Review, 138(9):3418–3433, 2010. ISSN 0027-0644. doi: 10.1175/2010MWR3380.1. URL http://journals.ametsoc.org/doi/abs/10.1175/2010MWR3380.1.J. S. Whitaker, G. P. Compo, and J.-N. The´paut. A Comparison of Varia-tional and Ensemble-Based Data Assimilation Systems for Reanalysis ofSparse Observations. Monthly Weather Review, 137(6):1991–1999, 2009.ISSN 0027-0644. doi: 10.1175/2008MWR2781.1.177BibliographyT. White, J. Wolf, F. Anslow, A. Werner, and R. Creative. Indicators ofClimate Change for British Columbia: 2016 Update. Technical report,Ministry of Environment - BC, 2016.D. Wilks. Statistical Methods in the Atmospheric Sciences, volume100. Elsevier-International geophysics series, 3rd edition, 2011. ISBN9780123850225. doi: 10.1016/B978-0-12-385022-5.00022-1.D. Wilks. The Stippling Shows Statistically Significant Grid Points: HowResearch Results are Routinely Overstated and Overinterpreted, andWhat to Do about It. Bulletin of the American Meteorological Society,97:2263–2274, 2016. doi: 10.1175/BAMS-D-00267.1.P. Xie and P. a. Arkin. Global Precipitation: A 17-Year Monthly Anal-ysis Based on Gauge Observations, Satellite Estimates, and NumericalModel Outputs. Bulletin of the American Meteorological Society, 78(11):2539–2558, 1997. ISSN 00030007. doi: 10.1175/1520-0477(1997)078〈2539:GPAYMA〉2.0.CO;2.P. Xie, M. Chen, S. Yang, A. Yatagai, T. Hayasaka, Y. Fukushima, andC. Liu. A Gauge-Based Analysis of Daily Precipitation over East Asia.Journal of Hydrometeorology, 8(3):607–626, 2007. ISSN 1525-755X. doi:10.1175/JHM583.1.S. Yue, P. Pilon, B. Phinney, and G. Cavadias. The influence of autocor-relation on the ability to detect trend in hydrological series. HydrologicalProcesses, 16(9):1807–1829, 2002. ISSN 08856087. doi: 10.1002/hyp.1095.X. Zhang, L. a. Vincent, W. D. Hogg, and A. Niitsoo. Temperature andprecipitation trends in Canada during the 20th century. Atmosphere -Ocean, 38(3):395–429, 2000. ISSN 07055900. doi: 10.1080/07055900.2000.9649654.X. Zhang, W. D. Hogg, and E´. Mekis. Spatial and Temporal Character-istics of Heavy Precipitation Events over Canada. Journal of Climate,14(9):1923–1936, 2001. ISSN 0894-8755. doi: 10.1175/1520-0442(2001)014〈1923:SATCOH〉2.0.CO;2.X. Zhou, Y. Zhu, D. Hou, Y. Luo, J. Peng, and R. Wobus. Performance of theNew NCEP Global Ensemble Forecast System in a Parallel Experiment.Weather and Forecasting, pages WAF–D–17–0023.1, 2017. ISSN 0882-8156. doi: 10.1175/WAF-D-17-0023.1. URL http://journals.ametsoc.org/doi/10.1175/WAF-D-17-0023.1.178H. Zhu, M. C. Wheeler, A. H. Sobel, D. Hudson, H. Zhu, M. C. Wheeler,A. H. Sobel, and D. Hudson. Seamless Precipitation Prediction Skill inthe Tropics and Extratropics from a Global Model. Monthly WeatherReview, 142(4):1556–1569, apr 2014. ISSN 0027-0644. doi: 10.1175/MWR-D-13-00222.1.Y. Zhu and Y. Luo. Precipitation Calibration Based on the Frequency-Matching Method. Weather and Forecasting, 30(5):1109–1124, 2015. ISSN0882-8156. doi: 10.1175/WAF-D-13-00049.1.E. Zsoter. Recent developments in extreme weather forecasting. ECMWFNewsletter, Spring(107):8–17, 2006.F. W. Zwiers and V. V. Kharin. Estimating extremes in transient climatechange simulations. Journal of Climate, (18):43, 2005. URL http://www.cccma.ec.gc.ca/papers/fzwiers/kz{_}2003.pdf.V. V. K. F. W. Zwiers and X. Z. M. Wehner. Changes in temperature andprecipitation extremes in the CMIP5 ensemble. Climate Change, 119:345–357, 2013. doi: 10.1007/s10584-013-0705-8.179Appendix ASurface weather stationsTable A.1: Surface weather stations description. Station name and abbre-viation, longitude and latitude in degrees, elevation in meters, variable aseither 2-m temperature (T) or 1-day accumulated precipitation (p), networkand whether it is a SYNOP station.Station Abbrev Lon Lat Elev Var Network1 Abbotsford YXX -122.4 49.0 59 T, p ECCC(S)2 Agassiz WZA -121.8 49.2 15 T, p ECCC(S)3 Alouette ALU -122.5 49.3 125 T, p BCH4 Atlin ATL -133.70 59.6 74 T, p ECCC5 Barkerville BAR -121.5 53.1 128 p ECCC6 Bella Coola YBD -126.6 52.4 36 T, p ECCC(S)7 Blind Channel BCH -125.5 50.4 23 T ECCC8 Blue River YCP -119.3 52.1 68 T, p ECCC(S)9 Castlegar YCG -117.6 49.3 496 T, p ECCC(S)10 Chatham Point WFM -125.4 50.3 23 p ECCC11 Clowhom Falls CLO -123.5 49.7 10 T, p BCH12 Comox YQQ -124.9 49.7 26 T, p ECCC(S)13 Coquitlam CQM -122.8 49.5 290 T, p BCH14 Cortes Island CIT -124.9 50.1 15 p ECCCTiber Bay15 Cranbrook YXC -115.8 49.6 939 T, p ECCC(S)16 Creston WJR 116.5 49.1 597 p ECCC17 Darfield DAR -120.2 51.3 412 T, p ECCC18 Dawson’s Creek YDQ -120.2 55.8 655 T ECCC(S)19 Dease Lake WDL -130.0 58.4 807 T ECCC(S)20 Dryad Point DAU -128.1 52.2 4 p ECCC21 Duncan Kelvin DKC -123.7 48.7 103 p ECCCCreek22 Estevan Point WEB -126.6 49.4 7 T, p ECCC23 Fernie FER -115.1 49.5 1001 p ECCCContinued on next page180Appendix A. Surface weather stationsTable A.1 – continued from previous pageStation Abbrev Lon Lat Elev Var Network24 Fort Nelson YYE -122.6 58.8 382 T, p ECCC(S)25 Fort St James VFS -124.3 54.5 686 T, p ECCC26 Fort St John’s YXJ -120.7 56.2 695 T, p ECCC(S)27 Germansen YGS -124.7 55.8 766 p ECCCLanding28 Glacier GLA -117.5 51.3 1323 T,p ECCC29 Gold Creek GOC -122.5 49.45 794 T, p BCH30 Golden YGE -117.0 51.3 785 T, p ECCC31 Grand Forks GRF -118.5 49.0 532 T ECCC32 Kamloops YKA -120.4 50.7 345 T, p ECCC(S)33 Kaslo KAS -116.9 49.9 591 T ECCC34 Kelowna YLW -119.6 49.8 417 T ECCC(S)35 Kelowna KQG -119.6 49.8 417 p ECCCQuails Gate36 Laidlaw LAI -121.6 49.4 27 p ECCC37 Langara WLA -133.1 54.3 41 T, p ECCC38 Little Qualicum VOQ -124.5 49.4 30 p ECCCHatchery39 McInnes Island WMS -128.7 52.3 26 T , p ECCC(S)40 Merritt VME -120.8 50.1 609 T, p ECCC41 Mica Dam MCD -118.6 52.1 579 T, p BCH42 Nanaimo City YCD -124.0 49.2 114 p ECCCYard43 Nass Camp NAC -129.0 55.2 290 p ECCC44 Oliver OLI -119.5 49.2 297 T, p ECCC45 Pachena Point PAP -125.1 48.7 37 p ECCC46 Penticton YYF -119.6 49.5 344 T, p ECCC(S)47 Port Alice POA -127.5 50.4 21 p ECCC48 Port Hardy YZT -127.4 50.7 22 T, p ECCC(S)49 Prince George YXS -122.7 53.9 691 T, p ECCC(S)50 Prince Rupert YPR -130.5 54.3 35 T ECCC(S)51 Princeton YDC -120.5 49.5 700 T, p ECCC(S)52 Quesnel YQZ -122.5 53.0 545 T ECCC(S)53 Quatsino WIF -127.7 50.5 8 T, p ECCC54 Quinsam River QUI -125.3 50 46 T, p ECCCHatcheryContinued on next page181Appendix A. Surface weather stationsTable A.1 – continued from previous pageStation Abbrev Lon Lat Elev Var Network55 Saanichton SAN -123.4 48.6 61 p ECCC56 Salmon Arm WSL -119.2 50.7 527 T, p ECCC57 Sandspit YZP -131.8 53.3 6 T, p ECCC(S)58 Shawnigan Lake SHL -123.6 48.6 138 T, p ECCC59 Smithers YYD -127.2 54.8 522 T, p ECCC(S)60 Stave Ridge STV -122.4 49.6 930 T, p BCHUpper61 Stewart ZST -130 55.9 7 T, p ECCC(S)62 Tatlayoko XTL -124.4 51.7 870 T ECCCLake63 Terrace YXT -128.6 54.5 7 T, p ECCC(S)64 Tofino YAZ -125.8 49.1 24 p ECCC(S)65 Ucluelet UKC -125.5 48.9 30 p ECCCKennedy Camp66 Vancouver YVR -123.2 49.2 4 T, p ECCC(S)67 Cheakamus CMU -123.1 50.1 880 T, p BCHUpper CMU68 Vavenby VAV -119.8 51.6 445 T, p ECCC69 Vernon WJV -119.3 50.3 427 T, p ECCC(S)70 Victoria YYJ -123.4 48.6 19 T, p ECCC(S)71 Wahleach Jones WAH -121.6 49.2 641 T, p BCHReservoir72 Wasa WAS -115.6 49.8 930 p ECCC73 Wistaria WIS -126.2 53.8 863 T ECCC74 William’s Lake WLK -122.1 52.2 940 T, p ECCC(S)75 Wolf Ridge WOL -125.7 49.7 1490 T, p BCHUpper182AppendixA.SurfaceweatherstationsFigure A.1: Elevation of 2-m temperature and 1-day accumulated precipitation stations.183Appendix BTrain and test stationsTable B.1: Train surface weather stations description. Station name andabbreviation, longitude and latitude in degrees, elevation in meters, vari-able as either 2-m temperature (T) or 1-day accumulated precipitation (p),network station.Station Abbrev Lon Lat Elev Var Network1 Alouette ALU -122.5 49.3 125 T BCH2 Barkerville BAR -121.5 53.1 128 p ECCC3 Bella Coola YBD -126.6 52.4 36 T, p ECCC4 Blue River YCP -119.3 52.1 68 p ECCC5 Chatham Point WFM -125.4 50.3 23 p ECCC6 Coquitlam CQM -122.8 49.5 290 T, p BCH7 Cortes Island CIT -124.9 50.1 15 p ECCCTiber Bay8 Creston WJR 116.5 49.1 597 T ECCC8 Darfield DAR -120.2 51.3 412 T, p ECCC10 Dawson’s Creek YDQ -120.2 55.8 655 T ECCC11 Duncan Kelvin DKC -123.7 48.7 103 p ECCCCreek12 Fernie FER -115.1 49.5 1001 p ECCC13 Fort Nelson YYE -122.6 58.8 382 p ECCC14 Fort St James VFS -124.3 54.5 686 T, p ECCC15 Fort St John’s YXJ -120.7 56.2 695 T, p ECCC16 Germansen YGS -124.7 55.8 766 T ECCCLanding17 Glacier GLA -117.5 51.3 1323 T ECCC18 Gold Creek GOC -122.5 49.45 794 T, p BCH19 Golden YGE -117.0 51.3 785 p ECCC19 Ingenika ING -125.1 56.73 711 T BCH20 Kamloops YKA -120.4 50.7 345 T, p ECCC21 Kaslo KAS -116.9 49.9 591 T ECCCContinued on next page184Appendix B. Train and test stationsTable B.1 – continued from previous pageStation Abbrev Lon Lat Elev Var Network22 Kelowna YLW -119.6 49.8 417 T ECCC23 Langara WLA -133.1 54.3 41 T, p ECCC24 McInnes Island WMS -128.7 52.3 26 T , p ECCC25 Mica Dam MCD -118.6 52.1 579 p BCH26 Nanaimo City YCD -124.0 49.2 114 p ECCCYard27 Oliver OLI -119.5 49.2 297 T ECCC28 Pachena Point PAP -125.1 48.7 37 p ECCC29 Penticton YYF -119.6 49.5 344 T ECCC30 Port Alice POA -127.5 50.4 21 p ECCC31 Port Hardy YZT -127.4 50.7 22 T ECCC32 Prince Rupert YPR -130.5 54.3 35 T ECCC33 Princeton YDC -120.5 49.5 700 T ECCC34 Quesnel YQZ -122.5 53.0 545 T ECCC35 Quatsino WIF -127.7 50.5 8 p ECCC36 Quinsam River QUI -125.3 50 46 T ECCCHatchery37 Saanichton SAN -123.4 48.6 61 p ECCC38 Salmon Arm WSL -119.2 50.7 527 T ECCC39 Sandspit YZP -131.8 53.3 6 p ECCC40 Shawnigan Lake SHL -123.6 48.6 138 p ECCC41 Smithers YYD -127.2 54.8 522 T, p ECCC42 Stave Ridge STV -122.4 49.6 930 T, p BCHUpper43 Stewart ZST -130 55.9 7 p ECCC44 Tatlayoko XTL -124.4 51.7 870 T ECCCLake45 Ucluelet UKC -125.5 48.9 30 p ECCCKennedy Camp46 Vancouver YVR -123.2 49.2 4 T ECCC47 Cheakamus CMU -123.1 50.1 880 p BCHUpper CMU48 Vavenby VAV -119.8 51.6 445 T, p ECCC49 Vernon WJV -119.3 50.3 427 T, p ECCC50 Victoria YYJ -123.4 48.6 19 T, p ECCC51 Wasa WAS -115.6 49.8 930 p ECCCContinued on next page185Appendix B. Train and test stationsTable B.1 – continued from previous pageStation Abbrev Lon Lat Elev Var Network52 Wistaria WIS -126.2 53.8 863 T ECCC53 William’s Lake WLK -122.1 52.2 940 T, p ECCCTable B.2: Test surface weather stations description. Station name and ab-breviation, longitude and latitude in degrees, elevation in meters, variable aseither 2-m temperature (T) or 1-day accumulated precipitation (p), networkstation.Station Abbrev Lon Lat Elev Var Network1 Agassiz WZA -121.8 49.2 15 T,p ECCC2 Barnes Creek BRN -118.35 50.01 1620 T BCH3 Blind Channel BCH -125.5 50.4 23 T ECCC4 Blue River YCP -119.3 52.1 68 T ECCC5 Cheakamus Creek CHK -123.03 50.08 640 T,p BCH6 Comox YQQ -124.9 49.7 26 T,p ECCC7 Cranbrook YXC -115.8 49.6 939 T ECCC8 Creston WJR 116.5 49.1 597 p ECCC9 Dease Lake WDL -130.0 58.4 807 T ECCC10 Duncan Kelvin DKC -123.7 48.7 103 T,p ECCCCreek11 Estevan Point WEB -126.6 49.4 7 T,p ECCC12 Fort Nelson YYE -122.6 58.8 382 T ECCC13 Golden YGE -117.0 51.3 785 T ECCC14 Goldstream GOL -118.6 51.67 600 p BCH15 Germansen YGS -124.7 55.8 766 p ECCCLanding16 Glacier GLA -117.5 51.3 1323 T,p ECCC17 Grand Forks GRF -118.5 49.0 532 T ECCC18 Kemano KEM -127.9 53.6 87 p ECCC19 Kelowna KQG -119.6 49.8 417 p ECCC20 Laidlaw LAI -121.6 49.4 27 p ECCC21 Little Qualicum VOQ -124.5 49.4 30 p ECCCHatchery22 Little Qualicum VOQ -124.5 49.4 30 p ECCCContinued on next page186Appendix B. Train and test stationsTable B.2 – continued from previous pageStation Abbrev Lon Lat Elev Var NetworkHatchery23 Oliver OLI -119.5 49.2 297 p ECCC24 Ootsa VSL -126 53.8 861 p ECCC25 Pack Lake PAK -123.04 55 675 p ECCC26 Penticton YYF -119.6 49.5 344 p ECCC27 Port Hardy YZT -127.4 50.7 22 p ECCC28 Prince George YXS -122.7 53.9 691 T ECCC29 Princeton YDC -120.5 49.5 700 p ECCC30 Quatsino WIF -127.7 50.5 8 T ECCC31 Quinsam River QUI -125.3 50 46 p ECCCHatchery32 Sandspit YZP -131.8 53.3 6 T ECCC33 Shawnigan Lake SHL -123.6 48.6 138 T ECCC34 Stave Ridge STV -122.4 49.6 930 p BCH35 Stewart ZST -130 55.9 7 T, ECCC36 Terrace YXT -128.6 54.5 7 T, p ECCC37 Tofino YAZ -125.8 49.1 24 p ECCC38 Vancouver YVR -123.2 49.2 4 p ECCC187

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0378419/manifest

Comment

Related Items