Dendrohydrological reconstruction and hydroclimatic variability in southwestern British Columbia, Canada

(1)

Dendrohydrological reconstruction and hydroclimatic variability in

southwestern British Columbia, Canada

by

Bryan J. Mood

B.Sc., Mount Allison University, 2013 M.Sc., University of Victoria, 2015

A Dissertation Submitted in Partial Fulfillment of the Requirements for the Degree of

DOCTOR OF PHILOSOPHY

in the Department of Geography

(2)

ii

Supervisory Committee

Dendrohydrological reconstruction and hydroclimatic variability in southwestern British Columbia, Canada

by Bryan J. Mood

B.Sc., Mount Allison University, 2013 M.Sc., University of Victoria, 2015

Supervisory Committee

Dr. Dan J. Smith (Department of Geography)

Supervisor

Dr. David Atkinson (Department of Geography)

Departmental Member

Dr. Tobi Gardner (Victoria Capital Regional District)

(3)

iii

Abstract

Supervisory Committee

Dr. Dan J. Smith (Department of Geography)

Supervisor

Dr. David Atkinson (Department of Geography)

Departmental Member

Dr. Tobi Gardner (Victoria Capital Regional District)

Outside Member

The hydrology of southwestern British Columbia is influenced by the region’s mountainous topography and climate oscillations generated from the Pacific Ocean. While much of the region is characterized as a temperate rainforest, recent summers are defined by record-breaking droughts that focus attention on the threat to regional water supply security likely to accompany future climate changes.

The limited length and distribution of hydrological records in southwestern British Columbia provide poor context for resource managers tasked with developing policy and water management strategies. The purpose of the dissertation was to describe long-term variability in several key hydroclimatic variables and hydroecological

interactions that may be used in updated water resource policy and management strategies. Specifically, the research focused on developing long-term proxy records of April 1 snow water equivalent (SWE), summer streamflow, spring lake levels, and salmon abundance from tree ring records. A secondary goal of the dissertation was to identify the role and influence of several key climate oscillations on regional long-term hydroclimatic and ecological variability.

Freshet contributions from melting snow are critical for sustained summer

streamflow in southwestern British Columbia. Even so, few manual snow survey stations exist within the region are of sufficient length to understand the full range of natural SWE variability. Long-term April 1 SWE records were constructed by establishing statistical relationships with the radial growth of high-elevation trees and April 1 SWE records. Explaining 51% and 73% of the total variance in the instrumental SWE records in coastal and continental settings, the reconstructions provide high-resolution

descriptions of April 1 SWE over the past three centuries and help position the remainder of the dissertation. Negative phases of the Pacific Decadal Oscillation (PDO) and El Niño-Southern Oscillation (ENSO) were shown to strongly influence April 1 SWE totals. Both reconstructions illustrate repeated step-changes in April 1 SWE during the last 300 years and show that coastal areas may be more sensitive to annual variability than snowpack that accumulates in more continental locations.

Water shortages in the Metro Vancouver area in recent summers are linked to low total winter snowpack and early spring melt. Dendrohydrological analysis of dry-season

(4)

iv

streamflow was conducted to determine if the instrumental range has been

underestimated over the past several centuries. A regionalized record of July-August streamflow for the Capilano and Seymour watersheds, which supply the Metro

Vancouver area, was modelled from present to 1711using tree-rings. Explaining 54% of total variance over the instrumental period, the models show that below-average

streamflow events are becoming more frequent. When compared to those characterizing the past 300 years, streamflow totals from 1977 to present have consistently fallen well-below the average long-term discharge. Further analyses indicated that negative ENSO and PDO conditions strongly influenced July-August runoff trends since 1711, as have climate regimes related to the Pacific North American pattern (PNA).

The increased frequency in recent years of reduced summer runoff in

southwestern British Columbia has led many communities to rely on natural and dammed reservoirs to supplement their water needs. Where communities rely on natural lakes, this dependence may have socioeconomic consequences if lake levels fall below those

necessary to supply built infrastructure. Unfortunately, there are few lake level records in southwestern British Columbia and none of sufficient duration to understand the full range of variability in natural lake systems. Harrison Lake is the only natural lake with a lake level record exceeding 50 years. Using the average April water level dataset, a dendrohydrological model was constructed that explained 49.5% of total variance. The model was used to reconstruct a proxy record of April water levels spanning the interval from 1711 to 1980. Averaging 9.37 m in depth, lake levels in Harrison Lake ranged from 8.9 to 10.0 m over the past 300 years. These variations were shown to be statistically associated with negative and positive phases of ENSO and positive phases of PDO. April water levels in Harrison Lake have been, on average, 0.13 m lower since the mid-1930s compared to the previous 200 years. This reduction in storage capacity amounts to a loss of almost 300-million litres of stored water since the start of instrumental records.

Salmon play a vital economic, cultural, and social role in many southwestern British Columbia communities. There is concern that salmon populations in the region are under threat, as changing climates alter and impact their spawning habitat. While various lines of research have sought to determine the response of salmon to these changing conditions, population records that extend only to 1951 hinder a complete understanding of the impacts. Two dendroecological models were constructed to provide a longer-term perspective of regional salmon-climate relationships. Explaining 48.2% and 48.9% of variance in observed Chinook and Coho salmon abundance since 1951, the models were used to construct proxy escapement records extending to the 1700s.

Spectral analysis revealed that the reconstructions account for generational life histories and that low-frequency climate variability was associated with fluctuations in abundance. Both the Chinook and the Coho reconstructions show phase dependent relationships to climate oscillations generated from the Pacific Ocean. The Coho record is strongly linked to negative winter and spring ENSO, while the Chinook record was shown to be

associated with negative PDO conditions. The identified relationships to teleconnections generated in the Pacific Ocean to our record indicates that both species are sensitive to oceanic interactions prior to entering natal habitats. Taken together, the reconstructions illustrate that the observational record encompasses a period of lower-than-average

(5)

v

abundance and that neither accounts for the full range of variability in annual abundance when considered over the past three centuries.

The proxy tree-ring records presented in this dissertation provide new information about climate-water resource relationships in southwestern British Columbia. Significant phase-dependent associations, especially to negative phases of the PDO and ENSO, were shown to exert long-term influences on the state of several critical hydroclimatic

variables over the last 300 years. Additionally, the research illustrates that over the instrumental period, both streamflow and lake volumes in the region have consistently remained below those characterizing the previous two to three centuries. These findings are of direct use to resource managers tasked with developing new policy and strategies under present and future climate change, in that they offer singular insights into the full range of natural hydroclimatic variability in southwestern British Columbia.

(6)

vi

List of Tables

Table 2.1: Descriptive information of manual snow survey stations used in this study.

Records were collected from the River Forecast Centre (2017) website. Two regions were targeted for reconstruction: (1) continental and (2) coastal. Continental locations are located east of the Pacific Ranges and characterized by lower average snowpack. Coastal stations are located in the Pacific Ranges and characterized by much greater snowpack averages (229 vs. 1367 SWE). Station ID is associated with River Forecast Centre designations; coordinates (latitude and longitude) are in decimal degrees; elevation is rounded to the nearest 10 m above sea level; mean April 1 SWE is calculated across the whole time series; span is the length of continuous April 1 SWE measurements from each station; length is the total number of years available for analysis; PCA loading represents the explanatory power of the first principal component from analysis. .……….. 19

