FORECASTING SPRING FRESHET EVENTS IN THE KISKATINAW RIVER BASIN,
BRITISH COLUMBIA
by
Hunter E. Gleason
B.Sc., Colorado State University, 2014

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENT FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES AND ENVIRONMENTAL STUDIES

UNIVERSITY OF NORTHERN BRITISH COLUMBIA
April 2018
© Hunter Gleason, 2018

FORECASTING SPRING RUNOFF IN THE KISKATINAW RIVER
Abstract
Incorporating climate information into hydrologic streamflow forecasts has allowed for
significant advancement in the ability to predict seasonal streamflow. The City of Dawson
Creek (CDC), BC, has depended on the Kiskatinaw River (KR) as its sole source of
municipal water for over 60 years. Hydro-meteorological changes in the KR along with
increasing population and growing industry have put stress on the CDC water supply. In this
study regional surface climate observations aggregated over the winter accumulation period
(15 November–25 March) integrated with global circulation indices were input into a series
of regression models providing spring runoff predictions in the KR. The surface climate
observations, indices of global circulation and snow cover provided good predictability of
both cumulative streamflow timing and volume in the KR. This study provides the CDC with
a tool for better informed releases and withdrawals from the KR during the spring freshet.

ii

iii
Table of Contents
Abstract ..................................................................................................................................... ii
Table of Contents ..................................................................................................................... iii
List of Figures ...........................................................................................................................vi
List of Tables ......................................................................................................................... xiii
Acknowledgments ..................................................................................................................xvi
Chapter 1:

Introduction.......................................................................................................... 1

1.1

Overview ..................................................................................................................... 1

1.2

Uncertainty in Hydrology ........................................................................................... 4

1.3

Hydrologic Response to Global Climate Change ....................................................... 5

1.4

City of Dawson Creek and the Kiskatinaw Watershed ............................................... 7

1.5

Hydro-Meteorological Changes in the Kiskatinaw Watershed................................... 8

1.6

Turbidity ..................................................................................................................... 9

1.7

Research Objectives .................................................................................................. 11

Chapter 2:

Study Area ......................................................................................................... 12

2.1

General Description .................................................................................................. 12

2.2

Climate ...................................................................................................................... 15

2.3

Hydrology ................................................................................................................. 16

2.4

Water Consumption .................................................................................................. 22

Chapter 3:

Data Sources and Assimilation .......................................................................... 24

iv
3.1

Synoptic Weather Patterns and Surface Climate Variability in British

Columbia ............................................................................................................................. 24
3.2

Data Sources ............................................................................................................. 24

Environment and Climate Change Canada historical climate data. ............................... 25
Manual snow survey observations. ................................................................................ 28
Satellite-derived snow covered fraction. ........................................................................ 29
Discharge data. ............................................................................................................... 34
Chapter 4:

Methods ............................................................................................................. 35

4.1

Timing of Cumulative Discharge and Total Runoff Volume ................................... 35

4.2

Spatial Correlations ................................................................................................... 35

4.3

Predicting Runoff Timing and Volume .................................................................... 37

Resampling. .................................................................................................................... 37
Model Selection.............................................................................................................. 38
Forecast Validation ............................................................................................................. 39
4.4

Predicting Turbidity Exceedance .............................................................................. 40

Chapter 5:

Results ............................................................................................................... 41

5.1

Average Timing and Volume of Cumulative Streamflow ........................................ 41

5.2

Streamflow Timing and Snow Ablation ................................................................... 41

5.3

Correlation with Global Circulation Indices ............................................................. 44

5.4

Spatial Correlations ................................................................................................... 44

v
5.4.1

Timing of spring runoff. .................................................................................... 45

5.4.2

Total runoff volume. .......................................................................................... 47

5.5

Runoff Timing .......................................................................................................... 48

Early runoff. ................................................................................................................... 49
Mid-season runoff. ......................................................................................................... 49
Late runoff. ..................................................................................................................... 51
5.6

Runoff Volume ......................................................................................................... 51

5.7

Model Validation ...................................................................................................... 53

5.8

Predicting Turbidity Exceedance .............................................................................. 56

Chapter 6:

Discussion .......................................................................................................... 59

6.1

Teleconnections ........................................................................................................ 59

6.2

Predictive Modelling vs. Causal Attribution ............................................................ 62

6.3

Decision Support Systems ........................................................................................ 63

6.4

Future Research and Recommendations ................................................................... 64

Chapter 7:

Conclusions ....................................................................................................... 66

References................................................................................................................................ 68
Appendix A.............................................................................................................................. 77
Appendix B. ............................................................................................................................. 79
Appendix C. ............................................................................................................................. 96
Appendix D.............................................................................................................................. 97

vi
List of Figures
Figure 1.1. Contour plot of monthly change in snow water equivalent (SWE) and
runoff as a function of elevation for the Peace River Basin and the Upper
Columbia River, British Columbia based on the Intergovernmental Panel on
Climate Change A1B scenario ensemble median. Figure reproduced from
Schnorbus, Werner, & Bennett (2014), used with permission from the American
Meteorological Society. ................................................................................................... 6
Figure 1.2. Average daily discharge at Water Survey Canada (WSC) 07FD001
Farmington stream gauge for two 26 year periods 1966–1991 and 1991–2016.
Peak spring flow during the latter period is reduced by 14% and occurs 13 days
earlier. ............................................................................................................................. 10
Figure 2.1. Digital elevation map (DEM) of the Kiskatinaw River Basin (KRB)
delineated at the Arras weir (circle) where the City of Dawson Creek (CDC)
withdraws water and delineated at the Water Survey of Canada (WSC) 07FD001
Farmington stream gauge (triangle). .............................................................................. 13
Figure 2.2. Relative height versus relative area in the Kiskatinaw River Basin defined
at Arras where a is the area of basin above height h, A is the total basin area, h is
the height above the outlet, and H is the total relief of the basin. .................................. 14
Figure 2.3. Map depicting the major and sub-basins of the Kiskatinaw River Basin
delineated at the Arras weir where the City of Dawson Creek withdraws water for
distribution and treatment............................................................................................... 17
Figure 2.4. Minimum (Min Q), mean (Mean Q) and maximum (Max Q) discharge at
Water Survey Canada (WSC) 07FD001 Farmington stream gauge, 1966–2016.
The black dotted line is the mean snow covered fraction of the Kiskatinaw River
Basin (KRB) derived from satellite imagery over 2002-2015. A significant
increase in streamflow during spring at Farmington corresponds with onset of
snowmelt in the KRB. .................................................................................................... 19
Figure 2.5. Hydrograph for Water Survey Canada hydrometric gauge at Farmington in
2016. During this year a peak discharge in response to snowmelt is seen
occurring March–May; however, a peak of greater magnitude is also observed in

vii
response to an intense cyclone that brought copious amounts of rain to the Peace
River Region, resulting in severe flooding over 10–16 June.∙ ....................................... 20
Figure 2.6. Analysis of wet (above zero) and dry (below zero) years based on mean
Water Survey Canada Farmington discharge, 1965–2016. Mean discharge for the
period is 10.7 m3 ·s-1. ...................................................................................................... 21
Figure 2.7. Amount of currently licensed water allocations from the Kiskatinaw River
by purpose represented as the percent of the total annual allocated volume
(24,200,000 m3) as of 15 April 2018. Four categories dominate water usage, the
remaining four make up less than 0.5% of the total volume. Data were obtained
from the online Government of British Columbia Water License Query tool
(Government of British Columbia, 2018). ..................................................................... 23
Figure 3.1. Map displaying the location of active automatic weather stations from both
Environment and Climate Change Canada and the BC Ministry of Transportation
and Infrastructure. Also included are the locations of manual surveys of snow
water equivalent obtained from the BC Ministry of Forests, Lands, Natural
Resource Operations, and Rural Development. The black circle denotes the 1200km synoptic radius encompassing the Kiskatinaw River Basin. .................................... 26
Figure 4.1: Diagram depicting the model structure implemented for making predictions
of spring runoff in the Kiskatinaw River. Climate observations obtained at
weather stations and manual snow survey locations were aggregated over the
winter accumulation period (15 November-25 March). Cumulative discharge at
the Farmington stream gauge (maroon curve) was determined for each runoff
period (25 March-1 July) from the raw discharge data (dashed line). The day of
year that discharge reached specified percentages of total runoff volume was
determined in 5% intervals. For instance, the R50% above (centre vertical red
line) would be the day of year that 50% of the total runoff volume was reached
for that year (2003 in this case). The R25%-R75% depicts the interquartile range.
The aggregated winter climate observations were input into a regression analysis
where candidate models were identified for each RX% and the total runoff
volume. The selected models were then used to predict the cumulative runoff for
that runoff period in discrete 5% intervals. .................................................................... 36

viii
Figure 5.1. Distribution of annual spring runoff timing and volume (25 March–1 July)
in the Kiskatinaw River since hydrologic year 1965 (calendar year 1966).The left
axis identifies the total runoff volume of water in cubic meters, whereas the right
axis shows the timing of cumulative runoff as a Julian day of year for the R5%,
R25%, R50%, R75% and R95% determined over the runoff period 25 March-1
July. ................................................................................................................................ 42
Figure 5.2: The observed (black solid line) and predicted (dashed maroon line)
cumulative discharge at Water Survey of Canada Farmington for hydrologic year
2002, over the runoff period 25 March (hydrologic day of year 207) to 1 (HDOY
304) July. The original predicted cumulative values are represented by the teal
X’s. A 10-day moving average is applied to the predictions. The various dots
represent the interquartile range of each model’s residuals, added to the
prediction, in different combinations of timing t, and volume v. For example, the
P25t is the predicted timing plus the 25th percentile of that model’s residuals.
The P75v is the predicted volume plus the 75th percentile of the residuals for the
total runoff volume model. ............................................................................................. 54
Figure 5.3. Predicted hydrograph at the Water Survey of Canada Farmington stream
gauge in the Kiskatinaw River for hydrologic year 2002. The blue line is the
predicted discharge and the black line is the observed discharge (linearly
interpolated from 5% cumulative intervals). See Figure 5.1 for legend
description. The Nash-Sutcliffe Efficiency Coefficient was determined to be 0.70.
The prediction is made on 1 April, and a 10-day moving average was applied to
the prediction (Figure 5.2). The abrupt changes in discharge are the result of
negative slopes in the predicted cumulative discharge, and can be filtered out............. 55
Figure 5.4. Daily average turbidity measured at the Arras weir as a function of daily
average discharge at Water Survey of Canada Farmington stream gauge since
2007 during spring runoff periods, 25 March–1 July. The horizontal blue line
marks the 500 nephelometric turbidity units (NTU) threshold value at which
water abstraction must cease. ......................................................................................... 57
Figure B.1. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume

ix
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed average
winter (15 November–25 March) temperature at Environment and Climate
Change Canada and Ministry of Transportation and Infrastructure climate
stations............................................................................................................................ 79
Figure B.2. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed average
winter (15 November–25 March) temperature at Environment and Climate
Change Canada and Ministry of Transportation and Infrastructure climate
stations............................................................................................................................ 80
Figure B.3. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed average
winter (15 November–25 March) temperature at Environment and Climate
Change Canada and Ministry of Transportation and Infrastructure climate
stations............................................................................................................................ 81
Figure B.4. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 95% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed average
winter (15 November–25 March) temperature at Environment and Climate
Change Canada and Ministry of Transportation and Infrastructure climate
stations............................................................................................................................ 82
Figure B.5. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed total winter
(15 November- 25 March) precipitation at Environment and Climate Change
Canada and Ministry of Transportation and Infrastructure climate stations. ................. 83
Figure B.6. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed total winter

x
(15 November–25 March) precipitation at Environment and Climate Change
Canada and Ministry of Transportation and Infrastructure climate stations. ................. 84
Figure B.7. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed total winter
(15 November–25 March) precipitation at Environment and Climate Change
Canada and Ministry of Transportation and Infrastructure climate stations. ................. 85
Figure B.8. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and surveyed 1 April
snow water equivalent from the Ministry of Forests, Lands, Natural Resource
Operations & Rural Development. ................................................................................. 86
Figure B.9. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and surveyed 1 April
snow water equivalent from the Ministry of Forests, Lands, Natural Resource
Operations & Rural Development. ................................................................................. 87
Figure B.10. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and surveyed 1 April
snow water equivalent from the Ministry of Forests, Lands, Natural Resource
Operations & Rural Development. ................................................................................. 88
Figure B.11. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed average
winter (15 November–25 March) atmospheric surface pressure at Environment
and Climate Change Canada and Ministry of Transportation and Infrastructure
climate stations. .............................................................................................................. 89
Figure B.12. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed average

xi
winter (15 November–25 March) atmospheric surface pressure at Environment
and Climate Change Canada and Ministry of Transportation and Infrastructure
climate stations. .............................................................................................................. 90
Figure B.13. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume
(25 March–1 July) occurs in the Kiskatinaw River, BC, and observed average
winter (15 November–25 March) atmospheric surface pressure at Environment
and Climate Change Canada and Ministry of Transportation and Infrastructure
climate stations. .............................................................................................................. 91
Figure B.14. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the
Kiskatinaw River, BC, and observed average winter (15 November–25 March)
temperature at Environment and Climate Change Canada and Ministry of
Transportation and Infrastructure climate stations. ........................................................ 92
Figure B.15. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the
Kiskatinaw River, BC, and observed total winter (15 November–25 March)
precipitation at Environment and Climate Change Canada and Ministry of
Transportation and Infrastructure climate stations. ........................................................ 93
Figure B.16. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the
Kiskatinaw River, BC, and surveyed 1 April snow water equivalent from the
Ministry of Forests, Lands, Natural Resource Operations & Rural Development. ....... 94
Figure B.17. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the
Kiskatinaw River, BC, and observed average winter (15 November-25 March)
atmospheric surface pressure at Ministry of Transportation and Infrastructure
climate stations. .............................................................................................................. 95
Figure C.1. Winter temperature (°C) departure from the 1951–2000 normal over
Canada in response to a positive phase of El Niño Southern Oscillation. Figure

xii
reproduced from Shabbar (2006), used with permission under the Creative
Commons Attribution 4.0 License. ................................................................................ 96
Figure C.2. Winter precipitation (mm d-1) departure from the 1951–2000 normal over
Canada in response to a positive phase El Niño Southern Oscillation. Figure
reproduced from Shabbar (2006), used with permission under the Creative
Commons Attribution 4.0 License. ................................................................................ 96

xiii
List of Tables
Table 1.1: Typical percentage of streamflow inputs (snowmelt and rain) known in
snowmelt dominated watersheds (northern hemisphere) with increasing time.
Table recreated from USDA (1972). ................................................................................ 4
Table 2.1: Percent of total area of the Kiskatinaw River Basin defined at Arras for
existing land cover/land use categories in the watershed as of 2010. Recreated
from Paul (2013). ........................................................................................................... 14
Table 2.2: Average daily minimum (Tmin) and maximum (Tmax) temperature and the
total seasonal precipitation (Total Ppt) over the Kiskatinaw River Basin defined
at the Arras weir (Figure 2.1) for 1950-2010. Values were calculated by Aseem
Sharma using the R statistical programming language and climate data provided
by the Australian National University (2017). ............................................................... 15
Table 3.1: Summary of datasets utilized for predicting spring runoff in the Kiskatinaw
River along with host organization and observation details. ......................................... 27
Table 5.1: Statistics describing the observed timing of cumulative runoff (RX%) and
total runoff volume (TRV) in the Kiskatinaw River, BC, defined at the Water
Survey of Canada Farmington stream gauge, 1965–2016. Timing is the
hydrologic day of year (HDOY, days since 1 September) that a given RX% is
reached, and TRV is in cubic meters.............................................................................. 43
Table 5.2: Table showing slope coefficients (β) and the jackknifed unbiased estimates
of the slope coefficients (βj) in estimating the hydrologic day of year that
cumulative discharge reaches 25%, 50% and 75% of total runoff volume in the
Kiskatinaw River. Bias is provided, and also shown as an absolute percentage of
the biased coefficient. MT stands for mean temperature, CPPT stands for
cumulative precipitation, and numbers below reference agency stations number.
Agency MOTI is Ministry of Transportation and Infrastructure, EC is
Environment and Climate Change Canada. ................................................................... 50
Table 5.5: Table showing slope coefficients (β) and the jackknifed unbiased estimates
of the slope coefficients (βj) in estimating the total runoff volume (ha m) in the
Kiskatinaw River. Bias is provided, and also shown as an absolute percentage of
the original coefficient. MT stands for mean temperature, CPPT stands for

xiv
cumulative precipitation, and numbers below reference agency stations number.
Agency MOTI is Ministry of Transportation and Infrastructure, EC is
Environment and Climate Change Canada. ................................................................... 52
Table 5.6: Logistic regression results predicting the probability of exceeding 500 NTU
at Arras from daily average discharge at the Water Survey of Canada Farmington
stream gauge, on the Kiskatinaw River. Model based on data since 2007 over the
spring runoff period 25 March-1 July. ........................................................................... 58
Table 5.7: Classification table showing results for logistic regression predicting
turbidity exceedance at the Arras weir in the Kiskatinaw River. A threshold of
500 Nephelometric Turbidity Units (NTU) was used as the classification
threshold. ........................................................................................................................ 58
Table A.1: Spearman rank correlation coefficients between the North Pacific (NP),
Oceanic Niño Index (ONI), Southern Oscillation Index (SOI) and Pacific Decadal
Oscillation indices and timing of incremented fractional runoff volume (RX%)
and total runoff volume (TRV). Intensity of blue (red) shading corresponds to
magnitude of the negative (positive) rank correlation. Note that at 48 degrees of
freedom a correlation of 0.28 is significant at a probability of 0.05.All seasonal
averages were calculated for the previous year leading up to the runoff period (25
March-1 July). Blank cells indicate no significant correlation exists. ........................... 77
Table D.1: Multiple regression results predicting the R5%. ................................................... 97
Table D.2: Multiple regression results predicting the R10%. ................................................. 98
Table D.3: Multiple regression results predicting the R15%. ................................................. 98
Table D.4: Multiple regression results predicting the R20%. ................................................. 98
Table D.5: Multiple regression results predicting the R25%. ................................................. 99
Table D.6: Multiple regression results predicting the R30%. ................................................. 99
Table D.7: Multiple regression results predicting the R35%. ................................................. 99
Table D.8: Multiple regression results predicting the R40%. ............................................... 100
Table D.9: Multiple regression results predicting the R45%. ............................................... 100
Table D.10: Multiple regression results predicting the R50%. ............................................. 100
Table D.11: Multiple regression results predicting the R55%. ............................................. 101
Table D.12: Multiple regression results predicting the R60%. ............................................. 101