Table 2.2: Times series information. Species/type are bold and italics represents time

series used as predictors in the reconstruction. Tree-ring site numbers in brackets are ITRDB codes. The Pacific silver rir regional index (PSFR) was developed by combining

the five site-level chronologies used a bi-weight robust mean method (Briffa and Melvin, 2011). Climate oscillation indices were collected from NOAA (2017) for teleconnections that are known to influence the overall hydroclimate of western North America. RBAR is the average value across whole index; length is the span of years used for reconstruction. For tree-ring series, the length used for reconstruction is used only for where EPS > 0.85. ………... 23

Table 2.3: Pearson correlations between model parameters used in the reconstruction. All

values shown are significant (p<0.01)………. 24

Table 2.4: Reconstruction, cross-validation, and descriptive statistics. D-W =

Durbin-Watson Statistic; VIF = variance inflation factor; SE = standard error; RE = reduction of error; RMSE = root mean squared error; CV = coefficient of variance. ………. 24

Table 2.5: Difference-of-correlations test results for Coastal and Continental SWEPC. Climate oscillation indices were grouped into negative and positive values then correlated to both the calculated and modeled Coastal and Continental SWEPC. Values shown are

p<0.05 while bold values indicate p<0.01. Where only one of two teleconnection phases

have a significant relationship, it illustrates the possibility of a non-stationary response.27

Table 3.1: Instrumental data information. Streamflow data was collected from the Water

Survey of Canada (2017); climate information was collected from Environment and Climate Change Canada (2017), and; snow survey information was collected from River Forecast Centre (2017) ..………..………. 42

Table 3.2: Times series information. Species/type are italicized and bold represents time

series used as predictors in the reconstruction. Mountain hemlock PC (PCMH) was

(10)

x

chronologies listed. Pacific silver fir did not have a sufficient number of sites to conduct PCA (<5) and site-level chronologies were combined into regional chronologies using a bi-weight robust mean method (Briffa and Melvin, 2011). Climate oscillation indices were collected from NOAA (2017) for teleconnections that are known to influence the overall hydroclimate of western North America. a = correlations among tree-ring series collected from each site; b = only the length of the tree-ring index where the expressed population signal (EPS) was ≥0.85 is documented; c_{= RBAR is the average value across}

whole index where EPS was >0.85 and PC loading is the total explained variance from PCA analysis. ………... 53

Table 3.3: Correlations between series used in the reconstruction model and snow water

equivalent (SWE) from the closest proximity station and regionalised streamflow. All values are significant (p<0.01). ………...…… 53

Durbin-Watson Statistic; VIF = variance inflation factor; SE = standard error; RE = reduction of error; RMSE = root mean squared error; CV = coefficient of variance. ………. 55

Table 3.5: Difference-of-correlations tests for measured ENSO, PDO, PNA values

against measured and reconstructed regional streamflow. Bold values indicate p<0.05;

Bold, underlined indicates p<0.01. Where only one teleconnection phase is significantly

related to QJA, it illustrates the possibility of a non-stationary response. ……….55 Table 3.6: Fifth- and ninety-fifth percentile flows from the reconstructed QJA record.

Bold represent flows during the instrumental record. Z-scores are the calculated average

for both reconstructed basins. ...…...……….. 61

Table 4.1: Instrumental data information retrieved from Environmental and Climate

Change Canada (2017), Water Survey of Canada (2017), and River Forecast Centre (2017). ……….. 81

Table 4.2: Time series information. Species/type area italicized and bold represents time

series used as predictors in the reconstruction. Mountain hemlock PC1 was developed by conducting a principal component analysis (PCA) on the five site-level chronologies listed. Climate oscillation indices were collected from NOAA (2017) for teleconnections that are known to influence the overall hydroclimate of western North America. a = correlations among tree-ring series collected from each site; b = only the length of the tree-ring index where the expressed population signal (EPS) was ≥0.85 is documented; c

= RBAR is the average value across whole index where EPS was >0.85 and PC loading indicates the variance explained by the first component of the PCA. ……… 83

Table 4.3: Correlations between model parameters, April Harrison Lake water levels,

and local snow survey stations (April 1 SWE values). Bold-underlined indicates p<0.01. ……….………..……… 86

(11)

xi

Durbin-Watson Statistic; VIF = variance inflation factor; SE = standard error; RE = reduction of error; RMSE = root mean squared error. ………. 86

Table 4.5: Difference-of-correlations tests for measured ENSO, PDO, PNA values

against modelled and measured April Harrison Lake Level. Bold indicates p<0.05. Where only one of the two teleconnection phases is significant, it illustrates a possible non-stationary response. ……….…. 88

Table 5.1: Times series information. Species/type are italicized and bold represents time

series used as predictors in the reconstruction. Mountain hemlock PC1 and PC2 were developed by conducting a principal component analysis (PCA) on the five site-level chronologies listed. Other species did not have a sufficient number of sites to conduct PCA (<5) and were combined into regional chronologies using a bi-weight robust mean method. Climate oscillation indices were collected from NOAA (2017) for

teleconnections that are known to influence the overall hydroclimate of western North America. a = correlations among tree-ring series collected from each site; b = only the length of the tree-ring index where the expressed population signal (EPS) was ≥0.85 is documented; c = RBAR is the average value across whole index where EPS was >0.85 and PC loading is the variance explained by PCA. ………..………. 106

Table 5.2: Watersheds (bold) and stream location names of Coho and Chinook salmon

escapement records collected from NuSEDS database suitable for use. Coordinates are not provided for each stream as they were not documented in the NuSEDS database. x = record used for Coho and/or Chinook principal component analysis. ………... 108

Table 5.3: Correlation values for salmon escapement and tree-ring records to

seasonalized climate oscillation indices. For salmon escapement records, only the most significant correlation was documented and others may exist (see Appendix A for full list of salmon correlations). DJF = December (of the previous year), January, and February; MAM = March, April, and May; JJA = June, July, and August; A = Annual. September, October, November was excluded from the Table 5.3 as no significant correlations were documented. ………..……… 110

Durbin-Watson Statistic; VIF = variance inflation factor; SE = standard error; RE = reduction of error; RMSE = root mean squared error; CV = coefficient of variance. ………... 114

Table 5.5: Difference-of-correlations test results for Chinook PC1 and Coho PC1.

Climate oscillation indices were grouped into negative and positive values then correlated to both the calculated and reconstructed Chinook and Coho PC1 values. Values shown are p<0.05 while bold values indicate p<0.01. Where only one of two teleconnection phases is significantly correlated, it illustrates the possibility of a non-stationary response. ………. 116

(12)

xii

List of Figures

Figure 2.1: Study area map showing all the tree-ring chronologies (white triangles) and

manual snow survey stations (continental = white square with centre black square; coastal = white circle with cross). Note that coastal snow survey stations are located in regions with high orographic precipitation (western areas) while continental snow survey stations are in the Coast Mountains rain shadow (east). ………..…. 12