xv
Table D.13: Multiple regression results predicting the R65%. ............................................. 101
Table D.14: Multiple regression results predicting the R70%. ............................................. 102
Table D.15: Multiple regression results predicting the R75%. ............................................. 102
Table D.16: Multiple regression results predicting the R80%. ............................................. 102
Table D.17: Multiple regression results predicting the R85%. ............................................. 103
Table D.18: Multiple regression results predicting the R90%. ............................................. 103
Table D.19: Multiple regression results predicting the R95%. ............................................. 103
Table D.20: Multiple regression results predicting total spring runoff volume (ha m). ....... 104

xvi
Acknowledgments
The completion of this thesis was made possible only through the dedicated help of
many people, and they deserve recognition. I am certainly grateful to have had the
opportunity to work under the guidance of Dr. Stephen Déry. The breadth of understanding
in the field of hydro-meteorology he demonstrates and his eagerness to teach others are
qualities I sincerely admire. Most of all, I am inspired by him as a mentor. I would like to
thank my supervisory committee members Dr. Mike Gillingham and Dr. John Rex who both
provided constructive feedback and support throughout this study. I would like to thank
Vanessa Foord and Richard Kabzems (FLNRORD) for all their support and for their
contributions the collection of valuable climate data in the Kiskatinaw watershed. I also wish
to thank Dr. Suzan Lapp (BC Oil and Gas Commission) for her insightful input and careful
review of this thesis.
I would like to thank the members of the Dawson Creek Watershed Stewardship
Program for their many contributions to this study and their important work in the
Kiskatinaw. I would also like to thank the members of the Bulkley Valley Research Centre
for their financial assistance, and providing me the opportunity to share my work in Smithers,
BC. Special thanks are due to the City of Dawson Creek for financial support.
I am so very grateful to my family and friends who have provided so much love and
guidance. Thank you.

FORECASTING SPRING RUNOFF IN THE KISKATINAW RIVER
Chapter 1:
1.1

Introduction

Overview
Water-supply forecasting is the science of predicting the volume of water that will

flow past a given point over a specified time period (USDA, 1990). Water supply is an
important consideration for both local and regional water-resource planners. Depending on
the region, the amount and timing of available water can affect municipal supply, ecological
needs, and agricultural as well as industrial uses. In addition, water quality and permitted
polluted discharges often depend on water supply, affecting river ecology and recreational
opportunities (Dunne & Leopold, 1998). Accurate hydrologic forecasts are particularly
important in municipal water supply and hydropower production because they are highly
sensitive to changes in climate (Alemu, Palmer, Polebitski, & Meaker, 2011).
Measurements that determine the amount of water stored in the snowpack, or snow
water equivalent (SWE), are frequently used to estimate and forecast the total amount of
runoff derived from snowmelt (USDA, 1972). Typically, seasonal water-supply forecasts are
generated monthly from January–June for most locations. In some cases, mid-month or more
frequent forecasts are generated by requests from users. Forecasts are updated monthly to
account for changes in weather and hydrologic conditions.
Models implemented for the purpose of generating prediction of streamflow can
typically be categorized as either deterministic or stochastic and spatially-lumped or
distributed (Beven, 2012). Deterministic models output a single prediction as a function of a
set of input variables, whereas stochastic models output a range of outcomes based on a
statistical distribution of input variables. Most continuous rainfall-runoff models historically
are deterministic. Stochastic models, however, are often applied for the purpose of predicting

1

2
spring streamflow totals using snowpack and precipitation measurements as inputs with
reasonable results (USDA, 1990).
Two approaches to modelling in hydrology are widely utilized. The first takes an
inductive approach relating inputs to outputs without adhering to the physics of the system.
The other, a deductive approach, requires that the model must adhere as close as possible to
our understanding of the physical processes at play, providing the ability to predict with
confidence values not within the range of the input data. Although the inductive approach is
often criticized for revealing little about the behaviour of the system, in many cases a databased analysis will lead to different conclusions about the behaviour of a system than a
theoretical counterpart may suggest. An inductive approach often suggests that our
theoretical understanding of the system may be incomplete. This quality makes the inductive
approach a suitable start for any rainfall-runoff modelling study (Beven, 2012). The most
commonly applied approach for making operational seasonal streamflow predictions is
through regression techniques. In the past, forecasts were usually made by a single
“optimized” regression and required the assumption of average future conditions. More
recently, different equations that change with time are utilized, and use only inputs available
at the time of the forecast. This approach is more difficult, but consistently provides
improved forecasts (USDA, 1990).
By incorporating climate information into hydrologic streamflow forecasts,
significant advancement in the ability to predict streamflow at seasonal leads (~4 months) has
been made over the past decades. The three main contributors to improved forecast skill are:
accurate forecasting of meteorological anomalies, the quantification of the water held within
the snowpack, and knowledge of soil-moisture conditions (Mahanama, Livneh, Koster,
Lettenmaier, & Reichle, 2012). The role of large-scale climate states such as the El Niño-

3
Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) on regional
streamflow variability in British Columbia (BC) and elsewhere has been affirmed in
numerous studies (e.g., Gobena, Weber, & Fleming, 2013; Grantz, Rajagopalan, Clark, &
Zagona, 2005). Researchers often utilize links between these large-scale atmospheric
circulation features (e.g., anomalies in 700 hPa geopotential heights, sea-surface
temperatures, etc.) and hydrologic conditions for forecasting streamflow at seasonal lead
times. Despite the strong links between these teleconnections and hydro-climatology in
north-western North America, there is significant spatiotemporal variation in the strength and
duration of effects on annual cycles of precipitation and temperature (Fleming & Whitfield,
2010). In addition, slight changes in large-scale atmospheric patterns commonly result in
major differences in surface climate suggesting that the standard indices of these atmospheric
patterns may not be optimal for describing variability in hydro-meteorology at a given
location and nonstandard indices might increase predictive performance (Regonda,
Rajagopalan, Clark, & Zagona, 2006).
In this study: (i) I explore the ability of regional surface climate inputs aggregated
over winter in combination with global circulation indices to make predictions of cumulative
streamflow in the Kiskatinaw River (KR) during spring runoff. (ii) I present spatial
correlations between climate metrics at the surface and runoff in the KR. (iii) I implement a
logistic regression model to predict the probability of exceeding the pumping threshold of
500 nephelometric turbidity units (NTU) based on the forecasted discharge. (iv) I discuss the
spatial structure of the correlations between regional surface climate and streamflow in the
KR, and relate them to established global circulation patterns. (v) I present a series of
regressions for predicting cumulative streamflow in the KR from aggregated surface climate
metrics.

4
1.2

Uncertainty in Hydrology
Models in hydrology serve as a sound method through which we can explore the

effects of alternative water-resource policies and systems, or predict future conditions.
Because all hydrologic models are simplifications of reality, and many driving factors are
only known so far into the future, they can only serve as an informed estimate of future
conditions (Loucks & Beek, 2005). There are three forms of uncertainty inherent in
hydrology: informational, model, and numerical uncertainties. Informational uncertainty
includes the spatial and temporal variability of hydrologic input data, as well as errors
associated with measuring both independent and dependent variables. Model uncertainty
includes incomplete conceptualization of model structure, uncertainty in parameter values,
and issues with inconsistency along temporal and spatial scales. Numerical uncertainty stems
from the wide use of algorithms that suffer from errors associated with truncation when
performing approximations, such as Taylor approximations (Loucks & Beek, 2005). A type
of uncertainty specific to operational streamflow forecasting is that associated with future
weather. Examples include abnormally intense storms and warm spring rains during
snowmelt; both would affect the timing and volume of runoff considerably (USDA, 1990).
As the melt season progresses, uncertainty associated with weather diminishes (Table 1.1).
Table 1.1: Typical percentage of streamflow inputs (snowmelt and rain) known in snowmelt
dominated watersheds (northern hemisphere) with increasing time. Table recreated from
USDA (1972).
Spring and
Factors Affecting
Snow Water Content
Date
Summer Rainfall
Streamflow Volumes
(% Known)
(% Known)
(% Known)
1 Jan
~50
~0
~25
1 Mar
~75
~0
~50
1 Apr
~100
~25
~75
1 May
~100
~50
~90

5
1.3

Hydrologic Response to Global Climate Change
Operational seasonal streamflow forecasting is complicated by the variability in the

water balance that is introduced as anthropogenic warming and land cover change continue to
affect global climate and hydrology (Koster et al., 2017). Snowmelt-dominated regions of the
globe, generally occurring at latitudes greater than about 45° north and south, are
significantly impacted by the effects of global warming (Barnett, Adam, & Lettenmaier,
2005). For instance, earlier onset of springtime snowmelt and streamflow has been observed
for most of western North America in response to warmer winter and spring temperatures
and increased precipitation as rainfall (Stewart, Cayan, & Dettinger, 2004). Throughout BC,
the projected impacts of climate change on hydrologic regimes vary. In the Nechako River
Basin, climate change is projected to result in warmer, wetter conditions, and have important
implications on releases at the Nechako Reservoir (Sharma & Déry, 2015). Snow water
equivalent (SWE) is projected to decline throughout the Peace and Campbell River Basins
and at lower elevations within the Columbia River Basin for the A1B scenario ensemble
median. At higher elevations (~1800 m) within the Columbia River Basin, SWE is projected
to increase in response to more winter precipitation. In these three watersheds, timing shifts
in runoff are expected (Schnorbus, Werner, & Bennett, 2014) (Figure 1.1). Coastal rain and
snowmelt-dominated regimes are the most sensitive to these changes, shifting to more
rainfall-dominated systems by the middle of the century. The two interior basins are
projected to remain as nival regimes by mid-century; however the timing of streamflow is
projected to shift to earlier freshet onsets due to increased mid-winter rainfall and snowmelt
(Schnorbus, Werner, & Bennett, 2014). In the Fraser River Basin, precipitation is projected
to increase slightly; however the fraction of precipitation falling as snow is projected to
decrease by nearly 50% by the 2050’s mainly over low to mid elevations. The onset of

6

Figure 1.1. Contour plot of monthly change in snow water equivalent (SWE) and runoff as a
function of elevation for the Peace River Basin and the Upper Columbia River, British
Columbia based on the Intergovernmental Panel on Climate Change A1B scenario ensemble
median. Figure reproduced from Schnorbus, Werner, & Bennett (2014), used with permission
from the American Meteorological Society.

7
snowmelt is projected to be 25 days earlier than the historical average, resulting in increased
winter and spring runoff and reduced low-flows during summer (Islam & Déry, 2017). The
changes in climate imposed by anthropogenic warming make it imperative that the
assumption of stationarity be carefully examined when employing statistical techniques in
hydrology (Montanari & Koutsoyiannis, 2014).
1.4

City of Dawson Creek and the Kiskatinaw Watershed
The KRB is defined by the Farmington Water Survey of Canada (WSC) hydrometric

station (ID 07FD001). The CDC and surrounding areas have utilized the KRB as a sole water
supply since the mid 1940’s, withdrawing water at the Arras pump house. The KR
headwaters originate in Bearhole Lake, ~60 km southwest of the CDC. During the Second
World War, the US Army established the water system for the CDC to aid in the construction
of the Alaska Highway. Since then there have been a number of enhancements to the system
to improve both water quality and quantity. Upwards of 20,000 people rely on the water from
the KR above the intake at the Arras weir. There is a single intake on the KR located at Arras
16 km west of the city limits. An intake pipe withdraws water at the Arras pumping station.
The pumping capacity of the station is 7,570 m3 per day. A weir to raise water levels during
periods of low flow was constructed 100 m downstream in 1992. Between the Arras pump
house and the water treatment plant (16 km) there is a booster pump and five reservoirs.
There is a total storage capacity of 704,087 m3 within the system. The KR is a naturally
turbid river; the CDC is forced to shut down pumping when turbidity levels exceed 500
nephelometric turbidity units (NTU) to avoid damage to pumping infrastructure. Periods of
high turbidity tend to be correlated with the spring freshet in the KR, and the times when the
CDC is licensed to pump. During periods exceeding 500 NTU, the CDC depends on stored

8
water. It can take considerable time to replenish what is lost during these periods (Dobson
Engineering Ltd. & Urban Systems, 2003).
Peak water demand typically occurs during June with a volume about 18% higher
than the rest of the year. Although little information is available regarding the timing of
withdrawn water by other uses, it can be estimated that the majority (~80%) of this volume
would be abstracted either during the spring freshet or during the summer. Additionally, the
natural gas industry has substantial permitted annual abstractions, but the actual amounts and
timing of these abstractions is poorly recorded, if at all. The potential issue for water
managers at the City lies in the seasonality of water use by the natural gas industry, and the
assumption that the bulk of water withdrawn by industry occurs during periods of high flow
(Forest Practices Board, 2011). There is potential for water shortage during periods of low
flows if this is not true.
1.5

Hydro-Meteorological Changes in the Kiskatinaw Watershed
There is a large body of evidence identifying global warming (Pachauri & Meyer,

2014). In general, studies have shown that the effect of climate change in the Peace River
region by 2050 will result in higher discharges during winter and an earlier freshet with
reduced flows during summer and fall in response to region-wide warming and increased
precipitation in the fall, winter and spring (e.g., Murdock, Sobie, Werner, & Shresthra, 2014;
Schnorbus, Werner, & Bennett, 2014). Several studies have reported changes in hydrology in
response to a changing climate in the KRB. These shifts in the hydrologic regime of the KRB
were confirmed in a cumulative effects assessment carried out by the Forest Practices Board
in 2011. One of the findings of the study was that peak discharge at WSC 07FD001
Farmington gauge for the period 1987–2007 declined by 45% from the prior period of 1966–
1986. However, these finding are likely influenced by the phase shift of the PDO from

9
negative to positive around 1978. More recently, peak discharge in response to snowmelt
occurs 13 days earlier, and a 14% reduction in peak flows exists between periods 1966-1991
and 1991-2016 (Figure 1.2).
The reduced flows during the latter period are likely in response to warming air
temperatures, particularly in winter and spring, leading to reduced winter snow accumulation
and earlier snowmelt. In the study Climate change induced precipitation effects on water
resources in the Peace Region of British Columbia, Canada, the Gridded Surface Subsurface
Hydrologic Analysis (GSSHA) modelling system was used to study the impacts of future
climate change induced precipitation on runoff in the KRB under the A2 and B1 greenhouse
gas emissions scenarios for 2020–2040 from the Intergovernmental Panel on Climate Change
(IPCC) (Saha, 2015). Projections from this study resulted in 5.5% and 3.5% increases in
annual precipitation for the A2 and B2 scenarios, respectively, compared to the reference
period (2000–2011). Annual mean temperature increases of 0.76°C and 0.57°C for the A2
and B2 scenarios were also projected, respectively. The projected result of the increasing
precipitation and air temperatures under climate change is a mean annual streamflow increase
of 15.5% and 12.1% for the A2 and B2 scenarios, with the greatest increase occurring during
spring for both scenarios (Saha, 2015).
1.6

Turbidity
Turbidity is an expression of the optical properties of a water sample that causes light

rays to be scattered and absorbed rather than transmitted in straight lines through the sample
and can be expressed as NTU (Anderson, 2005). Turbidity results from the presence of
suspended and dissolved matter such as clay, silt, finely divided organic matter, plankton,
organic acids, and dyes (American Water Works Association, 2005). The magnitude of

10

Figure 1.2. Average daily discharge at Water Survey Canada (WSC) 07FD001 Farmington
stream gauge for two 26 year periods 1966–1991 and 1991–2016. Peak spring flow during
the latter period is reduced by 14% and occurs 13 days earlier.

11
turbidity in streams and lakes is often proportional to suspended
sediment concentration, and can be related with regressions (Rasmussen, Gray, Glysson, &
Ziegler, 2009). Suspended sediment can negatively affect water supply, recreation, aquatic
life, flood control, transportation, fisheries, reclamation, and navigation (Angino & O'Brien,
1967).
1.7

Research Objectives

This study aims to provide the CDC with a decision support system (DSS) that will enable:
1) Spring runoff predictions of turbidity levels, volume and timing for the KR utilizing
readily available climate and turbidity data.
2) Identifying significant relationships between surface regional climate observations and
streamflow characteristics in the KR.
3) Building a stochastic model framework that may be extendable to other watersheds /
municipalities.
The first two chapters of this thesis provide context regarding hydrologic forecasting,
including limitations, and introduces the KRB and its importance to the CDC and
surrounding areas. Chapters 3 and 4 describe the datasets and methods, respectively, used for
developing the regression models implemented in predicting spring runoff in the KR. The
resulting correlations and models are presented in Chapter 5. Finally, Chapter 6 provides a
brief discussion of the results in the context of global circulation patterns, emphasizes the
difference between causal and predictive modelling, and suggests areas of future research.
Finally, the main conclusions of the study are summarised.

12
Chapter 2:
2.1

Study Area

General Description
The Kiskatinaw River Basin (KRB) delineated at the Arras weir is a 2,800 km2

watershed located southwest of the City of Dawson Creek (CDC), BC. This delineation
defines the land upstream of the CDC’s municipal intake at Arras, 78% of the total KRB
basin area (3,600 km2). The KRB extends roughly from 55°N-56°N latitude and
120°W-121°W longitude with elevation ranging from 579 m ASL at the outlet to 1428 m
ASL above Bearhole Lake (Figure 2.1). The KRB lies within the Boreal White and Black
Spruce (BWBS) bio-geoclimatic (BGC) zone (Meidinger & Pojar, 2017). The BWBS zone
contains several upland forest types including mixed stands of trembling aspen (Populus
tremuloides) and white spruce (Picea glauca) and mixed stands of lodgepole pine (Pinus
contorta) and black spruce (Picea mariana). Long, cold winters and short, but relatively
warm summers characterize climate in this BGC zone. The BWBS zone receives relatively
little precipitation and the least amount of snowfall of all northern BGC zones. Land cover in
the KRB consists primarily of several forest types, pasture and cropland (Table 2.1). The
KRB also provides important habitat for the declining woodland caribou (Rangifer tarandus)
(Committee on the Status of Endangered Wildlife in Canada, 2002). The topography of the
BWBS zone varies, with the KRB lying within the rolling Rocky Mountain foothills in the
northwest region of the zone (Meidinger & Pojar, 2017). The hypsometric distribution for the
KRB defined at Arras reveals that the geomorphology of the KRB is more eroded (Figure
2.2). The concave form suggests that the basin shows less total runoff with a higher fraction
of surface response than a convex counterpart primarily due to shallower water-table

13

Figure 2.1. Digital elevation map (DEM) of the Kiskatinaw River Basin (KRB) delineated at
the Arras weir (circle) where the City of Dawson Creek (CDC) withdraws water and
delineated at the Water Survey of Canada (WSC) 07FD001 Farmington stream gauge
(triangle).

14
Table 2.1: Percent of total area of the Kiskatinaw River Basin defined at Arras for existing
land cover/land use categories in the watershed as of 2010. Recreated from Paul (2013).
Land Cover Type
Total Area (%)
Cropland
0.7
Coniferous Forest
39.1
Deciduous Forest
28.8
Mixed Forest
12.9
Planted Forest
5.5
Forest Fire
0.9
Cut Block
2.1
Pasture
0.7
Water
6.5
Wetland
1.7
Built-up Area
1.2

Figure 2.2. Relative height versus relative area in the Kiskatinaw River Basin defined at
Arras where a is the area of basin above height h, A is the total basin area, h is the height
above the outlet, and H is the total relief of the basin.

15
positions and lower surface gradients (Vivoni, Benedetto, Grimaldi, & Eltahir, 2008). The
bedrock geology of the KRB is composed predominantly of interbedded cretaceous shales
sandstones while the surficial geology consists mainly of lacustrine deposits with alluvium
(Dobson Engineering Ltd. & Urban Systems, 2003).
2.2

Climate
Average air temperatures for the KRB are coldest during winter (-11.7°C) and

warmest during summer (13.1°C). Spring remains relatively dry, accounting for 18% of the
total annual precipitation, whereas the summer is relatively wet accounting for 39% of the
total annual precipitation (Table 2.2).
Table 2.2: Average daily minimum (Tmin) and maximum (Tmax) temperature and the total
seasonal precipitation (Total Ppt) over the Kiskatinaw River Basin defined at the Arras weir
(Figure 2.1) for 1950-2010. Values were calculated by Aseem Sharma using the R statistical
programming language and climate data provided by the Australian National University
(2017).
Variable
Statistic
Annual Winter Spring Summer Fall
Minimum
303.0
40.1
35.1
208.2
41.3
Maximum
864.8
230.3
177.6
419.7
218.6
Total Ppt (mm)
Mean
530.4
105.4
96.7
208.2
119.8
SD
100.6
40.7
29.7
72.2
39.4
Minimum
-8.1
-28.0
-9.2
4.7
-8.0
Maximum
-1.7
-9.8
-1.2
8.5
0.0
Tmin (°C)
Mean
-4.5
-17.1
-4.6
6.8
-3.4
SD
1.3
4.0
1.7
0.8
1.7
Minimum
4.5
-15.5
2.3
17.2
1.3
Maximum
9.3
-0.3
10.8
22.6
11.3
Tmax (°C)
Mean
7.0
-6.2
7.4
19.4
7.1
SD
1.2
3.5
1.8
1.2
2.0
On average, about 36% of the total annual precipitation falls as snow (Dobson
Engineering Ltd. & Urban Systems, 2003). Snowfall has important social, economic and
ecological implications in cool temperate and subarctic regions. The amount of water held

16
within the snowpack is of particular interest to water managers because its melt replenishes
water storage infrastructure and groundwater, and dictates the amount of water available for
irrigation, municipal use, power generation and ecological needs (Dunne & Leopold, 1998).
Poleward of 45° north and south in latitude, the hydrologic cycle is strongly influenced by
the accumulation and ablation of snow (Barnett, Adam, & Lettenmaier, 2005). Studies in
these snowmelt-dominated regions have provided evidence that, particularly in western North
America, as surface air temperatures rise and precipitation increases in response to
anthropogenic warming, more precipitation falls as rainfall and snow ablation onset is earlier
(e.g., Barnett, Adam, & Lettenmaier, 2005; Kang, Shi, Gao, & Déry, 2014; Stewart, Cayan,
& Dettinger, 2004). Warmer winter temperatures, in combination with a shift in the
precipitation regime, directly influence the timing of streamflow in response to snowmelt
(Kang, Gao, Shi, Islam, & Déry, 2016). These shifts in both climate and hydrology have
important impacts on regional water supplies (Barnett, Adam, & Lettenmaier, 2005).
2.3

Hydrology
Five major sub-basins exist in the KRB delineated at Arras. The East Kiskatinaw of

the KRB defines one of these basins and has its headwaters in Bearhole Lake, accounting for
1009.7 km2. Bearhole Lake is a controlled reservoir and represents an important source of
municipal water supply for the CDC contributing about 396,000 m3 during winter months.
Groundwater recharge likely represents a substantial input to the total water supply within
this reservoir with an estimated net groundwater influx of 209,000 m3 (Scruton, 2017). The
other major headwater basin is the West Kiskatinaw including the West Headwaters and
Jackpine covering 1004.1 km2, the largest of the basins. The West and East Kiskatinaw drain
into the Mainstem (430.1 km2) along with the smaller Halfmoon-Oetata and Brassy basins

17

Figure 2.3. Map depicting the major and sub-basins of the Kiskatinaw River Basin delineated
at the Arras weir where the City of Dawson Creek withdraws water for distribution and
treatment.

18
(Figure 2.3). Based on the annual hydrograph the lowest flows occur from September–
March, likely consisting primarily of groundwater discharge. This is followed by a swift
increase in streamflow in response to spring snowmelt lasting into mid-May when sporadic
summer rains augment flows well into July. By mid to late August a period of low flows
persists into winter until the next spring snowmelt (Figure 2.4). Fluctuations from year to
year are largely a result of synoptic-scale climate variations. Highest flows are typically
associated with spring snowmelt; however, intense summer rainstorms can cause high peak
flows comparable to those generated by the onset of snowmelt. An example of such an event
is 15 June 2016, when a low pressure system stalled in the region and upslope flow along the
northern Rockies caused extensive flooding in the Peace River region including Dawson
Creek (Figure 2.5). There is a strong link between the onset of snowmelt indicated by the
rapid decline in snow covered area within the KRB and a seasonal peak in discharge
occurring near 8 May, emphasizing the importance of snowmelt on both timing and volume
of spring and summer discharge in the KR (Figure 2.4).
The long-term annual discharge for WSC Farmington averages 10.7 m3·s-1 for the
period 1966–2016. Using this value, the deviations of annual means from the long term
mean-annual discharge are plotted identifying relatively wet and dry years in the KRB
(Figure 2.6). The wettest year on record is 1997, and the driest year on record is 1992 (Figure
2.6). Although seasonal low flows typically occur during the winter months, the lowest flows
on record occurred in response to drought conditions during August of 1992 (Dobson
Engineering Ltd. & Urban Systems, 2003). There is evidence in historical streamflow records
that there may be a trend to higher peak flows during the late summer months (Forest
Practices Board, 2011). This could result in the KR experiencing bi-modal seasonal peaks,

19

Figure 2.4. Minimum (Min Q), mean (Mean Q) and maximum (Max Q) discharge at Water
Survey Canada (WSC) 07FD001 Farmington stream gauge, 1966–2016. The black dotted
line is the mean snow covered fraction of the Kiskatinaw River Basin (KRB) derived from
satellite imagery over 2002-2015. A significant increase in streamflow during spring at
Farmington corresponds with onset of snowmelt in the KRB.

20

Figure 2.5. Hydrograph for Water Survey Canada hydrometric gauge at Farmington in 2016.
During this year a peak discharge in response to snowmelt is seen occurring March–May;
however, a peak of greater magnitude is also observed in response to an intense cyclone that
brought copious amounts of rain to the Peace River Region, resulting in severe flooding over
10–16 June.

21

Figure 2.6. Analysis of wet (above zero) and dry (below zero) years based on mean Water
Survey Canada Farmington discharge, 1965–2016. Mean discharge for the period is 10.7 m3
·s-1.

22
one in the spring and one in the late summer. This would likely result in higher suspended
sediment loads during the summer (e.g., Figure 2.5).
2.4

Water Consumption
There are currently 15 licenses granted by the BC Ministry of Environment under Part

2 of the Water Act that allow a maximum abstraction of water from the KR totaling
24,200,000 m3 annually (Figure 2.7) (Government of British Columbia, 2017). This does not
include non-stream abstractions. Water withdrawals from the KR tend to be seasonal,
occurring mostly during the peak flow or transition period (Forest Practices Board, 2011).
Outside of the records kept by the Water and Environmental Department of the CDC
regarding abstracted volumes at the Arras pump house there is little available information
about the actual amount of water abstracted from the KR. Water use by the CDC has steadily
increased since 1960 when abstractions were 0.034 m3·s-1. In 2007, usage increased to 0.077
m3·s-1. Maximum usage occurs during June with an 18% higher than average volume than for
the remainder of the year (Forest Practices Board, 2011). As mentioned in Section 1.4, little
is known about the actual amounts of the total annual allocated water abstracted from the
KRB.

23

Figure 2.7. Amount of currently licensed water allocations from the Kiskatinaw River by
purpose represented as the percent of the total annual allocated volume (24,200,000 m3) as of
15 April 2018. Four categories dominate water usage, the remaining four make up less than
0.5% of the total volume. Data were obtained from the online Government of British
Columbia Water License Query tool (Government of British Columbia, 2018).

24
Chapter 3:

Data Sources and Assimilation

This chapter describes the data sources that were assimilated for the purpose of
predicting spring discharge in the KR and explains why they were considered for model
input. A brief synopsis of synoptic-scale weather patterns in BC and their links to regional
climate anomalies is provided. The second section details specific data sources and where
they were obtained, including basic statistics performed and quality assurance/quality control
(QA/QC) measures taken.
3.1

Synoptic Weather Patterns and Surface Climate Variability in British Columbia
Typical data inputs utilized for the purpose of water supply forecasts in snowmelt-

dominated watersheds include SWE, precipitation, and antecedent streamflow. Of the inputs
affecting watershed runoff, melted snow and rain are most directly related to seasonal
streamflow (USDA, 1990). Variability of these inputs has been linked to variation in synoptic
weather processes (Moore & McKendry, 1996). Many studies have identified correlation
between major Pacific climate indices and surface climate in BC (e.g., Moore & McKendry,
1996; Shabbar & Khandekar, 1996; Shabbar, Bonsal, & Khandekar, 1997). Furthermore,
associations between the frequency of certain synoptic circulation types and surface climate
anomalies over BC have been identified in relation to major Pacific climate indices (Stahl,
Moore, & McKendry, 2006). Snow-covered fraction (SCF) at the basin scale has also been
shown to be a critical predictor of spring runoff response (Tong, Déry, & Jackson, 2009).
3.2

Data Sources
To determine the ability of surface-weather observations to capture the variability in

winter climate driving snow accumulation and precipitation in the KRB, several data sources
were assimilated. These sources included active weather station data within 1200 km
(features of ~1000 km are said to be of synoptic scale, or roughly the size of a mid-latitude

25
cyclone) of the geometric centroid of the KRB, manual observations of SWE on 1 March and
1 April at established snow courses, several indices capturing ocean-atmosphere-streamflow
teleconnections and the fraction of basin area covered in snow (SCF) for the KRB defined at
WSC Farmington obtained from satellite imagery (Figure 3.1). There were 78 climate
stations and 68 snow survey sites obtained within 1200 km of the KRB. Discharge data were
retrieved from Environment and Climate Change Canada’s (ECCC) WSC 07FD001
Farmington stream gauge because of the long period of record, and the reliability of its
measurements (Figure 2.1). The specific datasets and their providing agencies are detailed in
Table 3.1.
Environment and Climate Change Canada historical climate data. The ECCC
maintains a network of manual and automated weather stations throughout Canada. The
climate variables recorded at each station vary, but typically include air temperature,
atmospheric pressure, precipitation, solar radiation, wind speed, and wind direction. Past data
from the ECCC stations were obtained from the ECCC Historical Climate Data webpage
(Environment and Climate Change Canada, 2017a). To identify the optimal set of climate
stations to include, the same physiographic regions identified in Moore & McKendry (1996)
were utilized. Only active climate stations with a period of record longer than ten years were
retained. For each region the active stations were sorted by elevation, obtaining suitable low,
medium and high stations. The period of record of the stations was also taken into account,
giving preference to stations with longer records. A daily time step was implemented for
analysis. For each station the entire period of record was downloaded back to 1965 if data
were available. Temperature was averaged and precipitation was totaled for the period 15
November–25 March, following Moore & McKendry (1996). Additionally, snow depth at
each climate station was averaged for the period 24–25 March, capturing snow

26

Figure 3.1. Map displaying the location of active automatic weather stations from both
Environment and Climate Change Canada and the BC Ministry of Transportation and
Infrastructure. Also included are the locations of manual surveys of snow water equivalent
obtained from the BC Ministry of Forests, Lands, Natural Resource Operations, and Rural
Development. The black circle denotes the 1200-km synoptic radius encompassing the
Kiskatinaw River Basin.

27
Table 3.1: Summary of datasets utilized for predicting spring runoff in the Kiskatinaw River
along with host organization and observation details.
Organization
Dataset
Observation Type
Environment and Climate
Past weather and climate
Surface temperature,
Change Canada, 2017a
precipitation, and snow
depth (daily)
Environment and Climate
Change Canada, 2017b

Water Survey Canada

Historical discharge (daily)

BC Ministry of
Transportation and
Infrastructure, 2017

Avalanche and weather

Surface temperature,
precipitation, pressure and
snow depth (daily)

Ministry of Forests, Lands,
Natural Resource
Operations, and Rural
Development, 2017

Manual snow survey
observations

Manually surveyed snow
water equivalent

National Snow and Ice Data
Center, 2017

MOD10A2

Remotely sensed eight day
snow cover

National Oceanic and
Atmospheric
Administration, 2017a

NCEI Pacific Decadal
Oscillation Index

Based on NOAA’s extended
reconstruction of SSTs
(ERSST version 4)

National Oceanic and
Atmospheric
Administration, 2017b

Oceanic Niño Index (ONI)

National Oceanic and
Atmospheric
Administration, 2017c

Southern Oscillation Index
(SOI)

Three-month running mean
of ERSST.v4 SST
anomalies in the Niño 3.4
region (5°N–5°S, 120°W–
170° W)
Standardized index based on
the observed differences in
sea level pressure between
Tahiti and Darwin, Australia

National Center for
North Pacific (NP) Index by
Atmospheric Research, 2016 Trenberth and Hurrell
(1994)

Area-weighted mean sea
level pressure over the
region (30°N–65°N and
160°W–140°W)

28
accumulation near the beginning of ablation, much like peak SWE. For each hydrologic year
(HY) associated with the start of the snow accumulation period starting on 1 September,
average temperature, atmospheric pressure, and total precipitation were calculated only if
95% of total possible observations were present over the period 15 November-25 March to
avoid introducing biased observations.
BC Ministry of Transportation and Infrastructure: Avalanche and Weather
climate data. The BC Ministry of Transportation and Infrastructure (MOTI) weather station
network is intended primarily for winter use and is maintained from October–May. Some
stations remain operational over the summer, but many are decommissioned due to the risk of
lightning strikes (BC Ministry of Transportation and Infrastructure, 2017). Data are recorded
hourly and include temperature, atmospheric pressure, humidity, precipitation, snow depth,
wind speed, and wind direction. All active stations within BC were considered for model
input (Figure 2.1). Similar to ECCC climate data, MOTI stations were selected based on their
elevation and period of record for each physiographic region within BC. Atmospheric
pressure and temperature were averaged and precipitation totaled annually over the period 15
November–25 March for each station. Additionally, snow depth at each location was
averaged for the period 24–25 March. For each HY, average temperature, atmospheric
pressure, and total precipitation were calculated only if 95% of total possible observations
were present over the period 15 November-25 March to avoid introducing biased
observations.
Manual snow survey observations. One of the services provided by the BC Ministry
of Forests, Lands, Natural Resource Operations, and Rural Development (FLNRORD) is
gathering data about snowpack conditions through manual snow surveys. Snow-survey
measurements have been conducted by FLNRORD since 1935 in BC. Because quantifying