Figure 2.2: General climate information for the study region. (A) Climate normals

(1981-2010) for coastal (Grouse Mountain, station ID: 1105658) and continental

(Lillooet, station ID: 1114627) regions (ECCC, 2017); (B) time series of average coastal and continental April 1 SWE records (River Forecast Centre, 2017). ...……….. 14

Figure 2.3: Time plot of calculated April 1 SWE principal components (thin black line),

modelled (thick black line), and cross-validated (thick grey line) records of (A) coastal and (B) continental. ……….. 25

Figure 2.4: Time plot of reconstructed (A) coastal and (B) continental SWE plotted as

departures from average (z-scores; thin black line) with a 5-year running mean (thick black line) and error (grey area) calculated from the cross-validation RMSE. Black-filled areas illustrate long-term changes in the mean value as identified by intervention

analysis. Black-filled areas above 0 indicate long-term above average conditions while below indicate below average conditions. Step changes were detected using a two-sample t-test on the previous and future 15 years for year t. Long-term averages were calculated as the mean value between intervention years. ……….... 26

Figure 2.5: Morlet wavelet power spectrums of the coastal (A) and continental (B)

SWEPC reconstructions. The y-axis represents Fourier periods while the x-axis represents

time. White enclosed areas represent a 95% confidence interval where significant wavelengths in the series are detected. Black lines represent power spectrum ridges and highlight nonstationary, significant frequencies across the time series. The lower contrast, or faded, areas on the left and right extremes of each figure represent areas outside of the analysis that is susceptible to zero padding effects. Note that A and B figures have different x-axes due to differing reconstruction lengths between the coastal and

continental models. ……….. 28

Figure 3.1: Study site map showing the location of the Capilano, Seymour and

Coquitlam watersheds. Squares with X’s indicate locations of hydrometric stations; circles represent Pacific silver fir tree-ring chronology locations; triangles represent mountain hemlock tree-ring chronology locations. .……… 41

Figure 3.2: (A) Climate Normal (1981-2010) for Grouse Mountain (Station ID:

1105658); monthly streamflow averages for the (B) Capilano (08GA010) and (C) Seymour (08GA013) watersheds; and, (D) combined July-August streamflow record of the Capilano and Seymour watersheds (Sources: ECCC, 2017; WSC, 2017). …...……. 44

(13)

xiii

Figure 3.3: Time plot of the measured (blue), reconstructed (thick black), and validation

(thin black) QJA records from the model calibration. ...……… 56 Figure 3.4: Time plot of reconstructed regional streamflow (z-scores; thin black line)

with a 10-year running mean (thick black line) and gauged data (blue line). Error

associated with the reconstruction is presented as the grey area surrounding the data and calculated using RMSE from the validation. The black histograms (top) represent 5th

-percentile QJA for the reconstruction. The grey-filled area below the reconstruction

represents changes in the long-term mean. Grey-filled (bottom) areas above 0 indicate long-term above average conditions while below indicate below average conditions. Step changes from one long-term mean to another are where significant intervention years were detected using a two-sample t-test between the previous and future 15-years at year

t. Intervention averages were calculated as the mean value between intervention years. 56

Figure 3.5: (A) Wavelet power spectrum of the proxy streamflow record. The vertical

axis shows Fourier periods while the horizontal axis are years. White-enclosed regions along the time series where significant frequencies are detected. The faded regions on the left and right of the image represent areas where sample depth is too low to describe low-frequency variability and susceptible to zero-padding effects. (B) Average wavelet power as represented by the average power value of rows, or Fourier periods, in (A). Blue and red points represent statistically significant wavelet power averages (p<0.10). ...…... 57

Figure 3.6: Scatterplot of consecutive year July-August below-average streamflow

defined as recurring discharge values below the full-reconstructed mean (0) over greater than one year. Values lower in the y-axis indicate higher severity lower-than-average streamflow during July and August while longer duration values along the x-axis illustrate prolonged intervals of sustained below-average flows. Grey circles represent values based on the reconstructed proxy record prior to the instrumental record. Grey circles with black borders represent consecutive-year below-average July-August

streamflow within the instrumental period (1929-present). …...……….. 62

Figure 3.7: Comparison of the (A) this study to streamflow reconstructions of the (B)

Skeena, (C) Atnarko, and (D) Chilko watersheds. Thick black lines indicate a moving 10-year average. ……… 66

Figure 4.1: Study site map showing the location of Harrison Lake. Circles with X’s

indicate tree-ring chronology locations; the square with X at Agassiz indicates the location of the climate station used in this study, the black circle with white outline at Harrison Mills indicates the location of the water level station on Harrison Lake; the triangles with a black circle indicates the location of the manual snow survey stations used in the study (Source: WSC, 2017). ……….. 74