29
the entire snowpack over the extent of a watershed is an especially difficult task, snow
surveys carried out in the same location from one season to the next serve to provide a good
index of relative snowpack conditions. This information is valuable in water supply and flood
forecasting (Ministry of Forests, Lands, Natural Resource Operations, and Rural
Development, 2017).
In this study, active snow surveys within BC were obtained online from the Ministry of
Environment Manual Snow Survey Observation Data Archive (Ministry of Environment,
2017). Average SWE occurring on 1 March and 1 April were obtained for each year only if
there were no survey issues indicated by the survey flag code. The 1 March and 1 April snow
surveys were used in this study because these dates generally coincide with peak SWE
depending on the region. Survey locations were selected based on physiographic regions and
elevation in the same manner as the weather station data.
Satellite-derived snow covered fraction. The Earth Observing System (EOS) is one
of the main contributions to the NASA initiated Mission to Planet Earth now known as
NASA’s Earth-Sun Systems Missions. As part of this program two key satellites were
launched into polar orbits: Terra in 1999 and Aqua in 2002. The focus of these missions was
to provide scientific understanding of the earth-sun system and its response to human induced
changes. Of the five sensors included on the Terra platform, one consists of the Moderate
Resolution Imaging Spectro-Radiometer (MODIS). The MODIS sensor provides a two-day
repeat global coverage at spatial resolutions of 250, 500 or 1000 m, depending on the
wavelength and has 36 spectral bands with 12-bit radiometric sensitivity. The total field of
view of MODIS is ±55°, resulting in a swath width of 2330 km. These properties make it an
ideal sensor for various atmospheric, oceanic, and land surface tasks (Lillesand, Kiefer, &

30
Chipman, 2007). One of these products is the snow cover eight day L3 Global 500 m Grid
(MOD10A1) derived from Terra. The spectral bands on MODIS allow it to determine snow
cover on the landscape via the Normalized Difference Snow Index (NDSI). Snow on land is
typically highly reflective in the visible bands and poorly reflective in the shortwave infrared
bands, with the magnitude of this difference being exploited by the NDSI. Maximum snow
cover extent is determined by reading eight days of 500 m resolution MOD10A1 tiles and if
snow is observed in a cell on any day it is classified as snow covered for the eight day period.
If no snow is present it is classified by the clear view observation that occurred most often
(Hall & Riggs, 2016).
For this study, all MOD10A1 tiles back to 1999 were obtained from the National Snow
and Ice Data Center. The MOD10A1 tiles were gridded in the MODIS Sinusoidal Tile Grid,
comprised of 460 non-fill tiles each covering 10° × 10° at the equator. The H11V03 grid cell
is used for this study as the KRB is contained within it. Cloud coverage adversely affects the
overall accuracy of MODIS snow products to detect snow because visible and mid-infrared
radiation cannot penetrate. A spatial-filter method developed in Tong, Déry, & Jackson
(2009) reduces the error introduced by cloud cover through a reclassification process using a
moving window that determines cell identity based on the eight cells surrounding it,
improving the overall accuracy of snow cover detection. This method was applied to all of
the H11V03 tiles. Using a 10 m × 10 m resolution digital elevation map, the extent of the
KRB at Farmington was determined (Figure 2.1). The raster was resampled to the resolution
of the H11V03 tiles. The snow covered fraction (SCF) of the KRB for each tile was obtained
by dividing the number of pixels classified as snow by the total number of pixels within the

31
KRB. The date that SCF reached 50% during melt was correlated with streamflow timing and
volume during runoff in the KRB.
Ocean-Atmosphere teleconnections. Climate variability tends to manifest itself in
certain patterns that allow improved hydrologic forecasting at seasonal lead times (Garbrecht
& Piechota, 2006). Improving understanding of the interactions between atmospheric
circulation and variation in the hydrology and climate at the surface permits a more complete
interpretation of that variability in the context of climate change (Trenberth, 1990). The
influence of the El Niño-Southern Oscillation (ENSO), Pacific Decadal Oscillation (PDO),
and variation in the intensity of the Pacific North American (PNA) circulation pattern have
been linked to variations in temperature, precipitation, snowpack storage and seasonal
streamflow in western Canada (e.g., Gobena, Weber, & Fleming, 2013; Moore & McKendry,
1996; Woo & Thorne, 2003). Because of the influence of these teleconnections on surface
hydrology in BC, correlations between the timing and volume of spring runoff in the KR and
major Pacific teleconnections were explored.
El Niño/Southern Oscillation. The El Niño-Southern Oscillation (ENSO) is the
phenomenon of anomalous warming of sea-surface temperatures (SSTs) over the eastern
equatorial Pacific Ocean with corresponding changes in the atmospheric sea-level pressure in
the tropical Pacific (Garbrecht & Piechota, 2006). The warm phase of ENSO is referred to as
El Niño, while the cool phase is known as La Niña. Typically ENSO has a periodicity of 2–7
years. Normally equatorial trade winds blow towards the west and push along warm surface
water towards the western Pacific, resulting in sea-surface levels and temperatures that are
higher in Indonesia than off the coast of Ecuador. The SSTs are warmer in Indonesia because
of upwelling of cold water off the coast of South America. Under these conditions the warm
water in the west results in rising air and concurrent precipitation, leaving the eastern Pacific

32
relatively dry. Under El Niño conditions, the trade winds weaken causing the thermocline to
fall in the eastern Pacific. The reduction in the thermocline in the eastern Pacific initiates
eastward flow of warm surface waters, which in turn results in changes in the global
atmospheric circulation. Conversely, under La Niña conditions the thermocline falls in the
western Pacific and inclines in the eastern Pacific, resulting in warmer SSTs in the western
Pacific, and overall drier conditions. Western (including BC), Northwestern and Central
Canada experience anomalously mild winters and springs during El Niño. El Niño has been
shown to affect streamflow volumes and patterns in streamflow timing because of its
influence on both winter temperature and precipitation fields over BC (Fleming & Whitfield,
2010).
To capture the influence of ENSO on spring runoff characteristics in the KR, the
Oceanic Niño Index (ONI) was obtained online and used as a potential predictor of spring
runoff in the KR. This index is based on a three-month running mean of Extended
Reconstructed Sea Surface Temperature (ERSST.v4) anomalies in the Niño 3.4 region
(5°N–5°S, 120°–170°W) (Huang et al., 2014). In addition, the Southern Oscillation Index
(SOI), a standardized index based on the observed differences in sea-level pressure between
Tahiti and Darwin, Australia, was included as a possible predictor of spring runoff in the KR.
This index captures the large-scale fluctuations in air pressure occurring between the western
and eastern Pacific during ENSO episodes. The SOI index is highly correlated with changes
in ocean temperatures across the eastern tropical Pacific (National Oceanic and Atmospheric
Administration, 2017c).
Pacific Decadal Oscillation. The Pacific Decadal Oscillation (PDO) is another Pacific
atmospheric pattern that has been associated with variations in the Pacific Basin and North
American climate (Zhang, Wallace, & Battisti, 1997). The PDO has a warm and a cool phase

33
similar to ENSO, and is defined by SSTs in the northeast and subtropical Pacific Ocean.
When SSTs are anomalously cool in the central North Pacific and warm along the eastern
Pacific coast, and sea-level pressures are correspondingly below average over the North
Pacific, the PDO is in its warm phase. Conversely, when SSTs are warm in the central North
Pacific and cool along the North American coast or above average sea-level pressures exist
over the North Pacific, the PDO is in its cool phase (Trenberth & Hurrell, 1994). There is
evidence of only two full cycles of the PDO occurring over the past century. The PDO has
expressed its cool phase for 1890–1924, and again from 1947–1976. The PDO has reflected
its warm phase for periods 1925–1946, and again from 1977 to around 1995 (Minobe, 1997).
Since 1995, the PDO has primarily reflected its cool phase. The influence of the PDO on the
spring hydrologic response in the KR is explored through use of NOAA’s PDO index that is
based on the ERSST.v4 dataset integrated with the Mantua Index (National Oceanic and
Atmospheric Administration, 2017a). The PDO index values were averaged seasonally
(three-month periods) for the year leading up to the runoff period.
Pacific North American pattern. The PNA pattern is a predominant mode of low
frequency climate variability in the Northern Hemisphere extratropics acting primarily during
winter (Wallace & Gutzler, 1981). The pattern materializes as anomalies in the 500 hPa
geopotential height fields. Regionally, 500 hPa field anomalies south of the Aleutian Islands
and over the southeastern United States correspond to anomalies of similar sign over the
Aleutian climatological low. Anomalies of opposite sign occur near Hawaii and over central
Canada during fall and winter (National Oceanic and Atmospheric Administration, 2017d).
The pattern has four centres of action: near Hawaii (20°N, 160°W), over the North Pacific
Ocean (45°N, 165°W), over Alberta (55°N, 115°W) and over the Gulf Coast of the United
States (30°N, 85°W). The PNA index is defined as a linear combination of the normalized

34
height anomalies at the four pattern centres (Wallace & Gutzler, 1981). An alternative index
to the PNA is the North Pacific (NP) index developed by Trenberth & Hurrell (1994)
computed as the area-weighted mean sea-level pressure over the region 30°N–65°N and
160°W–140°W. This index is generally more robust than the PNA and better captures the
intensity of the Aleutian low during winter (Trenberth & Hurrell, 1994). In this study, the NP
anomaly index (National Center for Atmospheric Research, 2016) based on the 1008.9 hPa
sea-level pressure mean over 1925–1989 averaged from November through March was
correlated with runoff metrics in the KR. Seasonal averages of the NP index were also
correlated with runoff metrics in the KR.
Discharge data. There are two sources of discharge data available on the KR. Because
the CDC withdraws water at the Arras weir, discharge at Arras would be the optimal choice
for a source of historical discharge data in the KR. However, the logging system used to
record stage at the Arras weir has been plagued by a communication error resulting in
erroneous spikes in the recorded stage. In addition, the period of record with available
discharge data at Arras is relatively brief (January 2013–2018) and is insufficient for multiple
regression analysis. As an alternative, the WSC Farmington (07FD001) stream gauge has
data beginning in 1944, and is reliable from 1965 to present day. Because it is downstream of
Arras, it adds 800 km2 to the total basin area (Figure 2.1). Using a simple non-linear
regression between concurrent discharge at Arras and Farmington, however, should provide
adequate means of estimating flows at Arras from simultaneous daily average flows at
Farmington. Therefore, discharge data at Farmington were used to calculate incremented
annual cumulative discharge in the KR, as well as total runoff volume during the runoff
period.

35
Chapter 4:

Methods

This chapter describes the derivation of streamflow metrics regarding volume and
runoff timing from the Farmington discharge data. Methods for determining spatial
correlations between the streamflow metrics and aggregated climate inputs are discussed.
The regression modelling approach for predicting spring runoff in the KR is explained.
Finally, the use of logistic regression for predicting turbidity exceedances in the KR is
outlined.
4.1

Timing of Cumulative Discharge and Total Runoff Volume
Discharge data at Farmington were used to calculate cumulative discharge for each

hydrologic year (HY) during the spring runoff period, defined here as 25 March–1 July, since
1965. The hydrologic day of year (HDOY), or days since 1 September, that cumulative
discharge reached 5% intervals (RX%, where X refers to the interval) of total runoff volume
(TRV) for the spring runoff period (25 March–1 July) was also determined (Figure 4.1). The
HDOY that the incremented RX% occur, and the TRV in the KR are what this study aimed
to predict from the aggregated winter climate data, and global circulation indices described in
Chapter 3.
4.2

Spatial Correlations
Correlation analysis was used to determine the spatial relationships among the

surface climate variables described in Chapter 3 with the incremented RX% and TRV in the
KR. To quantify the degree of association, the Spearman rank correlation was used because it
does not require assumptions about the distribution of the input climate data (Woo & Thorne,
2003). For average winter temperature, atmospheric pressure, total 15 November-25 March
winter precipitation and 1 April SWE at each station or survey location, correlations with the
R25%, R50%, R75% and TRV were calculated. The 1 March SWE and snow depth were not

36

Figure 4.1: Diagram depicting the model structure implemented for making predictions of
spring runoff in the Kiskatinaw River. Climate observations obtained at weather stations and
manual snow survey locations were aggregated over the winter accumulation period (15
November-25 March). Cumulative discharge at the Farmington stream gauge (maroon curve)
was determined for each runoff period (25 March-1 July) from the raw discharge data
(dashed line). The day of year that discharge reached specified percentages of total runoff
volume was determined in 5% intervals. For instance, the R50% above (centre vertical red
line) would be the day of year that 50% of the total runoff volume was reached for that year
(2003 in this case). The R25%-R75% depicts the interquartile range. The aggregated winter
climate observations were input into a regression analysis where candidate models were
identified for each RX% and the total runoff volume. The selected models were then used to
predict the cumulative runoff for that runoff period in discrete 5% intervals.

37
included in the correlation analysis, but were still included as potential model input. The
correlations were then plotted spatially and contoured from an inverse distance weighted
interpolation. The contour maps of the Spearman rank correlations between the surface
observations and streamflow were then used for identifying any spatial patterns present.
4.3

Predicting Runoff Timing and Volume
Multiple regression analysis was used to make predictions of the RX% and TRV in

the KR from the input climate variables described in Chapter 3. Because this study aimed
only to make predictions of streamflow and not necessarily attribute causation, no initial
hypothesis specifying relationships between timing and volume of discharge in the KR to the
various predictors was proposed. To this end, the objective of the regressions was predictive
in nature. For a dataset with one response variable and n predictors there are (2n-1) possible
regression models. Given these many predictors, there was a large number of possible
regressions, even with modern computing capabilities. These ranged from models with just
one predictor, to the fully saturated model (n = 396). This raised the dilemma as how to
choose the best model. To avoid choosing models biased by the data, a resampling technique
(described below) was applied to determine unbiased coefficients, and improve out of sample
predictive reliability.
Resampling. To reduce predictive error and quantify bias in the regression
coefficient (β) estimates, a jackknife technique was applied when estimating the β
coefficients. This technique provided a first-order unbiased estimator for the β coefficients in
the multiple regressions (Tukey, 1958). To illustrate, let 𝛽̂ be the value of a β when estimated
from the entire dataset, and let 𝛽̂(𝑗) be the estimate of β with the jth observation omitted,

38
where j=1,…,N. The pseudo-value 𝛽̂𝑗 ∗ for the jtth observation is determined by Equation
(4.1):

̂ 𝒋∗ = 𝑵𝜷 − (𝑵 − 𝟏)𝜷
̂ (𝒋) .
𝜷

(4.1)

The mean of the pseudo-values is the alternative, first-order unbiased estimator of the
β coefficient, and the standard error around the mean of the pseudo-values is an estimator of
the standard error of 𝛽̂ .The bias of an estimator is found as the difference between the
expected value of the estimator, 𝛽̂ , and the unbiased jackknife estimator 𝛽̂𝑗 . In this study, the
decision to determine bias of the model coefficients through resampling was an effort to
improve out of sample predictive performance, and to confirm that the estimated coefficients
were not overly biased by the modelling approach, or by measurement error of the climate
variables. The jackknifed standard error was used in model selection to establish confidence
intervals for each coefficient, rejecting models with insignificant (p > 0.05) coefficients. The
final adj-r2 was calculated from the unbiased model coefficients and used in model selection.
Model Selection. In multiple regression analysis the objective is to choose a model
that has sufficient variables to properly account for the variation in the data and minimize the
residual sum of the squares, but is not simply fit to random noise in the data (Gotelli &
Ellison, 2013). In this study the performance of a model was measured using the adj.-r2 value
as defined in Equation (4.2):

𝑨𝒅𝒋𝒖𝒔𝒕𝒆𝒅 𝒓𝟐 = 𝟏 − (

𝒏−𝟏
𝒏−𝒌

) (𝟏 − 𝒓𝟐 ),

(4.2)

where n is the sample size, k is the number of parameters in the model, and r2 is the
coefficient of determination (Gotelli & Ellison, 2013). Information criterion statistics balance
reduced performance (sum of the squares) with the addition of parameters. The adj-r2 is
essentially the coefficient of determination penalized for additional predictors. There were a

39
number of model selection strategies including forward selection, backward elimination, and
stepwise methods. In this study, all permutations of one, two or three input variables were
iterated using Java and the Apache Commons (2017) Math Version 4.0 statistical package in
Eclipse Neon. For each permutation, the adj.-r2 was determined from the resampled β
estimates. Only models with an adj.-r2 ≥ 0.85 (based on jackknifed coefficients), and n ≥ 15
were retained. For each RX% and the TRV models with the highest adj.-r2were chosen if
they did not violate the assumption of linear regression. These assumptions include linearity,
homogeneity of variance, absence of outliers, and independence. Linear relationships were
confirmed through plots of residuals for each model. To test for heteroscedasticity, the
Breusch-Pagan / Cook-Weisberg test was applied to the final candidate models (Breusch &
Pagan, 1979). To test for outliers, studentized residuals ≥ 2 were flagged for each model. In
addition, possibly influential observations were detected using the DFFIT command (Belsley,
𝑘

Kuh, & Roy, 1980). DFFIT values that exceeded the threshold 2√𝑛 were assessed for
validity. Independence of variables was determined via the variance inflation factor (VIF).
The VIF provides an index to determine how significantly the variance of an estimated
regression coefficient increases due to collinearity (Belsley, Kuh, & Roy, 1980). Models with
any VIF score > 10 were discarded.
Forecast Validation. Once a model was selected for each RX% and the TRV, a
prediction of the cumulative discharge curve for HY 2002 was carried out to demonstrate the
forecasting procedure. To account for uncertainty in the predicted RX% and the TRV, the
residuals for each model were simultaneously predicted and stored. The minimum, first
quartile, third quartile and maximum value for the set of residuals was determined for each
model capturing the uncertainty inherent around each prediction given past observations. To

40
determine the accuracy of the discharge predictions, the Nash-Sutcliffe Efficiency (NSE)
coefficient was calculated from the resulting forecast to provide a metric of the model’s
ability to accurately predict discharge (World Meteorological Organization, 1986).
4.4

Predicting Turbidity Exceedance
The CDC is not only concerned with the timing and volume of runoff, but the

associated turbidity levels. Although there are several methods for predicting sediment
transport in streams, this study utilized a logistic regression model for determining the
probability of exceeding 500 NTU in the KR (the limit at which pumps at the Arras weir
must be shut down to avoid damage). Logistic regression cannot be fit using least squares
methods, but is instead fit using maximum log likelihood estimation. The log likelihood (LL)
was calculated with iteration using maximum likelihood estimation (MLE). The LL serves as
the primary test statistic for logistic regression (Gotelli & Ellison, 2013). To determine the
accuracy of the logistic model, the area under the receiver operating characteristic (ROC)
curve was determined. The ROC is an estimate of predictive accuracy of logistic models that
is not biased by threshold probabilities and is the preferred measure when presence and true
absence data are available (Fielding & Bell, 1997). The ROC score is calculated by plotting
the occasions when turbidity exceedance is correctly predicted divided by the total number of
exceedances (sensitivity), against the fraction of incorrect cases where exceedance is
predicted (1-specicifity) across thresholds (Psyllakis, 2006). ROC values of 0.5 indicate that
the explanatory variables do not improve discrimination beyond random assignment and a
score of 1.0 indicates perfect discrimination (Manel, Williams, & Ormerod, 2001).

41
Chapter 5:

Results

In this Chapter, the HDOY for which incremented RX% are reached during the runoff
period and the TRV are summarized. A significant relationship between snow covered area
and spring runoff response exists in the KRB during the early stage of the runoff period. The
correlations between the global atmospheric circulation indices and incremented timing RX%
and TRV in the KR indicate that streamflow early on is significantly correlated with both the
ONI and SOI index. The PDO also demonstrated significant correlations with early runoff
response in the KR. The resulting correlation fields between the aggregate winter climate
inputs and the R25%, R50%, R75% and TRV demonstrated that regional surface climate is
significantly correlated with runoff response in the KR. Several multiple regression models
were identified that provide good predictive skill of runoff volume and timing, with inputs
shown to be significantly correlated with runoff timing. An example prediction is made
validating the forecast procedure. The results of the logistic regression analysis predicting
probability of exceeding 500 NTU at Arras as a function of discharge reveal that turbidity
exceedance can be predicted with moderate sensitivity.
5.1

Average Timing and Volume of Cumulative Streamflow
Descriptive statistics summarizing the minimum, mean, interquartile range (IQR) and

maximum HDOY for each incremented RX% and the TRV are found in Table 5.1. Based on
the IQR, the greatest variability in streamflow timing occurs over the R60%, with less
variability during early and late runoff. The distribution of annual RX% and TRV over time
since HY 1965 is seen in Figure 5.1.
5.2

Streamflow Timing and Snow Ablation
A significant correlation between the date that SCF reaches 50% (SCF50) and the

R25% was identified in the KRB (rs = 0.73, n = 16, p < 0.002).

42

Figure 5.1. Distribution of annual spring runoff timing and volume (25 March–1 July) in the
Kiskatinaw River since hydrologic year 1965 (calendar year 1966).The left axis identifies the
total runoff volume of water in cubic meters, whereas the right axis shows the timing of
cumulative runoff as a Julian day of year for the R5%, R25%, R50%, R75% and R95%
determined over the runoff period 25 March-1 July.

43

Table 5.1: Statistics describing the observed timing of cumulative runoff (RX%) and total
runoff volume (TRV) in the Kiskatinaw River, BC, defined at the Water Survey of Canada
Farmington stream gauge, 1965–2016. Timing is the hydrologic day of year (HDOY, days
since 1 September) that a given RX% is reached, and TRV is in cubic meters.
Statistic R5%
R10% R15% R20%
R25% R30% R35% R40%
Min
208.6
210.8
212.7
214.6
217.0
219.8
222.5
225.1
Mean
229.8
234.0
237.3
240.5
243.4
246.2
249.1
252.3
IQR
13.4
13.6
14.5
14.8
14.8
14.4
14.0
13.2
Max
244.3
247.6
252.0
258.0
267.2
270.5
274.6
282.1
Statistic R45% R50% R55% R60%
R65% R70% R75% R80%
Min
228.5
233.1
236.9
239.3
242.1
244.9
247.5
251.6
Mean
255.4
258.3
261.6
264.9
268.5
272.2
275.9
279.7
IQR
13.2
14.0
15.3
18.9
20.6
18.2
19.6
17.8
Max
287.1
288.4
289.1
289.7
290.3
293.5
296.0
299.9
Statistic R85% R90% R95% TRV
Min
251.6
257.7
266.7
4.91 × 107
Mean
279.7
283.3
288.0
2.16 × 108
IQR
17.8
14.8
13.0
1.86 × 108
Max
299.9
300.9
302.1
6.26 × 108

44
The mean SCF50 was HDOY 227, with an interquartile range of 17 days. The correlation
with early timing confirms the importance of input from snowmelt on streamflow response in
the KR, which is well captured by the reduction in snow covered fraction over the KRB.
5.3

Correlation with Global Circulation Indices
Correlations between the incremented RX% and TRV in the KR with the various

global circulation indices described in Chapter 3 are summarized in Table A.1. Generally, all
global circulation indices are more strongly correlated with early streamflow timing (R25%),
than later during the runoff period (R75%). The SOI and ONI indices are the most highly
correlated with streamflow timing in the KR. The SOI is positively correlated with the timing
of runoff and the ONI is negatively correlated with the timing of runoff in the KR.
The PDO is generally negatively correlated with the timing of runoff in the KR,
especially PDO-SON early on in the runoff season. It becomes slightly positively correlated
with the timing of runoff later in the season (R95%).
The NP index is correlated with the SOI and ONI indices, and generally positively
correlated with the timing of streamflow in the KR, but to a lesser extent than the ENSO
indices. There also exists a significant positive correlation between the NP over JJA and the
TRV in the KR (rs = 0.29, n = 50, p < 0.04). The SOI over MAM is also correlated with TRV
(rs = 0.24, n = 50, p < 0.1).
5.4

Spatial Correlations
For the various climate observations, significant regional Spearman rank correlations

with streamflow timing and volume exist. The various patterns identified are described in the
following sections.

45
5.4.1 Timing of spring runoff.
Mean winter temperature. Temperature in the BC Central Interior and southernmost
NE BC region displays strong negative correlations with the R25% (Figure B.1). Local
temperature in the KR is negatively correlated with the R25%. At the MOTI Braden Rd.
station within the KRB there is a strong negative correlation with temperature and the R25%
(rs = -0.74, n = 13, p < 0.004). The widespread negative correlation with temperature and the
R25% is the result of warmer temperatures leading to less snow accumulation and an earlier
onset of snowmelt, resulting in earlier runoff. Average winter temperatures in the low
elevation areas of the BC Northern Interior also display strong negative correlations with the
R25%. Temperature remains negatively correlated with the R50%, with the most significant
correlations again occurring in the BC Central Interior, but are generally weaker than
correlations with R25% (Figure B.2). Temperature is only weakly correlated with streamflow
timing mid to late in the runoff season (R75%). A region of slightly negative correlation
exists in southern BC. The low lying regions of NW BC have slight positive correlations with
the timing of the R75%. Interestingly, temperature becomes generally more highly correlated
with timing at the end of the runoff period (R95%) (Figure B.4). A region of negatively
correlated stations exists east of the Rockies in Alberta (AB), with stronger correlations in the
south near the Kootenays. A region of general positive correlations remains in NW BC, but is
enhanced relative to the R75%.
Total winter precipitation. The timing of early (R25%) streamflow in the KR is
generally less correlated with observed cumulative winter precipitation than temperature
(Figure B.5). Still, a trough of significant negative correlations is seen in SW BC, notably
MOTI Red Bluffs (rs = -0.67, n = 18, p < 0.003) and MOTI Squilax (rs = -0.69, n = 18, p <
0.002). The region of NW BC also has several stations with positive correlations. Later

46
during the runoff season (R50%), the same trough seen in the correlations with the R25%
exists in the Thompson-Okanagan region but has expanded into the BC Southern/Central
Interior and the South Coast region (Figure B.6). Cumulative winter precipitation is only
weakly correlated with the R75%, and the spatial distribution is convoluted. The stations
within NW and SW BC are generally negatively correlated with the R75%. While a
significant positive correlation exists in the Kootenay-Columbia region at ECCC
Kickinghorse (rs = 0.79, n = 7, p < 0.04), it is based on few observations. Stations in the BC
Central Interior exhibit modest positive correlations with the R75% (Figure B.7).
1 April snow water equivalent. Generally, SWE on 1 April is positively correlated
with the timing of early runoff in the KR (Figure B.8). In the BC Central Interior near
Williams Lake and Prince George, there are a few sites with significant positive correlations
with the R25%, and to a lesser degree the R50% in the KR (Figure B.9). Notably, the snow
survey site 1C33A (rs = 0.92, n = 9, p < 0.001). The BC Northeast, Central Interior and
locations in the southern half of the Coast Mountains have significant negative correlations
with the R50%; i.e., site 4A21 (rs = -0.40, n = 38, p < 0.02) in the NE and 1B06 in the
Central Interior (rs = -0.46, n = 26, p < 0.02). Snow water equivalent on 1 April is primarily
no longer positively correlated with runoff timing later during the runoff period (R75%).
Regions of moderate negative correlation with the R75% exist over the southern half of the
Coast Mountains (Figure B.10). Regions of moderate negative correlation also occur in the
BC Central Interior, and NE BC, including the KRB.
Atmospheric surface pressure. Average surface pressure during the accumulation
period at and near ECCC West Twin Creek, BC, is negatively correlated (rs = -0.59, n = 12,
p < 0.05) with the R25% (Figure B.11). Atmospheric surface pressure at stations within NE
BC is positively correlated with the R25% in the KR. This pattern remains correlated with

47
the R50%, but weakens slightly in both regions (Figure B.12). Stations within the SE corner
of BC up into the Kootenay and southern Columbia Mountains are negatively correlated with
the R75%. Generally, in NE BC, the correlation with streamflow timing in the KR becomes
slightly negative, as with the R75% (Figure B.13).
5.4.2 Total runoff volume.
Mean winter temperature. Average winter temperature over three prominent regions
is correlated with TRV in the KR (Figure B.14). Winter temperatures east of the continental
divide in the prairie regions of AB and Saskatchewan (SK) into the Northwest Territories
(NT) are negatively correlated with runoff volume in the KR, notably, ECCC High Level (rs
= -0.40, n = 43, p < 0.008), and ECCC Craigmyle (rs = -0.39, n = 41, p < 0.013). Generally,
stations in the NE region of BC and SW region of the Yukon (YT) have no coherent
correlation with TRV in the KR.
Total winter precipitation. Correlations between TRV and total winter precipitation
are less spatially coherent than temperature (Figure B.15). Total precipitation is generally
positively correlated with TRV. Stations near and local to the KR have significant positive
correlations with TRV, notably the MOTI 73 Mile station (rs = 0.57, n = 13, p < 0.05), and
ECCC Chetwynd (rs = 0.51, n = 31, p < 0.004). Stations within the Central Interior of BC
have significant positive correlations with the TRV, the highest being at MOTI Lee’s Hill
station (rs = 0.54, n = 18, p < 0.02). Stations in the Central Plains of AB reveal strong
positive correlations with TRV, i.e., ECCC Entwistle (rs = 0.50, n = 20, p < 0.03), ECCC
Cold Lake (rs = 0.47, n = 48, p < 0.001), and ECCC Craigmyle (rs = 0.49, n = 41, p <
0.002). The correlations here indicate a response to the strong influence of El Niño on winter
precipitation in this region (Figure C.2).

48
Snow water equivalent. Generally, snow water equivalent on 1 April is positively
correlated with TRV over BC (Figure B.16). Stations in the BC Central Interior have strong
positive correlations with the TRV, with the highest occurring at site 1A23 (rs = 0.52, n = 23,
p < 0.02). Stations occurring near the KRB also have strong positive correlations with the
TRV, as more snow accumulation leads to greater runoff volumes in response to melt
including site 4A25 (rs = 0.42, n = 40, p < 0.008). Stations in the BC Southern Interior and
SW BC have moderately strong positive correlations with TRV, notably site 1C23 (rs = 0.40,
n = 39, p < 0.02). Surveys in the Cariboo and Kootenay Mountains reveal positive
correlations with TRV, but surveys are sparse. The 2C15 survey west of Calgary has the
highest correlation with TRV (rs = 0.50, n = 44, p < 0.001) indicating that more surveys in
this region should be included as both precipitation and temperature are highly influenced by
El Niño in this region (Figure C.1 and C.2).
Atmospheric surface pressure. Atmospheric surface pressure is generally positively
correlated with TRV (Figure B.17). There is a general gradient of weaker correlations in
southern BC, and stronger correlations in the BC Central and Northern Interior. The positive
correlations are likely the result of higher pressures during winter correlating with colder
temperatures, and slightly more precipitation in the KRB. Mean pressure at Braden Rd. also
has a significant positive correlation with the TRV (rs = 0.56, n = 13, p < 0.05), as with
stations nearby. Stations in the Thompson-Okanagan region also display strong positive
correlations with TRV.
5.5

Runoff Timing
The winter climate inputs and global circulation indices were shown to provide good

predictions of streamflow response in the KR. Prediction of early runoff response was

49
typically more accurate than prediction of later runoff response, and this was also indicated
by the erosion of the contour fields of the Spearman rank correlations with the climate
observations from the R25% to the R75%. Many of the climate inputs in the selected
regression models were located in the regions shown to display coherent correlations with
runoff timing based on visual interpretation of the contoured Spearman rank correlations.
Multiple regression results for streamflow timing are presented for the three primary stages
of runoff: early, mid and late season(i.e., R25%, R50% and R75%). The detailed regression
models for each incremented RX% and TRV can be found in Appendix D.
Early runoff. The model selected for predicting the R25% in the KR included total
cumulative winter precipitation (CPPT) at the MOTI Red Bluffs climate station, mean winter
temperature (MT) at the MOTI Bella Coola station, and the SOI index over JJA (Table D.5).
Some evidence of bias in the regression coefficients were identified (Table 5.2). The SOI
index over JJA shows a slight upward bias, about 8% of the original coefficient. The
combination of these three aggregated winter climate inputs provided good predictability of
the R25% (adj-r2= 0.94, F(3,16) = 58.37, p < 0.001). No major violations regarding
assumptions of linear regression were identified based on regression diagnostics described in
Chapter 4. Hydrologic years 1997, 1998 and 2015 were identified as influential; however,
1997 was the highest discharge on record for the KRB. Year 2015 has the largest
interquartile range (IQR), and 1998 also has a substantial historical IQR (Figure 5.1).
Mid-season runoff. The optimal model identified for predicting the R50% in the KR
included CPPT at the MOTI Red Bluffs climate station, MT at the ECCC Chetwynd station,
and MT at the ECCC Smithers station (Table D.10). Little evidence of considerable (Percent
𝛽̂ > 5) bias in the regression results was identified (Table 5.2). The combination of these
three aggregated winter climate observations provided the most accurate predictions of

50
Table 5.2: Table showing slope coefficients (β̂) and the jackknifed unbiased estimates of the
slope coefficients (β̂j ) in estimating the hydrologic day of year that cumulative discharge
reaches 25%, 50% and 75% of total runoff volume in the Kiskatinaw River. Bias is provided,
and also shown as an absolute percentage of the biased coefficient. MT stands for mean
temperature, CPPT stands for cumulative precipitation, and numbers below reference agency
stations number. Agency MOTI is Ministry of Transportation and Infrastructure, EC is
Environment and Climate Change Canada.
̂𝒋
̂
̂
R25%
𝜷
Bias
𝜷
Percent 𝜷
CPPT @ MOTI
-0.1223
-0.1232
0.0009
0.7185
(24225)
MT @ MOTI
-5.9987
-6.0578
0.0591
0.9848
(47121)
Avg. SOI-5.0373
-4.6149
-0.4224
8.3854
Summer
Const.

249.6216

249.4138

0.2078

0.0832

R50%
CPPT @ MOTI
(24225)
MT @ EC
(1181508)
MT @ EC
(1077500)

̂
𝜷

̂ 𝒋.
𝜷

Bias

̂
Percent 𝜷

-0.1148

-0.1171

0.0023

1.9835

5.5404

5.3789

0.1615

2.9155

-11.8854

-12.1634

0.2780

2.3387

Const.

270.9989

269.3354

1.6635

0.6138

R75%
MSD @ MOTI
(47091)
1 April SWE @
(1E03A)
1 April SWE @
(4C05)

̂
𝜷

̂ 𝒋.
𝜷

Bias

̂
Percent 𝜷

0.4422

0.4391

0.0031

0.7083

-0.0399

-0.0444

0.0045

11.2905

-0.2217

-0.2214

-0.0003

0.1509

Const.

313.6378

315.8872

-2.2494

0.7172

51
the R50% when input into a multiple linear regression framework (adj-r2 = 0.86, F(3,15) =
44.66, p < 0.001). Cumulative precipitation at MOTI Red Bluffs showed up as input in many
of the candidate regressions not selected. No major violations regarding assumptions of
linear regression were identified based on regression diagnostics described in Chapter 4.
Hydrologic years 2007 and 2010 were identified as influential; both years had an abnormally
tardy R50% (Figure 5.1).
Late runoff. The model identified for predicting the R75% in the KR includes snow
depth averaged for 24–25 March (MSD) at the MOTI Lee’s Hill station, 1 April SWE at
Trophy Mountain, and 1 April SWE at Fort Nelson Airport (Table D.15). Some evidence of
bias in the regression results were identified (Table 5.2). Notably, 1 April SWE at Fort
Nelson shows a moderate downward bias, about 11% of the original coefficient. The
combination of these three inputs were able to predict the R75% well, but with less
uncertainty then the timing of runoff early during the runoff period (adj-r2 = 0.82, F(3,15) =
30.38, p < 0.001). No major violations regarding assumptions of linear regression were
identified based on regression diagnostics described in Chapter 4. Hydrologic year 2010 was
identified as an outlier and highly influential observation, as the R75% was quite late that
year (Figure 5.1).
5.6

Runoff Volume
Streamflow volume (ha m) in the KR is predicted by CPPT at ECCC Cold Lake

Airport station, MT at ECCC Chetwynd Airport station, MT at ECCC Pleasant Camp station,
and CPPT at ECCC Mayo Rd. (Table D.20). The Mayo Rd. shows slight upward bias (Table
5.5). A high adj- r2 was observed suggesting good predictive skill (adj-r2 = 0.89, F(4,23) =
33.64, p < 0.001). Years 1989, 1994, 1996, 2005, 2006 and 2010 were identified as
influential observations. Year 1996 had the highest recorded TRV in the KRB, 1989 and