Figure 4.2: (A) Climate normal (1981-2010) for Agassiz climate station (Station ID:

(14)

xiv

with squares indicate average monthly temperatures. (B) Harrison Lake average water levels (08MG012; 1933-2016) for each month of the year. .………... 75

Figure 4.3: Time plot of measured, modelled, and cross-validated records from the

model calibration. ………...………...………... 87

Figure 4.4: Time plot of reconstructed April Harrison Lake water level (thin black line)

with a 10-year running mean (thick black line) and gauged data (dashed line). Straight dashed line indicates the modelled mean lake level; solid black line indicates the measured mean value of April lake level. Error associated with the reconstruction is represented as the grey area surrounding the model. The thin horizontal dashed line represents the 300-year reconstructed mean April lake level while the thicker solid

horizontal line from 1933-2018 represents the gauged record mean. …….………. 87

Figure 5.1: Map of study area showing watershed boundaries. White circles with crosses

represent locations of interests (e.g., towns or mountains), triangles are tree-ring chronologies used in this study, and the outlined/faded regions are the three primary watersheds with salmon escapement information. ……….. 98

Figure 5.2: Climate Normal information for southwestern British Columbia (1981-2010)

for coastal (Grouse Mountain; Station ID: 1105658) and continental (Lillooet; Station ID: 1114619) regions. (Source: ECCC, 2017). ……….. 99

Figure 5.3: Full reconstructed record of (A) Chinook PC1 and (B) Coho PC1.The

thin-black line represents actual modeled values, thick-thin-black lines are a 5-year running mean, and grey areas are error calculated by the cross-validation RMSE. The grey-filled areas below each reconstruction show changes in the long-term mean. Grey-filled areas above 0 indicate long-term above average conditions while below indicate below average conditions. Step changes from one long-term mean to another are where significant intervention years were detected using a two-sample t-test on the previous and future 15 years for year t. long-term averages were calculated as the mean value between

intervention years. ………..……….... 113

Figure 5.4: Time plot of (A) Coho PC1 and (B) Chinook PC1 calculated (thin black),

reconstructed (thick black), and cross-validated (grey) for the calibration period (1951-1992). ………. 115

Figure 5.5: Spectral analysis of Coho (A, C) and Chinook (B, D) PC1 reconstructions:

Morlet wavelet power spectrum of the proxy Coho (A) and Chinook (B) records. Black enclosed areas represent a 95% confidence interval based using a white-noise background spectrum. The faded areas represent periods and frequencies of the analysis that is

susceptible to zero padding effects due to sample depth. Multitaper method (MTM) spectral analysis of Coho (C) and Chinook (D) proxy records. Red curves represent 90%, 95%, and 99% (bottom-top) significance. Significant (p ≤ 0.10) power exists at labelled frequencies. ………...………. 118

(15)

xv

Acknowledgements

This dissertation is the culmination of countless complaints and copious cups of coffee. I wish to thank those that have endured me since my arrival in Victoria. I have cherished my time at the University of Victoria and the people that made it great.

Dan, thanks for always being around to help or lend an ear to my insatiable ranting. You took a true Maritimer and turned him into a west coaster – for better or worse. You’ve spoiled me with trips to remote locations throughout the Coast Mountains and beyond. I will never forget spending 19 days in the backcountry eating ‘spaghetti wraps’ with ‘nachos for texture’ because most of our food was buried on day two...

Thank you, David and Tobi, for your help with getting this dissertation put together. Your thoughts and perspectives were appreciated as I ventured towards the completion of my research.

To the UVTRL, thanks for having me. I know that I was distracting, inappropriate, and – well – full of shit 90% of the time but it was a great few years. I’ve forged friendships with so many people in the lab and know they will continue. A special thanks to Bethany, Cedar, Jill, and Jodi. You may not know it but I look up to each of you as role models and mentors. Each of you is a badass. Lauren, Lee, Ale, and Tavy, thanks for the help in the field and company! You’re all gems!

To my futsal, seven-a-side, and intramural soccer teams, thanks so much for having me! Even with my constantly sprained ankle! I made a great group of friends while playing with you all. A special thanks to Aiden, ‘baps’, Ben, Britt, Byron, Evan, Kyle, and Shannon for being a great bunch of people and making me part of your amazing crew. Cody, Devin, and Henry – thanks for always listening or being generally happy to see me on my rare visits to Nova Scotia. You’re my best friends and I couldn’t ask for better. Although you really didn’t understand what I was going through, you always supported me.

Mom, Dad, Julie, and Haley, thanks for being a great supportive family! Never once did you ask those classic PhD questions ‘when are you done’ or ‘when are you getting a job.’ You trusted the process and me. I am forever grateful. Bodie, you can’t read this, but you’re cool too. Another thanks to Aris, you made sure I kept up with exercise and my mental health in check. Thanks, my ‘Ariba’. You’re adorable… most of the time.

Finally, thanks to my wonderful wife, Kaitlin. I’m here because I decided the best way to get to know you was studying. I would have flunked out of university if it weren’t for you. You challenge and bring the best out of me… Even when I hate it! You are the best. You are my rock. I love you.

(16)

xvi

Dedication

(17)

1

Chapter 1: Introduction

1.1 Hydrology of southwestern British Columbia

The hydrologic regime of southwestern British Columbia (BC) is influenced by its dynamic topography and climate oscillations originating in the Pacific Ocean

(Rodenhuis et al., 2009; Church and Ryder, 2010; Mote et al., 2018). Positioned adjacent to the Pacific Ocean and characterized by the Pacific Ranges of the BC Coast Mountains, moisture is delivered to the region primarily through orographic processes during the winter months (Church and Ryder, 2010). While much of the region is a temperate rainforest, recent summers were characterized by record-breaking, severe, drought conditions (River Forecast Centre, 2017). Attempts to understand the dynamics of these events in southwestern BC and to forecast the likelihood of future events are hindered by a relatively sparse network of instrumental records that rarely extend beyond 50 years (Coulthard et al., 2016). Limited records provide poor context for resource managers tasked with developing sound policy and management strategies for the region’s water supply.

Elsewhere in western North America, proxy reconstructions of critical

hydrological variables and their association to climate oscillations is an area of focused research (e.g., Meko et al., 2001; Pederson et al., 2011; Ballinger et al., 2018; Pathak et al., 2018; Welsh et al., 2019). The collective findings of this research illustrate that instrumental records underestimate worst-case scenario drought conditions and that current conditions are increasingly outside the range of observed variability (Belmechari et al., 2015; Coulthard et al., 2016). Of particular consequence, these studies reinforce the

(18)

2

strong regional hydroclimatic teleconnection to sea surface temperature variability in the Pacific Ocean (Rodenhuis et al., 2009; Ballinger et al., 2018; Pathak et al., 2018). While research has shown that related climate oscillations influence the hydroclimatic character of southwestern BC (Rodenhuis et al., 2009; Spry et al., 2014), directed research focused on examining and expanding upon the nature of these relationships over the long-term is sparse (e.g., Coulthard et al., 2016).

1.2 Hydroclimate relationships

Over the last 50 years, climate oscillations originating in the Pacific Ocean have been repeatedly linked to hydrological variability in southwestern BC (Bonsal and

Shabbar, 2008; Rodenhuis et al., 2009). While research shows large-scale teleconnections are responsible for upwards of 30% of the annual hydroclimate variability in snowpack accumulation, there is only limited appreciation for how these relationships influence regional hydrology. In western North America, recent research highlights the role of non-stationary relationships between hydrological variability and climate oscillations that have yet to be documented in southwestern BC (Marcinkowski et al., 2013; Coulthard et al., 2016; Mote et al., 2018; Litzow et al., 2018). What is appreciated is that three large-scale climate oscillations play a part in defining regional hydroclimate dynamics: (1) the El Niño-Southern Oscillation (ENSO); (2) the Pacific Decadal Oscillation (PDO), and; (3) the Pacific North American pattern (PNA) (McCabe et al., 2004; Rodenhuis et al., 2009; Abatzoglou, 2011; Starheim et al., 2013a; Coulthard et al., 2016; Mote et al., 2018).

(19)

3

ENSO is characterized by temperature and pressure fluctuations in the equatorial ocean surface off the western coast of South America that lead to ‘seesaw’ interactions with the southeastern tropical Pacific and Australian-Indonesian region (McCabe et al., 2004; Bridgman and Oliver, 2006;). During warm (positive; El Niño) ENSO events, the jet stream splits producing both high and low latitude storm tracks that rarely influence southwestern BC (Shabbar et al., 1997). The Pacific storm track migrates equatorward and downstream during El Niño events in response to local enhancement of the Hadley circulation off the coast of western North America (Chang et al., 2002). Winters that follow the onset of El Niño conditions during the previous spring, tend to be warmer and drier than normal (Stahl et al., 2006). In contrast, during La Niña (cool, negative ENSO) years winters tend to be cooler in response to a northward shift in the Pacific jet stream that brings wetter conditions to southwestern BC and the Pacific Northwest (Stahl et al., 2006).

The PDO is a long-lived ENSO-like pattern of Pacific climate variability defined by sea surface temperatures (SST) (Mantua et al., 1997; Bridgman and Oliver, 2006). The PDO describes a forced response from the North Pacific to atmospheric variability

associated with the Aleutian Low (AL) (Newman et al., 2003). While its behaviour is similar to ENSO, PDO tends to persist for decades rather than months with dominant periodicities at 15-25 and 50-70 years (Minobe, 1997; Mantua and Hare, 2002). Warm (positive) phases of PDO are associated with negative SST anomalies in the southwest Pacific that lead to positive winter temperature anomalies and higher-than-average precipitation totals in southwestern British Columbia and the Pacific Northwest. Cool PDO phases are associated with lower-than-average rainfall totals during the winter

(20)

4

months (Mantua et al., 1997; Chang et al., 2002; Bonsal and Shabbar, 2008; Abatzoglou, 2011).

The PNA pattern is an internal mode of atmospheric variability with alternating pressure patterns in the central Pacific Ocean and centres of action over western Canada and the southeastern United States (Latif and Barnett, 1994). It is characterized by out-of-phase atmospheric flows between the west coast of North America and eastern

Pacific/southeastern USA (Yu and Zwiers, 2007). Its development is associated with episodes of atmospheric blocking introduced by high-pressure ridges in the eastern Pacific Ocean (Wallace and Gutzler, 1981; Bridgeman and Oliver, 2006). Along with the North Atlantic Oscillation (NAO), the PNA plays a crucial role in setting up winter temperatures in Pacific North America by creating strong pressure gradients off the coast of British Columbia (Wallace and Gutzler, 1981). Positive phases of the PNA are

characterized by a strong AL, with southerly airflow along the west coast of North America and a high-pressure ridge over the Rocky Mountain cordillera (Wise et al., 2015). Negative phases of the PNA are associated with a weak AL with westerly zonal flows. While the Pacific storm track shifts southward during positive PNA phases, it shifts northward shifts during weaker phases (Yu and Zwiers, 2007).

There are clear relationships between climate oscillations and moisture delivery to southwestern British Columbia but our current understanding does not extend beyond the mid-20th century (Rodenhuis et al., 2009). While research highlights general changes between positive (warm) and negative (cool) phases of climate oscillations, most records are too short to understand how hydrological variables change over time and across phases. Longer-term understanding of spatial and temporal relationships between

(21)

5

dominant variables and hydroclimate are essential for the development of realistic policy and management strategies.

1.3 Research motivation

Present and future climate change poses a significant challenge for water resource management in western North America (Doll et al., 2015). Recent events in California offer an illustration of the economic, social, and ecological impacts of sustained droughts (Cook, 2019). Following five years of drought, reservoir storage in California decreased by upwards of 60% and many regions were left without running water (Cook, 2019). Hydroelectric power generation capacity was significantly reduced, necessitating an increased reliance on fossil fuels to supplement energy needs (Cook, 2019). Over 100-million trees died following drought-induced stress (Howitt et al., 2015), leading to a state-wide mortality event that continues to impact regional carbon and water cycles. While many ecosystem-services were also adversely affected, the related social impacts were psychological in nature and continue to be expressed by lingering health issues (Cochrane et al., 2004; Ding et al., 2011; Padhy et al., 2015; Kjellstrom et al., 2016).

While ongoing research focuses on past, present, and future trends of snowpack and streamflow in the western United States, there is only a limited understanding of the long-term hydroclimate state of southwestern BC. In a region where the water supply of many communities is reliant upon snowpack meltwater released during the spring and summer months, water shortages are certain if the winter snowpack is low or melts earlier in the year (Cayan et al., 2001; McCabe et al., 2004; Moore et al., 2010; Belmechari et al., 2015; Mote et al., 2018).

(22)

6

The downstream consequences of record low snowpack totals in 2014 and 2015 heightened water supply concerns in southwestern BC (River Forecast Centre, 2017). Uncertainty surrounds our understanding of how future climate change will impact the region’s water supply, especially with expected reductions in annual snowpack totals (Mote et al., 2018). Developing an informed understanding of how the regional hydroclimate has varied in the past will provide the insights water resource managers require to improve current strategies and policies.

The purpose of this dissertation was first to examine how southwestern BC’s hydroclimate varied prior to the instrumental record. Secondly, the research was intended to assess the long-term influence of low-frequency climate oscillations on that variability. The research objectives were to:

1. Develop multi-century reconstructions of snowpack and other hydroclimate-related variables in southwestern BC and determine their natural variability over the past few centuries;

2. Determine spatial and temporal relationships, synchrony, and variability of snowpack across southwestern BC, and;

3. Determine relationships between hydroclimate-related variables and climate oscillations generated from the Pacific Ocean and place our present understanding in a longer-term context.

1.4 Organization of dissertation

Following this chapter, Chapters 2, 3, 4, and 5 present the primary results of the dissertation. They are presented as individual manuscripts written and formatted for

(23)

7

refereed journal submission. Chapter 2 presents reconstructions of April 1 SWE over the past 300 years for continental and coastal regions experiencing different hydroclimate regimes in southwestern BC. Chapter 3 presents a regional 300-year dry-season (July-August) streamflow reconstruction for the Capilano and Seymour watersheds that provide 66% of the water supply for the 2.6-million residents of Greater Vancouver Regional District. Chapter 4 presents a long-term model of pre-instrumental April water level fluctuations in Harrison Lake over the past 300 years. Chapter 5 presents a regional-scale 300-year paleorecord of Chinook and Coho salmon escapement history within the Lower Fraser River watershed, and associates their associations to low-frequency Pacific Ocean teleconnections. The dissertation concludes with Chapter 6, where the main contributions and linkages between the different elements of the dissertation research are summarized, and ends with suggestions for future research.

(24)

8

Chapter 2: Snow Water Equivalent Variability Over Last 300 Years in

Southwestern British Columbia Driven by Teleconnections

2.1 Article Information

Chapter 2 consists of a manuscript prepared for submission. The text and figures from the manuscript have been renumbered and reformatted for consistency within the dissertation.

2.1.1 Author names and affiliations

Bryan J. Mood1* and Dan J. Smith1,

University of Victoria Tree-Ring Laboratory, Department of Geography, University of Victoria, British Columbia V8W 3R4, Canada.

*Corresponding Author Email: bjmood@uvic.ca

2.1.2 Author contributions

Mood developed the study, hypothesis, conducted and led field and laboratory work, statistical testing, wrote the manuscript, and produced all tables and figures. Smith provided funding for the research, provided guidance in formatting the study design, reviewed and edited the manuscript.

2.2 Abstract

Freshet contributions from melting seasonal snowpack are critical during hot, dry summer months for streamflow supply in southwestern British Columbia (BC), Canada. Recent below-average winter snow totals have generated cascading socioeconomic and ecological impacts that draw attention to the impending consequences of ongoing climate change. Within this region, knowledge of the year-to-year winter snowfall variability is

(25)

9

largely derived from a sparse network of short-duration (≤50 year) snow survey stations. In this paper, two long-term April 1 snow water equivalent (SWE) records were

developed from living tree-ring chronologies. The dendrohydrological models accurately reconstruct April 1 SWE for coastal and continental regions to AD 1710 and 1725. The models demonstrate annual April 1 SWE dynamics in this region are strongly associated with negative (cold) phases of the Pacific Decadal Oscillation and the El Niño-Southern Oscillation, and show that coastal SWE dynamics demonstrate more year-to-year variability than in more continental settings. Continental SWE is also shown to be strongly associated with positive phases of the Pacific North American Oscillation while coastal SWE is not. The coastal SWE reconstruction identifies eleven significant step changes and continental SWE record illustrates eight step chances since the early 1700s. Both reconstructions contain generally synchronous step changes in the scale of annual SWE magnitudes. Two periods, between 1804-1821 and 1887-1927, suggest that coastal regions are more sensitive to minor changes in SWE when compared to continental areas on the leeward side of the BC Coast Mountains. The two models provide the first high-resolution description of April 1 SWE over the past 300 years in southwestern BC and offer significant insights for water resource policy makers and planners.

2.2 Introduction

Seasonal snowpack is a critical input to regional hydrologic dynamics in many regions of western North America. Storing winter precipitation that is released as

meltwater during the spring and early summer freshet, these contributions provide runoff supplements at a time when the ecological, social, and industrial demands are highest

(26)

10

(Mote al., 2005, 2018; Rodenhuis et al., 2009). Snow is, however, vulnerable to changes in climate that influence both cover and depth (Abatzoglou, 2011; Mote et al., 2018). There is growing concern that the recent declines in April 1 SWE in the British Columbia (BC) Coast Mountains, Canada, are an outcome of shifting climates that are likely to persist (Rodenhuis et al. 2007; Déry et al. 2009, 2012; Mote et al., 2018).

In southwestern BC the water supply for many communities is dependent upon direct snowmelt contributions to streamflow and reservoir storage, as well as indirectly through contributions to groundwater recharge and throughflow (Eaton and Moore, 2010; Beaulieu et al., 2012; Olmstead, 2014; Koop et al., 2017). Record low April 1 SWE values in 2014 and 2015 contributed to summer-long water shortages in the Metro Vancouver area and emphasize the substantial water management challenges that

potentially lie ahead for many communities in southwestern BC. With further reductions in mountain snowpack certain to result in increasingly severe and frequent summer streamflow droughts, focused research is required to describe the full range of natural variability in April 1 SWE and to understand the underlying relationships to long-term climate oscillations generated from the Pacific Ocean.

The instrumental record of winter snowfall variability in southwestern BC is derived from a sparse network of short-duration (≤50 year) snow survey stations (Rodenhuis et al., 2009; Mote et al., 2018). As short instrumental records provide only limited insight into long-term April 1 SWE variability, attempting to forecast future trends from these data for water supply management purposes is problematic (Rodenhuis et al., 2009). Short records lead to challenges in understanding whether recent low April 1 SWE values are ‘extreme’ relative to what has occurred in the past and whether the

(27)

11

assumed climatic relationships are stable over time. Most watersheds in the region have limited storage capacity and variability in their annual runoff depends primarily on broad-scale climate oscillations (Yu and Zwiers, 2007; Bonsal and Shabbar, 2008; Abatzoglou, 2011). Together several distinct climate oscillations act to enhance or diminish snow delivery during the winter (Rodenhuis et al., 2009) and influence snow availability into the summer.

In this chapter, April 1 SWE records extending over several centuries are

reconstructed for southwestern BC from snow-sensitive Pacific silver fir (Abies amabilis Douglas ex J. Forbes) tree-ring chronologies. I use April 1 SWE records from coastal and continental regions to assess the influence of large-scale climate oscillations originating in the Pacific Ocean. The reconstructions provide the first high resolution descriptions of April 1 SWE over the past 300 years in southwestern BC, and are of immediate use to water resource managers charged with developing the strategies and policies required for adaptation to changing mountain snowpack.

2.3 Study Area and Research Background

The study area includes the Metro Vancouver area and the southern Pacific Ranges extending from Joffre Lakes to the Lower Fraser River within the southwestern BC Coast Mountains (Figure 2.1). Coastal settings in this region are moderated by

proximity to the Pacific Ocean and experience short, cool summers and long, wet winters. Average air temperatures above 1000 m above sea level (asl) remain below 0°C for 0 to 5 months of the year and above 10°C for 1 to 3 months of the year (Pojar et al., 1991; Kottek et al., 1996). The study area includes BC’s wettest ecological zone, where

(28)

12

Figure 2.1: Study area map showing the location of the tree-ring chronologies (white

triangles) and snow survey stations (continental = white square with centre black square; coastal = white circle with cross) incorporated into this research. Note that coastal snow survey stations are located in regions with high orographic precipitation (western areas), while continental snow survey stations are in the Coast Mountains rain shadow (east).

(29)

13

precipitation totals range between 1700 to 5000 mm/yr and 70% falls as snow and/or rain during the winter months (Pojar et al., 1991; Spry et al., 2014) (Figure 2.2A). Most precipitation falls as a result of orographic interactions on the windward slopes of the Pacific Ranges. Precipitation totals decrease substantially on the eastern slopes, where persistent rain shadow conditions prevail (Pojar et al., 1991; Kottek et al., 2006; Church and Ryder, 2010; Moore et al., 2010) (Figures 2.2A and 2.2B).

The impact of decreased winter snow delivery and storage in southwestern BC is demonstrated by below-average April 1 SWE totals in recent years. In 2015, the 2014-2015 winter snowpack in the region was 49% of normal in January and reached a record low of 0% by June 1st (River Forecast Centre, 2017). Low April 1 SWE resulted in a severe summer streamflow drought and led to depleted water supplies necessitating use restrictions for the 2.6-million people residing in Metro Vancouver. Highlighting the importance of mountain snowpack for sustaining summer streamflow, the depleted runoff had cascading socioeconomic and ecological consequences that drew attention to the impending consequences of ongoing climatic changes in this region (Cook et al., 2018, 2019; Lledo et al., 2018).

Previous research within the study area identified a causal relationship between inter-annual/-decadal climate variability generated by atmosphere-ocean interactions originating in the Pacific Ocean and April 1 SWE variability (Sellars et al., 2008; Spry et al., 2014). These interactions include those characterized by the Pacific Decadal

Oscillation (PDO), the El Niño-Southern Oscillation (ENSO), and the Pacific North American (PNA) pattern (Rodenhuis et al., 2009; Abatzoglou, 2011; Mote et al., 2018). Warm/cool phase relationships between April 1 SWE and ENSO show that an average

(30)

14

Figure 2.2: General climate information for the study region. (A) Climate normals

(1981-2010) for coastal (Grouse Mountain, station ID: 1105658) and continental

(Lillooet, station ID: 1114627) regions (ECCC, 2017) stations; (B) time series of average coastal and continental April 1 SWE records (River Forecast Centre, 2017).

(31)

15

variability of greater than 30% characterized the interval between 1951 and 2007 (Rodenhuis et al., 2009).

Relationships between April 1 SWE and PDO variability are less pronounced over the instrumental period, with only a 15% change in snowpack noted between warm/cool phases (Rodenhuis et al., 2009). The PDO is a large-scale climate system that influences the surface climate and hydrology of western North America (Whitfield et al., 2010). Phase changes in the PDO have a significant influence on the overall hydroclimate of southwestern BC (Moore and McKendry, 1996; Rodenhuis et al., 2009). The PDO is typically coupled with ENSO when describing changes in temperature and precipitation regimes associated with their variability (Whitfield et al., 2010). Temperature

fluctuations are readily apparent between warm and cool phases of the PDO (Kiffney et al., 2002; Whitfield et al., 2010), and the PDO shift in 1976-1977 resulted in an abrupt change in average winter temperatures (Hartman and Wendler, 2005) throughout BC (Stahl et al., 2006; Fleming and Whitefield, 2010). In the BC South Coast region,

precipitation is greater during cool phases of PDO and less during warm phases (Kiffney et al., 2002; Stahl et al., 2006). In southwestern BC, April 1 SWE is also greatly

influenced by PDO (Moore and McKendry, 1996; McCabe and Dettinger, 2002). During cold phases of PDO the greatest seasonal SWE values are characteristically recorded on April 1st, while during warm phases lower than average SWE values are common (Moore and McKendry, 1996).

No statistical relationships between April 1 SWE and PNA are confirmed in the BC Coast Mountains, although it has been associated with the seasonal snow line elevation (Abatzoglou, 2011). The PNA does affect the overall hydroclimate of

(32)

16

southwestern BC through changes in winter storm intensity and frequency (Rodenhuis et al., 2009). Positive PNA results in anomalously high freezing levels and reduced snow coverage (Abatzoglou, 2011). April 1 SWE values are lower in continental regions when compared to coastal settings. I associate the different values with orographic precipitation patterns between coastal (windward; higher overall precipitation) and continental

(leeward; rain shadow region) regions. This association suggests that continental locations experience more significant influences from the PNA, as snowfalls in these locations are more reliant on freezing level relationships than are coastal settings.

2.4 Methods

Pacific silver fir (PSF) is a conifer tree common to maritime climate regions in Pacific North America (Crawford and Oliver, 1990). Normally found at sites with deep, well-drained soils, stands of PSF are often found in association with other high-elevation tree species including western hemlock (Tsuga heterophylla (Raf.) Sarg.), mountain hemlock (Tsuga mertensiana (Bong.) Carr.), yellow-cedar (Chamaecyparis nootkatensis D. Don), and western red cedar (Thuja plicata Donn ex D. Don) (Pojar et al., 1991).

Snowpacks that persist into the summer months truncate the growing season of many tree species located close to the altitudinal treeline (Ettinger et al., 2011). At high elevation sites, PSF tends to be found where precipitation totals range from 750 to over 6500 mm/yr (Crawford and Oliver, 1990). Near the upper limit of its altitudinal range, the radial growth of PSF is dependent upon soil moisture derived from melting seasonal snow (Crawford and Oliver, 1990), and is also known to be limited by late-lying

(33)

17

trees shares similar environmental requirements and, due to their strong, negative relationship to April 1 SWE, were previously used to reconstruct SWE dynamics in Pacific North America (e.g., Peterson and Peterson, 2001; Pederson et al., 2011;

Coulthard et al., 2016; Welsh et al., 2019). By comparison, PSF has seen limited use as an environmental proxy, although it demonstrates strong statistical relationships to April 1 SWE at high elevation in the area (Ettinger et al., 2011).

2.4.1 Tree-ring data

PSF tree ring samples were collected with 5.2-mm increment borers at high-elevation sites in southwestern BC in the summer of 2016 using standard

dendrochronological techniques (Fritts, 1976; Speer, 2010). After air drying and

processing, the annual ring widths were measured at the University of Victoria Tree-Ring Laboratory (UVTRL) using a WinDendro™ digital measurement system (v. 2016a; Guay et al., 1992). Supplementary PSF tree-ring chronologies were downloaded from the International Tree-Ring Data Bank (ITRDB; NOAA, 2017). Site-specific tree-ring chronologies were constructed by converting the ring width measurements to

standardised indices using the R package dplR to remove growth related trends (Bunn, 2008) with a 50-year cubic smoothing spline. The site-specific chronologies were subsequently combined with the supplementary chronologies using the bi-weight robust mean method to create a regional PSF chronology (Briffa and Melvin, 2011).

2.4.2 Snow water equivalent data

April 1 SWE data from manual snow survey stations were accessed using the online portal managed by the BC River Forecast Centre (River Forecast Centre, 2017). April 1 SWE data exceeding 45 years in length were downloaded and quality checked.

(34)

18

Where missing data was identified, the long-term mean was used to replace the null value. Data were separated into coastal and continental SWE based geographic location and linear associations between snow survey sites. Individual coastal and continental SWE sites were combined using a principal component analysis (PCA), with only principal components loadings >10% used for further analysis, to extract the underlying regional variance (Table 2.1). Coastal SWEPC only used the first component for analysis,

and successfully explained 89% of the variance in the time series. Continental SWEPC

was similar, using the first component of the analysis which explained 69% of the underlying variance (Table 2.1). The data could not be combined into a larger PCA because the two groups were not linearly associated.

2.4.3 Climate oscillation data

Proxy PDO and PNA records were downloaded from the NOAA Paleoclimate website (NOAA, 2017). The datasets were developed from tree-ring records that have repeatedly been linked to winter hydroclimate variability in western North America (D’Arrigo and Wilson, 2006; Trouet and Taylor, 2009; Asong et al., 2018). The PDO reconstruction provides an estimate of the Asian expression of the phenomenon, while the PNA reconstruction provides a winter period record (D’Arrigo and Wilson, 2006; Trouet and Taylor, 2009). While the Asian expression of the PDO is similar to observed trends in Pacific North America, it has different ENSO-related modulating

characteristics.

2.4.4 Model estimation

Coastal and continental April 1 SWE records (coastal and continental SWEPC)

(35)

19

Table 2.1: Descriptive information of manual snow survey stations used in this study.

Records were collected from the River Forecast Centre (2017) website. Two regions were targeted for reconstruction: (1) continental and (2) coastal. Continental locations are located east of the Pacific Ranges and characterized by lower average snowpack. Coastal stations are located in the Pacific Ranges and characterized by much greater snowpack averages (229 vs. 1367 SWE). Station ID is associated with River Forecast Centre designations; coordinates (latitude and longitude) are in decimal degrees; elevation is rounded to the nearest 10 m above sea level; mean April 1 SWE is calculated across the whole time series; span is the length of continuous April 1 SWE measurements from each station; length is the total number of years available for analysis; PCA loading represents the explanatory power of the first principal component from analysis.

Station ID Latitude (DD) Longitude (DD) Elevation (m asl) Mean April 1 SWE Span (years AD) Length

(years) PCA Loading

Continental 1C05 50.70 -122.62 1720 609 1952-2014 62 1C06 50.91 -121.82 1214 44 1955-2014 59 1C09A 50.50 -120.98 1550 98 1958-2014 56 1C14 50.78 -122.79 1400 167 1963-2014 51 1C19 50.46 -121.05 1640 123 1967-2014 47 2G06 49.50 -120.80 1490 334 1960-2014 54 Continental SWEPC 89% Coastal 1D08 49.58 -122.32 1195 1529 1968-2014 46 1D10 49.83 -122.06 1555 1360 1968-2014 46 3A01 49.38 -122.08 1130 1231 1936-2014 78 3A09 49.46 -123.03 880 1479 1946-2014 68 3A10 49.37 -122.96 1010 1235 1945-2014 69 Coastal SWEPC 68%

(36)

20

determine initial relationships between the variables. A pool of predictors was developed for the both coastal and continental April 1 SWE records and input through a multiple linear regression to provide a best-fit model. Model selection was based on adjusted R2, multicollinearity (variance inflation factor or VIF), autocorrelation at lag 1 (Durbin-Watson or DW), and cross-validation statistics (Reduction of Error or RE).

2.4.5 Analysis of the reconstruction

The coastal and continental SWEPC models selected for reconstruction were

compared to measured winter ENSO, PDO, and PNA climate oscillations (October of the previous year to March) to reveal low-frequency influences on snowpack over the

reconstructed and instrumental periods. A difference-of-correlations test was conducted on both the measured coastal and continental SWEPC records, as well as on the

reconstructed values, against positive and negative phases of the selected climate

oscillations. Where only one phase illustrates a statistical relationship with the models, it indicates non-stationary (alterations in significance over time) responses to the specific oscillation.

A Morlet wavelet analysis using the R-package WaveletComp (Roesch and Schmidbauer, 2015) was used to illustrate fluctuations in power over the length of the full reconstruction period and to determine the average variability. An intervention analysis was conducted to determine significant long-term changes in reconstructed averages. The intervention analysis used a 30-year moving window to highlight where significant differences occurred between the preceding and proceeding 15 years of year t using a two-sample t-test. Periods of change in the long-term mean were calculated as departures from the average SWEPC value over the full reconstruction.

(37)

21

2.5 Results

Five PSF chronologies were selected for use as coastal and continental SWEPC

proxies (Table 2.2). All the chronologies were significantly correlated to each other and were combined into a single index using regional curve standardisation (PSFR).

Consolidating the site-level tree-ring records enhanced the explanatory power and possible reconstruction period to the interval from 1710-2015. The indices had RBARs ranging from 0.350 to 0.491 and were significantly correlated to coastal and continental SWEPC (Table 2.3). The selected climate oscillation reconstructions (PDO and PNA)

span from 1565-1988 and 1725-1999, respectively (D’Arrigo and Wilson, 2006; Trouet and Taylor, 2009). The PDO reconstruction was significantly correlated with coastal SWEPC while both PDO and PNA reconstructions were related to continental SWEPC

(Table 2.3).

Two models were developed from the coastal and continental SWEPC records.

The first principal component for coastal and continental SWEPC were identified for

reconstruction and used PSFR, PDO, and PNA as predictors:

Coastal SWEPC = -1.75*PSFR + -0.34*PDO + 1.67

Continental SWEPC = -2.41*PSFR + -0.59*PDO + 0.22*PNA + 2.56

where coastal SWEPC is the first component April 1 SWE from coastal manual snow

survey stations; PSFR is the regional Pacific silver fir tree-ring width chronology; PDO is

a Pacific Decadal Oscillation reconstruction from D’Arrigo and Wilson (2006); continental SWEPC is the first component April 1 SWE from continental manual snow

(38)

22

survey stations; and PNA is a Pacific-North American pattern tree-ring reconstruction from Trouet and Taylor (2009).

The Coastal SWEPC reconstruction extends from 1710 to 1988 and explains

50.7% of variance, accounting for lost degrees of freedom, in the first principal component. The model has a DW statistic of 1.86, VIF of 1.62, and F-Ratio of 13.3 (Table 2.4). It was successfully cross-validated using the LOO method (RE = 0.41) indicating the model is robust. Regression and cross-validation statistics are shown in Table 2.4 and a time plot of the actual, modelled, and cross-validated coastal SWEPC is

shown in Figure 2.3A. The full reconstruction is illustrated in Figure 2.4A.

The continental SWEPC reconstruction extends from 1725 to 1988 and explains

73.1% of variance, accounting for lost degrees of freedom, in the first principal

component. The continental SWEPC model has a DW statistic of 2.27, VIF of 1.095, and

F-Ratio of 23.63 (Table 2.4). The model was successfully cross-validated using the LOO method (RE = 0.57) suggesting that the model is robust. Regression and cross-validation statistics are shown in Table 2.4 and a time plot of the actual, modelled, and cross-validated continental SWEPC are shown in Figure 2.3B. The full reconstruction is shown

in Figure 2.4B.

Both reconstruction models and calculated SWEPC were tested against large-scale

climate teleconnections originating in the Pacific Ocean. A difference-of-correlations test between ENSO, PDO, and PNA reveal significant relationships between Continental SWEPC and negative ENSO, positive PDO, and positive PNA while Coastal SWEPC was

(39)

23

Table 2.2: Times series information. Species/type are bold and italics represents time

series used as predictors in the reconstruction. Tree-ring site numbers in brackets are ITRDB codes. The Pacific silver rir regional index (PSFR) was developed by combining

the five site-level chronologies used a bi-weight robust mean method (Briffa and Melvin, 2011). Climate oscillation indices were collected from NOAA (2017) for teleconnections that are known to influence the overall hydroclimate of western North America. RBAR is the average value across whole index; length is the span of years used for reconstruction. For tree-ring series, the length used for reconstruction is used only for where EPS > 0.85.

Study Site (Site No.) Source

Lat (DD)

Lon (DD)

Elev.

(m asl) RBAR Length

Pacific Silver Fir

Deek’s Lake This study 49.52 -123.21 1050 0.350 1696-2015 Callaghan Lake This study 50.18 -123.13 975 0.360 1696-2015 Seymour (CANA107)

Dobry et al.

(1996) 49.52 -123.07 1000 0.491 1686-1992 Hurricane Ridge (WA081)

Schweingruber

et al. (1991) 46.15 -122.15 1200 0.372 1698-1983 Mt. St. Helens (WA082)

Schweingruber

et al. (1991) 47.98 -123.47 1200 0.437 1648-1980

Pacific Silver Fir Regional (PSFR) 0.540 1648-2015

Teleconnections

Pacific Decadal Oscillation

(PDODW) D'Arrigo and Wilson (2006) 1565-1988

Pacific-North American

(40)

24

Table 2.3: Pearson correlations between model parameters used in the reconstruction. All

values shown are significant (p<0.01).

Time Series Coastal SWEPC Continental SWEPC Pacific Silver Fir Regional (PSFR) -0.67 -0.68 Pacific Decadal Oscillation (PDODW) -0.67 -0.78 Pacific North American (PNATT) -0.66

Durbin-Watson Statistic; VIF = variance inflation factor; SE = standard error; RE = reduction of error; RMSE = root mean squared error; CV = coefficient of variance.

Reconstruction R2 _{Adj. R}2 _D-W _VIF _SE _F-Ratio Coastal 0.55 0.51 1.87 1.57 0.70 13.30 Continental 0.76 0.73 2.27 1.10 0.59 23.60

Cross Validation RE RMSE Coastal 0.41 0.66 Continental 0.57 0.54

(41)

25

Figure 2.3: Time plot of calculated April 1 SWE principal components (thin black line),

modelled (thick black line), and cross-validated (thick grey line) records of (A) coastal and (B) continental.

Dendrohydrological reconstruction and hydroclimatic variability in southwestern British Columbia, Canada