52
Table 5.3: Table showing slope coefficients (β̂) and the jackknifed unbiased estimates of the
slope coefficients (β̂j ) in estimating the total runoff volume (ha m) in the Kiskatinaw River.
Bias is provided, and also shown as an absolute percentage of the original coefficient. MT
stands for mean temperature, CPPT stands for cumulative precipitation, and numbers below
reference agency stations number. Agency MOTI is Ministry of Transportation and
Infrastructure, EC is Environment and Climate Change Canada.
̂ 𝒋.
̂
̂
TRV (ha m)
𝜷
Bias
𝜷
Percent 𝜷
CPPT @ EC
540.7233
562.4713
-21.7480
0.0402
(3081680)
MT @ EC
-7039.12
-7088.6034
49.4834
0.0070
(1181508)
MT @ EC
5820.866
5825.4768
-4.6108
0.0008
(1206197)
CPPT @ EC
-156.0408
-171.5196
15.4788
0.0992
(2100709)
Const.

-28751.83

-29286.2270

534.3970

0.0186

53
2010 were also considerably high. Year 2005 had one of the lowest TRV in the KRB (Figure
5.1).
5.7

Model Validation
An example forecast for hydrologic year 2002 is provided in this section as validation

of the resulting regression models’ ability to predict a continuous cumulative discharge
curve. Because the different data inputs specified in Chapter 2 are not always available for a
given year, particularly before 2000, there are only a number of years with all inputs for
timing and volume present. Hydrologic year 2002 happens to have all required inputs
available that year. Hydrologic year 2002 is a relatively ‘normal’ year in regards to TRV
(18438.62 ha m) and the timing of the R50% (HDOY 251). Cumulative discharge for HY
2002 in the KR was predicted from the resulting regression models (Figure 5.2). The
predicted cumulative curve has been smoothed by a 10-day moving average filter; the
cumulative curve was differenced daily to compute an average daily discharge (Figure 5.2).
From the smoothed discharge, the Nash-Sutcliffe Efficiency Coefficient (NSE) was
calculated. The resulting NSE of 0.70 indicates good predictive ability.

54

Figure 5.2: The observed (black solid line) and predicted (dashed maroon line) cumulative
discharge at Water Survey of Canada Farmington for hydrologic year 2002, over the runoff
period 25 March (hydrologic day of year 207) to 1 (HDOY 304) July. The original predicted
cumulative values are represented by the teal X’s. A 10-day moving average is applied to the
predictions. The various dots represent the interquartile range of each model’s residuals,
added to the prediction, in different combinations of timing t, and volume v. For example, the
P25t is the predicted timing plus the 25th percentile of that model’s residuals. The P75v is the
predicted volume plus the 75th percentile of the residuals for the total runoff volume model.

55

Figure 5.3. Predicted hydrograph at the Water Survey of Canada Farmington stream gauge in
the Kiskatinaw River for hydrologic year 2002. The blue line is the predicted discharge and
the black line is the observed discharge (linearly interpolated from 5% cumulative intervals).
See Figure 5.1 for legend description. The Nash-Sutcliffe Efficiency Coefficient was
determined to be 0.70. The prediction is made on 1 April, and a 10-day moving average was
applied to the prediction (Figure 5.2). The abrupt changes in discharge are the result of
negative slopes in the predicted cumulative discharge, and can be filtered out.

56
5.8

Predicting Turbidity Exceedance
The rate of wash-load sediment transport is determined primarily by the rate of fine-

grained sediment supply from the drainage basin rather than the transporting capacity of the
flow. High suspended sediment concentrations are generally associated with periods of high
discharge, and empirical relationships between the two are often established, however
(Walling & Webb, 1992). To provide some indication of suspended sediment concentrations
associated with springtime runoff in the KR, logistic regression was utilized to provide
probabilities of exceeding 500 NTU (The NTU limit at which pumps must be shut down) at
Arras, as a function of average daily discharge. Average daily turbidity at Arras as a function
of average daily discharge at the WSC Farmington stream gauge for years 2007 and later
during the spring runoff period 25 March–30 June were plotted to determine the relationship
between these two variables in the KR. Generally turbidity increases with daily average
discharge (Figure 5.4).
Discharge at the WSC Farmington stream gauge is found to be a significant predictor
of turbidity exceedance of 500 NTU at Arras (χ2(1) = 254.46, n = 939, p < 0.001) (Table 5.4).
The model has moderate sensitivity, only predicting 45.5% of the observed exceedances;
however, it exhibits robust specificity by correctly predicting 97.7% of observations below
500 NTU at a threshold probability of 0.5. Overall, 91.6% of observations were correctly
classified (Table 5.5). The area below the ROC curve is 0.91, indicating that the model has
good discrimination ability (Manel, Williams, & Ormerod, 2001). However, because of the
model’s moderate sensitivity, other inputs should be explored to improve the ability to
predict exceedances of 500 NTU.

57

Figure 5.4. Daily average turbidity measured at the Arras weir as a function of daily average
discharge at Water Survey of Canada Farmington stream gauge since 2007 during spring
runoff periods, 25 March–1 July. The horizontal blue line marks the 500 nephelometric
turbidity units (NTU) threshold value at which water abstraction must cease.

58

Table 5.4: Logistic regression results predicting the probability of exceeding 500 NTU at
Arras from daily average discharge at the Water Survey of Canada Farmington stream gauge,
on the Kiskatinaw River. Model based on data since 2007 over the spring runoff period 25
March-1 July.
P(Excd. 500NTU) β
SE
z
P>|z|
95% Conf. Interval
Discharge (m3s-1)
0.0605
0.0050
12.04
<0.001
0.0507
0.0704
Const.
-4.1715 0.2494
-16.73
<0.001
-4.6603
-3.6827
Table 5.5: Classification table showing results for logistic regression predicting turbidity
exceedance at the Arras weir in the Kiskatinaw River. A threshold of 500 Nephelometric
Turbidity Units (NTU) was used as the classification threshold.
Classified
Obs. NTU≥500
Obs. NTU<500
Total
Pr(NTU≥500)≥0.5
50
19
69
Pr(NTU≥500)≤0.5
60
810
870
Total
110
829
939

59
Chapter 6:

Discussion

In this chapter relationships between synoptic circulation patterns and teleconnections
with streamflow timing and volume in the KR are discussed. The difference between the
causal and predictive modelling approach applied in this study is emphasized. The
application of this study as a decision support system is explored. Finally, areas of future
research are mentioned.
6.1

Teleconnections
Of the global circulation indices considered in this study, the ENSO is the most

highly correlated with variation in early runoff timing in the KR (Table A.1). The significant
positive (negative) correlations observed between the timing of streamflow in the KR,
particularly early on, and the SOI (ONI) index, respectively, indicate delayed flows in
response to more snow accumulation, less mid-winter melt, and gradually increasing spring
temperatures associated with La Niña conditions in the KR. The fall component of the SOI
and winter component of the ONI index are highly correlated with streamflow timing, as this
is when the ENSO signal is most prominent (Woo & Thorne, 2003).
The ENSO has a defined influence on temperature fields over Canada. An
intensification of the Canadian Ridge and a deepening of the Aleutian Low lead to strong
southwesterly flow along the west coast of Canada producing positive surface temperature
anomalies in BC, and tropospheric patterns consistent with a positive PNA pattern (Shabbar
& Khandekar, 1996). Conversely, during La Niña, abnormally low geopotential height
thicknesses over the Canadian Rockies and a positive anomaly centre over the NE Pacific
Ocean lead to northwesterly advection of cold air over western Canada. The result is a strong
influence on snowmelt-dominated rivers in western Canada. The PDO has been shown to
have a possible enhancing (attenuating) effect on ENSO depending on whether their phases

60
coincide (Gershunov & Barnett, 1998), and also explains much streamflow variation in the
KR. The positive PDO phase is also highly correlated with early streamflow timing,
particularly early on (Table A.1).
Of the surface observations, temperature was the most highly correlated with
streamflow timing in the KR (Figure B.1). It would be expected that prior winter
temperatures would yield strong influence on the timing of streamflow as snow ablation
magnitude and timing is largely a function of temperature. The regions negatively correlated
with the R25% and R50% in the KR are consistent with those that express positive
temperature anomalies in Figure 4 of Shabbar (2006), and Figure 2 of Fleming & Whitfield
(2010) in response to strong postive ENSO events (Figure C.1). In addition, pressure
anomalies correlated with timing in the KR correspond regionally with positive anomalies in
the 500-1000 hPa thickness field associated with El Niño found in Shabbar & Khandekar
(1996) (Figure B.12).
The ENSO is also linked to anomalies in precipitation fields over Canada.
Precipitation in southern Canada has a significant response to ENSO events, and is most
pronounced during winter (Figure C.2). These patterns are linked to anomalies in the 500 hPa
geopotential height circulation, which resemble the positive (negative) phase of the PNA
pattern that has been linked to ENSO events (Shabbar, Bonsal, & Khandekar, 1997).
Negative precipitation anomalies over the KRB in response to ENSO are seen in Shabbar
(2006) (Figure C.2). The correlation between the timing of streamflow in the KR and global
circulation patterns is less consistent for precipitation than temperature, particularly for
indices specific to ENSO. Precipitation at several SW BC stations was found to have
significant negative correlations with the R50% (Figure B.5). Past studies have shown that
typically, precipitation over southern BC into the Thompson-Okanagan and the Lower

61
Mainland has expressed negative anomalies during positive ENSO conditions (e.g., Fleming
& Whitfield, 2010; Shabbar, Bonsal, & Khandekar, 1997). In more recent studies, a small but
significant location of positive precipitation anomaly focused over Vancouver Island
eastward into the Coast Mountains, in response to El Niño has been identified (Figure C.2).
Generally, increases in precipitation near this region show negative correlations with timing
of streamflow in the KR (Figure B.5). There are stations shown to have significant positive
Spearman rank correlation with ONI over DJF and MAM occurring in this region and
include the ECCC Whistler station (rs = 0.55, n = 25, p < 0.06) and the MOTI Red Bluffs
station (rs = 0.55, n = 19, p < 0.02), respectively. The MOTI Red Bluffs station is
significantly correlated with the ECCC Whistler station (rs = 0.71, n = 12, p < 0.01). In the
resulting regression models, precipitation in this region, particularly MOTI Red Bluffs,
proves to be an important predictor of streamflow timing in the KR (see Appendix D). For
instance, to predict the R25%, precipitation at Red Bluffs along with mean temperature at
MOTI Bella Coola and the SOI over JJA are used as input. The models predicting early
streamflow response typically include some variable highly correlated with ENSO (such as
Red Bluffs), with the other variables that capture the remaining variability not associated
with ENSO. This study supports the finding that the response of surface meteorological
regimes to ENSO vary considerably over BC (Fleming & Whitfield, 2010), and play an
important role in annual hydrologic variability.
The NP index averaged over the previous winter (JJA) has a significant positive
correlation with TRV in the KR (Table A.1). This is in part because an enhanced PNA signal,
especially during winter, results in widespread reductions in snow accumulation in BC
(Moore & McKendry, 1996). Winter NP was also found to be positively correlated with total
precipitation during the spring runoff period in the KR based on the local MOTI Braden Rd.

62
station (rs = 0.42, n = 14, p < 0.2). Summer PNA conditions provide information about both
snowpack conditions come spring, in addition to summer convective patterns. This would
explain the strong correlation with TRV in the KR. The PDO, particularly during winter, is
also highly correlated with spring and summer precipitation at the MOTI Braden Rd. station
(rs = -0.56, n = 14, p < 0.04), and also frequently occurs as input to models not selected, but
that provided good predictability of TRV in the KR.
6.2

Predictive Modelling vs. Causal Attribution
Statistical modelling has traditionally been applied for determining causal explanation

in the sciences. One of the consequences of failing to include predictive modelling (data
mining) as an exploratory approach is losing the ability to test the relevance of existing
theories and to discover new causal mechanisms (Shmueli, 2010). Regression is a powerful
tool for testing causal hypotheses. The process of modelling to attribute causation and the
process of modelling to provide prediction are quite different. One difference between the
two approaches is that predictive modelling is forward-looking, whereas explanatory
modelling is retrospective as it is used to test an established hypothesis (Shmueli, 2010). As
there are numerous regression models in this study, it is important to note that the approach
taken in this study is in essence an exercise in data mining, and causation should not
necessarily be attributed. Take the regression model determined for the R50% for instance
(Table D.10). This model suggests that for every millimeter increase in total winter
precipitation at MOTI Red Bluffs, the R50% decreases by -0.12 days, if the other
independent variables remain constant. To imply that precipitation at Red Bluffs is causing
the R50% to be earlier in the KR would certainly be false.
Precipitation at Red Bluffs is significantly correlated with larger synoptic processes
that are driving hydrologic inputs in the KR during spring runoff, however, and therefore

63
provides good predictability of the R50%. Understanding the underlying physical processes
of these relationships will be an important step in reducing uncertainty around predictions,
and furthering knowledge of hydro-meteorology in the KR, however.
There is much to be lost if predictive modelling is omitted from scientific research,
however, and that predictive modelling allows us to learn from data, verifying past theory on
new data (Feelders, 2002). That said, I emphasise the modelling approach in this study, and
not so much the models themselves.
6.3

Decision Support Systems
Computer-based models combined with interactive interfaces are typically known as

a decision support system (DSS). The overall objective of a DSS is to provide timely
information that supports decision makers (Loucks, 1995). Water-resource systems constitute
many branches, including environmental, hydrological, social and administrative (Alemu,
Palmer, Polebitski, & Meaker, 2011). To be considered practical and of utility, a DSS must
provide required information in a digestible form to those making the decisions. Time is also
especially important, as decision makers require information when the window of
opportunity to use that information exists. A DSS consists of two major components: a
mathematical model to generate the data, and some visual medium to convey the results
(Loucks & Beek, 2005). Although this study explored possible atmospheric mechanisms
responsible for variation in runoff in the KR, providing the CDC a practical tool allowing
timely forecast of spring runoff in the KR was a key objective and influenced the modelling
approach taken in this study. As an extension to this research, a DSS has been provided to the
CDC, providing a tool for more informed water management based on available climate data.
This DSS provides both predictions of discharge at Farmington and the probability of

64
exceeding 500 NTU at Arras, along with a number of programs to aide in data assimilation
and aggregation.
6.4

Future Research and Recommendations
There are several areas of opportunity for future research as an extension to the

findings in this study. This effort provides evidence that aggregated winter climate metrics
(i.e., temperature, precipitation, pressure, SWE, snow cover / depth) in combination with
global circulation indices can offer accurate predictions of streamflow timing in the KR.
Studies have shown that non-standard indices of global circulation patterns can improve
forecast of winter snowpack and streamflow in western North America (Grantz, Rajagopalan,
Clark, & Zagona, 2005). The aggregated surface climate inputs seem to be achieving this
goal as they are manifestations of anomalies in synoptic weather patterns. Using clusters of
synoptic weather patterns derived from daily mean sea-level pressure grids, the frequency of
synoptic weather pattern types could also be provided as possible model input (Stahl, Moore,
& McKendry, 2006). These synoptic pattern frequencies could possibly improve predictions
beyond using surface observations alone, or in combination with surface observations.
Although streamflow variation was successfully correlated with large-scale climate
data in this study, local processes certainly influence the water balance of the KR. Therefore,
major landscape changes (e.g., from forest fires) would likely result in significant increases
in predictive uncertainty. Because the regression models are empirical, they are not robust to
major landscape changes at the basin scale. Therefore, in the case of a major disturbance,
research will have to be carried out as to what the effects are on the hydrograph, and how to
interpret the predicted hydrograph in the context of these changes.
As many of the regression models identified typically had relatively small sample
sizes, ~n ≤ 15, updating of the models on new data will be important and should lead to less

65
uncertainty in predictions. In addition, possible new model structures may be identified, and
the nature of any non-stationarities affecting predictions better understood.
Although discharge alone provides moderate predictive skill of turbidity exceedance
of 500 NTU at Arras, many exceedances were not predicted. The moderate sensitivity is
likely because the sediment supply from processes not associated with the carrying capacity
of the flow is not being captured, as these are associated with events such as rainfall and
upstream land use activities. This makes turbidity particularly difficult to predict. However,
adding variables such as maximum 24-hr rainfall and upstream turbidity would likely provide
improved predictions of turbidity up to several weeks ahead of time.
This study’s methods used to forecasting streamflow in response to snowmelt in the
KR could easily be applied to another basin. The only requirements for implementing the
methods defined in this study are: (a) there must be an adequate record of non-regulated
discharge data available in the basin of interest; and (b) the basin must be within synoptic
reach of an adequately dense weather and snow survey network. There is therefore the
opportunity to apply the methods in this study to other basins and compare performance of
the modelling approach when applied to other watersheds.

66
Chapter 7:

Conclusions

There is evidence in this study that climate observations aggregated during the snow
accumulation period within synoptic range of the KRB including surveys of SWE in
combination with large-scale atmospheric circulation indices, can provide accurate
predictions of cumulative spring runoff in the KR. In this study, incremented timing of
cumulative streamflow (RX%) during spring runoff in the KR was predicted using multiple
regression. The day of year that cumulative discharge reaches 50% (R50%) can be predicted
by a combination of mean winter temperature and total winter precipitation at specific
locations with good accuracy (adj-r2 = 0.86, F(3,15) = 44.66, p < 0.001). The optimal
regression model selected for predicting total runoff volume also relied on mean winter
temperature and total winter precipitation as inputs, and provided good accuracy (adj-r2 =
0.89, F(4,23) = 33.64, p < 0.001). The predicted discharge for (HY) 2002 demonstrated a
Nash-Sutcliffe Efficiency Coefficient of 0.70.
The probability of turbidity exceeding 500 NTU at the Arras weir was predicted using
logistic regression. Implementing average daily discharge at WSC Farmington as a predictor,
a ROC score of 0.91 was obtained, indicating good discrimination ability. However, only
moderate sensitivity was achieved, meriting the exploration of possible other inputs that
might describe basin processes responsible for suspended sediment supply not captured by
discharge alone.
The resulting regression models are supported by evidence of correlations between
the aggregated climate metrics at the surface, known global circulation indices, and
streamflow metrics in the KR. There is compelling evidence that they serve as non-standard
indices of variation in synoptic weather patterns during the winter into spring and early
summer, driving much of the observed variation in streamflow in the KR over the last 50

67
years. It will be important to update and test the models identified here on new data as they
become available, particularly as many of the sample sizes are currently relatively small. This
will lead to less uncertainty in the predictions. Additionally, the exploration of other, possibly
non-linear model structures may be implemented as more is understood about the underlying
physics responsible for the regional correlations identified between winter surface climate
and runoff response in the KR.
Ideally streamflow predictions provided from this work will aid the CDC in
optimizing the timing of water releases at Bearhole Lake and downstream water withdrawals
from the Arras weir so that it can be abstracted during periods of high flow and low turbidity,
allowing adequate storage reserves during the low flow period in late summer and fall. A
decision support system (DSS) can be easily implemented on any desktop computer. Benefits
of the modelling approach include the ability to easily update the regression models as
additional input data become available over time likely improving predictions, adapting the
model to changes in climate, and improving estimates of predictive uncertainty. Because the
only requirements for implementing the modelling approach in this study are an adequate
record of discharge data, and a relatively dense regional network of climate stations and snow
survey locations, the model can be applied to other basins, an opportunity for future research.

FORECASTING SPRING RUNOFF IN THE KISKATINAW RIVER
References
Alemu, E. T., Palmer, R. N., Polebitski, A., & Meaker, B. (2011). Decision support system
for optimizing reservoir operations using ensemble streamflow predictions. Water
Resources Planning and Management, 137(1), 72-82.
American Water Works Association. (2005). Standard methods for the examination of water
and waste-water (21st ed.). Washington, D.C.: American Public Health Association.
Anderson, C. W. (2005). 6.7 Turbidity (Version 2.1, 9/2005). In Field measurements: U.S.
Geological Survey techniques of water-resources investigations (pp. 5-52). U.S.
Geological Survey.
Angino, E. E., & O'Brien, W. J. (1967). Effects of suspended material on water quality.
Symposium on Geochemistry, Precipitation, Evaporation, Soil Moisture, Hydrometry,
General Assembly of Bern (pp. 78, 120-128). International Association of Science
Hydrology.
Apache Commons. (2017). Statistics. Retrieved from Math:
http://commons.apache.org/proper/commons-math/userguide/stat.html
Australian National University. (2017). ANUSPLIN Vrsn 4.4. Retrieved from
http://fennerschool.anu.edu.au/research/products/anusplin-vrsn-44
Barnett, T. P., Adam, J. C., & Lettenmaier, D. P. (2005). Potential impacts of a warming
climate on water availability in snow-dominated regions. Nature, 438(17), 303-309.
BC Ministry of Transportation and Infrastructure. (2017). Avalanche and Weather Programs.
Retrieved from https://prdoas3.pub-apps.th.gov.bc.ca/saw-paws/weatherstation
Belsley, D. A., Kuh, E., & Roy, E. W. (1980). Regression diagnostics: Identifying influential
data and sources of collinearity. Hoboken: John Wiley & Sons, Inc.
Beven, K. J. (2012). Rainfall-Runoff Modelling. Chichester, UK: John Wiley & Sons, Ltd.

68

69
Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and random
coefficient variation. Econometrica, 47(5), 1287-1294.
Committee on the Status of Endangered Wildlife in Canada. (2002). COSEWIC assessment
and update status report on the woodland caribou Rangifer tarandus caribou in
Canada. Ottawa, ON: Committee on the Status of Endangered Wildlife in Canada.
Dobson Engineering Ltd. & Urban Systems. (2003). Kiskatinaw River watershed
management plan. Dawson Creek.
Dunne, T., & Leopold, L. B. (1998). Water in Environmental Planning. New York: W.H.
Freeman and Company.
Environment and Climate Change Canada. (2017a). Past Weather and Climate. Retrieved
from Government of Canada:
http://climate.weather.gc.ca/historical_data/search_historic_data_e.html
Environment and Climate Change Canada. (2017b). Water Level and Flow. Retrieved from
Government of Canada: https://wateroffice.ec.gc.ca/search/historical_e.html
Feelders, A. (2002). Data mining in economic science: In dealing with the data flood.
Netherlands: STT/Beweton, The Hague.
Fielding, A. H., & Bell, J. F. (1997). A review of methods for the assessment of prediction
errors in conservation presence / absence models. Environmental Conservation, 24,
38-49.
Fleming, S. W., & Whitfield, P. H. (2010). Spatiotemporal mapping of ENSO and PDO
surface meteorological signals in British Columbia, Yukon, and Southeast Alaska.
Atmosphere-Ocean, 48(2), 122-131.
Forest Practices Board. (2011). Cumulative Effects Assessment: A Case Study for the
Kiskatinaw River Watershed Special Report. Victoria: Forest Practices Board.

70
Garbrecht, J., & Piechota, T. C. (2006). Climate Variations, Climate Change, and Water
Resources Engineering. American Society of Civil Engineers.
Gershunov, A., & Barnett, T. P. (1998). Interdecadal modulation of ENSO teleconnections.
Bulletin of the American Meteorological Society, 79, 2715-2725.
Gobena, A. K., Weber, F. A., & Fleming, W. S. (2013). The role of large-scale climate
modes in regional streamflow variability and implications for water supply
forecasting: A case study of the Canadian Columbia river basin. Atmosphere-Ocean,
51(4), 380-391.
Gotelli, N. J., & Ellison, A. M. (2013). A primer of ecological statistics. Sunderland, MA:
Sinauer Associates, Inc.
Government of British Columbia. (2018). Water Licences Query. Retrieved from
http://a100.gov.bc.ca/pub/wtrwhse/water_licences.output?p_Source_Name=Kiskatina
w+River&p_Licence_No=&p_Priority_Issue_Date=&p_POD_Purpose=&chk_Appur
tenant_Land=&p_POD_Qty_Equality=%3D&p_POD_Qty=&chk_Licence_Comment
s=&chk_POD_Qty_Flag_Desc=&chk_Date_Update
Grantz, K., Rajagopalan, B., Clark, M., & Zagona, E. (2005). A technique for incorporating
large-scale climate information in basin-scale ensemble streamflow forecast. Water
Resources Research, 41, W10410.
Hall, D. K., & Riggs, G. A. (2016). MODIS/Terra Snow Cover 8-Day L3 Global 0.05Deg
CMG, Version 6. Retrieved from National Snow & Ice Data Center:
http://dx.doi.org/10.5067/MODIS/MOD10C2.006.
Huang, B., Banzon, V. F., Freeman, E., Lawrimore, J., Liu, W., Peterson, T. C., . . . Zhang,
H. M. (2014). Extended reconstructed sea surface temperature version 4 (ERSST.v4).
Journal of Climate, 28, 931-951.

71
Islam, S. U., & Déry, S. J. (2017). Future climate change impacts on snow and water
resources of the Fraser River Basin, British Columbia. Journal of Hydrometeorology,
18, 474-496.
Kang, D. H., Gao, H., Shi, X., Islam, S. U., & Déry, S. J. (2016). Impacts of a rapidly
declining mountain snowpack on streamflow timing in Canada’s Fraser River Basin.
Scientific Reports, 6, 19299.
Koster, R. D., Betts, A. K., Dirmeyer, P. A., Bierkens, M., Bennett, K. E., Déry, S. J., . . .
Yuan, X. (2017). Hydroclimatic variability and predictability: A survey of recent
research. Hydrology and Earth System Sciences, 21, 3777-3798.
Lillesand, T. M., Kiefer, R. W., & Chipman, J. W. (2007). In Remote Sensing and Image
Interpretation (pp. 340-356). Hoboken: John Wiley & Sons, Inc.
Loucks, D. P. (1995). Developing and implementing decision support systems: A critique
and challenge. Water Resources Bulletin, 31, 572-582.
Loucks, D. P., & Beek, E. V. (2005). Model sensitivity and uncertainty analysis. In Water
resources systems planning and management (pp. 255-287). Italy: United Nations
Educational, Scientific and Cultural Organization.
Mahanama, S., Livneh, B., Koster, R., Lettenmaier, D., & Reichle, R. (2012). Soil moisture,
snow, and seasonal streamflow forecast in the United States. Journal of
Hydrometeorology, 13, 189-203.
Manel, S., Williams, C. H., & Ormerod, S. J. (2001). Evaluating presence-absence models in
ecology: the need to account for prevalence. Journal of Applied Ecology, 38, 921931.
Meidinger, D., & Pojar, J. (2017). Ecosystems of British Columbia. Retrieved 08 07, 2017,
from https://www.for.gov.bc.ca/hfd/pubs/docs/bro/BRO49.PDF

72
Ministry of Environment. (2017). Manual Snow Survey Observations Data Archive.
Retrieved from https://catalogue.data.gov.bc.ca/dataset/manual-snow-surveyobservations-data-archive
Ministry of Forests, Lands, Natural Resource Operations, and Rural Development. (2017).
Snow Survey Data. Retrieved from River Forecast Centre:
https://www2.gov.bc.ca/gov/content/environment/air-land-water/water/water-sciencedata/water-data-tools/snow-survey-data
Minobe, S. (1997). A 50-70 year climatic oscillation over the North Pacific and North
America. Geophysical Research Lettters, 24, 683-686.
Montanari, A., & Koutsoyiannis, D. (2014). Modeling and mitigating natural hazards:
Stationarity is immortal! Water Resources Research, 50(12), 9748-9756.
Moore, R. D., & McKendry, I. G. (1996). Spring snowpack anomaly patterns and winter
climatic variability, British Columbia, Canada. Water Resources Research, 32(3),
623-632.
Mottishaw, C. (2016). Kiskatinaw River 2015/2016 Update; Watershed Stewardship
Program. Retrieved from Dawson Creek, British Columbia:
http://www.dawsoncreek.ca/wordpress/wp-content/uploads/background-watershedmanagement-plans/WSP-Council-Update_April-2016.pdf
Murdock, T. Q., Sobie, S., Werner, A., & Shresthra, R. (2014). Climate change in Northeast
BC. Pacific Climate Impacts Consortium.
National Center for Atmospheric Research. (2016). NP index data provided by the Climate
Analysis Section. Boulder, CO, USA.

73
National Oceanic and Atmospheric Administration. (2017a). Pacific Decadal Oscillation
(PDO). Retrieved from National Centers for Environmental Information:
https://www.ncdc.noaa.gov/teleconnections/pdo/
National Oceanic and Atmospheric Administration. (2017b). Climate Prediction Center.
Retrieved from Historical El Nino:
http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.sht
ml
National Oceanic and Atmospheric Administration. (2017c). Southern Oscillation Index
(SOI). Retrieved from National Centers for Environmental Information:
https://www.ncdc.noaa.gov/teleconnections/enso/indicators/soi/
National Oceanic and Atmospheric Administration. (2017d). Pacific-North American (PNA).
Retrieved from National Centers for Environmental Information:
https://www.ncdc.noaa.gov/teleconnections/pna/
Pachauri, R. K., & Meyer, L. A. (2014). Climate change 2014: Synthesis report.
Contribution of working groups I, II and III to the fifth Assessment Report of the
Intergovernmental Panel on Climate Change. Geneva, Switzerland: IPCC.
Paul, S. S. (2013). In Analysis of Land Use and Land Cover Change in Kiskatinaw River
Watershed: A Remote Sensing, GIS & Modeling Approach (Thesis) (pp. 1-132).
Prince George: University of Northern British Columbia.
Psyllakis, J. M. (2006). A multi-scale analysis of forest structure and vertebrate diversity.
Prince George, BC: University of Northern British Columbia.
Rasmussen, P. P., Gray, J. R., Glysson, G. D., & Ziegler, A. C. (2009). Applications of
hydraulics: Sediment and erosion techniques. Reston, VR: U.S. Geological Survey.

74
Regonda, S. K., Rajagopalan, B., Clark, M., & Zagona, E. (2006). A multimodel ensemble
forecast framework: Application to spring seasonal flows in the Gunnison River
Basin. Water Resources Research, 42(9), W09404.
Román, M., Hall, D., & Riggs , G. (2017). MODIS/Terra snow cover 8-Day L3 global 500m
grid, version 6. Retrieved from National Snow and Ice Data Center:
http://nsidc.org/data/MOD10A2
Saha, G. C. (2015). Climate change induced precipitation effects on water resources in the
Peace Region of British Columbia, Canada. Climate, 3(2), 264-282.
Schnorbus, M., Werner, A., & Bennett, K. (2014). Impacts of climate change in three
hydrologic regimes in British Columbia, Canada. Hydrological Processes, 28(3),
1170-1189.
Scruton, M. (2017). Bearhole Lake water balance. Victoria, BC: Kerr Wood Leidal
Associates.
Shabbar, A. (2006). The impact of El Niño-Southern Oscillation on the Canadian climate.
Advances in Geosciences, 6, 149-153.
Shabbar, A., & Khandekar, M. (1996). Impact of El Niño Southern Oscillation on the
temperature field over Canada. Atmosphere Ocean, 34, 401-406.
Shabbar, A., Bonsal, B., & Khandekar, M. (1997). Canadian precipitation patterns associated
with the southern oscillation. Journal of Climate, 10, 3016-3027.
Sharma, A. R., & Déry, S. J. (2015). Climate change impacts on water resources in the
Nechako River Basin, BC. American Geophysical Union. San Fransisco, CA:
University of Northern British Columbia.
Shmueli, G. (2010). To explain or to predict? Statistical Science, 25(3), 289-310.

75
Spearman, C. (1904). The proof and measurement of association between two things. The
American Journal of Psychology, 15, 72-101.
Stahl, K., Moore, D. R., & McKendry, I. G. (2006). The role of synoptic-scale circulation in
the linkage between large-scale ocean-atmosphere indices and winter surface climate
in British Columbia, Canada. International Journal of Climatology, 26, 541-560.
Stewart, I. T., Cayan, D. R., & Dettinger, M. D. (2004). Changes toward earlier streamflow
timing across Western North America. Journal of Climate, 18, 1136-1155.
Tong, J., Déry, S. J., & Jackson, P. L. (2009). Interrelationships between MODIS/Terra
remotely sensed snow cover and the hydrometeorology of the Quesnel River Basin,
British Columbia, Canada. Hydrology and Earth System Sciences, 13, 1439-1452.
Trenberth, K. E. (1990). Recent observed interdecadal climate changes in the northern
hemisphere. Bulletin of the American Meteorological Society, 71, 988-993.
Trenberth, K. E., & Hurrell, J. W. (1994). Decadal atmosphere-ocean variations in the
Pacific. Climate Dynamics, 9, 303-319.
Tukey, J. W. (1958). Bias and confidence in not-quite large samples. Annals of Mathematical
Statistics, 29, 614.
USDA. (1972). Snow Survey and Water Supply Forecasting, Section 22. In SCS National
Engineering Handbook (pp. 1-35). Soil Conservation Service.
USDA. (1990). Water supply forecasts: A field office guide for interpreting streamflow
forecasts. Portland, Oregon, USA: West National Technical Centre.
Vivoni, E. R., Benedetto, F. D., Grimaldi, S., & Eltahir, E. A. (2008). Hypsometric control
on surface and subsurface runoff. Water Resources Research, 44, W12502.
Wallace, J. M., & Gutzler, D. S. (1981). Teleconnections in the geopotential height field
during the northern hemisphere winter. Monthly Weather Review, 109, 784-812.

76
Walling, D. E., & Webb, B. W. (1992). Monitoring Programmes. In The Rivers Handbook:
Hydrological and Ecological Principles (pp. 124-143). Oxford: Blackwell.
Woo, M., & Thorne, R. (2003). Comment on 'Detection of hydrologic trends and variability'
by Burn, D.H. and Hag Elnur, M.A., 2002. Journal of Hydrology, 277, 150-160.
World Meteorological Organization. (1986). Results of an intercomparison of models of
snowmelt runoff. Modelling Snowmelt-Induced Processes (pp. 103-112). Geneva,
Switzerland: IAHS.
Zhang, Y., Wallace, J. M., & Battisti, D. S. (1997). ENSO-like interdecadal variability: 190093. Journal of Climate, 10, 1004-1020.

FORECASTING SPRING RUNOFF IN THE KISKATINAW RIVER

77

Appendix A.
Table A.1: Spearman rank correlation coefficients between the North Pacific (NP),
Oceanic Niño Index (ONI), Southern Oscillation Index (SOI) and Pacific Decadal
Oscillation indices and timing of incremented fractional runoff volume (RX%) and total
runoff volume (TRV). Intensity of blue (red) shading corresponds to magnitude of the
negative (positive) rank correlation. Note that at 48 degrees of freedom a correlation of
0.28 is significant at a probability of 0.05.All seasonal averages were calculated for the
previous year leading up to the runoff period (25 March-1 July). Blank cells indicate no
significant correlation exists.
Index

R5%

R10%

R15%

R20%

0.30

R25%

R30%

R35%

R40%

R45%

R50%

NP_ANOM
NP_FALL

0.28

NP_WINT
NP_SPR
NP_SUM

0.35

0.29

0.30

ONI_FALL

-0.37

-0.4

-0.37 -0.35

-0.33

-0.33

-0.31

-0.32

-0.35

-0.34

ONI_WINT

-0.35 -0.37

-0.35 -0.33

-0.32

-0.31

-0.31

-0.33

-0.37

-0.36

ONI_SPR

-0.32 -0.36

-0.38 -0.39

-0.41

-0.41

-0.39

-0.38

-0.33

-0.28

ONI_SUM

-0.41 -0.42

-0.40 -0.39

-0.38

-0.38

-0.36

-0.37

-0.38

-0.36

SOI_FALL

0.42

0.45

0.44

0.44

0.42

0.42

0.4

0.41

0.42

0.41

SOI_WINT

0.41

0.42

0.38

0.37

0.34

0.34

0.31

0.32

0.36

0.37

0.30

0.33

0.36

0.37

0.38

0.36

0.37

0.33

0.29

0.45

0.39

0.34

0.31
-0.28

-0.30

-0.30

-0.28

SOI_SPR
SOI_SUM

0.45

PDO_FALL

-0.32 -0.34

-0.30 -0.28

PDO_WINT

-0.34 -0.36

-0.28

PDO_SPR
PDO_SUM

78

Table A.1 (Continued)
Index

R55%

R60%

R65%

R70%

R75%

R80%

R85%

R90%

R95%

TRV

NP_ANOM
PNA_FALL
PNA_WINT
PNA_SPR

0.29

PNA_SUM
ONI_FALL

-0.31

-0.31

-0.28

ONI_WINT

-0.34

-0.33

-0.30

ONI_SUM

-0.33

-0.32

-0.29

SOI_FALL

0.41

0.39

0.35

SOI_WINT

0.36

0.34

0.28

SOI_SPR

0.30

0.28

ONI_SPR

SOI_SUM
PDO_FALL
PDO_WINT
PDO_SPR
PDO_SUM

FORECASTING SPRING RUNOFF IN THE KISKATINAW RIVER
Appendix B.

Figure B.1. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed average winter (15
November–25 March) temperature at Environment and Climate Change Canada and
Ministry of Transportation and Infrastructure climate stations.

79

80

Figure B.2. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed average winter (15
November–25 March) temperature at Environment and Climate Change Canada and
Ministry of Transportation and Infrastructure climate stations.

81

Figure B.3. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed average winter (15
November–25 March) temperature at Environment and Climate Change Canada and
Ministry of Transportation and Infrastructure climate stations.

82

Figure B.4. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 95% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed average winter (15
November–25 March) temperature at Environment and Climate Change Canada and
Ministry of Transportation and Infrastructure climate stations.

83

Figure B.5. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed total winter (15
November- 25 March) precipitation at Environment and Climate Change Canada and
Ministry of Transportation and Infrastructure climate stations.

84

Figure B.6. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed total winter (15
November–25 March) precipitation at Environment and Climate Change Canada and
Ministry of Transportation and Infrastructure climate stations.

85

Figure B.7. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed total winter (15
November–25 March) precipitation at Environment and Climate Change Canada and
Ministry of Transportation and Infrastructure climate stations.

86

Figure B.8. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and surveyed 1 April snow water
equivalent from the Ministry of Forests, Lands, Natural Resource Operations & Rural
Development.

87

Figure B.9. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and surveyed 1 April snow water
equivalent from the Ministry of Forests, Lands, Natural Resource Operations & Rural
Development.

88

Figure B.10. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and surveyed 1 April snow water
equivalent from the Ministry of Forests, Lands, Natural Resource Operations & Rural
Development.

89

Figure B.11. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 25% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed average winter (15
November–25 March) atmospheric surface pressure at Environment and Climate Change
Canada and Ministry of Transportation and Infrastructure climate stations.

90

Figure B.12. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 50% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed average winter (15
November–25 March) atmospheric surface pressure at Environment and Climate Change
Canada and Ministry of Transportation and Infrastructure climate stations.

91

Figure B.13. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between hydrologic day of year that 75% of spring total runoff volume (25
March–1 July) occurs in the Kiskatinaw River, BC, and observed average winter (15
November–25 March) atmospheric surface pressure at Environment and Climate Change
Canada and Ministry of Transportation and Infrastructure climate stations.

92

Figure B.14. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the Kiskatinaw River,
BC, and observed average winter (15 November–25 March) temperature at Environment
and Climate Change Canada and Ministry of Transportation and Infrastructure climate
stations.

93

Figure B.15. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the Kiskatinaw River,
BC, and observed total winter (15 November–25 March) precipitation at Environment
and Climate Change Canada and Ministry of Transportation and Infrastructure climate
stations.

94

Figure B.16. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the Kiskatinaw River,
BC, and surveyed 1 April snow water equivalent from the Ministry of Forests, Lands,
Natural Resource Operations & Rural Development.

95

Figure B.17. Spatially interpolated Spearman’s rank correlation coefficients (black
contours) between total spring runoff volume (25 March–1 July) in the Kiskatinaw River,
BC, and observed average winter (15 November-25 March) atmospheric surface pressure
at Ministry of Transportation and Infrastructure climate stations.

FORECASTING SPRING RUNOFF IN THE KISKATINAW RIVER
Appendix C.

Figure C.1. Winter temperature (°C) departure from the 1951–2000 normal over Canada in
response to a positive phase of El Niño Southern Oscillation. Figure reproduced from
Shabbar (2006), used with permission under the Creative Commons Attribution 4.0 License.

Figure C.2. Winter precipitation (mm d-1) departure from the 1951–2000 normal over
Canada in response to a positive phase El Niño Southern Oscillation. Figure reproduced from
Shabbar (2006), used with permission under the Creative Commons Attribution 4.0 License.

96

FORECASTING SPRING RUNOFF IN THE KISKATINAW RIVER

97

Appendix D.
The following tables detail the resulting multiple regressions chosen for making
predictions of fractional streamflow timing (RX%) and total spring (25 March–1 July) runoff
volume (TRV) in the Kiskatinaw River (KR). Models with the highest adj.-r2 were selected if
no violation of assumptions were found, and if jackknifed (JK) coefficients were significant
at p = 0.05. See Chapter 4 for more information regarding methods. The following
abbreviations are used in these tables: (a) 24-25 March Snow Depth (MSD); (b) total
cumulative precipitation (CPPT); (c) mean temperature (MT); (d) snow water equivalent
(SWE); (e) Environment and Climate Change Canada (EC); (f) Ministry of Transportation
and Infrastructure (MOTI). Numbers refer to station or survey identifiers unique to each
agency. Jackknife replications used to determine coefficients in each case equal n.
Table D.1: Multiple regression results predicting the R5%.
n = 17
F3,16 = 26.77
p < 0.001
adj.-r2= 0.86
R5%
JK β
JK SE
t
p
MSD @ EC
-0.0693
0.0282
-2.6100
0.0190
(1173242)
CPPT @ MOTI
-0.1016
0.0183
-5.7100
<0.001
(22123)
MT @ EC
-2.5656
0.4076
-6.4400
<0.001
(1200560)
Const.

235.0992

4.2766

55.0100

<0.001

95% Conf. Interval
-0.1334

-0.0139

-0.1429

-0.0655

-3.4878

-1.7596

226.2073

244.3392

98
Table D.2: Multiple regression results predicting the R10%.
n = 17
F3,16 = 71.91
p < 0.001
adj.-r2 = 0.92
R10%
JK β
JK SE
t
p
MSD @ MOTI
0.0618
0.0085
7.1600
<0.001
(47121)
MT @ EC
-3.4850
0.9779
-4.0000
<0.001
(3053600)

95% Conf. Interval
0.0430

0.0792

-5.9801

-1.8340

Avg. PDO-Spring

-5.6196

0.7978

-6.8000

<0.001

-7.1147

-3.7323

Const.

197.4846

3.6155

54.3100

<0.001

188.6860

204.0152

Table D.3: Multiple regression results predicting the R15%.
n = 17
F3,16 = 34.54
p < 0.001
adj.-r2 = 0.90
R15%
JK β
JK SE
t
p
MSD @ MOTI
0.0538
0.0096
5.4600
<0.001
(47121)
MT @ EC
-3.7308
0.8450
-4.7800
<0.001
(3053600)

95% Conf. Interval
0.0319

0.0724

-5.8323

-2.2495

Avg. PDO-Spring

-5.2803

1.0226

-5.0600

<0.001

-7.3470

-3.0114

Const.

201.9501

4.3679

46.1200

<0.001

192.1940

210.7129

Table D.4: Multiple regression results predicting the R20%.
n = 17
F3,16 = 43.48
p < 0.001
Adj.-r2 = 0.93
R20%
JK β
JK SE
t
p
CPPT @ MOTI
-0.1138
0.0120
-9.2600
<0.001
(24225)
MT @ MOTI
-6.5967
0.8767
-7.4400
<0.001
(47121)

95% Conf. Interval
-0.1369

-0.0859

-8.3776

-4.6604

Avg. SOI-Summer

-4.7667

1.4598

-3.4700

0.0030

-8.1598

-1.9707

Const.

242.2735

3.6353

66.6300

<0.001

234.5239

249.9371

99
Table D.5: Multiple regression results predicting the R25%.
n = 17
F3,16 = 58.37
p < 0.001
adj.-r2 = 0.94
R25%
JK β
JK SE
t
p
CPPT @ MOTI
-0.1232
0.0129
-9.4500
<0.001
(24225)
MT @ MOTI
-6.0578
0.9118
-6.5800
<0.001
(47121)

95% Conf. Interval
-0.1498

-0.0949

-7.9317

-4.0657

Avg. SOI-Summer

-4.6149

2.2387

-2.2500

0.0390

-9.7832

-0.2914

Const.

249.4138

3.5725

69.8700

<0.001

242.0483

257.1949

Table D.6: Multiple regression results predicting the R30%.
n = 17
F3,16 = 71.74
p < 0.001
adj.-r2 = 0.94
R30%
JK β
JK SE
t
p
CPPT @ EC
-0.1917
0.0316
-6.0000
<0.001
(3081680)
CPPT @ MOTI
-0.1091
0.0126
-9.1100
<0.001
(24225)
MT @ MOTI
-5.2417
0.5756
-9.3800
<0.001
(46127)
Const.

258.2644

4.2974

60.1400

<0.001

Table D.7: Multiple regression results predicting the R35%.
n = 17
F3,16 = 51.22
p < 0.001
adj.-r2 = 0.92
R35%
JK β
JK SE
t
p
MSD @ MOTI
-0.2920
0.0612
-4.8500
<0.001
(46127)
CPPT @ MOTI
-0.1380
0.0148
-9.0700
<0.001
(24225)
MT @ EC
-11.1283 1.5566
-7.0800
<0.001
(1113581)
Const.

282.1729

3.2381

87.0000

<0.001

95% Conf. Interval
-0.2566

-0.1225

-0.1410

-0.0878

-6.6172

-4.1770

249.3445

267.5644

95% Conf. Interval
-0.4265

-0.1671

-0.1658

-0.1029

-14.3217

-7.7221

274.8418

288.5707

100
Table D.8: Multiple regression results predicting the R40%.
n = 16
F3,15 = 61.63
p < 0.001
adj.-r2 = 0.90
R40%
JK β
JK SE
t
p
CPPT @ MOTI
-0.1594
0.0137
-11.7400 <0.001
(24225)
MT @ EC
2.7793
0.5861
4.5600
<0.001
(2100709)
MT @ EC
-9.6778
1.5140
-6.3300
<0.001
(1077500)
Const.

281.6419

7.0406

39.8900

<0.001

Table D.9: Multiple regression results predicting the R45%.
n = 16
F3,15 = 36.57
p < 0.001
adj.-r2 = 0.89
R45%
JK β
JK SE
t
p
CPPT @ MOTI
-0.126
0.014
-8.910
<0.001
(24225)
MT @ EC
5.099
1.274
4.060
<0.001
(1181508)
MT @ EC
-12.526
2.259
-5.410
<0.001
(1077500)
Const.

264.663

8.449

31.510

<0.001

Table D.10: Multiple regression results predicting the R50%.
n = 16
F3,15 = 44.66
p < 0.001
adj.-r2 = 0.86
R50%
JK β
JK SE
t
p
CPPT @ MOTI
-0.1171
0.0140
-8.2100
<0.001
(24225)
MT @ EC
5.3789
1.3993
3.9600
<0.001
(1181508)
MT @ EC
-12.1634 2.3168
-5.1300
<0.001
(1077500)
Const.

269.3354

9.2389

29.3300

<0.001

95% Conf. Interval
-0.1897

-0.1314

1.4243

3.9226

-12.8177

-6.3636

265.8747

295.8881

95% Conf. Interval
-0.157

-0.096

2.451

7.881

-17.039

-7.407

248.244

284.263

95% Conf. Interval
-0.1446

-0.0850

2.5578

8.5230

-16.8235

-6.9474

251.3068

290.6911

101
Table D.11: Multiple regression results predicting the R55%.
n = 15
F3,14 = 123.11
p < 0.001
adj.-r2 = 0.91
R55%
JK β
JK SE
t
p
CPPT @ EC
-0.1846
0.0370
-5.0200
0.0000
(1181508)
CPPT @ MOTI
-0.0230
0.0098
-2.1200
0.0520
(38326)

95% Conf. Interval
-0.2648

-0.1063

-0.0420

0.0002

Avg. ONI-Summer

-12.7391

1.1759

-10.7100

0.0000

-15.1160

-10.0719

Const.

291.2632

7.5034

38.5900

0.0000

273.4660

305.6523

Table D.12: Multiple regression results predicting the R60%.
n = 16
F3,15 = 54.05
p < 0.001
adj.-r2 = 0.89
R60%
JK β
JK SE
t
p
MT @ MOTI
-5.1234
1.2009
-4.2700
0.0010
(250921)

95% Conf. Interval
-7.6830

-2.5637

Avg. SOI-Spring

15.6883

1.7109

9.1700

0.0000

12.0416

19.3351

1 April SWE @
(4A09)

-0.1011

0.0195

-5.1700

0.0000

-0.1428

-0.0594

Const.

280.0379

8.0251

34.9000

0.0000

262.9328

297.1429

Table D.13: Multiple regression results predicting the R65%.
n = 15
F3,14 = 25.74
p < 0.001
adj.-r2 = 0.85
R65%
JK β
JK SE
t
p
MSD @ EC
0.0412
0.1066
0.3900
0.7050
(3062451)
CPPT @ EC
-0.2071
0.0719
-2.8800
0.0120
(1181508)
12.1430
1.9064
6.3700
0.0000
Avg. SOI-Spring
Const.

283.0682

4.9382

57.3200

0.0000

95% Conf. Interval
-0.1874
0.2698
-0.3613

-0.0530

8.0542

16.2318

272.4767

293.6596

102
Table D.14: Multiple regression results predicting the R70%.
n = 15
F3,14 = 14.44
p < 0.001
adj.-r2 = 0.85
R70%
JK β
JK SE
t
p
MSD @ MOTI
0.8118
0.1409
5.9500
<0.001
(47091)
1 March SWE @
-0.0164
0.0045
-3.9500
0.0010
(4B01)
1 March SWE @
-0.4026
0.0636
-6.4100
<0.001
(4C05)
Const.

309.1992

7.1198

43.6200

<0.001

Table D.15: Multiple regression results predicting the R75%.
n = 16
F3,15 = 30.28
p < 0.001
adj.-r2 = 0.83
R75%
JK β
JK SE
t
p
MSD @ MOTI
0.4391
0.1447
3.0600
0.0080
(47091)
1 April SWE @
-0.0444
0.0112
-3.5500
0.0030
(1E03A)
1 April SWE @
-0.2214
0.0245
-9.0700
<0.001
(4C05)
Const.

315.8872

6.9284

45.2700

<0.001

Table D.16: Multiple regression results predicting the R80%.
n = 15
F3,14 = 36.22
p < 0.001
adj.-r2 = 0.83
R80%
JK β
JK SE
t
p
MT @ EC
-5.0914
0.9377
-5.2000
<0.001
(3060903)
1 April SWE @
-0.1735
0.0176
-9.9100
<0.001
(4C05)
1 March SWE @
-0.2150
0.0487
-4.0600
0.0010
(1A23)
Const.

283.5928

4.5164

62.8100

<0.001

95% Conf. Interval
0.5360

1.1406

-0.0275

-0.0081

-0.5444

-0.2715

295.2657

325.8064

95% Conf. Interval
0.1338

0.7507

-0.0639

-0.0159

-0.2739

-0.1696

298.8702

328.4054

95% Conf. Interval
-6.8840

-2.8617

-0.2123

-0.1367

-0.3024

-0.0934

273.9963

293.3696

103
Table D.17: Multiple regression results predicting the R85%.
n = 18
F3,17 = 120.69
p < 0.001
adj.-r2 = 0.85
R85%
JK β
JK SE
t
p
CPPT @ EC
0.0857
0.0250
3.5900
0.0020
(4031581)
MT @ EC
-5.0444
0.5619
-9.0800
<0.001
(1168520)
1 April SWE @
-0.0507
0.0079
-6.6100
<0.001
(1B07)
Const.

287.5202

5.1207

56.2500

<0.001

Table D.18: Multiple regression results predicting the R90%.
n = 17
F3,16 = 41.07
p < 0.001
adj.-r2 = 0.89
R90%
JK β
JK SE
t
p
MT @ EC
-9.0833
1.1311
-8.3000
<0.001
(1060841)
1 April SWE @
-0.0534
0.0072
-7.6400
<0.001
(1B07)
1 April SWE @
-0.00463 0.0010
-4.5400
<0.001
(3B01)
Const.

335.4224

4.8776

69.0000

<0.001

Table D.19: Multiple regression results predicting the R95%.
n = 17
F3,16 = 60.31
p < 0.001
adj.-r2 = 0.91
R95%
JK β
JK SE
t
p
MSD @ EC
-0.1128
0.0192
-5.6100
<0.001
(1175122)
MT @ EC
-6.3661
0.5749
-11.2400 <0.001
(1060841)
1 April SWE @
-0.0235
0.0041
-5.9500
<0.001
(1B07)
Const.

327.2625

2.8063

116.6700

<0.001

95% Conf. Interval
0.0369

0.1426

-6.2903

-3.9191

-0.0690

-0.0356

277.2355

298.8430

95% Conf. Interval
-11.7865

-6.9911

-0.0699

-0.0395

-0.0070

-0.0025

326.2047

346.8847

95% Conf. Interval
-0.1481

-0.0668

-7.6825

-5.2450

-0.0335

-0.0159

321.4474

333.3454

104

Table D.20: Multiple regression results predicting total spring runoff volume (ha m).
n = 24
F4,23 = 33.64
p < 0.001
adj.-r2= 0.89
Total
Volume
JK β
JK SE
t
p
95% Conf. Interval
(ha m)
CPPT @
EC
540.7233
72.7654
7.4300
<0.001
390.1965
691.2500
(3081680)
MT @ EC
1068.9820
-6.5800
<0.001
(1181508) 7039.1200
9250.4780 4827.7620
MT @ EC
5820.8660 1495.6820
3.8900
0.0010
2726.8110 8914.9210
(1206197)
CPPT @
EC
-156.0408
51.4442
-3.0300
0.0060
-262.4612
-49.6205
(2100701)
Const.
28751.830 5272.1310
-5.4500
<0.001
39658.060 17845.600
0
0
0