Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical simulation of snow hydrology for management purposes Woo, Ming-Ko 1972

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

Item Metadata

Download

Media
831-UBC_1972_A1 W66.pdf [ 9.4MB ]
Metadata
JSON: 831-1.0101517.json
JSON-LD: 831-1.0101517-ld.json
RDF/XML (Pretty): 831-1.0101517-rdf.xml
RDF/JSON: 831-1.0101517-rdf.json
Turtle: 831-1.0101517-turtle.txt
N-Triples: 831-1.0101517-rdf-ntriples.txt
Original Record: 831-1.0101517-source.json
Full Text
831-1.0101517-fulltext.txt
Citation
831-1.0101517.ris

Full Text

NUMERICAL SIMULATION OF SNOW HYDROLOGY FOR MANAGEMENT PURPOSES by MING-KO WOO,. M.A., University of Hong Kong, I 9 6 7 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of Geography We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA February, 1972 In presenting this thesis in partial fulfilment of the requirements for an advanced, degree at the University of British Columbia, I agree that the Library shall make i t freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of G e o g r a p h y The University of British Columbia Vancouver 8, Canada Date 30th M a r c h 1972 ABSTRACT i i The objective of th i s study i s to provide a method of estimating snow storage and melt-water release for small coastal mountain basins i n the temperate forest zone of southwestern B r i t i s h Columbia. It adopts a numerical simulat-ion approach, using d a i l y p r e c i p i t a t i o n and temperature as input parameters. The data are usually obtainable from 'base-stations* located i n the v a l l e y s . Most stations have a short history so that long-term prediction of snow d i s t r i b u t i o n i n small basins requires both temporal and s p a t i a l extensions of data. Temporal extension involves treating p r e c i p i t a t i o n and temperature f i r s t as random variables, and then as inputs to a deterministic evaluation of snow storage. In order to model p r e c i p i t a t i o n , a second-order Markov chain technique i s used to simulate the occurrence of wet and dry days, while the quantity of d a i l y p r e c i p i t a t i o n i s also generated for wet days. Daily temperature i s simulated by l i n e a r l y superimpos-ing various deterministic and random components of temperature time s e r i e s . S p a t i a l extension involves a recognition of a l t i t u d -i n a l and vegetation cover zones. The extension of p r e c i p i t a t -ion and temperature records from a base-station to other a l t i t u d i n a l zones and to forest sectors i s achieved by empiri-c a l functions derived from observations at the University of B r i t i s h Columbia Research Forest at Maple Ridge, B.C. i i i With an understanding of p r e c i p i t a t i o n and temper-ature, i t i s possible to simulate snow events. For snowfall, a p r o b a b i l i t y matrix based on minimum d a i l y temperature i s constructed from data available at the University of B r i t i s h Columbia Research Forest, For snowmelt, a d i s t i n c t i o n i s made between clear-weather and r a i n melts, and d a i l y maximum temperature and d a i l y r a i n f a l l are used to estimate melt-water release. A simulation model i s developed to extend d a i l y temperature and p r e c i p i t a t i o n from short-term records. Using the simulated data as inputs, a second model simulates snow-pack storage as the balance between incoming p r e c i p i t a t i o n and outgoing snowmelt i n both the forested and the open sectors of each of the several a l t i t u d i n a l zones. The model i s tested by comparing one hundred years of simulated temperature, p r e c i p i t a t i o n and snow storage data with h i s t o r i c a l records, and the correspondence i s s a t i s -factory. The model i s useful because even when prediction i s d i f f i c u l t owing to meagre input data, i t i s possible with this method to produce long-term estimates of snow resources i n small coastal basins. IV CONTENTS Page Abstract 1 1 Contents 1 V L i s t of tables ' v i i i L i s t of figures i x L i s t of symbols Acknowledgment X V 1 I INTRODUCTION 1 1-1 Research methodology 1 1-2 Statement of a problem 6 1-2.1 Objective of the study 8 1-3 Choice of an approach 9 1-3.1 Temporal model 9 I- 3.2 S p a t i a l model 12 I- 4 Presentation of thesis 15 II STUDY AREA 1? I I - l Hydrology 17 II-2 Topography 22 I I - 3 Geology 27 II - 3 . 1 Lithology 27 II-3.2 Structure 29 II-3.3 S u r f i c i a l deposits 30 Vegetation 32 III INSTRUMENTATION 35 I I I - l Data a v a i l a b i l i t y 35 III-2 P r e c i p i t a t i o n measurements 36 III-3 Temperature measurements 38 III-4 Measurement of snowmelt 39 III-5 F i e l d programme 42 Page PRECIPITATION 45 IV-1 Physical model 45 IV-2 Choice of a temporal model 46 IV-2,1 Mathematicalmodel for temporal d i s t r i b u t i o n of p r e c i p i t a t i o n 51 IV-2.2.1 P r e c i p i t a t i o n trend and season-a l i t y 51 IV-2.2.2 Wet and dry s p e l l s 53 IV-2.2.3 Quantity of d a i l y p r e c i p i t a t i o n 54 IV-2.3 Result of analysis of h i s t o r i c a l data ' 56 IV-2.3.1 Trend 56 IV-2.3.2 Spells 59 IV-2.3.3 Daily p r e c i p i t a t i o n 6 l IV-3 S p a t i a l variations of p r e c i p i t a t i o n 6 l IV-3.1.1 Physiographic considerations -physical model 63 IV-3.1.2 Physiographic considerations -mathematical model and s t a t i s t i c a l analysis 65 IV-3.2 E f f e c t s of vegetation on p r e c i p i t a t -ion 66 IV-3.2,1 Forest canopy and r a i n f a l l -physical model 69 IV-3.2.2 Forest canopy and r a i n f a l l -s t a t i s t i c a l analysis 72 IV-3.2.3 Forest canopy and snowfall -physical model 73 IV- 3.2.4 Forest canopy and snowfall -mathematical model and r e s i i l t of analysis 76 IV- 4 Discussion 81 TEMPERATURE 83 V- 1 Physical model 83 V-2 Temporal variations of temperature 84 V- 2.1 Theoretical model 85 v i Page V-2.1.1 Trend component 88 V-2.1.2 O s c i l l a t o r y component 89 V-2.1.3 Markov component 91 V-2.1.4 Random component 92 V-2.1.5 Diurnal temperature ranges 92 V-2.2 Results of analysis 93 V-2.2.1 Fourier component 93 V-2.2.2 Cosine wave component 96 V-2.2.3 Markov component 98 V-2,2.4 Random component 98 V-2,2.5 Diurnal temperature ranges 99 V-3 Spa t i a l v a r i a t i o n of temperature 99 V-3.1 Physiographic considerations -physical model and r e s u l t of analysis 99 V- 3.2 Forest cover and temperature - physi-c a l model and r e s u l t of analysis 105 V- 4 Summary 108 VI SNOWFALL AND SNOWMELT EVENTS 112 VI- 1 Snowfall events - s t a t i s t i c a l model 112 VI-2 Snowmelt events 113 VI- 2.1.1 Multiple regression approach 113 VI-2,1.2 Energy balance approach 115 VI-2.1.3 Temperature index approach 116 VI-2.2 Measurement of snowmelt 117 VI-2.2,1 Indirect measurement of snowmelt through streamflow analysis 11? VI-2.2,2 Indirect measurement by snow survey 117 VI-2.2.3 Laboratory measurement of snow-melt 118 VI-2,2.^ Lysimeter measurement of snowmelt 118 VI-2.3 S t a t i s t i c a l model for snowmelt 121 VI-2.3.1 Clear-weather melt 121 VI-2.3.2 Rain melt 122 v i i Page VI-3 Spatial extension of snowmelt equations 124 VI- 4 Summary and discussion 126 VII SIMULATION MODEL 127 VII- 1 Algorithms 127 VII-1.1 Sub-model I 129 VII-1.2 Sub-model II 129 VII-1.3 Computer time 131 VII-2 Model te s t i n g 131 VII-2.1 P r e c i p i t a t i o n data 134 VII-2.2 Temperature data 134 VII-2.3 Snowpack data 137 VII-3 Analysis of simulated snowpack data 137 VII-3.1 Magnitude of snowpack 137 VII-3.2 Snowpack persistence 140 VII- 4 Model l i m i t a t i o n s 142 VIII DISCUSSION 145 VIII- 1 Model performance 145 VIII-2 Model and management applications 146 VIII-3 Model extension 146 References 149 Appendix LIST OF TABLES v m Page 1.1 Examples of predictive techniques used i n streamflow analysis 7 II.1 Morphometric and vegetation c h a r a c t e r i s t i c s of Blaney Creek basin and Jacob Creek basin 25 111.1 Summary of p r e c i p i t a t i o n and temperature recording stations 37 111.2 Comparison of lysimeter readings 4 l IV.1 Summary of selected references on stochastic studies of p r e c i p i t a t i o n 49 IV.2 Transition matrices for the occurrence of wet days, and mean d a i l y p r e c i p i t a t i o n by month 60 IV.3 Two-sample Kolmogorov-Smirnov tests on the d i s t r i b u t i o n of s p e l l lengths 62 IV.4 Least-square f i t : d a i l y p r e c i p i t a t i o n vs al t i t u d e , University of B r i t i s h Columbia Research Forest 67 IV.5 Covariance analysis of throughfall (dependent) vs r a i n f a l l i n open (independent variable) 74 IV.6 Selected references on snow interception and forest snow accumulation studies 77 V . l Parameters from various components of temper-ature time-series (temperature records of U.B. C. Admin) 95 V.2 Least-square f i t * maximum temperature vs al t i t u d e , U.B.C. Forest 106 V.3 Least-square f i t : minimum temperature vs al t i t u d e , U.B.C. Forest 107 LIST OF FIGURES ix Page 11.1 Location and physiographic map of the Univer-s i t y of B r i t i s h Columbia Research Forest. Maple Ridge, B.G. 18 11.2 Mean discharge-basin area r e l a t i o n s h i p for eleven selected coastal mountain streams 19 11.3 Flood frequency curves for three selected coastal mountain streams 20 11.4 Variations i n monthly discharge for two coastal mountain streams 21 11.5 Daily discharge at three hydrometric stations, University of B r i t i s h Columbia Research Forest 23 11.6 Frequency of slope directions i n Blaney and Jacob Creek basins. Inset shows frequency of fa u l t and joint s t r i k e directions 2k 11.7 Geological map of the University of B r i t i s h Columbia Research Forest 28 11.8 D i s t r i b u t i o n of forested and open areas i n the study area Jk 111.1 Automatic snow lysimeter and manual snow lysimeter (inset) ko 111.2 Hydrometric and meteorologic stations at the University of B r i t i s h Columbia Research Forest kk IV.1 Surface pressure pattern at 0600 PST, 31st Dec. 1968 showing the continental A r c t i c and mari-time A r c t i c fronts and the accompanying p r e c i -p i t a t i o n belt k7 IV.2 Flow-chart for p r e c i p i t a t i o n analysis 57 IV.3 Theoretical and empirical p r o b a b i l i t y d i s t r i -butions of d a i l y p r e c i p i t a t i o n by month 58 IV.k P r e c i p i t a t i o n - a l t i t u d e r e l a t i o n s h i p , Univer-s i t y of B r i t i s h Columbia Research Forest, 1967-68 68 X Page IV,5 T h r o u g h f a l l - r a i n f a l l r e l a t i o n s h i p from two s i t e s at Seymour Basin, B.C. 71 IV.6 T h r o u g h f a l l - r a i n f a l l r e l a t i o n s h i p , University of B r i t i s h Columbia Research Forest 75 IV.7 Throughfall-snowfall r e l a t i o n s h i p , University of B r i t i s h Columbia Research Forest 80 V . l Flow-chart for the decomposition of temperat-ure time-series 94 V.2 Theoretical and empirical d i s t r i b u t i o n s of residual mean wave temperatures and wave-lengths 97 V.3 P r o b a b i l i t i e s of temperature ranges as condit-ioned by d a i l y means 100 V.4 Decomposition of temperature time-series: f i t t i n g Fourier series to mean dai l y temperat-ures 101 V.5 Decomposition of temperature time-series: f i t t i n g cosine waves to r e s i d u a l temperatures 102 V.6 Decomposition of temperature time-series: f i t t i n g f i r s t order Markov component to r e s i -dual temperatures 103 V.7 Minimum d a i l y temperature: the forest vs the open, Placid Lake, 1969-70 water year 109 V.8 Maximum d a i l y temperature: the forest vs the open, Placid Lake, 1969-70 water year 110 VI.1 Probability of snowfall as conditioned by minimum dai l y temperatures 114 VI.2 R a i n f a l l , temperature, and snowmelt records, Placid Lake s i t e , February to A p r i l , 1971 120 VI.3 Snowmelt vs maximum d a i l y temperature, Placid Lake s i t e 123 VI.4 Snowpack outflow vs da i l y r a i n f a l l and maximum dai l y temperature, Placid Lake 125 VII.1 Flow-chart for the simulation of daily temper-ature, p r e c i p i t a t i o n , and snowpack 128 x i Page VII.2 Flow-chart for temperature and p r e c i p i t a t i o n simulation sub-model 130 VII.3 Flow-chart f o r snow storage water-balance subroutine 132 VII.4 Flow-chart for snow storage simulation sub-model 133 VII.5 D i s t r i b u t i o n of maximum (upper band) and minimum (lower band) d a i l y temperatures based on 41 years of records from Stave F a l l s , B.C. 135 VII.6 D i s t r i b u t i o n of maximum (upper band) and minimum (lower band) d a i l y temperatures based on 100 years of simulated record 136 VII.7 Snowpack water-equivalent exceedance probabi-l i t y for various a l t i t u d e s on 1st A p r i l 138 VII.8 Monthly changes i n snow water-equivalent exceedance p r o b a b i l i t y at 1200 m., based on 100 years of simulated data 139 VII.9 Medians and i n t e r - q u a r t i l e ranges of snow storage for various a l t i t u d i n a l zones, based on 100 years of simulated data l 4 l VII.10 Means and variations i n the persistence of snowpack i n the opening, based on 100 years of simulated data 143 VII.11 Means and variations i n the persistence of snowpack i n the forest, based on 100 years of simulated data 144 VIII. 1 Schematic diagram showing possible areas of extension for the present model 148 x i i LIST OF SYMBOLS Chapter a IV, v, VI empirical constant a i V c o e f f i c i e n t for the i-th harmonic i n a Fourier series A, A' II map and surface area (L ) a II slope angle a V amplitude of (temperature) wave b IV, v, VI empirical constant b I index, r e f e r r i n g to base-station b. 1 V c o e f f i c i e n t for the i-th barmonic i n a Fourier series P V an angle (Ji I Borel f i e l d Cov V covariance 6 E I , I V an error term an event E [ ] IV, V expected value of f IV, V subscript, denoting forest s i t e f IV, V a function of, or density function of F IV d i s t r i b u t i o n function g I, IV a function of h I, IV, V altitude (L) H I, IV variable related to topography H V seasonal component in temperature time series Dimension of some variables are given i n the force (F), length (L), time (t) system. Chapter i I, I I , IV, V an index I IV snow interception (L) j I, IV, V an index k I, IV, V an index L V h a l f period i n harmonic analysis 3 -1 L I rate of water loss (L t ) X V wave-length (\ IV mean da i l y p r e c i p i t a t i o n ( L t - 1 ) Mc VI clear-weather melt (L) M R VI outflow from snowpack, associated with rain-on-snow events (L) n IV, V, VI an index, or frequency of events N V, VI t o t a l number of events 5 I, V a random variable or component o IV, V subscript, denoting open s i t e w V an o s c i l l a t o r y component i n temperature time series Si I sample space P 3 -1 I rate of p r e c i p i t a t i o n ( I r t ) P V persistence component i n temperature time series Pr I, IV, V, VI pr o b a b i l i t y of i> I pr o b a b i l i t y measure 0 V phase angle i n harmonic analysis Dimension of some variables are given in the force (F), length (L), time (t) system. x i v Chapter Q I discharge ( I r t ~ ) Q E VI latent heat flux (FLt" 1) Q H VI sensible heat flux ( F L t - 1 ) Q-^  VI long-wave radiati o n (FLt-"'") Q M VI energy for snowmelt ( F L t - 1 ) Q N VI net all-wave radiation (FLt"" 1) Q s VI short-wave r a d i a t i o n (FLt x) r ^ IV, V k-th order autocorrelation c o e f f i c i e n t R IV, VI r a i n f a l l (L) V k-th order autocorrelation c o e f f i c i e n t (population parameter) s IV, V number of lag periods S I a storage term (L- ) S IV snowfall (L) I a set of events o* V standard deviation t I, IV, V a time parameter t a V standard normal deviate at pro b a b i l i t y l e v e l a T IV t o t a l time period % V trend component i n temperature time series 0 V, VI temperature <5 V mean temperature 6* V temperature range Dimension of some variables are given i n the force (F), length (L), time (t) system. XV Chapter e't 6" V r e s i d u a l temperature 9 V a temperature sequence U I transfer function V I variable related to vegetation cover Var IV, V variance W I, IV state of a multi-stage process or hydrologic system X I, IV, V a variable Z V a deterministic component ACKNOWLEDGMENT xvi I wish to thank Dr. Olav Slaymaker for his encour-agement and help during a l l stages of my thesis work, and to thank Drs. Bertram Goodell and Gary Gates for the many f r u i t f u l discussions. I am indebted to Professors John Stager and Richard Copley for th e i r useful advice, and to Professor Coulthard for his valuable comments and the permission for me to quote his data from the M i l i t z a basin. Dr. Rolf K e l l e r h a l s , now at the University of Alberta, has provided me with streamflow records he gathered at various s i t e s i n the University of B r i t i s h Columbia Research Forest, and Mr. John Walters, the Director of the Research Forest, has supplied me with maps and a i r photographs of the Research Forest, together with the permission to stay at the Loon Lake camp during the winter f i e l d seasons. To them, I wish to express my appreciation, I am also g r a t e f u l to Messr. John Zeman and Robert Gi l b e r t who offered me much help i n the construction and ins-t a l l a t i o n of f i e l d stations, I g r a t e f u l l y acknowledge f i n a n c i a l assistance at various times from the following sources! National Research Council, Nation Advisory Committee on Water Resources Research and the University of B r i t i s h Columbia Research Fund. A fellowship from the University of B r i t i s h Columbia, Institute of Animal Resource Ecology enabled me to concentrate on my research project during the period 1 9 6 9 - 7 1 . 1 I INTRODUCTION Management and planning of water resources require information regarding the quantity, r e l i a b i l i t y and d i s t r i -bution of water supplies. Accurate assessment of available water resources derives i d e a l l y from f i e l d data? but thi s often proves to be uneconomical, p a r t i c u l a r l y for areas of poor a c c e s s i b i l i t y . In mountainous t e r r a i n with a humid, temperate forest environment, snow constitutes an important source of water supply. While s p a t i a l variations of the mountain snowpacks necessitate the maintenance of a fine net-work of snow courses, snow surveys are expensive operations. It i s therefore desirable to develop some i n d i r e c t approaches which would enable the evaluation of snow storage given other forms of hydrologic or meteorologic records which are r e a d i l y a v a i l a b l e . The quest for such i n d i r e c t assessment techniques i n place of extensive f i e l d programmes i s a s i g n i f i c a n t research enterprise, and i s the primary objective of t h i s thesis. Although several techniques may reach equivalent solutions, usually one technique i s more e f f i c i e n t than the others. In terms of the present problem of snow hydrology i n mountainous areas, research methodology therefore merits some consideration, 1-1.0 Research methodology The aims of most resources studies may be defined as either explanatory or predictive. The former emphasizes 2 the processes, the l a t t e r i s more concerned with the response to such processes. Although better understanding of processes leads to more accurate prediction, incomplete knowledge of processes and the lack of relevant data commonly leads research workers to analyse responses f i r s t of a l l . Another consideration i s the scope of the research project. Most water resources studies f a l l into the categories of component and o v e r a l l catchment models 1 (Dawdy and O'Donnell, I 9 6 5 ) . For explanatory purposes, a model should take into account as large a number of elements as i s required to produce a f a i t h f u l r e p l i c a of i t s prototype. Unfortunately, errors inherent i n the data and errors a r i s i n g from approximation of the physical system (Dawdy and Bergmann, 1 9 6 9 , p. 9 5 8 ) often lead to the following paradoxical s i t u a t i o n . When more v a r i a -bles are introduced, a model may appear more r e a l i s t i c , but the " f i t t e d parameters may f i t r e a l i t y less and less i n the numerical value" (Dawdy and O'Donnell, 1 9 6 5 ) . Hence, most predictive models attempt to st r i k e a balance between theoret-i c a l soundness and s t a t i s t i c a l goodness of f i t . For e f f i c i e n t 2 modelling , assumptions are made about the behaviour of the 1 . In the present context, model i s defined as ideal i z e d representation of r e a l i t y . 2 . Popov ( I 9 6 8 , p. 4 5 3 ) has provided several guidelines for modelling. These include: (a) p h y s i c a l l y correct idea about the major processes (b) minimum number of variables and parameters necessary for a s u f f i c i e n t l y complete description of the processes (c) clear physical meaning of the parameters and the p o s s i b i l i t y of determining th e i r numerical values (d) minimum assumptions l i m i t i n g the use of the model (e) p o s s i b i l i t y of further development and improvement of the model on the basis of new information. 3 less c r u c i a l elements i n the model (Dooge, 1 9 ° 8 , p. 6 3 ) . The problem i s 1 . to i d e n t i f y the p r i n c i p a l elements 2 . to select model equations, and 3 . to evaluate the parameters for such equations (Vemuri and Vemuri, 1 9 7 0 , p. 2 5 ) . Closely related to the question of model complexity i s the choice of scale (Riley et a l . , I 9 6 8 , p. 3 2 ) . As with most geographical enquiries, studies may range from s i t e scale to basin or regional scales (Slaymaker, 1 9 6 9 ) and from a time base of minutes to years. Phenomena which are important at one scale may be subsumed at another. Thus d a i l y fluctuations i n snowmelt cannot be detected from the annual water budget, nor can effects of latitude on temperature variations i n a basin be f e l t i f the basin i s small. The choice of scale therefore influences the selection of pertinent components for the models. When model components are i d e n t i f i e d , solutions to a model may be obtained with a va r i e t y of techniques. With s u f f i c i e n t and accurate information, analytic solutions often provide s a t i s f a c t o r y r e s u l t s that are exact and p h y s i c a l l y meaningful. This unfortunately seldom applies to complex models or to areas where data are scarce. Operational cons-t r a i n t s i n terms of time and resources often r e s u l t i n incom-plete data c o l l e c t i o n networks. This c a l l s for some predictive techniques such as systems investigations (Amorocho and Hart, 1 9 6 4 ) . A system comprises inputs and outputs and i n t e r -4 r e l a t i o n s amongst various components which are mathematically described by transfer functions (Vemuri and Vemuri, 1 9 7 0 ) , With the transfer functions known, outputs can be predicted. Advancement i n computer technology enables the development of numerical simulation techniques. Simulation has been defined as a "model of some s i t u a t i o n i n which the elements of the s i t u a t i o n are represented by arithmetic and l o g i c a l processes that can be executed on a computer to predict the dynamic properties of the s i t u a t i o n " (Emshoff and Sisson, 1 9 7 0 ) . Problems of great complexity which are not amenable to analytic solutions may be tackled with simulation models. As applied to hydrologic investigations, there are examples such as the Harvard Water Program (Hufschmidt and Fi e r i n g , 1 9 6 6 ) , Stanford Water Model IV (Crawford and Linsley, 1 9 6 6 ) , or the Streamflow Synthesis and Reservoir Regulation model of the U.S. Army (Rockwood, I 9 6 I: Schermerhorn and Kuehl, I 9 6 8 ) . A d i s t i n c t i o n i s sometimes made between deterministic r and p r o b a b i l i s t i c models (Chow, 1 9 6 4 ) , In the case of deter-ministic models, most mathematical functions r e l a t i n g hydro-logic processes and variables (or transfer functions) are known. Numerical values for the transfer functions are derived from physical p r i n c i p l e s , from least-square f i t t i n g of empirical observations or by parameter optimization (Lichty et a l , , 1 9 6 8 ) , Thus, these models require h i s t o r i c a l data inputs. In the absence of s u f f i c i e n t records or when many of the transfer functions are unknown, p r o b a b i l i s t i c models provide an al t e r n a t i v e . Such models make use of the mathematics of random processes. 5 To define a random process, consider a set of out-comes which depend on parameter t . As t i s allowed to vary, a Borel f i e l d 1 of events i s induced with a correspond-ence between the elementary events and the values of t . Then, the random process w i l l be defined as the parameter set, the Borel f i e l d and the correspondence between them. In pa r t i c u l a r , i f t i s a time parameter, the process i s known as stochastic. The Borel f i e l d induced by varying parameter t corresponds to a family of random variables defined on the p r o b a b i l i t y space (fl, <8 , ) where fl i s the sample space encomp-assing a l l possible elementary events <B i s the c o l l e c t i o n of a l l events i n the sample space f i i s a p r o b a b i l i t y measure defined on the sets of & . When the behaviour of such random variables i s known, i t 1, <3 i s a Borel f i e l d when the following requirements are s a t i s f i e d * (1) the elements of & are subsets of some set $ (2) the set # i s an element of <B (3) i f and Eg are subsets of $ contained i n (6 as elements, then the sets E-^  , E 2 , E2. u E2 a n d E-jflEg are also elements of <B . This implies that the empty set i s also an element of £ (4) for any denumerable sets of subsets E-^ , Eg, ... of $ , which are a l l elements of © , t h e i r sum and product also belong to (B . 2. I f fl i s a sample space with a pr o b a b i l i t y measure -p ; i f X i s a real-valued function defined over the points of fi , then X i s c a l l e d a random vari a b l e . 6 w i l l be possible to simulate data which are s t a t i s t i c a l l y indistinguishable from h i s t o r i c a l records. Depending on the nature and scope of the problem and the operational constraints, various techniques have been developed for hydrologic studies, A comprehensive survey i s beyond the present scope, but several examples are summarized i n table 1 ,1 , These examples i l l u s t r a t e the range of e x i s t i n g techniques to handle problems of d i v e r s i f i e d nature. On the whole, short-term forecasts are commonly handled determin-i s t i c a l l y , drawing on h i s t o r i c a l data available f o r c e r t a i n s i t e s (e.g. Anderson and Rockwood, 1970; Crawford and Linsley, 1966) . Long-term predictions, however, are usually achieved with the aid of a stochastic approach, making use of the s t a t i s t i c a l properties of the random variables involved (e.g. Chow and Ramaseshan, 1965s F i e r i n g , 1964) . A t h i r d category i s the assessment of hydrologic responses by varying the parameters within r e a l i z a b l e l i m i t s . In th i s case, both deterministic or stochastic models have been applied (e.g. Jeng and Yevjevich, I9665 Shen, I 9 6 5 ) . Mathematical models have to be developed i n accord-ance with the physical models because research methodology cannot be separated from the problems under investigation. The following sections w i l l describe the research problem together with an outline of the research approach, 1-2.0 Statement of a problem The small mountainous r i v e r basins around Vancouver Table 1 . 1 Examples of predictive techniques used i n streamflow analysis Author Objective Time scale S p a t i a l scale Technique Anderson & Rockwood, 1970 snowmelt runoff synthesis day basin deterministic Beard, 1 9 6 7 streamflow prediction month s i t e Markov Beard, 1 9 6 8 design flood hour basin deterministic Carlson et a l . , 1970 streamflow description year basin time series Chow & Ramaseshan, 1965 r a i n f a l l & runoff synthesis hour s i t e Markov Crawford & Linsley, I 9 6 6 runoff synthesis hour basin deterministic F i e r i n g , 1 9 6 4 low-flow prediction day s i t e Markov Harms & Campbell, 1 9 6 7 streamflow prediction month s i t e Markov Hufschmidt & Fiering, I 9 6 6 r i v e r system study month basin Markov Huggins & Monke, I 9 6 8 runoff synthesis minute basin deterministic Jeng & Yevjevich, I 9 6 6 lake outflow prediction year s i t e Monte Carlo Langbein, 1 9 5 8 reservoir storage prediction month s i t e queuing theory Lichty et a l . , 1968 peak-flow synthesis storm s i t e optimization Payne et a l . , 1968 streamflow synthesis day s i t e Markov Ri l e y et a l . , 1968 streamflow synthesis minute day ,si t e deterministic Schermerhorn & Kuehl, 1968 streamflow synthesis and reservoir regulation day basin deterministic Shen, 1 9 6 5 flood analysis basin deterministic Thomas & Fier i n g , 1962 flood and streamflow synthesis hour, month basin, s i t e Markov Young et a l . , I969 reservoir regulation day basin Markov 8 do not provide s u f f i c i e n t data on the basis of which water resources can be managed. The most r e a d i l y available types of data are d a i l y temperatures and p r e c i p i t a t i o n , but most of the stations are confined to the glaciated valley f l o o r s . As the demands for domestic and i n d u s t r i a l water supplies increase, there i s a growing need for assessing the d i s t r i b u t -i o n a l aspects of water resources. The bulk of water supply comes as winter p r e c i p i t a t i o n which reaches the basins i n the forms of r a i n or snow, depending on the al t i t u d e of the freezing l e v e l . Consequently, simultaneous accumulation and depletion of snowpacks within the same basin i s a common pheno-menon, and d a i l y streamflow from these basins becomes highly e r r a t i c . Drainage basins l y i n g below the t r e e - l i n e are part-i c u l a r l y prone to such alternating accretion and melt processes. The majority of these small basins, however, have not been gauged, nor are any snow data available for such basins. The prediction of snow resources i n such basins poses a p r a c t i c a l problem. 1-2.1 Objective of the study In the l i g h t of the problems stated above, the object-ive of the present study i s : "to develop a model for the prediction of snow storage i n small drainage basins'*" i n a humid, temperate forest environment, using commonly available data of temperature and p r e c i p i t a t i o n " . 1 . Following Chow ( 1 9 6 2 ) , "a small drainage basin may be defined as one that i s so small that i t s s e n s i t i v i t y to high-intensity r a i n f a l l s of short duration and to land-use i s not suppressed by the channel c h a r a c t e r i s t i c s " . 9 1 - 3 . 0 Choice of an approach With a limited amount of supporting data, i t i s d i f f i c u l t to adopt an a n a l y t i c a l approach. On the other hand, numerical simulation offers an appropriate tool i n handling the complex problem outlined above. Daily p r e c i p i t a t i o n and temperature from valley stations serve as the basic input data to the model. Although some of the v a l l e y stations have been i n operation for over 3 0 years, others have a much shorter history. It becomes necessary to r e p l i c a t e , then to extend, the records of these 'base-stations* or 'index-plots' i n order to acquire s u f f i c i e n t l y lengthy input data for the long-term prediction of snow resources. While these records may be regarded as point-measurements i n space, an assessment of basin snow storage requires the knowledge of d i s t r i b u t i o n s over an area. Hence, data from base-stations have to be broadened from a micro-to a meso-scale (Slaymaker, 1 9 6 9 ) . In other words, the present study w i l l concentrate on both the temporal and s p a t i a l d i s t r i b u t i o n s of snow resource„ 1 - 3 . 1 Temporal model The study of temporal variations i n basin storage 10 w i l l centre around the water-balance approach 1 which hinges on the l i n e a r equation: t t t £ Qdt = S pdt - £ Ldt + S ( 1 . 1 ) o o o where Q i s runoff (IPT""1) 3 -1 p i s p r e c i p i t a t i o n input (LJ T ) L i s water loss (L T ) S i s storage (L^) The equation i s integrated between time 0 and time t . Runoff i s represented as the amount of p r e c i p i t a t i o n , less evapotranspiration loss, plus (or minus) the quantity of water freed from (or added to) basin storage. Ground water storage i s minimized for areas underlain by impervious bedrock with th i n s u r f i c i a l deposits, such as most basins i n the- coastal region. Also, evapotranspiration i s much reduced i n winter, so that runoff may be e f f e c t i v e l y expressed i n terms of pre-c i p i t a t i o n and conditions of the snowpack. 1 . The purely s t a t i s t i c a l approach of r e l a t i n g discharge to other possible co-varying elements i s more suitable for predicting runoff on a regional scale (e.g. Mustonen, 1967? Quick, I 9 6 55 Sopper and L u l l , 1 9 6 5 ) . For small basins, these predictive equations may vary from one period to another and from one basin to another so that they are neither stationary i n time nor are they s p a t i a l l y trans-f e r a b l e (see for instance, Wilson's comments on snow course-discharge r e l a t i o n s h i p , I 9 6 6 ) . On the other hand, the energy-balance approach provides a more r a t i o n a l approach (Crawford and Linsley, 1 9 6 6 ; Erickson and McCorquodale, I 9 6 7 ; U.S. Army, 1 9 5 6 ) . To be operative, however, f a i r l y exact knowledge i s assumed of many meteorological elements and snowpack conditions. In view of climatic data constraints for the coastal areas of B r i t i s h Columbia, th i s approach i s not adopted. 11 Conceptually, water c y c l i n g through the basin may be considered as an endogenous variable which undergoes changes i n states from p r e c i p i t a t i o n to snow i n the pack to melt-water that ultimately drains out of the basin. As the water element i s transformed through i n t e r a c t i o n with some exogenous variables (such as temperature), i t may be said to have undergone some multi-stage process (Bellman and Kalaba, 1 9 6 5 ) . I f an input of an i n i t i a l state W(0) i s transformed by some mathematical function U , the output at the end of stage one w i l l be of the state W(l) = U 1[W (0 ) ] ( 1 . 2 ) At stage (j+ 1 ) , the water element w i l l be at a state W(j+1) = u . + 1 [ w ( j ) j = u j + 1 [ u j [ w ( ( i - i ) ] ] = uj+i["« U 1[W (0 ) ] ... J ( 1 . 2 . a ) Hence, given an i n i t i a l input of p r e c i p i t a t i o n and a l l the transfer functions , the quantity of snow storage w i l l be expressions of several states i n the multi-stage process. Complete knowledge of the hydrologic system, however, i s seldom available, so that i t i s appropriate to incorporate an element of uncertainty into each stage, and equation I . 2 . a becomes W(j+1) = u j + 1 [ W(j), l(j+1) ] ( 1 . 3 ) where ? i s a time-dependent random variable. The inclusion of a random variable into the equation advances the problem into the realm of stochastic processes , One major asset of the stochastic approach i s i t s capacity to generate events from combinations of data input and system parameters. The present study i s primarily concerned with long-term properties of snow hydrology, so that i n d i v i d u a l events are s i g n i f i c a n t only i n t h e i r c o n t r i -bution to statements regarding the p r o b a b i l i t i e s of the mag-nitude, frequency and duration of the whole array of events. Owing to our i n a b i l i t y to sample through a broad expanse of time, the true nature of the sample space w i l l never be obtained. Furtunately, by synthesizing a long sequence of events through a stochastic multi-stage process, the sample s t a t i s t i c s would l i k e l y converge to the population character-i s t i c s under the law of large numbers (Lamperti, I 9 6 6 ) . One assumption i s e s s e n t i a l : that the behaviour of the processes, the input variable and the exogenous variables do not change with time. In the absence of a trend, i t w i l l be possible to extrapolate events s t o c h a s t i c a l l y on the basis of short h i s t o r i c a l records available at the base-station. 1-3.2 S p a t i a l model Larson (1965) distinguished basin hydrology into a 1. The l i n k between stochastic and deterministic processes i s as follows: I f for any fixed t , there exists a unique value ^ ( t ) such that Pr[x-6 < 5(t) < x+6] = 1 for any 6>0 , then the only possible outcome of X(t) w i l l be x , and the process i s deterministic. The deterministic process i s thus a special case of the stochastic process ( K a r t v e l i s h v e l i , 1967) . 13 land phase and a channel phase. Discharge at the lower course of a stream may be regarded as the sum of land phase runoff from various sectors of the basin. In other words, i t i s possible to fashion a basin into hydrologically uniform sectors"*" to study the water balance of each sector. This concept has been applied to r e l a t i v e l y f l a t t e r r a i n where a basin was subdivided into square c e l l s and the runoff routed through each c e l l (Huggins and Monke, 1 9 6 8 ) . Instead of using regular grids, i t i s more convenient to segment the densely vegetated and rugged basins on the bases of topography/ and vegetation cover. For small basins, water balances for each sector are not independent of one another because the influence of weather systems often transcends t h e i r boundaries. Such s p a t i a l autocorrelation may be described by a set of mathematical functions which relate conditions at the base-sta t i o n to other parts of the basin: W k(t) = f [ W b(t), v, H ] (1.4) where and are the states of the hydrologic system at time t f o r the base-station and for sector k respectively v i s a vegetation factor H i s a topographic factor. In the absence of s u f f i c i e n t data, a subdivision l o Uniformity i s a r e l a t i v e term depending on the scale of generalization and the order of magnitude at which hydro-logic events manifest themselves. Hence, the degree of refinement required of the solut i o n and the a v a i l a b i l i t y of data would play a deciding role i n determining the size of hydrologically uniform sectors i n a basin. 14 of the basin i s kept simple so that sectors are defined i n terms of a l t i t u d i n a l zones and the presence or absence of continuous forest cover 1. Equation 1.4 w i l l then be s i m p l i f i e d considerably. Given a small basin, i t i s assumed that the entire surface f a l l s under the same weather system during any given day. Then, hydro-meteorological readings for the open sector of a l t i t u d i n a l zone k w i l l be estimated from readings taken at a base-station as follows: w k ( t ) = g l [ w b ( t ) , (h k-h b) ] ( 1 . 5 ) where and W^  are hydro-meteorolo-g i c a l conditions at base-station and the open sector of zone k h^ and h^ are alt i t u d e s of the base-station and of zone k Introducing forest e f f e c t s , hydro-meteorological readings at the forest sector of a l t i t u d i n a l zone k w i l l be W k = g 2 ^ W k ( t ) ^ = g 2[ g j l W b(t), (h k-h b) ] ] ( 1 . 6 ) where represents hydro-meteorolo-g i c a l readings at the forest sector of a l t i t u d i n a l zone k . When transfer functions g 1 and g 2 are known, the input component i n the water balance equation can be estimated for a l l sectors based on point-measurements at the base-station. This enables subsequent evaluation of hydrologic outputs from 1. Coastal clear-cut logging makes i t p r a c t i c a l to divide the land surface into forested and open sectors. 15 each sector. One assumption underlying the s p a t i a l model i s that the c h a r a c t e r i s t i c s of random processes are transfera-ble from s t a t i o n to station, though the parameters would take up d i f f e r e n t numerical values. Accepting t h i s , the problem i s to obtain these functions which would suitably describe the patterns of s p a t i a l v a r i a t i o n s . 1-4.0 Presentation of thesis To operationalize the temporal and s p a t i a l models described above, the d i s s e r t a t i o n w i l l be presented i n the following steps: (1) the selection and description of an appropriate research area. (2) the evaluation of the transfer functions which enable an extension of the short-time-period p r e c i p i t a t i o n and temperature records. (3) the development of routing procedures to determine the patterns of s p a t i a l variations i n p r e c i p i t a t i o n and temperature, (4) the description of a snow storage simulation model and the presentation of simulation r e s u l t s . Unless otherwise stated, the metric system w i l l provide the standard units of measurement. In keeping with available meteorological and hydrologic records, data analyses and simulation use d a i l y i n t e r v a l s as basic time units . A water-year i s assumed to have 365 days and begins on 1st October so that winter seasons are not interrupted by the 16 beginning of a calendar year. A l l computations were made at the University of B r i t i s h Columbia Computing Centre with i t s IBM 3 6 0 / 6 7 . Computer programmes were written i n FORTRAN IV language except for those which are available as l i b r a r y programmes at the Computing Centre. 17 II STUDY AREA To operationalize the temporal and s p a t i a l models outlined by equations 1 . 3 ? I . 5 P and 1.6, several mathematical functions have to be derived and many c o e f f i c i e n t s have to be determined empirically. It i s assumed that intensive observations made i n one small basin over a period of three years w i l l provide an estimate of the relevant mathematical functions. The University of B r i t i s h Columbia Research Forest at Maple Ridge, B.C., was selected as the study area. It covers an area of 50 square kilometres, at approximately 4 9 ° 19* N, and 1 2 2 ° 3 4 ' W ( f i g . I I . 1 ) . The following para-graphs are designed to indicate the representativeness of the U.B.C. Research Forest within the southwestern B r i t i s h Colum-b i a Coast Mountain region. I I - 1 . 0 Hydrology Coastal mountain streams are characterized by highly i r r e g u l a r d a i l y discharge. Mean d a i l y discharge tends to increase with basin area ( f i g . I I . 2 ) , but i s by i t s e l f a poor indicator of streamflow regime because the amplitude of da i l y fluctuations i s exceedingly large. In general, winter i s the period of high flow and the peak may occur at any time between October and February. Most streams have a wide range i n flood discharge ( f i g . 1 1 , 3 ) . Even monthly discharges exhibit great differences. Figure I I . 4 shows that i n certa i n years, d i s -charge for some winter months may be rather low. In which F i g . I I . 1 L o c a t i o n and p h y s i o g r a p h i c map o f t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a R e s e a r c h F o r e s t , M a n i e R i d e e , B.C. 19 1— i — r GA22 o SQUAMISH Basin Area (A) in km F i g • I I . 2 Mean d i s c h a r g e - b a s i n a r e a r e l a t i o n s h i p f o r e l e v e n s e l e c t e d c o a s t a l m o u n t a i n s t r e a m s 20 •05 I l 3 -4 -5 -6 -7 -8 -9 -9? Cumulative Probability F i g . II. 3 F l o o d f r e q u e n c y c u r v e s f o r c o a s t a l m o u n t a i n s t r e a m s t h r e e s e l e c t e d 21 F i p . I I . 4 V a r i a t i o n s i n m o n t h l y d i s c h a r g e c o a s t a l m o u n t a i n s t r e a m s f o r two 22 case, higher spring runoff from snowpack compensates low flows of e a r l i e r months. This pattern becomes more marked for streams further inland or when the basins occupy subs-t a n t i a l areas i n the alpine zone. In these cases, spring-melt runoff becomes more pronounced. For small basins i n the temperate forest environ-ment, the Blaney and Jacob basins are t y p i c a l . They both l i e within the University of B r i t i s h Columbia Research Forest. Their channels are steep because they flow on surfaces pre-sculptured by Pleistocene g l a c i e r s . Log-jams are common along the channels, and burst occasionally under t o r r e n t i a l flows. Streamflow i s rapid and i s marked by numerous hydraulic jumps. Most streams respond very rapi d l y to i n d i v i d u a l r a i n -storms or to rain-on-snow events, and then recede almost as quickly ( f i g . II.5)• The only exceptions occur where sizeable lakes bring about storage e f f e c t s . II-2.0 Topography Pleistocene g l a c i a t i o n has l e f t the Research Forest with steep mountain slopes and glaciated v a l l e y s , many of which are occupied by lakes. The main valleys trend north-south and NNE-SSW, but many tri b u t a r y valleys run NW-SE. As a r e s u l t , slope directions have f a i r l y pronounced west- and southwest-facing components ( f i g . II.6). The presence of flat-bottomed valleys and r e l a t i v e l y f l a t ridge-tops i s r e f l e c t e d i n the d i s t r i b u t i o n of slope angles with altitude (table II.1) which shows that most steep 24 F i r - F r e q u e n c y o f s l o p e d i r e c t i o n s i n fllaney and J a c o b C r e e k b a s i n s . I n s e t shows f r e q u e n c y o f f a u l t and j o i n t s t r i k e d i r e c t i o n s 2 5 Table II.1 Morphometric and vegetation c h a r a c t e r i s t i c s of Blaney Creek basin and Jacob Creek basin Name: Blaney Creek basin Basin area: 8 . 9 0 sq. km. Height range: 2 2 - 8 0 8 m. Cumulative Cumulative Percent Mean slope % height % area forested angle ( i n ) 0 . 4 0 . 0 3 . 6 0 . 1 6 . 7 0 . 4 9 . 9 1 . 5 1 3 . 1 3 . 3 1 6 . 3 4 . 8 1 9 . 5 6 . 2 2 2 . 6 7 . 1 2 5 . 8 7 . 6 2 9 . 0 8 . 4 3 2 . 2 8 . 7 3 5 . ^ 9 . 6 3 8 . 5 1 2 . 0 41 . 7 2 7 . 2 4 4 . 9 40 .2 4 8 . 1 • 5 1 . 2 5 1 . 3 5 9 . 7 5k. 5 6 5 . 5 5 7 . 6 7 2 . 5 6 0 . 8 7 9 . 2 64 . 0 8 6 . 2 6 7 . 2 9 0 . 5 7 0 . 4 9 2 . 6 7 3 . 5 9 3 . 9 7 6 . 7 9 4 . 9 7 9 . 9 9 6 . 1 8 3 . 1 9 7 . 5 8 6 . 3 9 8 . 3 8 9 . 4 9 8 . 8 9 2 . 6 9 9 . 3 9 5 . 8 9 9 . 8 1 0 0 . 0 - -1 0 0 . 0 0.4 1 0 0 . 0 1 . 2 1 0 0 . 0 0 . 7 8 9 . 5 0 . 8 1 0 0 . 0 0 . 6 9 3 . 3 0 . 6 1 0 0 . 0 1 . 3 1 0 0 . 0 1 . 9 8 8 . 9 0 . 5 1 0 0 . 0 0 . 7 1 0 0 . 0 1 . 3 9 8 . 5 1 . 1 48 . 9 4 . 2 6 0 . 8 1 0 . 8 ?4.4 9 . 7 8 0 . 2 1 2 . 1 81.4 1 6 . 1 8 2 . 1 1 5 . 1 78.4 1 0 . 9 89.4 1 2 . 0 8 2 . 3 1 3 . 0 9 2 . 0 1 7 . 6 5k. 3 1 9 . 8 6 9 . 7 2 7 . 7 1 0 0 . 0 2 5 . 1 1 0 0 . 0 1 7 . 5 1 0 0 . 0 1 9 . 3 1 0 0 . 0 20.4 1 0 0 . 0 1 3 . ^ 1 0 0 . 0 1 1 . 3 Table II. 1 (continued) Namei Jacob Creek basin Basin area. 1 0 . 5 9 sq. km. Height rangei 3 0 8 - 1 0 3 9 m. Cumulative Cumulative % height % area 2 . 3 0 . 6 5 . 7 6 . 2 9 . 2 9 . 5 1 2 . 6 14 .2 1 6 . 0 2 0 . 2 1 9 . ^ 2 7 . 1 2 2 . 8 3 2 . 4 2 6 . 3 3 9 . ^ 2 9 . 7 4 9 . 9 3 3 . 1 5 7 . 5 3 6 . 5 64.1 3 9 . 9 7 0 . 6 4 3 . 4 7 5 . 8 46 . 8 7 9 . 0 5 0 . 2 8 3 . 9 5 3 . 6 8 7 . 8 5 7 . 0 8 8 . 8 6 0 . 5 9 1 . 9 6 3 . 9 9 3 . 7 6 7 . 3 9 4 . 8 7 0 . 7 9 5 . 5 7 ^ . 1 9 6 . 7 7 7 . 6 9 7 . 3 81 .0 9 7 . 5 84 . 4 9 8 . 4 8 7 . 8 9 9 . 2 9 1 . 2 9 9 . ^ 9 4 . 7 9 9 . 8 Percent Mean slopg forested angle ( i n ) 1 0 0 . 0 2 . 3 8 6 . 7 7 . 0 8 8 . 2 1 2 . 2 9 1 . 3 1 2 . 3 9 3 . 3 1 1 . 2 2 5 . 4 8 . 2 1 0 0 . 0 1 2 . 3 5 5 . 8 1 3 . 7 5 ^ . 5 1 5 . 1 64 .2 1 9 . 2 9 1 . 9 2 5 . 0 7 0 . 0 2 2 . 2 7 3 . 5 2 5 . 2 8 5 . 5 2 2 . 2 7 8 . 9 2 2 . 3 7 0 . 5 2 0 . 6 8 6 . 5 2 4 . 2 9 6 . 2 2 3 . 6 8 8 . 7 1 7 . 5 1 0 0 . 0 1 0 . 9 1 0 0 . 0 1 1 . 8 6 6 . 7 1 1 . 3 1 0 0 . 0 1 1 . 7 1 0 0 . 0 2 1 . 0 80 .0 1 7 . 7 5 0 . 0 1 7 . 8 1 0 0 . 0 6 . 1 5 0 . 0 1 2 . 5 2 7 slopes occur i n the middle a l t i t u d i n a l range of 5 0 0 - 7 0 0 m. One of the effects of the topography i s that basin map areas are less than t h e i r surface areas. For example, Blaney basin, with a map area of 8 . 7 sq. km., has a surface area of 8 . 9 sq. km. when corrected for slope e f f e c t s 1 . S i m i l a r l y , the map area of Jacob basin i s 1 0 . 1 sq. km., but the basin area i s 1 0 . 6 sq. km. Despite the small size of the basins, a l t i t u d e s range from 2 0 to over 1 , 0 0 0 m. The basins are exclusively i n the humid, temperate forest environment, and the alpine zone i s not represented. I I - 3 . 0 Geology Geology plays a s i g n i f i c a n t part i n influencing the land phase of small basin hydrology and w i l l therefore be described i n greater d e t a i l i n the following sections, I I - 3 . 1 Lithology The basins are underlain by plutonic rocks of the Coast intrusions, ranging from granite to d i o r i t e , with the more a c i d i c members outcropping l o n g i t u d i n a l l y along the central part of the Research Forest ( f i g . I I . 7 ) . Texture varies from fine to coarse, but the predominant c r y s t a l 1 . Correction was made by summing the surface over various a l t i t u d i n a l zones, with surface area A? of zone i given by: 1 A i = A i s e c a i where A^ i s map area of zone i a i i s average slope angle of zone i F i g . I I . 7 G e o l o g i c a l nap o f t h e U n i v e r s i t v o f B r i t i s h C o l u m b i a R e s e a r c h F o r e s t 2 9 diameter i s around 2 - 3 mm. Chemical weathering i s s l i g h t , so that the re s i d u a l s o i l mantle i s t h i n . The pluton i s traversed by numerous dykes whose composition ranges from granophyre to d o l e r i t e and whose trends are N-S and NW-SE. In the plutonic rock, inclusions are abundant 0 Size of inclusions vary from xenoliths several centimetres across to roof-pendants of older sediments. There are no sharp contacts between the plutonic rocks and the older (pre-Jurassic) sedimentaries. Strong metamorphism has altered the pre-Jurassic sedimentaries into meta-arkose, hornblende-granulite, amphibolite or gneiss, so that th e i r geo-hydrological properties are simi l a r to the plutonic rocks. The basins are therefore quite impervious except where f a u l t s or joints abound or where pockets of s u r f i c i a l deposits develop into great thickness. I I - 3 . 2 Structure Numerous fa u l t s c r i s s - c r o s s the basins. Although minor f a u l t s are more common, some fault-zones reach a width of 8 to 10 metres and are followed by topographic features such as valleys or axes and rims of depressions (such as Gwendoline Lake or the eastern margin of Eunice Lake). Pre-dominant s t r i k e directions of f a u l t s or joints are N-S and NW-SE, which are followed by many tributary valleys i n the basin ( f i g . I I . 6 ) , Most f a u l t - and joint-planes have high angles of dip. Even e x f o l i a t i o n joints are rather steeply i n c l i n e d partly because of steep topography (Bradley, 1963> Chapman, 1 9 5 8 ) . Joint-spacing i s uneven, becoming closer 30 around zones of disturbances such as f a u l t s or dykes. Open joints are sometimes found near the ground surface, presumably as a r e s u l t of block d i s i n t e g r a t i o n and the release of stress by unloading of over-burden. In some cases, joints may open to a width of about 10 cm. and are f i l l e d with s o i l washed from ground surface. Open joints and f a u l t s provide suitable passages i n massive rocks f o r water seepage and ground-water storage which sustains low flow during dry seasons. II-3.3 S u r f i c i a l deposits Gl a c i a t i o n has l e f t most of the upper slope surfaces bare of any unconsolidated material. In general, s u r f i c i a l deposits may be categorized according to the s i t e of 1 occurrence . (1) Deposits at ridge-tops and steep mountain slopes They are mainly slopewash deposits which are often u n s t r a t i f i e d and have organic materials l y i n g above a zone of fine sand. The contact between such deposits and bedrock i s quite d i s t i n c t , and marks the zone of rapid subsurface flow (interflow) when r a i n water or melt-water f a i l s to penetrate 1. Cf. Lacate's (1965) subdivision of the Research Forest into four land associations J (A) mountainous to strongly r o l l i n g , granitic-cored uplands (B) h i l l y to gently r o l l i n g , granitic-cored uplands and valleys (C) f l a t to gently r o l l i n g , complex of g l a c i o - f l u v i a l and s u b - s t r a t i f i e d d r i f t deposits (D) f l a t to gently r o l l i n g , complex of glacio-marine deposits and s u b - s t r a t i f i e d d r i f t . into bedrock. L o c a l l y , when the t e r r a i n f l a t t e n s (e.g. between Gwendoline and Eunice Lakes), pockets of thicker and s t r a t i f i e d colluvium may accumulate, with sub-angular gravels set i n a matrix of fin e sand. These deposits are patchy i n occurrence and are by no means as extensive as Roddick has id e a l i z e d i n his geological map (Roddick, 1 9 6 5 ) . (2) Valley f l o o r deposits Many parts of the main valleys are or have been occupied by lakes. Lacustrine deposits, however, are not abundant. On the other hand, where the valleys f l a t t e n i n the v i c i n i t y of some present lakes (such as the head of Marion Lake), thick layers of organic deposits accumulate. These deposits are invariably water-logged but are of minor hydrologic significance because of lim i t e d s p a t i a l d i s t r i -bution. (3) Deposits at lower a l t i t u d e s (below 400 m.) Southward and southwestward towards the f o o t h i l l s , s u r f i c i a l deposits increase i n thickness. G l a c i a l t i l l covers the higher f o o t h i l l areas. Its surface i s weathered and i s t h i n l y veneered by slopewash deposits. Unweathered t i l l possesses a bluish-grey s i l t y to clayey matrix which makes i t impervious. Under t h i s t i l l i s a thick layer of contorted (overridden) sand which i s exposed at a gravel p i t 1 km, northeast of the main gate. According to Armstrong ( 1 9 5 7 ) t t h i s i s the Surrey T i l l of the Vashon g l a c i a t i o n , the f i n a l ice-sheet to occupy th i s area. 32 Below 1 0 0 m., Whatcom glacio-marine d r i f t (Armstrong, 1 9 5 7 ) underlies gentle slopes. This consists of boulders embedded i n bluish-grey clay. The marine d r i f t was deposited during the s h e l f - and sea-ice stages of the r e t r e a t i n g ice-sheet (Armstrong and Brown, 1 9 5 ^ ) . In places, i t has been disturbed by slumping to produce folded structures, Interbedded with, and superimposed onto the marine d r i f t are gravel and sand from submarine slopewash and beach gravels. The Whatcom marine d r i f t i s overlain i n part by the d e l t a i c and f l u v i a l Capilano gravel deposits which are post-Vashon i n age (Armstrong, 1 9 5 7 ) . Although the gravel and sand deposits have high pote n t i a l for ground-water, they are not important for the present study because of l i m i t e d occurrences i n the basin, 1 1 - 4 , 0 Vegetation Most areas i n the Research Forest are covered with Douglas f i r (Pseudotsuga menziesii)' 1', western hemlock (Tsuga  heterophylla) and western red cedar (Thuja p l i c a t a ) which are t y p i c a l of forests of Southern P a c i f i c Coast section (Rowe, 1 9 5 9 ) . Hardwood such as alder, broadleaf maple and vine maple occupy some of the lower or wetter s i t e s , p a r t i c u -l a r l y on logged areas. V i r g i n stands of conifers occur mainly at the northern end of the Research Forest, with an average tree age of about 3 0 0 years. The areas south of Loon Lake 1 , The Puget Sound and Columbia sections of the Society of American Foresters retains the s p e c i f i c name of t a x i f o l i a rather than adopting menziesii. 33 have recovered from f i r e s started by s e t t l e r s i n 1868 and are now stocked with dense mature stands intersparsed with some older trees. Upper Blaney basin, however, was extensively harvested by r a i l r o a d logging between 1920 and 1931t but regrowth has again established a 30-year old forest . Since 1953t logging has been practised mainly over the northern parts of the Forest. As with most areas i n the Coast Moun-tains, c l e a r - c u t t i n g i s practised which leaves plots of completely open areas, making i t feasible to c l a s s i f y cover conditions as forested and open ( f i g . II.8), F i g . I I . 8 D i s t r i b u t i o n o f f o r e s t e d and open a r e a s i n t h e s t u d y a r e a 3 5 III INSTRUMENTATION The University of B r i t i s h Columbia Research Forest i s well served with gravel roads. In winter, main roads to Loon Lake camp and to Spur 17 are kept open to enable Forest personnel to record temperature and p r e c i p i t a t i o n data twice a day. Where the snow was not ploughed, snow-shoes had to be used to gain access to higher s i t e s . T r a v e l l i n g took up h a l f of the time i n the f i e l d , though accommodations made available at Loon Lake camp helped to conserve daylight hours during winter f i e l d seasons. Public access to the Research Forest i s limited (except to hik e r s ) . This i s an asset because f i e l d i n s t a l l -ations are less l i a b l e to serious disturbance. I I I - l Data a v a i l a b i l i t y The Research Forest has four o f f i c i a l Canadian Atmospheric Environment Service weather stations and other less well-maintained s i t e s , two hydrometric stations operated by the Inland Waters Branch, plus other discontinuous stream-flow measurements (Kellerhals, 1970? Zeman, I969) and an excellent set of contour maps at a scale of 1*4,800. For supplemental data, thermographs, r a i n gauges and snowmelt lys.imeters were set up to cover most parts of the basins. Following i s an account of instrumentation related to thi s study. 36 I I I - 2 P r e c i p i t a t i o n measurements (1) The present study draws on d a i l y p r e c i p i t a t i o n records published by the Atmospheric Environment Service for three stations inside the Research Forest and one s t a t i o n at Haney (table I I I . l ) , In addition, d a i l y p r e c i p i t a t i o n records from two stations at M i l i t z a basin were available for the water year 1 9 6 7 - 6 8 . These stations cover an a l t i t u -d i n a l range up to over 550 m. (2) To compare the amount of r a i n f a l l inside and outside the forest canopy, two s i t e s were selected at elevations of 230 m, (near a bridge across Blaney Creek) and 520 m. (near Placid Lake). One Forester's r a i n gauge (No. 9 x 2 0 0 , Western Fire Equipment Ltd.) was placed i n the opening at each of the s i t e s . Three p l a s t i c wedge gauges were placed inside the forest at 230 m. and four standard gauges located inside the forest at 520 m. Readings were taken weekly or biweekly through a year. To check the v a l i d i t y of these gauge r e s u l t s , eight home-made gauges using funnels and jars were located at each of the s i t e s for a period of one to three months, T-tests were performed on the mean throughfall values as obtained from the o r i g i n a l gauges and the eight additional gauges. Differences between the means were found to be s t a t i s t i c a l l y i n s i g n i f i c a n t , (3) To measure fresh snow deposition inside and outside the forest, three s i t e s were chosen i n the Research Forest at 360 m. (road junction s i t e ) , 520 m. (Placid Lake s i t e ) and Table I I I . l Summary of p r e c i p i t a t i o n and temperature recording stations Station Long. (W) o « Lat. ( N ) o « Elevation (m.) Source Haney 532 [P] 122 3? 49 15 3 Atmosph. Environ. Ser. UBC Admin [P &T] 122 34 49 1 7 1 4 6 Atmosph. Environ. Ser. UBC Forest [P&T] 122 35 49 1 8 183 Atmosph. Environ. Ser. Loon Lake [P&T] 122 35 49 1 8 373 Atmo s ph, Environ. Ser. Placid Lake [T] 122 3 4 49 19 510 Present study Lower M i l i t z a [p] 122 35 49 19 518 Prof. Coulthard Middle M i l i t z a [p] 122 3 5 49 19 564 Prof. Coulthard P p r e c i p i t a t i o n data only T temperature data only P&T both p r e c i p i t a t i o n and temperature data 3 8 535 m. (Gwendoline Lake s i t e ) . The number of sampling s i t e s was r e s t r i c t e d by a c c e s s i b i l i t y . One round of measure-ments i n mid-winter, for example, took approximately 8 to 10 hours. Each s i t e had paired sample plots, one i n the opening and the other i n the fores t . In the open plot, two to three sampling points were used while four to six points were used for the forest plot because of higher variance of snow deposited under the canopy. New snow deposited on a sample point was measured by weighing the column c o l l e c t e d above a one-square foot p l a s t i c sheet mounted on a wooden frame. The use of p l a s t i c sheets i n place of snow-boards preserves the albedo of the pack, though extra care was required to prevent puncturing the p l a s t i c when snow was scooped up for weighing. After each measurement, the frames were reset onto the surface of the pack. III-3 Temperature measurements The present study makes use of d a i l y temperature records published by the Atmospheric Environment Service for three s i t e s located i n the Research Forest (table I I I . l ) . In addition, an R. Fuess hygrothermograph was set up at an opening on the shore of Placid Lake. These four stations covered an a l t i t u d i n a l zone up to 510 m. To study the e f f e c t of forest canopy on a i r temperature, another R. Fuess hygrothermograph was set up i n the forest to pair with the Pla c i d Lake st a t i o n . From 39 time to time, both recorder readings were checked with mercury thermometers and the readings did not depart by more than 0 . 5°C. Temperature traces from these thermographs provided d a i l y temperature maxima and minima for a period of over one year, III-4 Measurement of snowmelt Before i n s t a l l i n g a f u l l y automatic snowmelt lysimeter"*", i t was necessary to evaluate the performance of melt-water c o l l e c t o r s . To do so, p l a s t i c garbage cans were used as containers of melt-water, and were buried flush with ground surface at Loon Lake and Placid Lake s i t e s . Holes were d r i l l e d through the c i r c u l a r wooden l i d s (diameter 1 7 . 5 " ) to enable rapid passage of melt-water, Por each lysimeter, a plexiglass access tube projected above the l i d to enable the reading of melt-water depths with water-test paper (CANLAB P 1 0 6 6 C ) attached to a long rod (inset of f i g , I I I . l ) . The accuracy of manual gauge readings was ± 5 mm. To assess the adequacy of a single lysimeter used i n open system (see section V I - 2 , 2 . 4 ) , three lysimeters were arranged i n a row at each s i t e so that t h e i r r e s u l t s could be compared. Table II I . 2 shows that discrepancies amongst three lysimeters were low. In view of serious boundary ef f e c t s which can 1 . Lysimeter i s defined as an apparatus f o r measuring the quantity of water percolating through the s o i l ; or consists of a container f i l l e d with the material under examination, with measuring devices to assess the amount of water and dissolved materials which have passed through (Monkhouse, 1 9 6 5 ) . i g . I I I . l A u t o m a t i c snow l y s i m e t e r and m a n u a l snow l y s i m e t e r ( i n s e t ) 41 Table 111 ,2 Comparison of lysimeter readings Date Placid N Lake : C lysimeters S Loon N Lake lysimeters C S 3 0 t h D e c , 1969 211 231 234 4 t h Jan., 1970 10 0 8 1 0 t h Jan., 1970 23 25 28 2 0 t h Jan., 1970 147 137 137 2 5 t h Jan., 1970 69 97 101 1 s t Feb., 1970 155 173 205 5 t h Feb., 1970 30 30 25 23 25 25 - a l l measurements are i n mm. - readings are not d i f f e r e n t i a t e d for.clear-weather or rain-on-snow melt Analysis of variance performed on Placid Lake lysimeter data using a completely randomized design Source of v a r i a t i o n Deg. of freedom Sum of squares Mean square F-ratio Among sample 2 889 4 4 4 . 5 . 0 5 3 6 Within sample 15 124448 8 2 9 6 . 6 Total 17 125337 At 5% significance l e v e l , F ( 2 , 1 5 ) = 3 , 6 8 > C a l c u l a t e d Or difference among lysimeters i s s t a t i s t i c a l l y i n s i g n i f i c a n t 42 arise from isolating the column of snow above the lysimeter, the open system principle was adopted for the construction of a self-recording lysimeter. The automatic lysimeter was installed at Placid Lake, not far from a Weatherwise tipping-bucket rain gauge and the Fuess hygrothermograph (see the last section) which supplied r a i n f a l l and temperature data against which snowmelt events could be correlated. Melt-water entered a plastic container which was kept three-quarters f u l l ( f i g . I I I . l ) . The reservoir in the container filtered out debris which otherwise could clog the outlet. Since water in the reservoir was always l e f t standing to the height of the outlet, any melt-water would drain down a buried hose to a storage vessel at the bottom of an insulated s t i l l i n g well. Water level in the storage vessel was then recorded by a Leupold-Stevens type-F recorder. When water in this vessel rose above the top of a siphon, automatic drainage was affected in about 8 hours 1. III-5 Field programme Most observations were taken during the period I969-71. Preliminary studies included biweekly snow surveys during winter and spring of 1969» and the selection of sites for snow, interception and snowmelt studies. Rainfall interception 1, The automatic siphon was designed by Mr, R, Gilbert who was inspired by a report of Goodell et a l . (1967). ^3 measurements were made during the summers of 1969 and 1970, while temperature relationships inside and outside the forest were observed during the 1969-70 water year. The winter of 1970 turned out to be too mild for snow studies, with most parts of the Research Forest devoid of snow before March. The following year, however, witnessed more snow-storms and provided ample opportunities for snow throughfall and snowmelt studies. During both winters (1969-71), snow surveys using Mount Rose snow sampler were conducted at 360 m, (road junct-ion s i t e ) , 370 m. (Loon Lake s i t e ) , 520 m. (Placid Lake site) and 535 m, (Gwendoline Lake si t e ) . Additional snow depth readings were taken with the aid of marked wooden dowels planted before the snow seasons. Figure III.2 summarizes the location of the instru-mentation sites. F i g . I I I . , ? H y d r o - m e t r i c and m e t e o r o l o g i c s t a t i o n s 31 t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a R e s e a r c h F o r e s t IV PRECIPITATION 45 P r e c i p i t a t i o n i n the forms of snow and r a i n c o n s t i t -utes the input to the water cycle. The d i s t r i b u t i o n of pre-c i p i t a t i o n i s highly variable both i n time and i n space. The problem of time d i s t r i b u t i o n of p r e c i p i t a t i o n w i l l be tackled by an analysis of p r e c i p i t a t i o n events at a 'base-station'. Subsequently, the s p a t i a l aspect w i l l be investigated by r e l a t i n g p r e c i p i t a t i o n events at the base-station to those at other parts of the basin. Temporal and s p a t i a l variations i n p r e c i p i t a t i o n are partly controlled by weather systems over the l i t t o r a l zone of southern B r i t i s h Columbia. A b r i e f description of a i r mass climatology of the west coast region f r e e l y abstracted from Kendrew and Kerr (1955) i s therefore included. IV-1 .0 Physical model Lying along the path of the westerlies, the coastal region of B r i t i s h Columbia receives the bulk of i t s p r e c i p i t a t -ion i n winter as various cyclonic fronts sweep inland from the P a c i f i c . From the south, polar front ( r e s u l t i n g from a meeting of the t r o p i c a l and polar a i r masses) sometimes moves up the coast to deposit p r e c i p i t a t i o n . In the north i s the maritime A r c t i c front which i s created between the maritime A r c t i c a i r mass and the more southerly maritime polar a i r mass. These maritime a i r masses produce most of the f r o n t a l p r e c i p i t a t i o n 46 events along the coast. Cold, dry s p e l l s , on the other hand, are caused by occasional breaks of the continental A r c t i c a i r mass whose source region i s northern Canada ( f i g . IV. 1 ) . In summer, cyclonic a c t i v i t i e s are much reduced. Shallow depressions, however, sometimes encroach the coastal areas as the polar front s h i f t s northward. Hot, dry s p e l l s a r r i v e with the northward spreading of continental t r o p i c a l a i r masses, though i n late summer and autumn, i n s t a b i l i t y may induce occasional thunderstorms. On the whole, dry periods are prolonged i n summer and heavy r a i n f a l l i s quite uncommon. IV - 2 . 0 Choice of a temporal model Various attempts have been made to model the stochastic nature of p r e c i p i t a t i o n at d i f f e r e n t time scales. Hourly p r e c i p i t a t i o n i s often given a Markov treatment"*" since the p r o b a b i l i t y of observing a p r e c i p i t a t i o n event i n the pre-sent hour i s conditioned by what happened i n the previous hour. Pattison ( 1964) went further to examine the conditional p r o b a b i l i t y structure as related to p r e c i p i t a t i o n conditions of the past s i x hours, and has thus generalized the approach to an n-th order Markov process. Gabriel and Neumann (1957 and 1962) have applied the Markov model to study wet and dry s p e l l s , using days as the basic time unit. Their d e f i n i t i o n 1. A stochastic process D ^ ] i s said to have the Markovian property i f P r [ x t + 1 = j [ x ^ k ^ x ^ ^ x t_ 1=k t_ 1,x t=i] = P r [ x t + 1 = j ) x t = i ] for a l l t = 0 , l , . . . and every sequence i t j t k Q, k-^ , k-t«i ( H i l l i e r and Lieberman, 1967 , p. 4 0 3 ) . F i g . I V . 1 S u r f a c e p r e s s u r e p a t t e r n a t 0600 PST, 3 1 s t D e c , 1968 s h o w i n g the c o n t i n e n t a l A r c t i c and m a r i t i m e A r c t i c f r o n t s and t h e a c c o m p a n y i n g p r e c i p i t a t i o n b e l t 48 of a cycle i s one which comprises several consecutive days that have recorded no r a i n (dry) or otherwise (wet). This approach was further pursued by Green (1964) who made use of the assumption that the lengths of wet and dry s p e l l s have exponential density functions. Other researchers have employed d i f f e r e n t types of density functions to characterize the s p e l l durations (table IV.1). Such d i s t r i b u t i o n s were either t h e o r e t i c a l l y derived from some assumptions about the nature of p r e c i p i t a t i o n (such as Todorovic and Yevjevich, 1969; Verschuren, 1968), or they were chosen because of f a i r agreement with observed data (Grace and Eagleson, I 9 6 7 1 Sariahmed and K i s i e l , 1968). On the whole, the choice of a p r e c i p i t a t i o n simulation model depends on: (1) the objectives of the study so that the s i m p l i c i t y of a model i s balanced against the demand for accuracy, (2) the basic time unit involved. The shorter the time scale, the more prevalent i s the e f f e c t of persistence, (3) physical interpretations of the model. S t a t i s t i c a l models with parameters not amenable to physical explanations are i n -f e r i o r models. In the present study, simulated p r e c i p i t a t i o n cons-t i t u t e s the input component. Any attempt to incorporate an over sophisticated component into the complex storage model i s not warranted i n terms of computing time. Besides, when the basic time unit used i n the study i s one day, natural periods of shorter durations such as cloud-bursts w i l i no Table IV.1 Summary of selected references on stochastic studies of p r e c i p i t a t i o n Author Technique Time unit L o c a l i t y On weather cycles Feyerherm & Bark, 1 9 6 7 Gabriel & Neumann, 1 9 5 7 Gabriel & Neumann, 1 9 6 2 Grace & Eagleson, 1 9 6 7 Green, 1 9 6 4 Green, 1 9 7 0 Hopkins & Robillard, 1 9 6 4 Longley, 1 9 5 3 Raudkivi & Lawgun, 1 9 7 0 Sariahmed & K i s i e l , I 9 6 8 Simpson & Henry, I 9 6 6 Second order Markov chain F i r s t order Markov chain day day F i r s t order Markov chain leading day to geometric d i s t r i b u t i o n Weibull d i s t r i b u t i o n for storm dur- minute ation and time between storms Alternating wet and dry runs follow- day ing renewal process Modified Markov chain Second order Markov chain F i r s t order Markov chain Weibull d i s t r i b u t i o n Weibull d i s t r i b u t i o n F i r s t order Markov chain day day day minute storm day N. U.S.A. Isr a e l I s r a e l NE U.S.A. Is r a e l England Canadian p r a i r i e s Canada, U.K. New Zealand SW U.S.A. Canadian p r a i r i e s Table IV.1 (continued) Author Technique Time Unit L o c a l i t y Todorovic & Yevjevich, I 9 6 9 Poisson d i s t r i b u t i o n for number of hour, day storms per period Verschuren, I968 Woo, 197 2 Poisson d i s t r i b u t i o n f o r number of storms per period day F i r s t and second order Markov chains day and exponential d i s t r i b u t i o n for s p e l l durations On quantity of p r e c i p i t a t i o n Bagley, 1964 Chow & Ramaseshan Grayman & Eagleson, 3967 Modified Poisson and geometric dis tributions Markov process Linear relationship with storm duration (Weibull) - day Gamma d i s t r i b u t i o n Sixth order Markov chain Ison et a l . , 1971 Pattison, 1 9 6 4 Todorovic & Yevjevich, I 9 6 9 Gamma d i s t r i b u t i o n Verschuren, I 9 6 8 Gamma d i s t r i b u t i o n hour storm storm hour storm storm C U.S.A. Colorado SW Canada W. U.S.A. SE U.S.A. NE U.S.A. C U.S.A. C a l i f o r n i a C U.S.A. Colorado 51 longer be preserved by the model. Some model more robust i n character has to be formulated from basic assumptions about the behaviour of d a i l y p r e c i p i t a t i o n . IV - 2 . 1 Mathematical model for temporal d i s t r i b u t i o n of  p r e c i p i t a t i o n Records of d a i l y p r e c i p i t a t i o n are t o t a l l y defined by the durations of wet and dry s p e l l s and the quantities of d a i l y p r e c i p i t a t i o n during a wet period. Over the year, these elements vary with the changing atmospheric c i r c u l a t i o n system. To develop a mathematical model for p r e c i p i t a t i o n , various factors such as s p e l l durations, p r e c i p i t a t i o n magnitudes and seasonality have to be considered. I V - 2 . 2 . 1 P r e c i p i t a t i o n trend and seasonality A trend i s defined as a "u n i d i r e c t i o n a l diminishing or increasing change i n the average value of a (hydrologic) variable" (Chow, 1964, p. 8 - 1 0 ) . This follows from the d e f i n i t i o n of s t a t i o n a r i t y which states that a stochastic process [ x(t), jt| < °°] i s considered to be ( s t r i c t l y ) stationary i f absolute time plays no role i n the sense that F t l t t 2 , t n ( x 1 , x 2, x n) ^tj+s, t 2+s, t n + s ( x l t x 2, x ) (IV.l) fo r a l l f i n i t e values of n and s ( K a r t v e l i s h v i l i , 1 9 6 9 , p. 5 6 ; Lamperti, 1966, p. 1 1 5 ) . When a trend e x i s t s , the 52 process i s said to be non-stationary or time-inhomogeneous 1. It i s not easy to v e r i f y that climatic or hydrologic processes are s t r i c t l y stationary, partly because of the lack of obser-vations over a s u f f i c i e n t l y extended period of time. On the other hand, process i s said to be weakly stationary i f there exists some value k , such that E [ x ( t ) J = constant E [ x 2 ( t ) ] < °° E [x(t )X(t+k)] = f(k) (IV.2) Should a trend component exist i n the p r e c i p i t a t i o n 2 record , i t i s necessary to calculate the magnitude of departure from the long-term mean conditions. For records of short duration, however, evidence i s i n s u f f i c i e n t to introduce the trend component and any short-term deviations from the long-term mean are regarded as merely fluctuations of a long climatic sequence (Curry, 1962). Seasonality, on the other hand, gives r i s e to a rhythmic v a r i a t i o n of p r e c i p i t a t i o n during the year. The annual march of monthly p r e c i p i t a t i o n has been investigated 1. Yevjevich and Jeng (1969) distinguished between inconsis-tency and non-homogeneity. The former appears as jumps i n the record such as when dams are b u i l t across a r i v e r . For non-homogeneous records, a trend occurs such as when increased water consumption leads to a depletion of water supplies. 2. Systematic deviations from long-term averages may be attributed to natural causes such as sunspot a c t i v i t i e s . Also, a trend may be induced by human a c t i v i t i e s such as weather modification or the increase of condensation nuclei near i n d u s t r i a l centres (Hobbs et a l , , 1970). 5 3 through harmonic analysis (Sabbagh and Bryson, 1 9 6 2 ) . Hence, to portray the p r e c i p i t a t i o n of a station, i t s long-term monthly mean or seasonal mean values can be used to characterize the e f f e c t of seasonality. I V - 2 . 2 . 2 Wet and dry s p e l l s For the southwestern B r i t i s h Columbian coastal environment, the occurrence of wet or dry days i s affected by past conditions so that a Markov chain model i s appropriate. A Markov chain i s characterized by i t s t r a n s i t i o n matrix and the i n i t i a l d i s t r i b u t i o n functions. In p a r t i c u l a r , a Markov chain with ergodic states w i l l be aperiodic and has f i n i t e recurrence time ( F e l l e r , 1 9 6 8 , p. 3 9 3 ) . With wet and dry periods, three assumptions are made: ( 1 ) The prob a b i l i t y of a period being wet or dry i s conditioned by the nature of the past period or a f i n i t e number of past-periods. This assumption establishes the Markov property. ( 2 ) The t r a n s i t i o n from a dry to a wet period (and vice versa) follows some probability which does not change with time except for the e f f e c t of seasonality. This assumption i s made to s a t i s f y s t a t i o n a r i t y i n the weak sense. ( 3 ) Neither wet nor dry periods w i l l extend to i n f i n i t e lengths. This and the above assumption establish the ergodic condition of the Markov chain. From h i s t o r i c a l records, a t r a n s i t i o n matrix can be constructed by noting the frequencies at which a wet day i s followed by a wet or dry day, and so on. The problem i s 54 how many previous days should be retraced before the e f f e c t of persistence can be ignored. The answer to thi s problem cannot be e a s i l y obtained from physical climatology. Rather, the problem i s best s e t t l e d by s c r u t i n i z i n g the s t a t i s t i c a l good* ness of f i t s between h i s t o r i c a l data and models of f i r s t and second order Markov chains, Markov chain models of higher order w i l l not be considered mainly because of the want of physical v e r i f i c a t i o n of a persistence e f f e c t that would l a s t for more than two days under a west coast environment. I V - 2 , 2 , 3 Quantity of d a i l y p r e c i p i t a t i o n By f i x i n g the time of p r e c i p i t a t i o n gauge reading, the basic time unit during which p r e c i p i t a t i o n i s recorded i s dis-synchronised with the int e r v a l s of continuous p r e c i p i -t a t i o n . Hence, the amount of p r e c i p i t a t i o n recorded per day i s d i s t i n c t from p r e c i p i t a t i o n amount per storm. This d i s -location of time bases for p r e c i p i t a t i o n record and for storm events i s contributive to the independence i n d a i l y p r e c i p i -t a t i o n . The time at which gauge-reading i s taken may occur during any part of the storm so that i n terms of the quantity of p r e c i p i t a t i o n gathered, the a r r i v a l of a gauge reading event i s t o t a l l y random. Once the gauge i s emptied, the quantity of p r e c i p i t a t i o n i s reset to zero and the amount col l e c t e d during the next period w i l l therefore be independent of the amount obtained i n the past. Let (X J " 1 be the mean quantity of p r e c i p i t a t i o n c o l l e c t e d when the gauge i s read, X w i l l be the mean 5 5 p r o b a b i l i t y of terminating the act of p r e c i p i t a t i o n c o l l e c t i o n ( i . e . the event of taking a gauge reading). Then, between the points when the gauge has gathered x and X + A X units of pr e c i p i t a t i o n , with A X being a small amount of incremental p r e c i p i t a t i o n , the p r o b a b i l i t y of i t s not being read w i l l be Pr Q(x+Ax) = P r Q ( x ) ( l - X pAx) (IV.3) where Pr Q(x) i s the pr o b a b i l i t y that the gauge was not read when i t had coll e c t e d up to x units of p r e c i p i t a t i o n . Re-arranging terms and l e t t i n g A X approach zero, we obtain the d i f f e r e n t i a l equation lim Pr 0(x+*x) - Pr Q(x) = p . ( x ) A X - * 0 A X = -X pPr 0(x) (IV.4) Since the gauge i s empty at the beginning of a c o l l e c t i o n period, the i n i t i a l condition w i l l be P r Q ( 0 ) = 1 , Then, equation (IV.4) w i l l have the solution Pr Q(x) = exp(-X px) (IV.5) I f X i s the cumulative quantity of p r e c i p i t a t i o n collected up to the a r r i v a l of a gauge-reading event, the pr o b a b i l i t y of the gauge not being read when i t has gathered up to x units of p r e c i p i t a t i o n (where x<X) w i l l be Pr Q(X>x) = Pr Q(x) = exp(-X px) , x>0 (IV.5.a) Therefore, the amount of p r e c i p i t a t i o n c o l l e c t e d up to the time when the reading i s taken w i l l have the d i s t r i b u t i o n function F(x) = 1 - exp(-X px) , x$0 (IV.6) 56 with i t s density function f(x) = F'(x) = X pexp(-X px) (IV.?) which shows that the amount of p r e c i p i t a t i o n c o l l e c t e d between two gauge-reading events follows an exponential d i s t r i b u t i o n . The f i r s t and second moments for the exponential d i s t r i b u t i o n are E [ x J = £ xX pexp(-X px)dx = U p ) ~ x (IV.8) and oo VAR = S (x - E [ x]) 2X pexp(-X px)dx = (-X p)~ 2 = ( E [ X ] ) 2 (IV . 9 ) where E [ X ] i s expectation VAR i s variance or second moment about the mean. The parameter X i s not constant throughout the year as i s witnessed by seasonal changes i n mean d a i l y p r e c i p i t a t i o n . Therefore, monthly values of X p w i l l be evaluated ( f i g . I V . 3 ) . IV - 2 . 3 Result of analysis of h i s t o r i c a l data Daily p r e c i p i t a t i o n records f o r the period 1957-70 were available for the weather sta t i o n near the Administration building of the University of B r i t i s h Columbia Research Forest. A computer programme was written i n FORTRAN IV language to analyse such data. The flow-chart for t h i s programme i s presented as figure IV.2, IV - 2 . 3 . 1 Trend Since t h i r t y years i s generally recognized by the Start ye&t counter / Input (by year) precipitation I n i t i a l i z e month counter Init tA.lt re day counter find weather conditions Ct-.| of day i-t Ci-2 ot day t-2 Sum pcecip into tn daily itation enth total tncre> courtte PPCM>|( went r for Increment counter for Increment year counter yes Increment tnonth counter I _L Increment dau counter No J e s . 57 yes Find mean daily precipitation 6y month Find conditional probabilities af occurrences ot wet days Output transition matrices And mean oldili/ precipitation i g . IV.2 F l o w - c h a r t f o r p r e c i p i t a t i o n a n a l y s i s DEC Ap -072 1 ~i—i—i—t r- r r—r—i i i "i I i i v | , | t f ! i i i •8-, O <V APR Ap .103 •8-r 0 l\k i — I — i — i ~ i — r JUL MAV Ap -125 XP .116 i * i i n i i i i i i v AUG Ap .120 I I I I I — 1 — I — I — I 20 AO 60 so D A U U precipitation (mm.) F i g . IV.3 T h e o r e t i c a l „f- , a n d e m p i r i c a l p r o b a b i l i t y d i s t r i b u t i o n s o i d a x l y p r e c i p i t a t i o n bv month S 59 Atmospheric Environment Service as the standard period on which c l i m a t i c s t a t i s t i c s should be based, i t was necessary to compare the Research Forest records with those from some nearby s t a t i o n with longer h i s t o r y . The station at Stave F a l l s , B.C. ( 1 2 2 ° 21* W; 49° 14' Nj elevation 55 m.) was chosen f o r the period 1 9 2 8 - 6 9 . The record length was s u f f i c -i e n t l y long to enable an analysis of p r e c i p i t a t i o n trend by autocorrelation of the t o t a l annual p r e c i p i t a t i o n s e r i e s . A lag of up to four periods was used i n the c a l c u l a t i o n of autocorrelation c o e f f i c i e n t s (equation V . 2 ) . None of these c o e f f i c i e n t s were found to be s t a t i s t i c a l l y s i g n i f i c a n t at the significance l e v e l (equation V,**-)"*". For p r a c t i c a l purposes, therefore, annual p r e c i p i t a t i o n amounts are considered to be independent of one another. IV - 2 . 3 . 2 Spells For the Markov chain models, t r a n s i t i o n matrices were obtained from h i s t o r i c a l data of Stave F a l l s . Based on these matrices (table I V . 2 ) , one hundred years of p r e c i p i t a t i o n data were simulated with a f i r s t order Markov chain model, and another hundred years of data with a second order Markov chain model. From these sets of data, the frequency d i s t r i b u t i o n s of wet and dry s p e l l s of various durations were abstracted. 1, The autocorrelation c o e f f i c i e n t s ( r k ) are: order k = 1 2 3 4 r k = . 2 ? 0 .122 . 2 3 1 . 0 9 3 Confidence l i m i t s f o r r-^ are - 0 . 3 3 9 and 0 , 2 8 9 at 9 5 $ confidence l e v e l . Table IV,2 Transition matrices for the occurrence of wet days, and mean d a i l y p r e c i p i t a t i o n by month Stave F a l l s 0 N D J F M A M J J A S p ( w j w ) P ( W i D ) P ( W J W W ) P ( W J W D ) P(W1DW) P ( W I D D ) . 7 2 7 .361 . 7 3 4 .708 .443 .318 .800 .395 .802 .793 . 4 7 0 .348 .829 . 4 3 4 .827 .838 . 5 ^ .352 .796 .333 .801 .772 .385 . 3 0 4 .781 .3^6 ,788 .759 . 4 5 0 .291 .763 .352 .767 .749 .426 .313 .718 .3^8 .717 . 7 2 3 .394 . 3 2 3 . 6 8 9 .264 .663 .746 . 3 1 0 .248 .671 . 2 7 4 . 6 7 7 .659 .3^1 .248 . 6 0 2 .146 . 5 3 4 .706 . 2 0 3 .136 . 6 3 0 .182 .618 .649 . 2 3 1 .172 .654 .231 . 6 3 0 . 6 9 9 .262 .221 Mean d a i l y precips (inm) 14.0 12 ,3 1 3 . 3 13 .9 12.3 10 .5 8.6 8 . 2 7.7 7.4 7.4 1 0 . 9 UBC Admin. 0 N D J F M A M J J A s p ( w l w ) P ( W ! D ) P ( W J W W ) P ( W l W D ) P ( W l D W ) P ( W l D D ) . 744 .363 .753 .717 . 4 0 7 .339 .820 .432 .846 .712 .510 .373 . 8 5 7 . 3 6 4 .849 .902 .342 .377 .811 . 3 8 9 .792 . 8 9 6 . 5 6 0 .284 . 7 8 7 . 3 7 5 . 8 0 2 . 7 3 1 . 4 7 2 .313 .772 .277 .763 .800 .373 .241 . 711 . 3 9 5 . 6 8 9 . 7 6 1 . 4 1 3 . 3 8 4 . 6 6 9 . 2 5 0 . 6 7 5 . 6 5 5 . 2 6 7 .244 . 6 0 8 . 2 5 7 . 6 0 2 .617 . 3 0 7 .240 .619 . 109 . 5 6 5 . 7 1 4 .189 . 0 9 8 . 6 6 2 .182 . 6 0 5 . 7 6 6 . 2 3 3 .172 . 6 5 0 . 2 3 2 . 6 5 4 .643 . 3 1 5 . 2 0 7 Mean d a i l y precip. (mm) 1 3 . 9 1 3 . 1 1 3 . 9 14 .6 1 2 . 3 1 0 . 7 9 . 7 8 . 1 6 . 6 8 . 6 8 . 3 1 1 . 5 61 A two-sample Kolmogorov-Smirnov test applied to each pair of simulated and h i s t o r i c a l data sets shows that a second order Markov chain model reproduces successfully the s p e l l durations for a l l four seasons (table I V . 3 ) . This model w i l l be adopted for the present study. I V - 2 . 3 . 3 Daily p r e c i p i t a t i o n To establish the independence of d a i l y p r e c i p i t a t i o n , wet days from the Research Forest for the 1968-70 water years were analysed for autocorrelation. Autocorrelation c o e f f i c i e n t for a lag of one day was found to be O . I 6 5 at 285 degrees of freedom. This i s much lower than the value quoted by Kotz and Neumann ( 1 9 5 9 ) . According to Brooks and Carruthers ( 1 9 5 3 , p. 3 2 4 ) , " i n random samples from uncorrelated populations, one out of 24 c o e f f i c i e n t s would be expected to exceed 0 . 2 " , The degree of co r r e l a t i o n between p r e c i p i t a t i o n of consecutive days at the Research Forest i s thus attributed to chance, IV-3 . 0 Spatial variations of p r e c i p i t a t i o n A gross pattern of p r e c i p i t a t i o n i n the Lower Fraser v a l l e y i s one with steep gradient r i s i n g towards the mountains i n the north. This s i m p l i f i e d isohyetal surface i s indented as the isohyets encounter the valleys and uneven mountain slopes (Wright, 1 9 6 6 ) . Detailed mapping of isohyets cannot be made because of i n s u f f i c i e n t weather stations to match the complexity of t e r r a i n and vegetation cover. Consequently, some in d i r e c t measure has to be found to estimate the incoming quantity of water to mountainous basins. Table IV 0 3 Two-sample Kolmogorov-Smirnov tests on the d i s t r i b u t i o n of s p e l l lengths D - s t a t i s t i c s for dry s p e l l s 1 s t order Markov chain 2nd order Markov chain Oct-Dec Jan-Mar Apr-Jun Jul-Sep Oct-Dec Jan-Mar Apr-Jun Jul-Sep 2nd order Markov chain . 0 7 6 * .081* . 0 3 6 . 0 3 4 Stave F a l l s . 0 7 8 * . 1 0 1 * . 0 5 8 . 0 3 6 .019 . 0 6 6 . 0 5 0 . 0 2 9 No. of samples 1272 1216 1335 1132 1250 1197 1329 1040 Number of samples for Stave F a l l s s Oct-Dec 507 Jan-Mar 501 Apr-Jun 562 Jul-Sep 464 D - s t a t i s t i c s for wet sp e l l s 1 s t order Markov chain 2nd order Markov chain Oct-Dec Jan-Mar Apr-Jun Jul-Sep Oct-Dec Jan-Mar Apr-Jun Jul-Sep 2nd order Markov chain .017 .024 . 0 1 5 . 0 2 7 Stave F a l l s . 0 3 5 . 0 2 2 . 0 5 0 . 0 4 7 . 0 3 1 . 0 2 7 . 0 6 5 . 0 2 0 No. of samples 1303 1205 1324 1123 1279 1194 1304 1137 Number of samples for Stave F a l l s : Oct-Dec Jan-Mar Apr-Jun Jul-Sep * D - s t a t i s t i c i s s i g n i f i c a n t at 5% l e v e l 525 494 551 464 63 J u l i a n et a l . ( 1967) have studied p r e c i p i t a t i o n variations i n the mountainous areas of western United States i n terms of meteorological, physiographic and human factors. To t h i s l i s t should he added the e f f e c t s of vegetation cover. For most of the remote basins, however, human influences are less prominent since effects of weather modification or the increase of i n d u s t r i a l pollutants (as condensation nuclei) have not been demonstrated for the entire study area. On the other hand, the influence of storm types on s p a t i a l d i s -t r i b u t i o n of p r e c i p i t a t i o n i s more evident. On f l a t t e r r a i n of I l l i n o i s , Huff and Shipp ( 1969) found that s p a t i a l c o r r e l -ation decays with distance between stations, and i s greatest i n thunderstorms, rain-showers and a i r mass storms, while minimum decay occurs with steady r a i n and the passage of low pressure centres. In more rugged t e r r a i n , physiography exerts a stronger influence through interactions with vegetation and weather conditions. Such interactions amongst several elements tend to produce a high degree of complexity i n s p a t i a l d i s t r i -bution of p r e c i p i t a t i o n . The following section w i l l discuss interactions between vegetation and various forms of p r e c i p i -t a t i o n . IV - 3 . 1 . 1 Physiographic considerations -physical model Elevation, slope angle, orientation and exposure are regarded as major elements which a f f e c t the d i s t r i b u t i o n of p r e c i p i t a t i o n (Linsley, 1 9 5 6 ; Schermerhorn, 1 9 6 7 ; Spreen, 1 9 ^ 7 ) . In r e l a t i o n to wind d i r e c t i o n at 500 mb l e v e l , for 64 example, one side of the v a l l e y may receive more r a i n than the slope across the valley, while the pattern i s reversed when upper wind moves i n the opposite d i r e c t i o n (Urfer-Henneberger, 1 9 7 0 ) . When r e l i e f i s not prominent, v a l l e y and ridge stations may have s i m i l a r quantity of catch (Dickison, 1 9 6 8 ; Sartz, 1 9 6 6 ) , otherwise, orographic e f f e c t i s s i g n i f i -cant (Sporns, 1 9 6 4 ) , The general increase of p r e c i p i t a t i o n with elevation does not hold for the positioning of the zone of maximum p r e c i p i t a t i o n whose pos i t i o n varies with storm types (Heigel, I 9 6 0 ; Hovind, 19655 Mink, i 9 6 0 ) . In northern Utah, storms associated with cold-centred lows at 500 mb l e v e l were found to y i e l d p r e c i p i t a t i o n f a i r l y evenly over d i f f e r e n t elevation zones. On the other hand, cold fronts or occluded fronts tend to produce marked precip-i t a t i o n contrasts between v a l l e y stations and ridge stations (Williams and Peck, 1 9 6 2 ) . Such complex interactions between physiography and storm types cannot be applied to the present study because of the a r b i t r a r y time base inherent i n the present study. Thus, a weaker physical model i s sought. It i s postulated that the amounts of p r e c i p i t a t i o n recorded at neighbouring stations within a small basin are i n t e r r e l a t e d i f a p a r t i c u l a r weather system passes over the basin within a c e r t a i n period. Then, i t would be appropriate to obtain some transfer functions with which p r e c i p i t a t i o n records from a base-station are extended to other parts of the basin. 65 IV-3.1.2 Physiographic considerations -mathematical model and s t a t i s t i c a l analysis Consider a storm W occurring over two nearby stations i and j during period ( t , t+M). The amount of p r e c i p i t a t i o n recorded at both stations may be written as x i(At) = f[w(At), H] x.(At) = f[w(At), HJ , t«T (IV.10) J where H i s a physiographic variable T i s t o t a l duration of storm W Since p r e c i p i t a t i o n at sta t i o n i r e f l e c t s the storm c h a r a c t e r i s t i c s and since p r e c i p i t a t i o n at s t a t i o n j i s derived from the same storm, i t i s possible to express the quantity recorded at j during time ( t , t+«vt) i n terms of p r e c i p i t a t i o n at sta t i o n i Xj(At) = g [ X i ( A t ) , H] (IV.11) In other words, i f the quantity of p r e c i p i t a t i o n for some base-station i i s known, the amount received by an unknown station j i n the same basin w i l l be estimated by the function g . To investigate thoroughly the e f f e c t of various topographic factors on p r e c i p i t a t i o n requires a very dense gauging network and i s much beyond the scope of the present study. Instead, elevation i s chosen as the major variable which influences d a i l y p r e c i p i t a t i o n i n the mountains at the Research Forest. Equation IV.11 i s modified to the form of: x-(At) = g 1 [ x i ( A t ) , h] (IV.11.a) where h i s elevation. 66 For the analysis, data obtained for the period 1967-68 from six stations at or near the Research Forest (see section I I I - 2 ) were used. Daily p r e c i p i t a t i o n was subdivided according to the amount recorded at the v a l l e y s t a t i o n . By least-square f i t t i n g , a family of equations of the following form was obtained t x. = a. + b,h. . x. > 0 (IV.12) where a and b are c o e f f i c i e n t s to be used when stat i o n i records p r e c i p i t a t i o n of magnitude k , h. . i s a l t i t u d i n a l difference between stations i and j , with i being the base-station. Results of the analysis (table IV.4 and f i g . IV.4) show l i t t l e c o r r e l a t i o n between p r e c i p i t a t i o n and elevation when the amount recorded at the valley s t a t i o n i s below 10 mm. As d a i l y p r e c i p i t a t i o n at the valley s t a t i o n increases, orographic e f f e c t becomes more marked. This i s i n accordance with Walkotten and Patric's (1967) finding near H o l l i s , Alaska. When the v a l l e y station r e g i s t e r s over 80 mm., however, the influence of elevation becomes reduced. It i s uncertain whether l i n e a r r e l a t i o n s h i p between elevation and p r e c i p i t a t i o n would hold for higher a l t i t u d e s . Caution must therefore be exercised when these equations are used f o r extrapolation. IV -3 .2 E f f e c t s of vegetation on p r e c i p i t a t i o n The presence of vegetation cover on the land surface tends to reduce the quantity of p r e c i p i t a t i o n that reaches the Table IV.4 Least-square f i t . d a i l y p r e c i p i t a t i o n vs a l t i t u d e , University of B r i t i s h Columbia Research Forest P r e c i p i t a t i o n at base-station (mm) Regression c o e f f i c i e n t Correlation c o e f f i c i e n t Degrees of freedom Standard error (mm) 0-10 .0018 .076 263 4 . 7 1 0 - 2 0 . 0 2 5 1 . 2 7 1 95 18 .3 20-30 . 0 3 7 7 .377 71 1 9 . 1 30-40 . 0 7 9 1 O608 29 2 1 . 8 40-50 .1240 . 9 6 0 5 9 . 0 5 0 - 7 0 .1465 . 8 8 9 5 19.8 9 0 - 1 0 0 . 0 6 9 1 . 8 7 5 5 9 . 5 120-130 .0757 .916 5 8 . 2 Data fromi Haney 532 UBC Admin UBC Forest (Marc) Loon Lake Lower M i l i t z a Middle M i l i t z a 1967-68 water year 68 0 100 200 X0 400 , ALTITUDE (IN M.) Fig. IV.4 P r e c i p i t a t i o n - a l t i t u d e relationship, Universit of B r i t i s h Columbia Research Forest, 1967-68 69 s o i l or snowpack surface. Vegetation influences on r a i n f a l l d i f f e r from that on snowfall so that they w i l l "be considered separately. Because the entire basin i s the basic unit of t h i s study, i t i s deemed s u f f i c i e n t to treat the basin surfaces as either forested or open rather than subdividing further i n accordance with forest density or tree age. I V - 3 . 2 . 1 Forest canopy and r a i n f a l l - physical model Rain f a l l i n g on a forest canopy i s partly stored on the surfaces of leaves and branches to be subsequently evapor-ated, and p a r t l y passed on to the forest f l o o r . The quantity of r a i n f a l l stored and evaporated i s known as interception loss while the amount which reaches the f l o o r d i r e c t l y or through dripping i s throughfall. The quantity of stemflow has been found to be n e g l i g i b l e f o r hardwoods (Rogerson and Byrnes, 19681 Rouse, I 9 6 5 ) as well as for pine (Rowe and Hendrix, 1 9 5 1 ; Thorud, 1 9 6 3 ) , Norway spruce (Eidmann, 1 9 5 9 ) or Douglas f i r stands (Rothacher, I 9 6 3 ) . Stemflow for conifers i n general i s estimated to range between 0 . 5 to 3*5$ of net r a i n f a l l (Thorud, 1 9 6 3 ) . so that for p r a c t i c a l purposes, i t i s excluded from the present discussion. The rate of interception during a storm i s not constant, A certain amount of r a i n f a l l i s required to saturate the l e a f and twig surfaces before dripping begins (Leonard, I 9 6 7 ) . A steady state condition i s reached when evaporative loss does not increase with increasing supply of moisture (Rowe and Hendrix, 1 9 5 1 ) . From then on, the proportion of 70 throughfall w i l l tend to be constant. Although transpiration rate i s reduced while the canopy i s wet (Harr, 1 9 6 6 ) , the reduction i s probably only moderately compensatory on interception l o s s . Seasonal variations i n throughfall are not apparent for hardwood forests i n northeastern parts of the United States (Coffay, 1962j Leonard, 1 9 6 1 ) . On the other hand, Rothacher (1963) obtained a higher r a t i o of winter throughfall i n Oregon because of i t s b r i e f inter-storm periods have been i n s u f f i c i e n t to dry out the foliage of Douglas f i r s so that the canopy i s maintained at a moist l e v e l . Pruning can also aff e c t the quantity of throughfall. Thorud's ( I 9 6 3 ) study, however, showed that pruning of pine i n Minnesota s i g n i f i c a n t l y increased throughfall only when t o t a l p r e c i p i t a t i o n i n the open was less than 150 mm. Nor i s the age of trees a good indicator of throughfall. This i s borne out by re-working Cartwright's ( 1970) throughfall data from Seymour basin near Vancouver, B.C., for a 50 year old Douglas fir-hemlock-cedar s i t e and a 9 year old replanted Douglas fir-hemlock s i t e ( f i g . IV.5) 1. Throughfall was adjusted according to r a i n f a l l 1, Covariance analysis of the two sets of data gives the following r e s u l t s J 50 year old s i t e R f = -0 .5 + 0.80Ro 9 year old s i t e R f = -0.2 + 0.81RQ Pooled data R f = -0 .3 + 0.81RQ corr. coef. i s . 9 9 8 at 3 7 d.f. R Q and R f are r a i n f a l l i n the open and the forest Neither the slope nor the intercept was s i g n i f i c a n t l y d i f f e r e n t . 71 120 F i g . I V . 5 T h r o u p h f a l l - r a i n f a l l r e l a t i o n s h i n f r o m two s i t e s a t Seymour B a s i n , B.C. 72 at opening near these two s i t e s . A t - t e s t performed on adjusted throughfall at these s i t e s shows that t h e i r difference i s s t a t i s t i c a l l y i n s i g n i f i c a n t . Hence, indices other than maturity of the forests are needed to estimate throughfall. From the above example, i t i s surmised that changes i n precipitation-vegetation interactions, rather than changes i n vegetation c h a r a c t e r i s t i c s alone, w i l l play the important r o l e i n throughfall determination. To look at the water balance of a forested area at another scale, measurements i n small forest clearings have shown some excess over measurements from open f i e l d s nearby (Fedorov and Burov, 1 9 6 7 ) . Higher p r e c i p i t a t i o n over the forested areas has been attributed to an increase i n dynamic roughness as the moisture-laden a i r flows over the fore s t . Krecmer ( 1967) has demonstrated a 100$ catch of r a i n f a l l i n a c i r c u l a r opening with a diameter of 0 . 7 times the average height of trees. These findings suggest that the measurement of incoming p r e c i p i t a t i o n should best be made i n forest openings, and that the distance of the gauge need not be more than h a l f a tree-height away from the margin of the openings. IV - 3 . 2 .2 Forest canopy and r a i n f a l l - s t a t i s t i c a l analysis The present study does not attempt to handle changes i n water balance as affected by forest types or tree si z e , so that i t i s reasonable to pool data from both s i t e s at 230 m. and 520 m, elevations (see section III-2) to obtain a more representative estimate of the variance. I n i t i a l l y , an assertion was made that seasonality may influence the r e s u l t s . 73 This assertion was not substantiated: covariance analysis showed that neither the slopes nor the intercepts of the least-square l i n e s f i t t e d to data pooled by three-month periods were s t a t i s t i c a l l y s i g n i f i c a n t (table IV.5 ) . Hence, a single equation of the following form was accepted as applicable throughout the year ( f i g . IV,6 ) : R f = a + bRQ, R f £ 0 (IV.13) where R f and R Q are r a i n f a l l inside and outside the forest a and b are constants. The throughfall r a t i o of b=0.74 i s i n keeping with findings of most other researchers (e.g. Eidmann, 1959? Rothacher, 1 9 6 3 ) . A positive intercept of a>0 , however, requires an awkward physical explanation that the forest receives more r a i n f a l l than the opening. On the other hand, i t i s more l i k e l y that measurement error and the lack of l i n e a r i t y towards the lower end of the equation have contributed to a p o s i t i v e intercept. In view of i t s small numerical value, t h i s term i s taken to be zero instead. I V - 3 . 2 . 3 Forest canopy and snowfall - physical model Unlike r a i n f a l l which stays i n the basin only for a b r i e f period, snowfall i s often stored i n the basin as part of the pack. Where forest canopy e x i s t s , the amount of snow p r e c i p i t a t i o n reaching the ground surface i s modified by snow interception. Intercepted snow i s a highly transient feature, part of which w i l l ultimately be evaporated, and the rest re-deposited onto the pack. Many studies have been conducted to evaluate the quantity, the time and s p a t i a l d i s t r i b u t i o n s Table IV.5 Covariance analysis of throughfall (dependent) vs r a i n f a l l i n open (independent variable) Season Slope Intercept Degrees of freedom Correlation c o e f f i c i e n t Standard error (mm) Oct-Dec .696 2 . 5 26 . 9 5 9 9 9 . 2 2 Jan-Mar ,748 0 . 5 13 . 9 8 9 9 4 . 3 0 Apr-Jun . 7 0 7 - 0 . 6 13 . 9 5 3 7 5 . 8 9 Jul-Sep . 8 0 5 - 0 . 9 28 .9711 6 . 7 1 Pooled .736 0 . 5 83 . 9 6 8 3 7 . 2 7 Test for significance of slopes F 3 / 7 6 = 5 2 9 Test for significance of l e v e l s F 3 / 7 9 = 0 , 7 7 1 Test for a l l data F 6 / 7 6 = 1 ' 1 5 8 None of these F-values suggest s t a t i s t i c a l s i g nificance 75 200 -§ z 120 -I a: rs o a: 40 ^ CD 8 .X o AUTUMN A VINTER + SPRING x SUMMER r — 1 r 40 80 120 ,60 RAINFALL IN OPEN (IN MM) 200 Fi,< I V . 6 T h r o u p h f a l l - r a i n f a l l r e l a t i o n s h i p , U n i v e r s i t y o f B r i t i s h C o l u m b i a R e s e a r c h F o r e s t l v e r s i t y 76 of intercepted snow. Their measurements either concern the snow on tree crowns or the snow which f i n a l l y reaches the forest f l o o r (table IV.6). To estimate the quantity of p r e c i p i t a t i o n entering a forest, i t i s more desirable to measure the pack underneath the canopy rather than to measure the intercepted snow on the crowns. I t was hypothesized that snow reaching the forest f l o o r i s rela t e d to snow-bearing storms of which the snowfall amount i s indicated by measurements made i n the opening nearby"*". The assumptions were that d r i f t i s n e g l i g i b l e i n a closed canopy and that l i t t l e melting occurred during storms. To ensure l i t t l e loss of fresh snow from melt, snow measure-ments were made soon aft e r storms. IV-3.2.4 Forest canopy and snowfall - mathematical model  and r e s u l t of analysis Studies on snow catch by coniferous crowns (Satterlund and Haupt, 196?; Watanabe and Ozeki, 1964) have shown crown interception to follow logarithmic 'growth' curves which tend asymptotically to steady state conditions as snowfall continues. Thus, one would expect the snow throughfall to increase 1. M i l l e r ' s (1966) review shows that the value of A A i s unstable, where A A i s difference i n snow measurement inside and outside the canopy. On the other hand, most of his A A ' s were calculated from snow depths and an assumed density of 0.1 and thus may be erroneous. Nevertheless, the present study i s not interested i n the values of A A since we are more concerned with actual quantities rather than the amounts of differences , Table IV.6 Selected references on snow interception and forest snow accumulation studies Author Description Interception approach Berndt & Fowler, 1969 NW U.S.A. Goodell, 1959 Colorado Hoover & Leaf, 1 9 ° 7 Colorado Japanese Government, 1952 Japan L u l l & Rushmore, 1 9 6 l a NE U.S.A. M i l l e r , 1964 review M i l l e r , 1966 review Miner & Trappe, 1957 NW U.S.A. Satterlund & Eschner, I 9 6 5 review Satterlund & Haupt, 1967 Idaho hoar fro s t and rime as negative forms of i n t e r -ception rapid evaporation of intercepted snow afte r storm strong wind may not dislodge a l l intercepted snow snow interception proceeds by bridging across adjacent branches hoar f r o s t formed on exposed intercepted snow surfaces on cold, calm nights snow accumulation re l a t e d to windspeed and geo-metry and closure of tree crowns snow interception r e l a t e d to a i r temperature and windspeed quantity of interception depends on tree type snow adhesion due to refreezing of snow meltwater so that wet snow i s more e a s i l y intercepted further accumulation due to cohesion of snow i n s i g n i f i c a n t evaporation during storms dense thicket as under-story of forest produces secondary interception differences i n energy and vapour balance between intercepted snow surface and snowpack surface storage capacity of intercepted snow d i f f e r s with species cumulative amount of intercepted snow follows some 'growth* curve Table IV,6 (continued) Author Description Satterlund & Haupt, 1970 Idaho Watanabe & Ozeki, 1964 Japan over 80% of i n i t i a l l y captured snow ultimately reaches the forest f l o o r age, pruning and thinning a f f e c t s the quantity of snow interception projected area of cedar changes during storm as intercepted snow increases Net accumulation (under canopy) approach Gary & Coltharp, 1967 Jeffrey, 1968 Kuz'min, 1963 L u l l & Rushmore, 1961b Meiman, 1968 M i l l e r , 1966 Seppanen, i960 New Mexico review steppes, U.S.S.R, models review review Finland - aspect and vegetation type influence the amount of snow accumulation - reviews the p r i n c i p a l processes which act upon the snowpack - snow d r i f t may penetrate into forest for a d i s -tance of 5 tree-heights from the margin - snow accumulation varies with species so that more i s deposited under hemlock than balsam, spruce or white pine - supplies a l i s t of annotated references on the influence of elevation, aspect and forest canopy on snow accumulation - differences i n snow accumulation between forest and opening vary with storms, forest types, density and age - uneven amount of snow accumulated around trees 79 gradually before crown interception has reached a steady state, a f t e r which throughfall should be d i r e c t l y proportional to snowfall: and dS. dS" dS, dT k l S o when d l dS. * 0 = k. when where d l d S o and S = 0 (IV.14) are snowfall inside f o and outside the forest i s snow interception. Integrating, one obtains S n = k-. S t a-i f 1 o 1 when d l ^ 0 and S f ~ k 2 S o + a2 when dl dS = 0 (IV.14.a) o The constants a^ and a 2 are dropped because throughfall cannot occur when there i s no snowfall recorded i n the open. From f i e l d measurements, c o e f f i c i e n t k^ was This resulted i n the following equation ( f i g . IV,7) 2 obtained , * for S f S Q (IV. 15) It was not possible to determine the moment when interception attains a steady state. For p r a c t i c a l purposes, when the computed value of exceeds the value of S Q , i t i s 1, The value of a 1 from least-square f i t t i n g was found to be 2.4 mm. which i s approximately zero, thus j u s t i f y i n g the exclusion of t h i s term i n equation IV.15. The standard error of estimate i s 5»81, while the co r r e l a t i o n c o e f f i c i e n t i s 0.95 at 19 degrees of freedom. 80 Snow in open (in mm) IV. 7 T h r o u g h f a l l - s n o w f a l l r e l a t i o n s h i p . U n i v e r s i t y o f R r i t i s h C o l u m b i a R e s e a r c h F o r e s t 81 assumed that steady state condition had been attained by the interception process, and S f = S Q (IV.15.a) It should be re-emphasized that these relationships hold for a f a i r l y closed canopy with l i t t l e d r i f t i n g , evaporation or melting of the fresh snow. Equations IV,15 and IV.15.a are therefore assumed applicable only to forests below the alpine zone. IV-4 Discussion A temporal and several s p a t i a l models have been proposed to extend p r e c i p i t a t i o n records from some v a l l e y s t a t i o n . Using the temporal model, short records from a base-station can be extended. With the s p a t i a l models, i t i s possible to budget t o t a l basin p r e c i p i t a t i o n input by routing p r e c i p i t a t i o n through various a l t i t u d i n a l and vegetation zones. On the other hand, the problem of snow storage cannot be tackled unless (1) we can determine whether p r e c i p i t a t i o n comes as r a i n or snow, and (2) the quantity of snowmelt i s known. Neither snowfall nor snowmelt processes can be inferred d i r e c t l y from p r e c i p i t a t i o n records, but are determined by available energy, an indicator of which i s a i r temperature. Although temperature i s an i n f e r i o r indicator compared v/ith r a d i a t i o n data, the ease of measurement has resulted i n far more temperature recording stations i n the coastal region. In 82 mountainous areas where clim a t i c elements change over short distances, i t i s more appropriate to employ an index with broader s p a t i a l coverage„ Adopting t h i s p r i n c i p l e , temperature w i l l be an important exogenous variable i n the present model. 83 V TEMPERATURE In the absence of data on energy fluxes, a i r temp-erature i s the only available variable to distinguish between r a i n and snow events, and to estimate the amount of snow melt. It i s therefore e s s e n t i a l to model d a i l y temperature i n d e t a i l to permit short-term budgeting of snow storage. Prior to the development of t h e o r e t i c a l models, there i s a b r i e f account of the climatic system which has a bearing upon the d i s t r i b u t i o n of temperature patterns. V-1.0 Physical model Being sheltered from the cold A r c t i c a i r masses by the Coast Mountains, the l i t t o r a l zone of B r i t i s h Columbia i s mild i n winter except at high a l t i t u d e s . In summer, oceanic influences keep the temperature lower than continantal regions of comparable l a t i t u d e . The re s u l t i s a rather small range i n mean monthly temperatures. Over the entire region, annual temperature range (amplitude) tends to increase with a l t i t u d e and distance from the sea, with the Coast Mountains demarcating the coastal b e l t of low temperature ranges from the i n t e r i o r sector which has higher annual ranges but lower annual means. Superimposed on t h i s annual cycle of temperature changes are the heat and cold waves which may l a s t from several hours to several days. In winter, temperature drops when the continental A r c t i c a i r mass breaks occasionally through the mountain barriers ( f i g . IV.1). Within a prolonged warm or 84 cold s p e l l , however, the day to day changes i n 'mean* temp-erature are rather gradual. It i s not uncommon to observe a steady increase i n d a i l y maxima accompanied by an increase i n d a i l y minima as a hot s p e l l continues. Such persistence e f f e c t i s p a r t i c u l a r l y marked i n summer when one type of a i r mass may dominate the weather for a long period. The t r a n s i t i o n from one a i r mass type to another often brings about large varia t i o n s i n temperature from one day to the next. However, th i s may be masked d i u r n a l l y by temperature range of consider-able magnitude. In the case of some v a l l e y stations, the diurnal temperature range may occasionally exceed seasonal temperature differences. On the other hand, the passage of low pressure centres i s frequently associated with l i t t l e diurnal temperature changes. A l l these modes of temperature changes suggest that several sources of climatic controls are operative at d i f f e r e n t time scales, variously described as seasonal cycle, heat and cold waves, day to day persistence and diurnal temperature ranges. I f i t i s assumed that t h e i r e f f ects are additive, then d a i l y temperature records w i l l bear resemblance to the algebraic sum of t h e i r e f f e c t s . V-2.0 Temporal variations of temperature Assuming that d a i l y temperature records manifest the combined e f f e c t s of several climatic processes, we attempt to decompose the records into appropriate components which corres-pond to the generating processes. This introduces the analysis of temperature time s e r i e s . 85 V - 2 , 1 Theoretical model A sequence of events 0 2 ' ' 6 n w h ^ - c h a r e sequentially arranged i n time w i l l constitute a time s e r i e s . These events may be discrete or continuous and are regarded as functions of time. In most cases, the events 6(t) are points recorded at c e r t a i n evenly spaced time in t e r v a l s or they are averages for the records within the i n t e r v a l At . For weakly stationary discrete time se r i e s , K n i s e l and Yevjevich ( I 9 6 7 ) found that information regarding the time series (measured by the variance unexplained by the deterministic element) decreases with an increase i n At , and with a decrease i n n where n i s length of the series. The time i n t e r v a l employed i n most hydrologic studies may be monthly (e.g. Roesner, I 965) or d a i l y (e.g. Quimpo and Yevjevich, 1 9 6 7 ) . The present study uses d a i l y i n t e r v a l s because published temperature records provide d a i l y maxima and minima. Daily temperature 'means' (computed as averages of da i l y maxima and minima) which constitute the time series are considered to comprise a deterministic Z(t) and a random component $(t) : 6(t) = Z(t) + \(±) (V.l) The deterministic component can be subdivided into several generating processes such as trends, o s c i l l a t i o n s or persistence. It should be emphasized that the decomposition of a time series does not imply that the series i s generated by the components operating as independent systems (Bhuiya, 1968). Rather, i t i s a device to characterize the r e a l system by assuming a d d i t i v i t y 8 6 of the components. To decompose the temperature time series, i t i s necessary to remove various deterministic elements successively from the o r i g i n a l data u n t i l the random component remains. A commonly used technique i s to test for s t a t i s t i c a l significance of the autocorrelation c o e f f i c i e n t s r k for sequences 6-^  and 6 2 separated by k time i n t e r v a l s , where k=0, 1, n (Matalas, 1963 J Matalas, 1 9 6 7 ) . 8 1 = 6 1 » V •••» 6 n 62 = 0k+l» ek+2» 6k+n k V[Var e + Var e t + k ] n S 9, Z 0 -k n-k n-k , -» t . , t+k 1 r 1 v fi o t=l t=l -1 " a(e t) t f(6 t + k) t n ^ / t V k — J (V.2) For c i r c u l a r time series, which i s one that has i t s l a s t value followed by the f i r s t value, autocorrelation c o e f f i c i e n t s are determined by r, = k n 9 n 2 £ 9 t 2 - " [4 e t ] ( Y . 3 ) which i s not greatly d i f f e r e n t from equation V.2 i f n i s large. I f , f o r a given significance l e v e l E [ r k ] = P k = 0 (v .4 ) where E [ J i s the expected value P i s autocorrelation c o e f f i c i e n t for the population :8? the series w i l l he considered as randomly generated. To decide whether a c e r t a i n c o r r e l a t i o n c o e f f i c i e n t shows s i g n i f i c a n t association between the sequences correlated, some confidence l i m i t s have to be set up. The available tech-niques are s t i l l not perfect, though the most commonly used one i s that of Anderson ( 1 9 4 2 ) . He showed that for a normal, random c i r c u l a r time series, the confidence l i m i t s for r, can be written as - 1 + t V(n-k-2) C.L.„ = — a n-k-1 (V .5) where t a i s standardized normal deviate at p r o b a b i l i t y l e v e l a for a two-tailed test The assumption of c i r c u l a r i t y becomes less r e s t r i c t i v e as n increases. When n tends to i n f i n i t y , the c i r c u l a r time series converges to an open time series. For graphical presentation, various values of r ^ are plotted against k and the r e s u l t i n g chart i s a correlogram. The maximum number of lags to be used i s usually taken as less than 10% of n , Lags exceeding t h i s number are not recommended because some o s c i l l a t i o n may be induced through spurious auto-c o r r e l a t i o n . It i s also obvious that for every increase of a l a g period, there i s an accompanying loss of one data point which leads to a reduction i n the degrees of freedom. With the aid of correlograms, generating processes are r e a d i l y detected. These processes have been broadly c l a s s -i f i e d as autoregressive, o s c i l l a t o r y or purely random (Chow, 88 1969» p.8). An autoregressive generating process exhibits decreasing values of r k as k increases, thus s i g n i f y i n g persistence e f f e c t which declines as lag periods increase. An o s c i l l a t o r y process i s characterized by rhythmic f l u c t u a t -ions of r^. which shows the existence of some wave components i n the record. Random processes, however, shows the absence of any deterministic element, so that r ^ fluctuates within the belt bounded by the confidence l i m i t s . In t h i s study and during each stage of time series decomposition, correlograms from residuals w i l l be examined to detect the types of processes remaining with the r e s i d u a l s . Based on our physical model, several possible generating pro-cesses have produced the following elements whose ef f e c t s are assumed add i t i v e : e(t) = ?(t) + #(t) + »(t) + P(t) + %(t) (V.6) where t i s a trend component ft a seasonal component (j an o s c i l l a t o r y component P an autoregressive com-ponent J* a random component. V-2.1.1 Trend component To demonstrate the existence of a trend requires a knowledge of the long-term (population) mean (see section IV-2.2.1). Markovic (1967) has described techniques with which differences i n mean values over two periods of time can be tested, given the values for a target period or values from a c l i m a t o l o g i c a l l y homogeneous sta t i o n . This target-control 89 approach requires a knowledge of the target mean or the mean value of a control station. It i s d i f f i c u l t to define a period which i s s u f f i c i e n t l y long to be representative of the population, and i t i s not easy to decide on a sta t i o n that i s cl i m a t o l o g i c a l l y homogeneous. For these reasons, annual means are assumed to be stable and th e i r second moments assumed f i n i t e (equation IV.2). V-2,1.2 O s c i l l a t o r y component Seasonality gives r i s e to very pronounced annual o s c i l l a t i o n s i n the temperature data. This component w i l l be represented by a Fourier series (Carson, 1963; Conrad and Pollak, 1 9 5 0 ) . e ( t ) = ^ + J ^ c o s i p + b . s i n i p ) ( v > 7 ) where L i s h a l f period. a, 4 f 6 ( t ) c o s ^ dt -L L 2L « Y 2 e(t) c o s ^ (V.?.a) L t=l L b. = i S e ( t ) s i n i ? i dt -L L Both a^ and b^ are more conveniently represented by the amplitude o f the harmonic o u and i t s phase angle 0^ , To obtain a and 0 from equation V.7, we make use of the trigonometrical i d e n t i t y sin(B + 0) = s i n 0 cos £ + s i n B cos 0 (V.8) 9 0 iiTt Putting p.=-T— , we obtain the identit y 1 ii ct.sin(p. + 0.) = a.sin 0. cos 0. + a.sin 0. cos 0. JL <1 = a icos £L + b^sin £^ (form v.?) Solving for and 0. a. = 4(al + b 2) ( V . 9 ) 0. = t a n " 1 ( a i / b i ) ( V . 1 0 ) Substituting back to equation V.7, 6(t) = B + E a.sin + 0. ) ( V . l l ) i= l 1 L 1 which i s more amenable to physical interpretation, with denoting the amplitude (or temperature range), and 0^ representing the time of a r r i v a l of the minimum temperature. To compute the amplitudes and the phase angles, mean monthly rather than d a i l y temperatures were used because of computational e f f i c i e n c y . For the following reasons, only the f i r s t harmonic (i=l) was calculated. ( 1 ) the variance accounted for by the f i r s t harmonic usually far exceeds that of other harmonics combined, (2) computation i s s i m p l i f i e d , (3) while the f i r s t harmonic conveys the interpretation of sea-sonality, physical meaning of higher harmonics are less apparent. Thus, to remove annual o s c i l l a t i o n from temperature data, the Fourier series w i l l be expanded far every day of the year, and such value w i l l be subtracted from the o r i g i n a l temperature record. A h e u r i s t i c approach was adopted to describe mathe-matically the heat and cold waves by a set of cosine waves. These cosine waves are represented by t h e i r wave-lengths X , amplitudes a , and the mean values of the waves 0 , 0*(t) = S ! + <x!!cos(t / x i ) (V.12) where the subscripts denote the i - t h wave occuring i n the year. F i t t e d cosine waves are smooth and symmetrical. By subtracting such smooth waves from the less regular o s c i l l a t i o n s from the o r i g i n a l data, i t i s l i k e l y to induce some minor fluctuations which need not be s e r i a l l y independent. The magnitude of such minor fluctuations, however, should be small so that f o r p r a c t i c a l purposes, they may be ignored. V-2.1.3 Markov component To model persistence e f f e c t , temperature r e s i d u a l of one day i s considered to be l i n e a r l y related to that of the previous day. This assumption enables the adoption of a f i r s t order Markov models e"(t) = P 10"(t-1) + 5(t) * r^ e - f t - i ) + *(t) (V.13) where P-^  i s the f i r s t order auto-c o r r e l a t i o n c o e f f i c i e n t , approximated by the sample f i r s t order autocorrelation c o e f f i c -ient r ^ $ i s an independent random var i a b l e . It i s not c e r t a i n whether the atmosphere retains a 'memory' of temperature for more than one day, so that higher order Markov 92 processes are discarded. S t a t i s t i c a l l y , there i s no s i g n i f i c a n t increase i n the c o e f f i c i e n t s of determination when higher order Markov processes are assumed. V-2.1,4 Random component Random component for the temperature time s e r i e s , ^ ( t ) , may be estimated from equation V.13 when the value of r ^ i s knowns * ( t ) = e"(t) - e'^t-Dr-L (v.13.a) Theoretically, i f a l l the deterministic components are removed from the time series, the random component should be an independent random variable which bears the properties of a Wiener-Levy process (Papoulis, 1965t P. 5 0 3 ) 1 . V-2.1.5 Diurnal temperature ranges Diurnal temperature ranges defined as d a i l y maxima minus minima, are related to 'mean' d a i l y temperatures through some conditional p r o b a b i l i t y structure. From h i s t o r i c a l records, i t i s possible to obtain frequencies at which diurnal ranges 6 are conditioned by d a i l y means e , By the law of large numbers (Lamperti, 1966) 1. Such process has the p r o p e r t i e s of (1) U O ) = 0 , ? ( t ) i s r e a l (2) E U ( t ) ] = 0 (3) i t has independent increments (4) t h e d e n s i t y o f [?(t+s) - | ( t ) ] depends on s and not on t , or, the process i s s t a t i o n a r y (5) VAR[5(t+s) - * ( t ) ] = f ( s ) (6) tat a l l s , t ^  0 , L5(t+s) - § ( t ) J i s normally d i s t r i b u t e d . 93 p - r f i * i s " ] - l i m n ( e " " l 3 ) where n(6*|5) i s the number of days with diurnal range 6 and dai l y mean 6" N i s the number of days with a mean d a i l y temperature § . A p r o b a b i l i t y matrix may then be constructed from d a i l y observations. This matrix enables the re-generation of d a i l y temperature maxima and minima given the d a i l y means (ef. the Monte Carlo technique used by Krumbein to generate stratigraphic sequences, 1 9 6 ? ) . V - 2 . 2 Results of analysis A computer programme was written i n FORTRAN IV language to decompose temperature time series ( f i g . V . l ) . For the present study, thirteen years of record for the st a t i o n at the Administration Building of the University of B r i t i s h Columbia Research Forest were analysed. Successive decomposition of the temperature data y i e l d s the following r e s u l t s . V - 2 . 2 . 1 Fourier component Table V . l l i s t s the annual Fourier parameters calculated from the records. Such parameters as annual means, amplitudes and phase angles were tested for autocorrelation which was found to be i n s i g n i f i c a n t at 9 5 ^ confidence l e v e l . Although h i s t o r i c a l data were not long enough to establish s t a t i s t i c a l d i s t r i b u t i o n s of these parameters, the fact that they are independent random variables enable the usage of the central l i m i t theorem (Lamperti, 1 9 6 6 ) . Thus, i f the period Initialize year counter Input tn«. <wd min. temperature 'main' temperature s £ (mM. + V\in) I find mentfily means and Fourier parameters (generate barter Component I Residual C I ) - mean temp. ~ fourier component Find cosine-vjave parameters < Gfencrate cosine-WAW component 0 94 Generate autoregressive component Final residual = residual U D -AUfci regressive component obtain probdbiUty-matrix for daily mean temp. vs find means, variances, frequency- distributions of parameters Increment year coittitec Residual CE) : residual CD -cosine-wave Component Test for trend Find BUcfotf parameter t 6 Output tempecature s statistics F l o w - c h a r t f o r t i m e - s e r i e s t h e d e c o m p o s 1 t i o n o f t e m p e r a t u r e 9 5 Table V . l Parameters from various components of temperature time series (temperature records of U.B.C. Admin) Fourier component Cosine wave Markov Water year 5 ttl r x for 5' r l 1 9 5 7 - 5 8 10 .67 7.44 2 . 6 9 0 . 6 5 2 . 5 0 7 1 9 5 8 - 5 9 9.04 7 . 3 5 2 . 6 0 6 . 5 1 2 . 4 9 4 1 9 5 9 - 6 0 9 . 3 9 7 . 7 8 2 . 5 6 1 . 7 2 9 .446 1 9 6 0 - 6 1 10.08 7 . 0 9 2 . 5 7 3 . 5 8 9 . 5 4 9 1 9 6 1 - 6 2 8.81 7.18 2 . 5 6 3 .680 . 6 0 2 1 9 6 2 - 6 3 9 . 8 8 7 . 2 8 2 . 4 7 0 . 6 1 1 . 5 2 7 1 9 6 3 - 6 4 9 . 1 2 6 . 5 3 2 . 4 1 6 . 7 0 1 . 4 3 5 1 9 6 4 - 6 5 9.04 8 . 0 3 2 . 5 7 ? . 6 0 8 . 4 9 4 1 9 6 5 - 6 6 9 . 5 6 7.08 2 . 4 3 8 . 7 0 0 .448 1 9 6 6 - 6 7 1 0 . 3 9 7 . 7 8 2.446 . 7 5 9 . 4 8 5 1 9 6 7 - 6 8 9 . 8 3 6 . 7 8 2 . 5 6 1 . 6 3 3 .481 1 9 6 8 - 6 9 8 . 7 4 8 . 9 4 2.641 . 7 3 2 . 4 9 2 1 9 6 9 - 7 0 9 . 7 6 6 . 6 9 2 . 5 2 4 . 6 7 7 . 4 9 8 Mean 9 . 5 6 7 . 3 8 2 . 5 4 4 . 6 6 0 . 4 9 7 Variance 0 . 3 7 0.41 0 , 0 0 7 . 0 0 5 . 0 0 2 i s mean annual temperature i n C i s mean temperature of cosine waves, i n °C i s annual temperature range i n °C i s phase angle i n radians i s f i r s t order autocorrelation c o e f f i c i e n t 96 of simulation, n , i s of s u f f i c i e n t length, these random variables w i l l assume normal d i s t r i b u t i o n s n ( 2 8 . - n B ) l i * F r [ _ l f l _ 6 - j i 6 7 = — - — £ e x p ( - t V 2 ) dt (V . 1 5 ) V(2TT)-» where the mean and standard deviation of the independent, i d e n t i c a l l y d i s t r i b u t e d , random variables 6-^9 © 2 , ,,,, 6^ are respectively 0" and o* . V - 2 . 2 . 2 Cosine wave component As for the cosine waves, i t was found that wave-lengths X are correlated with wave amplitudes a' t a' = a + bX (V . 1 6 ) For the Administration station, a = . 0 5 5 b = „383 c o r r e l a t i o n c o e f f i c i e n t = . 7 0 5 at 1 , 1 6 0 degrees of freedom standard error = . 7 3 Wave-lengths are undefined for periods of less than two days (i n view of the time base used i n the present study), so that a truncated s t a t i s t i c a l d i s t r i b u t i o n i s appropriate. Kolmogorov-Smirnov test shows that a truncated Poisson d i s -t r i b u t i o n gives very good f i t to the data ( f i g . V . 2 ) . The fact that wave-lengths are s e r i a l l y independent and that 98 t h e i r mean ( 4 . 0 1 ) i s close to the variance ( 3 . 5 9 ) are i n agreement with two basic c h a r a c t e r i s t i c s of Poisson d i s -t r i b u t i o n s . Unlike wave-lengths or wave amplitudes, mean wave temperatures are autocorrelated. In keeping with adopted procedure for time series decomposition, t h i s autoregressive element was removed. The r e s u l t i n g residuals constitute an independent, normally d i s t r i b u t e d random variable ( f i g . V . 2 ) . With s i m i l a r arguments as applied to the Fourier parameters, the annual autoregressive c o e f f i c i e n t s for mean wave temperat-ures would be normally d i s t r i b u t e d through the central l i m i t theorem. V - 2 . 2 . 3 Markov component Markov c o e f f i c i e n t s calculated for equation V . 13 are tabulated i n table V . l , These c o e f f i c i e n t s are s e r i a l l y independent, hence constituting some random variable whose d i s t r i b u t i o n i s normal when the simulation period i s s u f f i c -i e n t l y long (through the central l i m i t theorem). V - 2 . 2 . 4 Random component F i n a l temperature residuals from time series decomp-o s i t i o n are found to be normally d i s t r i b u t e d with a mean very close to zero and a variance of 2°C. These c o e f f i c i e n t s are not s e r i a l l y independent, but exhibit some form of o s c i l l a t i o n . The magnitude for this type of o s c i l l a t i o n i s small, so that f o r p r a c t i c a l purposes, t h i s element w i l l be ignored. 99 V - 2 . 2 . 5 Diurnal temperature ranges Figure V„3 i s a plot of the conditional p r o b a b i l i t y (z-axis) of diurnal temperature range (y-axis) as conditioned by d a i l y means. It was found that for days with r e l a t i v e l y high (over 15°C) mean values, diurnal ranges tend to increase. Small range or/curs when d a i l y means are moderate, while a tendency of increasing temperature range i s noticeable as dai l y mean drops below -5°C. Days with large diurnal ranges probably correspond to periods of clear skies and days with smaller ranges are related to the passage of low pressure systems. The various stages during which deterministic elements are removed are graphically presented as figures VA, V . 5 and V.6. V - 3 . 0 Spatial v a r i a t i o n of temperature The rationale of p r e c i p i t a t i o n d i s t r i b u t i o n over rugged t e r r a i n may be followed to conclude that d a i l y temp-erature of a- s t a t i o n should be a manifestation of interactions amongst climate, physiography and vegetation cover. It i s further assumed that d a i l y weather patterns are similar over the entire surface of a small drainage basin. It w i l l then be possible to express some re l a t i o n s h i p between temperatures of a base-station and those of the res t of the basin. V-3.1 Physiographic considerations - physical model and re s u l t of analysis Daily temperature records from mountain stations o o Fig. V.3 P r o b a b i l i t i e s of temperature ranges as conditioned by d a i l y means 30 'fflean' cfoify temperature Ifitted vdua (by Fourier secies) •2o 1— 0 T 60 - i 1 T 1 r 180 240 Vays since Ist0ct.i968 T F i p . V.A D e c o m p o s i t i o n o f t e m p e r a t u r e t i m e - s e r i e s : f i t t i n g F o u r i e r s e r i e s t o mean d a i l y t e m p e r a t u r e s 30 — l Xemperaturz residuals I Iritted i>Alu& (ty cosinz W&ves) 1 1 1 300 360 f i t t i n p c o s i n e waves 30 20 — c 10-V<xys since istOct.MS F i g . V.6 D e c o m p o s i t i o n o f t e m p e r a t u r e t i m e - s e r i e s : f i t t i n g f i r s t order M a r k o v component t o r e s i d u a l t e m p e r a t u r e s 10 if r e f l e c t the ambient temperature f i e l d at several le v e l s of generality, A conspicuous feature i n mountain t e r r a i n i s a l t i t u d i n a l temperature changes (Tanner, 1 9 6 3 ) induced by temperature lapse rates developed i n response to hydro-statics of the atmosphere (Hess, 1 9 5 9 ) . From data over B r i t i s h Columbia, Walker (I96I) has shown that lapse rate varies over the year, dropping to a low mean value i n winter. Superimposed on t h i s source of temperature differences i s a i r mass dynamics operating at a synoptic scale. Through turbulent heat and vapour transfer and with variat i o n s i n cloud cover, heat balance of a region i s altered by the passage of a i r masses. At a more l o c a l i z e d l e v e l , temperatures are modified by the d i s p o s i t i o n of ridges and va l l e y s (Lowry, 1 9 6 3 ) , aspects (Frank and Lee, 1 9 6 6 ; Rouse and Wilson, 1 9 6 9 ) and the l o c a t i o n of the station with respect to the path of cold a i r drainage (Bergen, 1 9 6 9 ? West, 1 9 6 1 ) . Radiation differences between ridges and valley bottoms or amongst slopes of d i f f e r e n t aspects often induce considerable temperature differences (MacHattie, 1 9 7 0 ) which give r i s e to convective c i r c u l a t i o n s i n the v a l l e y (MacHattie, 1 9 6 8 } Urfer-Henneberger, 1 9 6 4 ) . Thus, phenomena such as nocturnal temperature inversions are common. From an analysis of variance, Keen ( 1 9 5 6 ) demonstrated that interactions between elevation and aspect, elevation and location, and aspect and lo c a t i o n are highly s i g n i f i c a n t for the d i s t r i b u t i o n of low temperatures at an Oregon ponderosa pine fore s t . His study suggests that locations have i n d i v i d u a l c h a r a c t e r i s t i c s i n respect to temperature (p. 3 8 3 ) . Consideration of the entire range of physiographic io<5 factors i s too great a task for t h i s study; only the factor of a l t i t u d i n a l difference between stations w i l l be used. To investigate a l t i t u d i n a l e f f e c t on d a i l y temperature maxima and minima, records from four stations i n the Research Forest were used (see section I I I - 3 ) . On many occasions, noticeable temperature inversions occurred. The general tendency, however, i s a decrease of d a i l y maxima and minima with height. The lapse rate concept leads to a postulation of l i n e a r i t y so that equations of the following form were proposed and t h e i r c o e f f i c i e n t s tabulated i n tables V . 2 and V . 3 . \ = ak + V i J ( V ' 1 7 > where 6 . i s temperature of station • Jk when station i records temper-ature of magnitude k h.. i s a l t i t u d i n a l difference between stations i and j a and b are c o e f f i c i e n t s determined by least-square f i t . In p a r t i c u l a r , b i s the temperature lapse rate. V - 3 » 2 Forest cover and temperature - physical model and  re s u l t of analysis The presence of vegetation adds to the complexity of the temperature f i e l d . Forest cover performs the function of reducing turbulent heat transfer between land surface and the free atmosphere. Its presence also modifies the amount of incoming and outgoing r a d i a t i o n and increases heat loss through evapotranspiration (Rauner, 1 9 6 2 ) „ Maximum summer d a i l y temp-Table V .2 Least-square f i t i maximum temperature vs a l t i t u d e , U.B.C. Forest lemperature T ^ + _ _ + c 1 _ a Correlation Degrees of Standard at base-station Intercept Slope c o e f f i c i e n t freedom error ( bC) 0 - 3 3 . 0 - 0 . 0 0 6 9 - 0 . 5 8 8 19 1 .5 3 - 6 5 . 9 - 0 . 0 1 0 7 - 0 . 8 0 7 99 1 .2 6 - 9 8 . 7 - 0 . 0 0 8 8 - 0 . 5 9 6 203 1 .8 9-12 11 .6 - 0 . 0 0 7 8 -O.506 175 2 . 0 12-15 14 .6 - 0 . 0 0 9 5 - 0 . 6 5 1 195 1 .7 15-18 18 .2 - 0 . 0 1 0 8 - 0 . 7 1 3 107 1 .6 18-21 2 0 . 3 - 0 . 0 0 9 4 - 0 . 4 9 2 123 2 . 5 21-24 2 3 . 8 - 0 . 0 0 8 4 - 0 . 4 6 2 171 2 . 4 24-27 2 6 . 2 - 0 . 0 0 6 6 -0.440 87 2 . 0 2 7 - 3 0 3 0 . 1 - 0 . 0 0 6 1 - 0 . 3 2 8 51 2 . 7 Data fromj U.B.C. Admin, U.B.C. Forest (Marc) Loon Lake Placid Lake 1969-70 water year Table V .3 Least-square f i t : minimum temperature vs a l t i t u d e , U.B.C. Forest •Temperature at base-station ( i n °C) Intercept Slope Correlation c o e f f i c i e n t Degrees of freedom Standard error (°c) Below 3 - 3 A - 0 . 0 0 2 8 - 0 . 2 2 1 31 1 .9 3-6 - 0 . 3 -0.0047 - 0 . 3 9 1 103 ! • ? 6-9 2 . 3 - 0 . 0 0 5 8 - 0 0 4 5 8 271 1 .7 9-12 5 A - 0 . 0 0 7 2 - 0 . 5 5 1 243 1 .7 12-15 7 . 9 - 0 . 0 0 5 5 - 0 . 3 6 9 219 2 . 1 15-18 1 0 . 9 - 0 . 0 0 5 0 -0.406 263 1 .7 18-21 1 3 . 7 - 0 . 0 0 5 2 - 0 . 4 2 6 91 1 .7 21-24 1 5 . 8 -0.0077 -0.484 11 2 . 3 Data from: U.B.C. Admin. U.B.C. Forest (Marc) Loon Lake Placid Lake 1969-70 water year 108 erature i s not uniform (Geiger, I 9 6 5 ) 1 . We s h a l l not consider thi s source of v a r i a t i o n because snow accumulation and melt processes are more clos e l y related to the s t r a t a of a i r not far above the ground or snow surface. Therefore, only paired temperature observations were made, one i n the forest and the other i n the open nearby (at Placid Lake). Both d a i l y maxima and minima i n the forest were assumed to be l i n e a r l y r e l a t e d to those of the openi 6 f = a + b6 Q (V.18) Where the subscripts f and o re f e r to forest and open respectively a and b are c o e f f i c i e n t s determined by least-square f i t . The values of a and b as obtained by least-square f i t t i n g are as follows: Maximum temperature Minimum temperature Intercept (a) .40 1.61 Slope (b) .91 .96 Corr. c o e f f i c i e n t at 357 deg. of (r) freedom .99 .98 Standard error 1.08 .96 Relationships between temperatures i n the forest and temperat-ures i n the open are i l l u s t r a t e d by figures V .7 and V.8. V-4o0 Summary Daily temperature records have been considered as 1. Such temperature p r o f i l e s have been successfully simulated by Waggoner and Reifsneider ( I 9 6 8 ) . 109 F i g . V.7 Minimum d a i l y temperatures the forest vs the open, Placid Lake, 1969-70 water year 110 32 24-Q. E '4 E £ <6 8 -A o H o AUTUMN A WINTER + SPRING x SUMNER - 0 T" 16 T 0 8 i*> - 24 Maximum daily temp, in open ( cn °C ) 52 F i g . V.8 Maximum d a i l y temperature: the forest vs the open, Placid Lake, 1 9 6 9 - 7 0 water year I l l discrete-event time series comprising a deterministic and a random component. By decomposing the time s e r i e s , several deterministic elements were i d e n t i f i e d so that d a i l y temper-ature may be simulated for a longer period by recombining such components p r o b a b i l i s t i c a l l y . Temperature data from a base-station are extended to other parts of the same drainage basin, making use of several l i n e a r transfer functions which take into account a l t i t u d i n a l difference and the presence of a forest cover. This enables the simulation of temperature f i e l d s over the entire basin. The next section w i l l consider interactions between p r e c i p i t a t i o n and temperature which r e s u l t i n snowfall events or the melting of snowpacks. 112 VI SNOWFALL AND SNOWMELT EVENTS In the absence of cloud, wind or r a d i a t i o n data, temperature serves as a s i g n i f i c a n t indicator of snowfall and snowmelt events. In the discussions to follow, a i r temperature recorded i n Stevenson's screens w i l l be used because i t i s the most r e a d i l y available type of temperature data, VI-1 Snowfall events - s t a t i s t i c a l model Whether p r e c i p i t a t i o n reaches a basin as r a i n f a l l or snowfall depends on the nature of the storm-bearing clouds and conditions of the layers of a i r through which p r e c i p i t a t i o n penetrates. Thus, snowfall does not necessarily occur when screen temperature re g i s t e r s 0°C, though screen temperature r e f l e c t s c l o s e l y the heat conditions near the ground. To develop a p r o b a b i l i s t i c model for the occurrence of snowfall events given a ce r t a i n screen minimum temperature (U.S. Army, I 9 5 6 ) , we make use of the law of large numbers: Pr[snowlminimum temp.^0°c] = Jjf^ (n / N) (VI.l) where n i s the number of snow events s given a minimum d a i l y temperat-ure equal to or below 0°C N i s t o t a l number of p r e c i p i -t a t i o n events at a minimum d a i l y temperature equal to or below 0°C The University of B r i t i s h Columbia Research Forest has been recording snow and r a i n events since 1961 at UBC Admin, 113 Loon Lake and Spur 17 stations. These data were pooled to produce a snowfall p r o b a b i l i t y matrix as conditioned by minimum d a i l y temperatures ( f i g . V I . l ) . A category l a b e l l e d 'mixed events' shows that both r a i n and snow have been recorded on the same day. In subsequent simulations, the p r o b a b i l i t y of occurrence of t h i s class w i l l be s p l i t equally between r a i n and snow events to simplify computations. VI-2.0 Snowmelt events E x i s t i n g approaches to compute snowmelt f a l l under several broad categoriess VI-2.1.1 Multiple regression approach A highly empirical approach to the snowmelt problem i s to correlate melt to a l l types of available meteorological data (such as radiation, temperature, wind and humidity) to arrive at the 'best' predictive equation (Gartska et a l , , 1958; Pysklywec et a l . , 1968), Although for any l o c a l i t y at a given period of time, t h i s approach provides the best estimate, i t i s not certain i f the equations thus developed can be applied elsewhere (O'Neill and Gray, 1971). Also, the independent variables may be i n t e r r e l a t e d which gives r i s e to spurious correlations (Benson, 1965) that lead to physically meaningless r e l a t i o n s h i p s 1 . 1. e.g. Gartska et a l . (1958) found that the contribution of r a d i a t i o n to snowmelt equation was negative because temperature (another independent variable) i s in fact a function of r a d i a t i o n . 114 Probability of rft inUU 0 -25 -50 rj<j |.0 Probability ot snowfall F i g . VI.1 P r o b a b i l i t y o f s n o w f a l l as c o n d i t i o n e d by n i n i n u m d a i l y t e m p e r a t u r e s 115 VI-2.1.2 Energy balance approach The amount of energy available for transforming snow into water, Q M , i s calculated as the sum of net all-wave energy Q N , sensible heat flux Q H , latent heat flux Q E , advective energy flux Q A , and ground heat flux Q G (O'Neill and Gray, 1971; U.S. Army, 1 9 5 6 ) . Q M = Q N + Q H 1 Q E + Q A ± Q G ( V I < 2 ) with QN = QS* " QS* + Q l i " QL* (VI.2.a) where S and L r e f e r to short and long wave respectively, and upward and downward pointing arrows r e f e r to outgoing and incoming r a d i a t i o n respectively. Excellent estimates of point melt have been obtained with t h i s technique (Anderson, 1 9 6 8 ; Rockwood et a l . , 195*0. It i s b a s i c a l l y a model designed to compute snowmelt at a heavily instrumented s i t e though i t has been employed to simulate snow hydrology at a basin scale (Erickson and McCorquodale, I 9 6 6 ; Hendrick et a l . , 1971 ; Johnson and Boyer, 1 9 5 9 ) . Even though the energy balance approach i s superior to other models, i t requires an impractical amount of supporting meteorological data such as wind, vapour pressure and r a d i a t i o n from a base-station. A modified version makes use of several assumptions to simplify some of the parameters (Riley et a l . , 1 9 6 9 ) . Even the modified approach, however, requires a broader data base than most parts of the Coast Mountains can o f f e r . 116 VI-2.1.3 Temperature index approach A s i m p l i f i e d technique i s to consider snowmelt as a function of temperature. Degree-day''" i s the most commonly used indicator of available heat, though simpler indices l i k e mean or maximum temperatures were found to be equally s a t i s -factory. Most authors use l i n e a r equations to r e l a t e snowmelt to temperature indices. Such re l a t i o n s h i p s , however, i s time-dependent, with a general increase i n melt rate towards late spring (Anderson and Rockwood, 19?0; McKay, 1964). Going alongside t h i s phenomenon i s the r e l a t i o n s h i p between maximum da i l y temperature and degree-days as computed by summing hourly temperatures. Gartska et a l . (1958) presented l i n e a r equations which show that July maxima correspond with higher degree-days than do d a i l y maxima i n May. One reason i s that maximum da i l y temperature i s a one-dimensional indicator of magnitude while degree-day i s two-dimensional, being the product of duration and magnitude. As the number of hours with above-freezing temperatures increases towards late spring, degree-day values increase faster than maximum temperature values. In which case, a time-independent l i n e a r equation w i l l not adequately express such changing r e l a t i o n s h i p . Despite s i m p l i f i c a t i o n s , a temperature factor has been found to be "at least as good as, and i n many cases better than a combination of other factors used i n c o r r e l a t i o n analysis" 1. Gartska (1964) pointed out that the accurate computation of this index i s to sum up hourly exceedances above a reference temperature, say 0 ° C 117 (Gartska et a l . , 1958). This approach w i l l be adopted for the present study. In order to determine the necessary mathematical functions, snowmelt measurements were made at the University of B r i t i s h Columbia Research Forest. VI-2.2 Measurement of snowmelt The selection of measurement techniques hinges on the complexity of the model and the l e v e l of accuracy desired. Available techniques may be summed up as follows. VI-2.2,1 Indirect measurement of snowmelt through streamflow  analysis The quantity of snowmelt i s abstracted from stream-flow measured at the basin outlet (Gartska et a l . , 1958s Leaf, 1969). Thus, th i s technique computes t o t a l melt from a basin, but i s not amenable to the determination of melt from d i f f e r e n t a l t i t u d i n a l bands within the basin. For rain-on-snow areas with a fluc t u a t i n g freezing l e v e l , melting and accretion often take place simultaneously i n a basin. This r e s t r i c t s the predictive applications of such integrated melt measurement technique. VI-2.2.2 Indirect measurement by snow survey By computing the differences i n pack storage from day to day, i t should be possible to obtain a measure of snowmelt and snow accretion at a survey s i t e . Unfortunately, the error involved i n sampling and the choice of s i t e render t h i s technique inaccurate. Recent developments i n snow survey based on gamma ray attenuation through the pack (Smith 118 et a l . , 1967) has been able to improve the evaluation of snow storage, while snow pillows (Beaumont, 19&5J Peterson, I 9 6 8 ; Penton and Robertson, 1967) may provide continuous records of snowpacks. VI-2.2.3 Laboratory measurement of snowmelt One way of deriving mathematical function r e l a t i n g snowmelt to heat parameter i s to study melt under laboratory conditions (Quick, 1971). Using t h i s approach, external parameters such as temperature, humidity or r a d i a t i o n can be controlled, but the snow sample i s disturbed and boundary eff e c t s may pose a problem. It i s not c e r t a i n how d i f f e r e n t l y the sample behaves c©mpared with snow in'the pack which i s an open rather than a closed system. VI-2.2.k Lysimeter measurement of snowmelt If i t i s required to measure the amount of melt-water reaching the s o i l rather than the quantity that i s generated at the surface of the pack, then the snowmelt l y s i -meter provides the best solution. Several considerations are e s s e n t i a l when designing such lysimeters. One of these i s the size and shape of the container that receives melt-water. Generally speaking, large receivers minimize rim e f f e c t s , but are more costly to transport and to i n s t a l l . Several examples include the two triangular lysimeters of U.S. Army (1,300 sq. f t . and 600 sq. f t . surface area)5 the e q u i l a t e r a l triangular lysimeter of the University of New Brunswick (16 f t , sides) and the rectangular lysimeter of Intermountain Forest and Range, 119 Experiment Station ( 2 3 " by 17"). Another consideration i s the s e n s i t i v i t y of the melt-water c o l l e c t o r . Most lysimeters have perforated bottoms to ensure no delay i n transmitting melt-water to the recording device. Haupt ( 1 9 6 9 ) went further by scattering several centimetres of duff and needles over the meshed screen to r e p l i c a t e ground surface conditions. Perhaps the most s i g n i f i c a n t consideration i s the rationale regarding the column of snow whose melt-water one wishes to measure. This column above the lysimeter has been is o l a t e d as a closed system (Haupt, 1 9 6 9 ; Pysklywec et. a l . , 1 9 6 8 ) or has been l e f t undisturbed to interact with the r e s t of the pack as an open system (U.S. Army, 1 9 5 3 ) . As a closed system, undersirable rim e f f e c t cannot be avoided (Haupt, I 9 6 9 . p. 7 1 6 ) . As an open system, i t i s well known that melt-water moves l a t e r a l l y as well as downwards (Wakahama, I 9 6 8 ) . Unless a steady state i s attained, l a t e r a l flows may upset the water balance of the snow column d i r e c t l y above the lysimeter. On the other hand, i f the experiment aims at measuring the quantity of melt-water entering the ground at some point underneath the pack, the emphasis i s on the pack above and around the lysimeter rather than on the column of snow d i r e c t l y on top of i t . Taking these points into consideration, an economical snowmelt lysimeter operating on the open system p r i n c i p l e was designed for t h i s study (see section I I I - 4 ) . Based on point measurements from t h i s lysimeter ( f i g . VT . 2 ) , s t a t i s t i c a l models were developed to predict snowmelt. io H F i g . V I . 2 R a i n f a l l , t e m p e r a t u r e , anc! s n o w m e l t r e c o r d s , P l a c i d L a k e s i t e , F e b r u a r y t o A p r i l . 1971 1 2 1 VI-2>3- S t a t i s t i c a l model f o r snowmelt The following models make use of maximum a i r temperature and r a i n f a l l as indices of snowmelt. Following the U.S. Army ( 1 9 5 6 ) , a d i s t i n c t i o n was made between c l e a r -weather and r a i n melt. V I - 2 . 3 . 1 Clear-weather melt Most studies using degree-day indices express snow-melt as a l i n e a r function of d a i l y temperature. It i s doubtful i f the assumption of l i n e a r i t y i s v a l i d . Bruce ( 1 9 6 2 , p. 9 5 ) found that an equation of the form [melt ( i n inches) = O.O65 X temperature above 32°F] over-estimated snowmelt during days with low melt but under-estimated snowmelt during days when melt rate was high. From the Placid Lake lysimeter data, s i m i l a r undesirable r e s u l t was also produced when a l i n e a r equation was used. Instead of l i n e a r i t y , however, Placid Lake data strongly suggest c u r v i l i n e a r r e l a t i o n s h i p between melt and maximum d a i l y temperature. One reason for non-linearity i s that higher maximum temperature often accompanies days with prolonged hours of above-freezing temperatures. Since melt i s a function of the application of heat over an extended time onto the pack, a r i s e i n maximum dail y temperature generally increases towards the end of spring, which i s a period when snow durface albedo i s reduced and when the ripeness of the pack increases. In which case, higher melt per degree-day i s expected. For clear-weather melt, therefore, a power functional r e l a t i o n s h i p i s appropriates 122 = a + b6 2 max e max > 0 (VI.3) where a and b are c o e f f i c i e n t s NL i s clear-weather melt c 8 i s maximum da i l y temp-erature. Least-square f i t t i n g of data from Placid Lake lysimeter resulted i n ( f i g . VI.3) a = -0.73 b = 0.0724 corr. coeff. = 0.92 at 12 degrees of freedom standard error =2.2 Equation VI.3 does not generate posi t i v e melt for 6 max Physically, i t suggests the low l i k e l i h o o d of snowmelt below such a i r temperature, and for the present model, snowmelt w i l l be set to zero. VI-2.3.2 Rain melt same temperature as free a i r , U,S. Army ( I 9 5 6 , p. 180) has developed an expression of r a i n melt based upon heat transfer from r a i n water onto snowpacko Ov/ing to re-freezing of melt-water and rain-water during their passage through the pack, observations from Placid Lake lysimeter show that i t i s not possible to budget water input as the sum of r a i n f a l l and r a i n melt. Instead, the ef f e c t of melting and re-freezing are lumped i n the empirical equation which rel a t e s net outflow from the pack to meteorological factors of r a i n f a l l and maximum temperature„ Based on the assumption that raindrops carry the 123 Maximum daily tempewtuce (.°C) F i g . VI.3 Snowmelt v s P l a c i d L a k e maximum d a i l y t e m p e r a t u r e s i t e ' 12 It «R = a + t > < W > • % a x > 0 < « • * > where MR i s net outflow from the pack R i s d a i l y r a i n f a l l 6max m a x ^ I n u m d a i l y temper-ature . Least-square f i t t i n g gave the following r e s u l t ( f i g . VI.4). a = -1.23 b » 0.278 corr. coeff. = 0.95 at 13 degrees of freedom standard error = k The case when MR < 0 has no physical meaning, so that MR w i l l be considered as aero under such circumstances. When MR < R , the difference represents an increase i n the snowpack. When > R , the difference represents the quan-t i t y of melt-water which f i n a l l y escapes through the pack. Being empirical, the model has necessarily s i m p l i f i e d many of the int e r a c t i n g processes such as water retention and movements i n the pack (Rockwood et a l . , 195*0. In the absence of more elaborate data, however, equation VI.4 i s deemed adequate for the present study. VI-3.0 S p a t i a l extension of snowmelt equations To extend the clear-weather and r a i n melt equations sp a t i a l l j r , we make the following assumptions: (1) Terrain conditions such as aspects are not taken into consideration. I t i s assumed that the temperature and r a i n f a l l data used i n the computation are average values for a certa i n a l t i t u d i n a l zone, with minor r e l i e f elements at d i f f e r e n t parts of the zone compensating one another, (2) The equations are applicable to both forested and open environmentso Although short-wave ra d i a t i o n i n the forest i s 1 2 5 35 -, rUinfatl * Maximum daily temperature F i g . VI.4 Snowpack o u t f l o w v s d a i l y r a i n f a l l and maximum d a i l y t e m p e r a t u r e , P l a c i d L a k e 126 lower than the open, long-wave r a d i a t i o n i s emitted by the f o r e s t canopy to the snowpack, and albedo o f the pack i n the f o r e s t i s reduced by f a l l e n l e aves and twigs. A l l these tend to compensate the r e d u c t i o n of short-wave r a d i a t i o n . Indeed, r e l a t i v e l y t h i n packs underneath the canopy (below 600 m. l e v e l ) u s u a l l y d i s a p p e a r s e a r l i e r than adjacent packs i n the open. Under these assumptions, m e l t i n g o f snowpack at v a r i o u s a l t i t u d i n a l zones w i l l be computed by equations VI.3 and VI.4. VI-4.0 Summary and d i s c u s s i o n A s i m p l i f i e d temperature index approach was presented which enables the budgeting o f snowpack c o n d i t i o n s at v a r i o u s s e c t o r s w i t h i n a s m a l l b a s i n . Although the forms of the mathe-m a t i c a l f u n c t i o n s appear s i m i l a r to those used by other r e s e a r -c h e r s , n u m e r i c a l v a l u e s f o r the c o e f f i c i e n t s are d i f f e r e n t 1 . T h i s suggests a need f o r c o l l e c t i n g data from l o c a l i t i e s where snowpack s i m u l a t i o n s are to be performed. With t h i s l i m i t a t i o n , however, i t i s not d e s i r a b l e to apply the f o l l o w i n g s i m u l a t i o n programme to areas beyond the P a c i f i c l i t t o r a l zone of Canada. 1. Cf. s n o w f a l l event p r o b a b i l i t y matrix used by the U.S. Army (1956), or c o e f f i c i e n t s f o r the clear-weather melt equations as r e p o r t e d by Bruce (I962), Pysklywec e t a l . , (1968) and the U.S. Army (1956). 127 VII SIMULATION MODEL With the p r o p e r t i e s of v a r i a b l e s d e f i n e d and the nature o f t r a n s f e r f u n c t i o n s e s t a b l i s h e d i n the previous c h a p t e r s , a s i m u l a t i o n model was developed to handle two areas o f i n t e r e s t s (1) to extend a v a i l a b l e temperature and p r e c i p i t a t i o n r e c o r d s , (2) to evaluate snowpack storage c o n d i t i o n s . Depending on the purpose of s i m u l a t i o n , i t i s sometimes not necessary to invoke both s i m u l a t i o n procedures. Hence, two sub-models were i n t r o d u c e d i Part I - D a i l y temperature and p r e c i p i t a t i o n s i m u l a t i o n sub-model. Part I I - Snowpack and generated r u n o f f s i m u l a t i o n sub-model. A l l programmes were w r i t t e n i n FORTRAN IV language except when U.B.C. Computing Centre l i b r a r y s u b r o u t i n e s were c a l l e d . Magnetic tapes are used f o r d a t a - h a n d l i n g purposes because the sub-models generate a l o t of data output. T h i s serves the dual f u n c t i o n of keeping more 'durable' r e c o r d s f o r s t a t i s t i c a l a n a l y s i s , and o f p r o v i d i n g an i n t e r m e d i a t e break i n the s i m u l a t i o n sequence i n case only one of the s i m u l a t i o n procedures i s r e q u i r e d . More d e t a i l e d d e s c r i p t i o n s on the a l g o r i t h m s w i l l f o l l o w . A l i s t i n g o f programmes i s presented as the Appendix. VII-1,0 Algorithms F i g u r e V I I . l i s a f l o w - c h a r t i l l u s t r a t i n g the r e l a t i o n -s h i p between sub-models. Although shown i n sequence, Model II 128 Input •• fVriod «f simulation Tc»f«n»tUW l>At»mttcn Prodpit«t«M p*ra«et«r J Daily "Xtmpffatutet 0*»'y PretipitAtien SimuUtion for A t)«U output Stored >'n t»pc Output: Sample of te»perdtu»* Am* precipitation PART 1 Per'tei ai simulation Typt of *i«sal*tt<>* Simulation of SrtotvpMk *torA^ « And <j*oef«te4 rv»eff Runoff O u t f u t ' Sample <UU of InowpAtk condition* SnowpAck. Zonal fcnci-Ate* runoff 4 a t A output $t««d in tope conditions dot* output itored in tape PART H F i g . V I I . 1 F l o w - c h a r t f o r t h e s i m u l a t i o n o f d a i l y t e m p e r a t u r e , p r e c i p i t a t i o n , and snowpack 1 2 9 can be used without going through Model I i f h i s t o r i c a l records of temperature and p r e c i p i t a t i o n are used. V I I - 1 . 1 Sub-model I This sub-model simulates p r e c i p i t a t i o n ( i n mm.), maximum and minimum temperatures ( i n °C) at d a i l y i n t e r v a l s for a given s i t e . The approach i s e s s e n t i a l l y stochastic. For p r e c i p i t a t i o n , a second-order Markov chain simulates the occurrences of wet and dry days, while Monte Carlo tech-nique generates the amount of d a i l y p r e c i p i t a t i o n according to an exponential d i s t r i b u t i o n (see Chapter IV). For temp-erature, several components (seasonal temperature o s c i l l a t i o n s , heat and cold waves, day to day persistence and random f l u c -tuations of temperature) are superimposed l i n e a r l y to produce simulated data (see Chapter V). Details regarding the logic of i t s algorithm are represented by a flow-chart ( f i g . VII.2), VII-1.2 Sub-model II This sub-model simulates snow storage and snowmelt runoff for various a l t i t u d i n a l zones. Data input may be h i s t o r i c a l records from a base-station, or i t can take the form of outputs from sub-model I . The p o s s i b i l i t y of using either type of data increases i t s f l e x i b i l i t y . Other than the simulation of snowfall p r o b a b i l i t i e s , this sub-model i s deterministic. Temperature and p r e c i p i t a t i o n data from base-station are routed through a l l pre-specified zones (sections IV - 3 and V - 3 ) . With the water balance equation, snowpack storage and snowmelt release are generated 130 lnp»*t parameter* Input A&i), i*icntH »n4 yea* counter* Generate toKipsrfttuw Fourier coefficient* •fc' y * * r k-Generate Ami minimum tesnperatucej Generate daily precipitation bAjed m condition*! prob. f«r month "* No Inccement month yes 5'nnuUtuTi period cy»<ie4 Output cUilu mAXi^um 4nd M',m«u<n te<"pffAtti«Z» An* dA'Oy pTWip't*t><>n Plot t t w p c r A t u r e and precipitation d o t * N o F i £ . V I I . 2 F l o w - c h a r t f o r t e m p e r a t u r e and p r e c . i n i t n i o n s i m u l a t i o n s u b - m o d e l 1 3 1 for each zone ( f i g . VII. 3 ) so that a l l zones provide paired indices of these measurements, one for the open and one for the forested sector. This programme provides an option for data output: snowpack storage data or zonal generated runoff data. In "both cases, zonal data are stored onto magnetic tape. I f desired, the programme also y i e l d s graphical output showing da i l y changes i n snowpack storage by a l t i t u d i n a l zones. Details of the programme are shown by a flow-chart ( f i g . V I I . 4 ) . V I I - 1 . 3 Computer time The U.B.C. IBM 3 6 0 / 6 7 computer and the FORTRANH compiler require the following to execute the programmes: Model I Model II Core storage 3 9 , 4 9 4 bytes 4 3 , 0 3 2 bytes Programme compile time 18 sec. 2 0 sec. Execution time for 1 0 0 years of data app. 1 5 0 sec. app. 2 0 0 sec. for 2 5 zones V I I - 2 . 0 Model testing The purpose of the present simulation exercise i s to provide long-term pr o b a b i l i t y statements regarding the behaviour of snowpack storage and snowmelt runoff. It i s therefore necessary to examine the capacity of the model to reproduce data whose properties do not depart s i g n i f i c a n t l y from h i s t o r i c a l records. There i s 5 however, no exhaustive technique for model v e r i f i c a t i o n . One p o s s i b i l i t y i s to ensure that errors do not propagate unnoticed into other components. An important consideration here is that testing P tUuiy precipitat ion f o r «M *<mit RACK snon iterate for all lonm MMRUH tyhvt zof* c0i»trU>ufc»; Surface rt4«#ft RUKPFF surfaet runofl f B r o i l zone* < 0 Ort) <&y *)e» Clear wsafHet melt ( Minpk to N ) fi&xrun r 0 yes Runoff = melt I present zone No Reset t Runoff = Pack Pack = 0 yes .. — Enter: N, M»xr. MinpX, *Wx, p, f Ack 132 N toUl no. of elevation ron* MAXR »''eh«tz<w>« veceiviof ttinbil MfNPK loxMKt zone wift iww »tor*j« TMAX iwwtimiun (Utk) temperature fW All ZOrxfS -o Snow fsr Alt zone* Maxcun Miopk = MAXC+I > 0 B*th rain and jnaw in ba>i» yes R*tn below So«vglift«v Runoff = P •from zone 1 to Mftxr Rain Above Snowline Runoff = P from aone t +0 (Minpk-I) /UU Snow fo pAtk. from zone (MAxr+1) io N > Check for plotting (Miapk to frfcwrt AcW snow zone (MAxr+t) •fo M Runoff = 0 KA*run s present zone exec = 9-SS m«li AtW E io Po wk Runoff = Avail Runoff = malt Pack = 0 A v a i l Subtest (melt - P ) •frcw PACIV No / 4 ) c r t x "* < ? AvAtl Des g. VII.3 Flow-chart for snow storage water-balance subroutine Input • Length ol «i"uUtion period N o . c f M t ' W i n * l xonct Tijpo of ,i«nul*ti»n S n o w p * * plotting option Compute mid-rone A l t i t u d e * Xniti^lize DAU J And V«Af counter* Montfi counter Input : olAily t e m p e r A t u c e AMA p r c c i f ' t a t i o n o f Compute t c w p c c t t a r c o f v a r i o u s x o n M , W i 1 * in o p e n A n d i n f o t e * t CcaguVc p r e c l p i t A t i e n o f v A t i o u J z a n e j , o p e n AtA in f o c e * t . Update p A c k w i t h n e w 4 n o w , 'if A n y > Compute snowmelt, U p d a t e p a c k . ftnoV c o m p u t e , e e n e C A t e d runoff ( s e ? subroutine ZSNOW) Increment VeAC counter — ~ < i n c r e m e n t ^ N o n,ontti lounteC — < < I n c r e m e n t day counter Output Canto tape) snowpack condition) pa generated r u n o f f Check for plotting of snowpack V I I . 4 F l o w - c h a r t f o r snow s t o r a g e s i m u l a t i o n sub-modi 134 techniques should be independent of the way with which the model was derived. Owing to l i m i t a t i o n s on the amount of long-term records, only several key areas of the model can be examined. Despite the selective nature of these tests, they w i l l indicate the r e l a t i v e r e l i a b i l i t y of the simulated data. VII-2.1 P r e c i p i t a t i o n data Records from Stave F a l l s , B.C., were used to test against 100 years of simulated p r e c i p i t a t i o n data. The choice of t h i s s t a t i o n rests p a r t l y with i t s location ( l y i n g i n the Coast Mountain area of southern B.C.) and i t s history (41 years of record). Two-sampled Kolmogorov-Smirnov test shows that the duration of wet and dry s p e l l s from both sets of data did not d i f f e r s i g n i f i c a n t l y (table IV.3)> nor did Mann-Whitney U-tests detect any s i g n i f i c a n t differences between monthly p r e c i p i t a t -ion. The model for p r e c i p i t a t i o n i s therefore considered as s a t i s f a c t o r y . VII-2.2 Temperature data Data from Stave F a l l s , B.C., were again used to test against 100 years of simulated temperature records. In t h i s case, the median, upper and lower q u a r t i l e marks of both maximum and minimum temperatures for each calendar day were extracted ( f i g . VII.5 and VII.6). Mann-Whitney U-tests applied to these s t a t i s t i c s suggest no s i g n i f i c a n t differences and the simulation model i s thus considered adequate. T 1 1 1 r 1 r 50 100 t50 200 250 300 350 DAYS SINCE l>.fOCt F i g . V I I . 5 D i s t r i b u t i o n o f maximum ( u p p e r band) and minimum ( l o w e r band) d a i l y t e m p e r a t u r e s b a s e d on 41 y e a r s of r e c o r d s f r o m S t a v e F a l l s , B.C. ^ SimuUtcA -10 50 100 150 200 250 500 J5B CAYS SINCE Ist OCT. Fig. VII.6 Distribution of maximum (upper band) and minimum (lower band) daily temperatures based on 100 years of simulated record ON 137 V I I - 2 . 3 Snowpack data The only snowcourse with long history (since 1936) for the coastal area i s Grouse Mountain ( 1 2 3 ° 05* W; 4-9° 2 3 ' Nj elevation 1160 m.). Records for the beginning of A p r i l were compared with 100 years of simulated snowpack data for 1 s t A p r i l at the 1 , 2 0 0 m. zone ( f i g . VII. 7 ) . A Kolmogorov-Smirnov test shows;that the snow water-equivalent exceedance p r o b a b i l i t i e s for both sets of data are not s i g n i f i c a n t l y d i f f e r e n t . V I I - 3 . 0 Analysis of simulated snowpack data A l l simulated snowpack data are measured i n water-equivalent units. These data are average conditions for the a l t i t u d i n a l zones they represent. V11 - 3 . 1 Magnitu.de of snowpack Snowpack data for the f i r s t day of each month between November and July were extracted from the simulated records. From these, the p r o b a b i l i t y that snowpack storage exceeds a c e r t a i n magnitude was computed and plotted on a pr o b a b i l i t y scale. For a l l months, probability of higher snow storage increases with a l t i t u d e . Snow water-equivalents at most elevation bands tend to be normally dis t r i b u t e d , as evidenced by figures VII.? and VII.8 (nearly straight l i n e plots on p r o b a b i l i t y paper). Consider the changing snowpack conditions from November to July ( f i g . V I I . 8 ) . Between the altitudes of 700 and 1200 m., the p r o b a b i l i t y of reaching peak storage i s I 138 F i g . V I I . 7 Snowpack w a t e r - e q u i v a l e n t f o r v a r i o u s a l t i t u d e s on 1 e i r c e e d a n c e s t A p r i l p r o b a b i l i t y 139 MA? v I ~ ~| EXCEEDANCE PROBABILITY F i e V I I . 8 M o n t h l y c h a n g e s i n snow w a t e r - e q u i v a l e n t e x c e e d a n c e p r o b a b i l i t y a t 1200 m., b a s e d on 100 y e a r s o f s i m u l a t e d d a t a 140 highest i n May, while for lower a l t i t u d e s , peak storage often arrives e a r l i e r ( f i g . V I I . 9 ) . This corresponds with Schaerer's (1970) observation at Mount Seymour near Vancouver, B.C., that valle y s i t e s reach maximum snow-load i n January or early February, but higher s i t e s (about 1,200 m.) acquire maximum load i n May. The above phenomenon, together with the fact that lower zones generally have higher pr o b a b i l i t y of zero snow storage, suggest that throughout winter, both accumulation and depletion processes are active at lower a l t i t u d i n a l zones. In May and June, when melting i s acceler-ated, most lower zones are usually bare of snow, while higher grounds show greater variations i n pack conditions from year to year ( f i g . VII.8 and V I I . 9 ) . F i e l d observations bear wit-ness to t h i s phenomenon. In 1970, for instance, melting took place rather r a p i d l y and the pack was th i n . In 1971. however, low temperature, cloudiness and wetness maintained a heavy snowpack i n May and June at higher a l t i t u d e s . Snowpack i n the forest i s thinner than i t s counter-parts i n the open, and the year to year changes i n snow storage are smaller than those i n the open, p a r t i c u l a r l y for alt i t u d e s below 1 ,000 m0 ( f i g . V I I . 9 ) . This may be attributed to interception loss which reduces the p r o b a b i l i t y of creating thick snowpacks. VII -3c2 Snowpack persistence The duration of snowpack on the ground was examined. The mean number of days with snow cover increases with 750 F i g . V I I . 9 M e d i a n s and i n t e r - q u a r t i l e r a n g e s o f snow s t o r a g e f o r v a r i o u s a l t i t u d i n a l z o n e s , b a s e d on 100 y e a r s o f s i m u l a t e d d a t a 142 a l t i t u d e so that for elevations above 1,000 m., the ground i s usually covered with snow for over f i v e months ( f i g . VII. 10). Por lower a l t i t u d e s , January to March have the largest number of days with snow cover, the number being smaller for forested areas ( f i g . VII.11). A plot of the standard dev-iation s i n the persistence of snowpack shows that for lower elevations (below 200 m.), February and March are the months when snow duration varies considerably from year to year. This may be attributed to increasing uncertainties i n the occurrence of melt or accretion events. For higher grounds, January to May i s a period of snowpack accretion so that the ground i s constantly snow-bound. These areas have greater variations i n snow-covered days before the pack was f u l l y established (November and December) and during the period of spring melt (after March), This pattern applies to both the open sectors as well as to areas underneath the forest canopy. VII-4 Model l i m i t a t i o n s The simulation exercise i l l u s t r a t e s the p o s s i b i l i t y of extending temperature and p r e c i p i t a t i o n data and of pro-viding p r o b a b i l i s t i c statements regarding snowpack storage i n small basins. It should be emphasized that the transfer functions and the associated c o e f f i c i e n t s i n the simulation model were derived from the coastal mountain areas near Van-couver. Hence, the application of t h i s model should be res-t r i c t e d to small drainage basins not far beyond the coastal be l t , and at elevations below the alpine zone. fD 3 o w O cn 3 p."^ CO 3 (u C n JJ 31 a- < rt- rf C L 3* H -(U ft> QJ rt rt » O (J. 3 H -C5 3 cry a ra CO ro -3 a, ro o CO =1 H . CO t- rt ore 0 3 o ro CO fD S T A N D A R D D E V I A T I O N ! ( D A W ) MEAN SNOWPACK PERSISTENCE (DAHS) - to 3: o Z >H 144 MONTH F * B . V I I . l l Means and v a r i a t i o n s i n snowpack i n t h e f o r e s t , o f s i m u l a t e d d a t a t h e p e r s i s t e n c e o f b a s e d on 100 y e a r s 145 VIII DISCUSSION VIII-1 Model performance While most ex i s t i n g models related to snow hydrology-are deterministic, the present model incorporates a stochastic element which enables a simulation of snowpack data for areas where available records are scarce. By generating long sequences of climatic events which are s t a t i s t i c a l l y i n d i s -tinguishable from h i s t o r i c a l data, i t i s possible to explore a combination of a large number of events that r e s u l t i n hydrologic conditions l i k e l y to be experienced by the area under study. Thus, the model i s most suitable for investigat-ions of the snow hydrology of ungauged basins. The model i s intended for hydrologic studies involv-ing long-term rather than short-term conditions. Simulated data f o r in d i v i d u a l years are therefore of lim i t e d use. On the other hand, the aggregate contribution of i n d i v i d u a l data to the pr o b a b i l i t y statements regarding snow storage are use-f u l for planning purposes. In terms of s p a t i a l considerations, the model i s so constructed that the drainage basin forms the basic unit of investigation. Hence, i t disregards s i t e variations due to aspect, l o c a l ground roughness and other s i t e factors. The model then produces a rather generalized picture of the s p a t i a l d i s t r i b u t i o n s of snow storage. In order to reduce computer time without cutting down model s e n s i t i v i t y , i t i s advisable not to subdivide a l t i t u d i n a l zones to a fineness of less than 2 5 m. 146 i n t e r v a l s . On the whole, however, the algorithm i s rather e f f i c i e n t so that simulation can be performed expediently. VIII-2 Model and management applications The present model i s p a r t i c u l a r l y adaptable to areas of high r e l i e f where water resource d i s t r i b u t i o n exhibits far greater v a r i a b i l i t y than can be estimated by the e x i s t i n g information network. For resource management and planning, long-term p r o b a b i l i t y statements regarding the magnitude and duration of snowpacks over diverse t e r r a i n are more useful than short-range forecasts for c e r t a i n s i t e s . Possible areas of a p p l i c a t i o n of the model include the l o c a t i o n of s k i - or other winter-resort areas, or the assessment of requirements for road clearance in snow-prone zones. When coupled with streamflow routing procedures, the generated runoff model w i l l also enable the evaluation of floods and water storage i n small basins. Dual indices of snowpack conditions i n open and underneath forest canopies provide measures of differences between storage conditions i n the presence or absence of logging practice. In t h i s respect, the model i s useful i n exploring some hydrologic impacts due to changing land use. VIII-3 Model extension Although the number of hydro-meteorological stations i s growing rapid l y i n the coastal areas, t h e i r locations w i l l , for some time, be confined to valley areas. Hence, i n the near future„ the simulation of hydrologic phenomena w i l l prob-147 ably r e l y on information input from v a l l e y stations. The present model can expand i n two directions ( f i g . VIII.1). One of them i s to incorporate other components i n the water-balance equation. Another d i r e c t i o n of extension i s towards the formation of an a l l - b a s i n approach. In solving for other terms i n the water-balance equation, the present model w i l l have to include s o i l moisture, ground water storage and evapotranspiration. This w i l l r e s u l t i n an o v e r a l l water-balance model which encompasses a l l major aspects of water supply, water loss and water storage for coastal basins. On the other hand, the present model can be extended to other snow zones beyond the temperate forest sectors. One of the more immediate areas w i l l be the alpine and g l a c i e r belts which constitute the higher grounds of many Coast Mountain basins. The application of the model to regions outside the coastal areas i s possible i f the forms of the transfer functions (equations 1.3 and 1.4) are established from l o c a l data. Thus, an extension of the model can take the forms of adopting new regional data to the e x i s t i n g frame-work, or an integration of e x i s t i n g framework with other sub-models to y i e l d some comprehensive water-balance models for entire drainage basins. 148 Model Complexify Wrttor resource* 6j»C-*ttvtton Approach Z o n A l ApproAeW AU-lo*5tB AppfoAcW R e g i o n a l Approach Spat ia l dimension F i g . V111 .1 S c h e m a t i c d i a g r a m s h o w i n g p o s s i b l e a r e a s o f e x t e n s i o n f o r t h e p r e s e n t m o d e l 149 REFERENCES Amorocho, J. & Hart, W.E., 1964. A c r i t i q u e of current methods i n hydrologic systems investigation. Trans. Amer.  Geophys. U. 4 5 : 2 , pp. 3 0 7 - 3 2 1 . Anderson, E.A., 1 9 6 8 . Development and te s t i n g of snow pack energy balance equations. Water Resources Res. 4 : 1 , pp. 1 9 - 3 7 . Anderson, E .Ao & Rockwood, D.M., 1 9 7 0 . Runoff synthesis for rain-on-snow basin. 3 8 t h West).Snow Conf. pp. 82 -90. Anderson, R.L., 1 9 4 2 . D i s t r i b u t i o n of the s e r i a l c o r r e l a t i o n c o e f f i c i e n t s . Ann. Math. Stat. 1 3 , pp. 1 -13 . Armstrong, J.E., 1957. S u r f i c i a l geology of New Westminster map-area, B r i t i s h Columbia. Geol 0 Surv. Can. Paper 57-5» 25 P. Armstrong, J.E. & Brown, W.L., 195^. Late Wisconsin marine d r i f t and associated sediments of the lower Fraser v a l l e y , B r i t i s h Columbia, Canada. B u l l . Geol. Soc. Am. 6 5 : 4 , pp. 3 4 9 - 3 6 4 . Bagley, J.M., 1964. An application of stochastic process theory to the r a i n f a l l runoff process. Dept. of C i v i l Eng.  Stanford Univ. Tech. Rept. 3 5 , 118 p. Beard, L.R., 1967 . Hydrologic simulation i n water-yield analysis. Amer. Soc. C i v i l Eng. 93:IR1, pp. 33-42. Beard, L.R., 1 9 6 8 . Hypothetical flood computation for a stream system. Int. Assoc. Sci.' Hydr'ol. Pub. 81, pp. 2 5 8 -2 6 7 . Beaumont, R.T., 1 9 6 5 . Mount Hood pressure pillow snow gage. 3 3 r d West. Snow Conf. pp. 2 9 - 3 5 , Bellman, R.E. & Kalaba, R., I 9 6 5 . Dynamic programming and  modern control theory. Academic" Paperback, 112 p. Benson, M.A., I 9 6 5 . Spurious c o r r e l a t i o n i n hydraulics and hydrology. J. Hydraul. Div., Amer. Soc, C i v i l Eng. 91:HY4, pp. 3 5 - ^ 2 . Bergen, J.D., 1 9 6 9 . Nocturnal r a d i a t i o n loss estimates for a forested canopy. USDA. Forest Service, Res. Note RM-155 , ^ P. 150 Berndt, H.W. & Fowler, W.B., 1969. Rime and hoarfrost i n upper-slope forests of eastern Washington, J. Forest. 6?:2, pp. 92-95. Bhiuya, R.K., 1968. Analysis of periodic hydrologic time  se r i e s . Ph.D. Thesis, Colorado State Univ. Boe, K.N., 1970. Temperature, humidity and p r e c i p i t a t i o n at the Redwood Experimental Forest. USDA, Forest Service  Res. Note PSW-222, 11 p. Bradley, W.C., 1963. Large-scale e x f o l i a t i o n i n massive sandstone of the Colorado Plateau. B u l l . Geol. Soc. Am. 74, pp. 519-528. Brooks, C.E.P. & Carruthers, N.C., 1953. Handbook of  S t a t i s t i c a l methods i n meteorology. H.M. Stationery Office, London, 412 p. Bruce, J.P., I962. Snowmelt contributions to maximum flood. East n Snow Conf. Carlscn, R.F., MacCormick, A.J.A. & Watts, D.G., 1970. Application of lin e a r random models to four annual stream-flow s e r i e s . Water Resources Res. 6s4, pp, 1070-1078. Carson, J.E., 1963. Analysis of s o i l and a i r temperatures by Fourier techniques. J. Geophy. Res. 6818, pp. 2217-2232. Cartwright, C.A.C., 1970. Interception of r a i n f a l l i n immature  coniferous stands. B.S. Thesis, Forest. Faculty, Univ. of B r i t i s h Columbia, Chapman, C.A., 1958. Control of jo i n t i n g by topography. J. Geol. 66, pp. 552-558, Chow, V.T., 1962. Hydrologic determination of waterway areas for the design of drainage structures i n small drainage basins. Univ. of I l l i n o i s B u l l . 462, 104 p. Chow, V.T., 1964. S t a t i s t i c a l p r o b a b i l i t y analysis of hydrol-ogic data, part Is frequency analysis, i n Handbook of  applied hydrology (Ed. Chow), McGraw-Hill, Section 8-1. Chow, V.T., 1 9 6 9 . Stochastic analysis of hydrologic systems. Clearinghouse PB-189791. 34 p. Chow, V.T. & Ramaseshan, S,, I965. Sequential generation of r a i n f a l l and runoff data. J. Hydraul. Div,. Amer. Soc. C i v i l Eng. 9l:HY4, pp. 205 - 2 2 3 . Coffay, E.W., 1962. Throughfall i n a forest. Int. Assoc.  S c i . Hydrol. B u l l . 7»2, pp. 10-16. 151 Conrad, V, & Pollak, L.W,, 1 9 5 0 . Methods i n climatology. Harvard Univ. Pr. 4 5 9 p. Crawford, N.H. & Linsley, R.K., 1 9 6 6 . D i g i t a l simulation i n hydrologys Stanford Watershed Model IV. Dept. of  C i v i l Eng., Stanford Univ. Tech. Rept. 3 9 , 210 p. Curry, L., 1 9 6 2 . Climatic change as a random se r i e s . Ann.  Assoc. Amer. Geog. 5 2 , pp. 2 1 - 3 1 . Dawdy, D.R. & Bergmann, J.M., 1969. E f f e c t of r a i n f a l l v a r i a b i l i t y on streamflow simulation. Water Resources  Res. 5e5t PP. 9 5 8 - 9 6 6 . Dawdy, DeR. & O'Donnell, T., I 9 6 5 . Mathematical models of catchment behavior, J. Hydraul. Div., Amer. Soc. C i v i l  Eng. 91:HY4, pp. 1 2 3 - 1 3 7 . Dickison, R.B.B., 1 9 6 8 . The e f f e c t of minor r e l i e f on the d i s t r i b u t i o n of p r e c i p i t a t i o n . Dept. of Transport, Met.  Branch Tech. Memoir TEC 7 0 0 , 16 p. Dooge, J.C.I., I 9 6 8 . The hydrologic cycle as a closed system. Int. Assoc. S c i . Hydro!. B u l l . 13»lp pp. 5 8 - 6 8 . Eidmann, F.E., 1 9 5 9 . Die Interception i n Buchen- und Fichtenbestanden: Ergebnis mehrjahriger Untersuchungen im Rothaargebirge (Sauerland), Int. Assoc. S c i . Hydro!,.  Pub. 48 t l , pp. 5 - 2 5 . Emshoff, J.R. & Sisson, R.L., 1 9 7 0 . Design and use of com-puter simulation models. MacMillan Co., London, 302 p-. Erickson, D.M'. & McCorquodale, J.A. , 1 9 6 7 . Application of computers to the determination of snowmelt runoff, i n S t a t i s t i c a l methods i n hydrology. Symp0 No. 5 , Queen's Fri n t e r , Ottawa, pp. 3 6 1 - 3 7 8 . Federov, s.F. & Burov, A . S . , 1967* Influence of the forest on p r e c i p i t a t i o n . Soviet Hydrology Selected Paper No. 3 pp. 2 1 7 - 2 2 7 . F e l l e r , W., 1 9 6 8 . An introduction to pr o b a b i l i t y theory and  i t s applications, Vol. T, Wiley, 509 p. Feyerherm, A,M» & Bark, L.D., 1 9 6 7 . Goodness of f i t of a Markov chain model for sequences of wet and dry days, J. Applied'Met. 6 * 5 , pp 0 7 7 0 - 7 7 3 . F i e r i n g , M.B., 1 9 6 4 . A Markov model for low-flow analysis. Int„ Assoc. S c i . Hydro!. B u l l . 9 : 2 , pp. 3 7 - 4 7 . 152 Frank, E.C. & Lee, R., 1966. Potential solar beam i r r a d -i a t i o n on slopes. USDA Forest Service Res. Paper RM-18, 116 p. F r i t t s , H.C., 1961. Analysis of maximum summer temperatures inside and outside of a forest. Ecology 42, pp. 436-440. Gabriel, K.R. & Neumann, J., 1957. On a d i s t r i b u t i o n of wea-ther cycles by length. Quart. J. Roy. Met. Soc. 83, pp. 375-380. Gabriel, K.R. & Neumann, J., 1962. A Markov chain model for da i l y r a i n f a l l occurrence at Tel Aviv. Quart. J. Roy. Met. Soc. 88, pp 90-95. Gartska, W.U., 1964. Snow and snow survey, i n Handbook of applied hydrology (Ed. Chow). McGraw-Hill, Section 10. Gartska, W.U., Love, L.D., Goodell, B.C. & Bertie, F.A., 1958. Factors a f f e c t i n g snowmelt and streamflow. US Govt. Pri'nting Office", 189 p. Gary, H.L. & Coltharp, G.B., I967. Snow accumulation and disappearance by aspect and vegetation type i n the Sante Fe basin,. New Mexico. USDA Forest Service Res. Note RM-93, 11 P. Geiger, R., 1965. The climate near the ground. Harvard Univ. Pr., Cambridge, 6ll™p^ Goodell, B.C., 1959.. Management of forest stands i n western United States to influence the flow of snow-fed streams. Int. Assoc, S c i . Hydrol. Pub. 48, pp. 49-58. Goodell, B.C.* Watt, J.P.C. & Zorich* T.M., I967. Streamflow volumes and hydrographs by flouroscent dyes. Int. U.  Forest. Res. Organ., Iufro-Kongress. 1:01-02-11, pp. 3 2 5 - ^ : Grace, R.A, & Eagleson, P.S., I967. A model for generating synthetic sequences of short-time-interval r a i n f a l l depths. Proc. Int. Hydrol. Symp., Fort C o l l i n s , pp. 268-276. Grayman5 W.M. & Eagleson, P.S., 1969. Streamflow record length for modelling catchment dynamics. Hydrodynamic Lab. Rept. 114, Dept. of C i v i l Eng., M.I.T., 137 p. Green, J.R., 1964. A model for r a i n f a l l occurrences, J. Roy. Stat. Sec. Ser. B 26, pp. 345-353, Harms, A.A. & Campbell, T.H,, 196? 0 An extension of the Thomas-Fiering model for the sequential generation of streamflow. Water Resources Res. 3 : 3 , pp. 653-662. 153 Harr, D.R., 1966. Influence of intercepted water on evapo- transpiration losses from small potted trees. Ph.D. Thesis, Colorado State Univ. Haupt, H . P . , I969. A simple snowmelt lysimeter. Water  Resources Res. 5:3, pp. 714-718. Heigel, K . , i960. Die orographisch bedingte Veranderlichkeit des Niederschlags am Hohenpeissenberg. Int. Assoc. Sc i .  Hydrol. Pub. 53, pp. 278-284. Hendrick, R . L . , Fi lgate , B .D. , & Adams, W.M., 1971. Appl i -cation of environmental analysis to watershed snowmelt. J. Applied Met. 10«3 , pp. 418-429. Hess, S . L . , 1959. Introduction to theoretical meteorology. Hole, Rinehart & Winston, N .Y. , 362 p. H i l l i e r , F .S . & Lieberman, G.J., I967. Introduction to  operations research. Holden-Day, 639 p. Hobbs, P . V . , Radke, L . F . & Shumway, S . E . , 1970. Cloud condensation nuclei from industrial sources and their influences on precipitation in Washington State. J. Atmos, Sc i . 27sl, pp. 81-89, Hoover, M.D. & Leaf, C . F . , I967. Process and significance of interception in Colorado subalpine forest, in Proc. Symp. forest hydrology (Ed. Sopper & L u l l ) . Pergamon Pr. pp. 213-224T Hopkins, J.W. & Robil lard, P . , 1964. Some stat is t ics of daily r a i n f a l l occurrence for the Canadian prairie pro-vinces. J . Applied Met. 3, pp. 600-602. Hovind, E . L . , I965. Precipitation distribution around a windy mountain peak. J . Geophy. Res. 70, pp. 3271-3278. Huff, F .A. & Shipp, W.L. , 1968. Mesoscale spatial var iab i l i ty in mid-western precipitation. J . Applied Met. 7*5,'pp. 886-891. '• Hufschmidt, M.M. & Fiering, M.B. , 1966„ Simulation techniques for design of water resource systems. Harvard Univ. Pr. 212 p. " Ison, N . T . , Feyerherm, A.M. & Bark, L . D . , 1971. Wet period precipitation and the gamma distribution. J , Applied Met. 10s4, pp. 658-665. 1 5 4 Japan Government Forest Experiment Station, Lab, of Snow Damage, 1 9 5 2 . Study of the f a l l e n snow on the forest trees (snow crown). B u l l . Govt. Forest Expt. Station  Japan (Meguro) No. 54, pp. 115-164. Jeffrey, W.W., 1 9 6 8 . Snow hydrology i n the forest environ-ment, i n Snow hydrology workshop seminar. Can. Nat. Comm. Int. Hydrol. Decade, Federicton, pp. 1-19 # Jeng, R.I, & Yevjevich, V,M., 1 9 6 6 . Stochastic properties of lake outflows. Hydrology Pap. 14, Colorado State Univ. 2 1 p. Johnson, O.A. & Boyer, P.B., 1 9 5 9 - Application of snow hydrology to Columbia basin. J. Hydraul. Div., Amer.  Soc. C i v i l Eng. 85tKYl, p. 61-81. Jul i a n , R.W., Yevjevich, V. & Morel-Seytourx, H.J., I967. Prediction of water y i e l d i n high mountain watersheds based on physiography. Hydrology Pap. 2 2 , Colorado State  Univ. 2 0 p. K a r t v e l i s h v i l i , N.A., I969. Theory of stochastic processes  i n hydrology and r i v e r runoff regulation. Clearinghouse, T T 6 8 - 5 0 4 7 8 , 2 2 3 p. Keen,, F.P., 1 9 5 6 . Low temperature d i f f e r e n t i a l s i n one forest drainage. Ecology 3 7 * 2 , pp. 3 8 I - 3 8 3 . K e l l e r h a l s , R., 1 9 7 0 . Runoff concentration i n steep channel networks. Ph.D. Thesis, Univ. of B r i t . Columbia. Kendrew, W.G. & Kerr, D., 1 9 5 5 . The climate of B r i t i s h Columbia and the Yukon T e r r i t o r y . Queen's Printer, Ottawa 22~2 p. Knisel, W.G. J r . & Yevjevich, V,, 1 9 6 7 . The s t a t i s t i c a l measure of hydrologic time series„ Proc. Int. Hydrol.  Symp. y. 1. Colorado State Univ. pp. 3 0 6 - 3 1 3 . Kotz, S . & Neumann, J., I 9 6 3 . On the d i s t r i b u t i o n of pre-c i p i t a t i o n amounts for periods of increasing length. J.  Geophy. Res. 6 8 : 1 2 , pp. 3 6 3 5 - 3 6 4 0 . Krecmer, V,, 1 9 6 7 . Das Mikroklima der Kieferlochkahlschlage, IV. T e i l : v e r t i k a l e Niederschlage, Luftfeuchtigkeit. Wetter und Leben 1 9 , pp. 203-214. Krumbein, W.C, 1 9 6 7 . Fortran IV computer programs for Markov chain experiments in.geology. State Geol. Surv. Univ. of Kansas Tech. Rept. No. 6, AD 6 5 5 - 4 4 7 , 37 p. 155 Kuz'min, P.P., i960. Snow cover and snow reserves. Clear-inghouse, TT61-11467, 140 p. Lacate, D.S., I965. Forest land c l a s s i f i c a t i o n for the University of B r i t i s h Columbia Research Forest. Dept.  of Forest. Pub. No. 1107, 24 p. Lamperti, J., 1966. Probability. Benjamin Paperback, 150 p. Langbein, W.B., 1958. Queuing theory and water storage. J. Hydraul. Div., Amer. Soc. C i v i l Eng. 84:HY5» paper 1811. Larson, C.L., 19^5. A two phase approach to the prediction of peak rates and frequencies of runoff for small ungaged watersheds. Dept. C i v i l Eng., Stanford Univ. Tech. Rept. 5 3 , 109 p. Leaf, C.F., 1969. A e r i a l photographs for operational stream-flow forecasting i n the Colorado Rockies. 3?th West. Snow  Conf. pp. 19-28. Leonard, R.E., 1961. Interception of p r e c i p i t a t i o n by north-ern hardwoods. USDA, Forest Service, NE Forest Expt. Sta. Sta. Paper 159, 16 p. Leonard, R.E., 1967. Mathematical theory of interception, i n Proc. Symp. forest hydrology (Ed. Sopper & L u l l ) . Pergamon Pr. pp. 131-136^ Lichty, R.W. , Dawdy, D.R. & Bergmann, J.M. , I 9 6 8 . R a i n f a l l runoff model for small basin flood hydrograph simulation. Int. Assoc. S c i . Hydrol. Pub. 81:2, pp. 356-367. Linsley, R.K., 1956. The r e l a t i o n between r a i n f a l l i n t e n s i t y and topography i n northern C a l i f o r n i a . Dept. of C i v i l Eng., Stanford Univ. Res. Rept. No. 1, 18 p. Longley, R.W., 1953. The length of dry and wet periods. Quart. J. Roy. Met. Soc. 79, pp. 520-527. Lowry, W.P., 1963. Observations of atmospheric structure during summer i n a coastal mountain basin i n northwest Oregon. J. Applied Met. 2s6, pp. 713-721. L u l l , K,W. & Rushmore, F.M., 196la. Further observations of snow and fr o s t i n the Adirondacks. USDA, Forest Service, NE Forest Expt. Sta. Forest Res. Note 116, 4 p. L u l l , H.W. & Rushmore, F.M., 1961b,. Influence cf forest cover on snow and fros t i n the Adirondacks. 18th East. Snow Conf. pp. 71-79. 156 MacHattie, L.B., 1968 . Kananaskis v a l l e y winds i n summer. J. Applied Met. 7, pp. 3^8 - 3 5 2 . MacHattie, L.B., 1970 . Kananaskis v a l l e y temperatures i n summer. J. Applied Met. 9 : 4 , p p . 574-582 . Markovic, R.D., 1967. Discriminating the change i n means of hydrologic variables. Proc. Int. Hydrol. Symp. Colo- rado State Univ. p p . 581-588. Matalas, N.C., I 9 6 3 . Autocorrelation of r a i n f a l l and stream-flow minimums. U.S. Geol. Surv. Prof. Paper 434-B, 10 p. Matalas, N.C., 1967. Time series analysis. Water Resources Res. 313, p p . 817-829. McKay, G.A., 1964. Relationships between snow survey and cli m a t o l o g i c a l measurements. Int. Assoc. S c i . Hydrol. Pub. 6 3 , p p . 214-227. Meiman, J.R., 1968. Snow accumulation related to elevation, aspect and forest canopy. i n Snow Hydrology workshop  siminar, Can. Nat. Comm. Int. Hydrol, Decade, Fredericton, p p . 35 -^7 . M i l l e r , D.H,, 1964. Interception processes during snowstorms. USDA, Forest Service Res. Paper PSW-18, 24 p . M i l l e r , D.H., 1 9 6 6 . Transport of intercepted snow from trees during storms. USDA, Forest Service Res. Paper PSW-33, 30 p . Miner, N.H. & Trappe, J.M„, 1957. Snow interception, accu-mulation and melt i n lodgepole pine forests i n the Blue Mountains of eastern Oregon. USDA, P a c i f i c NW Forest &  Range Expt. Sta. Res. Note 153» 4 pT Mink, J.F., i960. D i s t r i b u t i o n patterns of r a i n f a l l i n the leeward Koolau Mountains, Oahu, Hawaii. J. Geophy. Res. 6 5 , p p . 2869-2876. Monkhouse, F.J„, I 9 6 5 . A dictionary of Geography. Arnold Ltd., London. 344 p . Moore, D.0.5 I968. Synthesising d a i l y discharge from r a i n -f a l l . J. Hydraul. Div., Amer. Soc. C i v i l Eng. 94:HY5, p p . 1283-1298. Mustonen, S .Eo, 1967. E f f e c t s of climatologic and basin c h a r a c t e r i s t i c s on annual runoff 0 Water Resources Res. 3 : 1 , p p . 123-130. 157 O'Neill, A.D.J. & Gray, D.M., 1971. Energy balance and melt theories, i n Runoff from snow and i c e . Hydrol, Symp, No. 8, Queen's Printer, Ottawa, pp. 31-58. Papoulis, A., 1965. Probability, random variables, and  stochastic processes. McGraw-Hill, 583 p. Pattison, A., 1964. Synthesis of r a i n f a l l data. Dept. of  C i v i l Eng., Stanford Univ. Tech. Rept. 40, 143 p. Payne, K., Neuman, W.R. & K e r r i , K.D., 1969 . Daily stream-flow simulation. J. Hydraul. Div., Amer. Soc. C i v i l Eng. 95»HY4, pp. I I 6 3 - I I 7 9 . Penton, V.E. & Robertson, A.C., I967. Experience with the pressure pillow as a snow measuring device. Water Resources  Res. 3s2, pp. 405-408. Peterson, N.R., 1968. Snow sensors i n C a l i f o r n i a - a progress report. 36th West. Snow Conf. pp. 99-105. Pierce, R.S., 1967, Evidence of overland flow on forested watersheds, i n Proc. Symp, forest hydrology (Ed. Sopper & L u l l ) . Pergamon Pr. pp. 247-254V Popov, E.G., I968. Pr i n c i p l e s of modelling and aspects of using electronic computers i n hydrologic forecasting. Int. Assoc. S c i , Hydrol. Pub. 81, pp. 448-455. Pysklywec, D.W., Davar, K.S. & Bray, D.I., I968. Snowmelt at an index plot. V/ater Resources Res. 4:5, pp. 937-946. Quick, M.C, I965. River flood flows, forecasts and proba-b i l i t i e s . J. Hydraul. Div., Amer. Soc, C i v i l Eng. 91:HY3, pp. 1-18. Quick, M.C, 1971. Experiments with physical snowmelt models, i n Runoff from snow and i c e . Hydrol, Symp. No. 8, Queen's Printer, Ottawa,- pp. 131-150. Quimpo, R.G. & Yevjevich, V., I967. Stochastic description of d a i l y r i v e r flows. Proc,_Int, Hydrol. Symp. v„ 1, Colorado State Univ. pp. 290-297. Raudkivi, A.J. & Lawgun, N., 1970. Synthesis of urban r a i n -f a l l * Water Resources Res. 6:2,, pp. 455-464, Rauner, Y.L., 1962. The heat balance of forests and i t s role i n the formation of the microclimate of wooded and treeless landscapes of the Moscow region 0 Soviet Geography 3, pp. 40-47, 1 5 8 Riley, J.P., Chadwick, D.G. & Eggleston, K.O., 1 9 6 9 . Snowmelt simulation. 3 7 t h West. Snow Conf. pp. 4 9 - 5 6 . Rockwood, D.M., 1961. Columbia basin streamflow routing by computer. Trans. Amer. Soc. C i v i l Eng. Paper 3H9. Rockwood, D.M., Boyer, P.B. & Hildebrand, C E . , 1 9 5 4 . Lysimeter studies of runoff from a deep snowpack. Int.  Assoc. S c i . Hydrol. Pub. 3 9 , pp, 1 3 7 - 1 6 5 , Roddick, J.A., I 9 6 5 . Vancouver north, Coquitlam, and P i t t Lake map-areas, B r i t i s h Columbia. Geol. Surv. Can. Mem. 3 3 5 , 2 7 6 p. Roesner, L.A., 1 9 6 5 . Analysis of time series of monthly  precipitation, and monthly r i v e r flows. M.S. Thesis, Colorado State Univ. Rogerson, T.L. & Byrnes, W.R., 1 9 6 8 . Net r a i n f a l l under hardwoods and red pine i n central Pennsylvania. Water Resources Res. 4 : 1 , pp. 5 5 - 5 7 . Rothacher, J., I .963. Net p r e c i p i t a t i o n under a Douglas f i r forest. Forest S c i . 9, pp. 4 2 3 - 4 2 9 . Rouse, W.R. & Wilson, R.G., I 9 6 9 . Time and space variations i n the radiant energy fluxes over sloping forested t e r r -ain and th e i r influence on seasonal heat and water balan-ces at a middle latitude s i t e , Geografiska Annaler 51:A3, pp. 1 6 0 - 1 7 5 . Rowe, J.S., 1 9 5 9 . Forest regions of Canada. Dept. of Northern A f f a i r s & Nat. Resources, Forest. B u l l . 1 2 3 , 71 p. Rowe, P.B. & Hendrix, T.M., 1 9 5 1 . Interception of r a i n and snow by second-growth ponderosa pine. Trans, Amer. Geophy. U. 3 2 : 6 , pp. 9 0 3 - 9 0 8 . Sabbagh, M„E 0 & Bryson, R.A., 1 9 6 2 . Aspects of the precip-i t a t i o n climatology of Canada investigated by the method of harmonic analysis, Ann. Assoc. Amer. Geog. 5 2 , pp. 4 2 6 - 4 4 7 . Sariahmed, A. & K i s i e l , C.C, 1 9 6 8 . Synthesis of sequences of summer thunderstorms volumes for the Atterbury water-shed i n the Tucson area. Int„ Assoc. S c i . Hvdrol. Pub. 81, pp 0 4 3 9 - 4 4 7 . Sartz, R.S., 1 9 6 6 , R a i n f a l l d i s t r i b u t i o n over dissected t e r r a i n i n southwestern Wisconsin,, Water Resouces Res. 2 : 4 , pp. 8 0 3 - 8 0 9 . 159 Satterlund, D.R. & Eschner, A.R., 1965. The surface geometry of a closed conifer forest i n r e l a t i o n to losses of intercepted snow. USDA Forest Service Res. Paper NE-3^> 16 p. Satterlund, D.R. & Haupt, H.F., 1967. Snow catch by conifer crowns. Water Resources Res. 6:2, pp. 649-652. Satterlund, D.R. & Haupt, H.F., 1970. The d i s p o s i t i o n of snow caught by conifer crowns. Water Resources Res. 6:2, pp, 649-652. Schaerer, P., 1970. Variation of ground snow loads i n B r i t i s h Columbia. 38th West. Snow Conf. pp. 44-48. Schermerhorn, V.P., 1967. Relations between topography and annual p r e c i p i t a t i o n i n western Oregon and Washington. Water Resources Res. 3:3, pp. 707-711. Schernerhorn, V.P. & Kuehl, D.W. 5 1968, Operational stream-flow forecasting with the SSARR model. Int. Assoc. S c i . Hydrol. Pub. 80, pp. 317-328. Seppanen, M,, i960. On the i n f l u e n c e of trees on accumulation of snow on pine dominated forest i n Finland. Int. Assoc.  S c i . Hydrol. Pub. 54, pp. 64-68. Shen, J.„ 1965. Use of analog models i n the analysis of flood runoffe U . S . Geol. Surv ; Prof. Paper 506-G, 24 p. Simpson, W. & Henry, CD., 1966. Dry and wet s p e l l s at Winnipeg. Dept. of Transport C i r . 4507. 15 p. Smith, J.L., Willen, D.W. & Owens, M.S., 1967. Isotope snow gages for determining hydrologic c h a r a c t e r i s t i c s of snow-packs. Geophy. Mono. No, 22 (Ed. Stout), Amer. Geophy. U. pp. 11-21. Slaymaker, H.O., 1969. Scale problems i n hydrology, i n Geography a t Aberytwyth. Univ. o f .Wales Pr,, Ca r d i f f , pp. 6" 8-867 '" ' Sopper, W.E. & L u l l , H.W., 1965. Streamflow c h a r a c t e r i s t i c s of physiographic units i n the n o r t h e a s t c Water Resources Res. 1:1, pp. 115-124. Sporns, U., 1964. The transposition of short-duration r a i n f a l l i n t e n s i t y data i n mountainous regions. Dept. of Transport Me!sJ^JLi£» ^ ° 3 2 » 6 p. Spreen, W.C, 1947. A determination of the ef f e c t of topo-graphy upon p r e c i p i t a t i o n . Trans. Amer. Geophy. U. 2:2, pp, 285-290. 160 Tanner, J.T., 1963. Mountain temperatures i n southeastern and southwestern United States during late spring and early summer. J. Applied Met. 2, pp. 473-483. Thomas, H.A. J r . & Fie r i n g , M.B., 196?. Mathematical syn-thesis of streamflow sequences for the analysis of r i v e r basins for simulation, i n Design of water-resource systems (Ed. Maass et a l . ) , Harvard Univ. Pr., pp. 459-493. Thorud, D.Be, 1963t E f f e c t s of pruning on r a i n f a l l i n t e r -ception i n a Minnesota red pine stand. Forest S c i . 9» pp. 452-455. Todorovic, P. & Yevjevich, V., 1969. Stochastic process of p r e c i p i t a t i o n . Hydrology Paper 35t Colorado State Univ. 61 p. U.S. Army, 1953^ Hydroraeteorological log of the Central  S i e r r a Snow.Laboratory, 1951-52 water year. Corps of Engineers, U.S. Army, 214 p. U.S. Army. 1956. Snow hydrology. Corps of Engineers, U.S. Army, 437 p. Urfer-Henneberger, C , 1964* Wind- und Temperaturverhaltnisse an ungestorten Schonwetterlagen im Dischmatal bei Davos, Mitt. Schv/eiz For. Anstalt 40:6, pp. 384-441. Urfer-Henneberger, C , 1970. Die Sommerniederschlage im Dischmatal bei Davos. Schweiz. Anstalt fur Forst. Ver. 46:2, pp. 69-110, Vemuri, V. & Vemuri, N., 1970. On the systems approach i n hydrology. Int. Assoc. S c i . Hydrol, B u l l , 15s2, pp, 17-38. Verschuren, J.P., 1968. A stochastic analysis of p r e c i p i t a t -i p n 3 Ph„D. Thesis, Colorado State Univ. Waggoner, P,E. & Reifsnyder, W.E., 1968. Simulation of the temperaturej humidity and evaporation p r o f i l e s i n a l e a f canopy. J. Applied Met. 7:3, pp, 400-409. Wakahama, G., 1968. Permeation of melt water into snow cover. Seppyo 30:6, pp. 175-188. Walker, E.F., 1961. A synoptic climatology for parts of the Western C o r d i l l e r a . McGill Univ. A r c t i c Met. Group, S c i .  Reptc No. 8, Pub, i n Met. No. 357 218 p e Walkotten, W.J, & Patr i c , J.H., 1967, Elevation e f f e c t s on r a i n f a l l near H o l l i s , Alaska. USDA Forest Service Res.  Paper PNW-53, 7 p. Watanabe, S. & oseki, Y,, 1964, Study o f fallen.snow on forest trees ( I I ) . Experiment on the nnow (• Town of the Japanese cedar. Japan Govt. Forest E x p t . Sba. B u l l , I69, pp 8 121-140. 161 West, A.J., 1961. Cold a i r drainage i n forest openings, USDA, Forest Service, P a c i f i c SW Forest & Range Expt.  Station Res. Note No. 180, 6 p. Williams, P. Jr. & Peck, E.L., I 9 6 2 . Terrain influences on p r e c i p i t a t i o n i n the intermountain west as r e l a t e d to synoptic s i t u a t i o n s . J. Applied Met. 8 , pp. 3 ^ 3 - 3 4 7 . Wilson, J.A., 1 9 6 6 . Determination and uses of best i n d i v i d -ual sample points on i n d i v i d u a l snow course. 3 4 t h West.  Snow Conf. pp. 82 -86. Woo, M.K., 1972. Comparative effectiveness of three stochas-t i c models of p r e c i p i t a t i o n . B.C. Geographical Ser. No. 1 6 , pp. 1 5 - 3 3 . Wright, J.B., 1 9 6 6 , P r e c i p i t a t i o n patterns over Vancouver c i t y and lower Fraser v a l l e y . Dept. of Transport, Met.  Branch C i r . 4 4 7 4 , 14 p. Yevjevich, V. & Jeng, R.I., 1 9 6 9 . Properties of non-homogeneous hydrologic s e r i e s . Hydrol. Paper 3 2 , Colorado State Univ. 33 P. Young, G.K., Somers, W.P., Pisano, W.C. & F i t c h , W.N., 1 9 6 9 . Assessing upland reservoirs using a d a i l y flow model. Water  Resources Res. 5«2, pp. 3 6 2 - 3 7 9 . Zeman, L.J., 1 9 6 9 . The r e l a t i o n s h i p of p r e c i p i t a t i o n to  runoff from a small westcoast watershed. M.S. Thesis, Univ. of B r i t i s h Columbia. A P P E N D I X A 1 c s i m u l a t i o n o f d a i l y t e m p e r a t u r e and p r e c i p i t a t i o n r a c o r d s c c c i n p u t s f r o m l o g i c a l u n i t 5 c c c o n t r o l c a r d : c c o l 1-6 no. o f y e a r s o f r e c o r d t o be s i m u l a t e d c c o l 7-12 no. o f f i l e s cn l o g i c a l u n i t 7 t o be s k i p p e d c b e f o r e e n c o u n t e r i n g t h a a p p r o p r i a t e f i l e s c c o l 13-18 i n t e r v a l ( i n y e a r s ) a t w h i c h s i m u l a t e d d a t a c a r e t o be p r e s e n t e d ( z s r o f o r no o u t p u t ) c c o l 19-30 a s h o r t t i t l e f o r t h e ba s e s t a t i o n w h i c h p r o -c v i d e s s i m u l a t i o n p a r a m e t e r s c c f o r m a t c a r d : c s p e c i f i e s t h e f o r m a t w i t h w h i c h t e m p e r a t u r e and c p r e c i p i t a t i o n p a r a m e t e r s a r e t o be r e a d i n c c i n p u t o f t e m p e r a t u r e p a r a m e t e r s : c 1. means and v a r i a n c e s o f f o u r i e r component c 2. wave component c 3. markov c o e f f i c i e n t a n d r e s i d u a l p a r a m e t e r s c 4. p r o b a b i l i t y m a t r i x o f temp, r a n g a vs mean d a i l y temp. c c i n p u t o f m o n t h l y p r e c i p . p a r a m e t e r s ( s t a r t i n g f r o m o c t . ) : c f i r s t mean d a i l y p r e c i p i t a t i o n by mo n t h s , t h e n c o n d i t i o n a l c p r o b . f o r t h e o c c u r r e n c e o f wet d a y s f o r 12 months, i n t h e c o r d e r o f p(w|ww), p(w|wd), p(w|dw) and p(w|dd) c c o u t p u t o f d a i l y d a t a ( u n f o r m a t t e d ) f r o m l o g i c a l u n i t 7: c e a c h y e a r i s one f i l e , f i r s t max. t e m p e r a t u r e f o r 365 d a y s c t h e n min. t e m p e r a t u r e f o r 365 d a y s , f i n a l l y p r e c i p i t a t i o n c f o r 365 d a y s c c INTEGER *2 LEN DIMENSION T I T L E ( 3 ) , FMT (20) DIMENSION H ( 2 , 2 , 1 2 ) , P { 1 2 ) , P P ( 3 6 5 ) DIMENSION T H A X ( 3 6 5 ) , THIN ( 3 6 5 ) , LOW ( 1 7 ) , PRAN (17,23) c 1 FORMAT (20A4) 2 FDR MAT (316, 3A4) 10 FORMAT (1B1, 5X, 'SUCCESSFUL SIMULATION BASED ON STATION:' 1, 3A4, / 5X, 'TEMP. S P R E C I P . OUTPUT FROM LOGICAL UNIT 7«, 1/ 5X, * STORED IN F I L E S ' , 1 7 , • TO', 17//) 11 FORMAT (// 5X, 'ERROR OCCURRED AT YEAR', 1 6 , 1*, DURING THE OUTPUT OF MAX. TEMP.'//) 12 FORMAT (// 5X, 'ERROR OCCURRED AT YEAR", 1 6 , 1', DURING THE OUTPUT OF MIN. TE3P.'//) 13 FORMAT (// 5X, 'ERROR OCCURRED AT YEAR', 1 6 , 1', DURING THE OUTPUT OF PRECIPITATION'//) 20 FDRMAT (1H1, // 31X, 'PARAMETERS INPUT *, // 17X, A 2 21 22 23 24 c c c c c c c c c c c c c c c c c 50 c c c c c 1'SIMULATION OF DAILY TEMPERATURE AND PRECIPITATION•, /// 15X, '80. OF YEARS TO BE SIMULATED', 7X, 14, // 5X, 1 'BASED ON DATA FROM', 12X, 3A4/I FORMAT (/5X,'TEMPERATURE PARAMETERS', / 51X, 'MEAN', 5X, 1'VARIANCE',//11X, 'FOURIER COMPONENT ANNUAL MEAN', 6X, 1F7.2, 4X, F7.2, 4X, '(DEG.C)', / 32X, * AMPLITUDE•, 8X, 1F7.2, 4X, F7.2, 4X, ' (DEG. C)', / 32X, 'PHASE ANGLE', 6X, 1F7.3, 4X, F7.3, 4X, '(RADIAN)'/) FORMAT (11X, 'WAVE COMPONENT AUTOCORRELATION', 2X, 1F7.3,/ 34X, 'RESIDUAL', 7X, F7.2, 4X, F7.2, 4X, ' (DEG. C) ' 1,/ 32X, 'WAVE LENGTH', 6X, F7.2, 15X, '(DAY)*, / 32X, 1'AMPL. VS LENGTH', / 34X, •INTERCEPT', 6X, F7.2, / 34X, 1'SLOPE', 10X, F7.2, / 34X, 'STAND. ERROR', 3X, F7.2, 15X, 1 • (DEG. C) » /) FORMAT (11X, 'MARKOV COEFFICEINT', 20X, F7.3, // 11X, 1 'RESIDUAL COMPONENT', 20X, F7.2, 4X, F7.2, 4X, * (DEG.C) •, 1// 11X, »PROB. MATRIX FOR TEMP. RANGE VS MEAN TEMP.:', / 1 20X, 'NOT PRESENTED', /// 5X, 'PRECIP. PARAMETERS:', // 18X, • MONTH DAILY PRECIP (MM)', 5X, 'P(W|WW) P(WJWD)', 14X, 'P (W|DW) P (WJDD) •/) FORMAT ( 9X, 12, 8X, F6-2, 13X, 4(F6.3, 5X)) read in control card READ (5,2) N YEAR, NFILE, NPRINT, TITLE read: variable format for input of temperature and precip. READ (5,1) FMT read in temperature time seri e s parameters fou r i e r component: average, amplitude, phase angle -their means and variances READ (5,FMT) FA, VFA, FAM, VFAM, FAN, VFAN read i n wave components mean wave temp, autocorrelation c o e f f i c i e n t , mean & variance of their residuals; wave length mean durat-ion i n days; wave amplitude's r e l a t i o n with length: intercept, slope and standard error READ (5,FMT) WREG, WAVE, WVAR, WLA, WAA, WAB, WSE markov coef. ; mean and variance of temp, residuals READ (5,FMT) RONE, WIEN, RANVAR proba b i l i t y matrix of temp, range vs mean daily temp, range from 0 to 22 deg. c, read i n as rows; mean temp from -10 to 22 deg c at 2 deg in t e r v a l s , read i n as calumns DO 50 J = 1, 23 READ (5,FMT) (PRAN(I,J), I = 1, 17) input p r e c i p i t a t i o n parameters mean daily p r e c i p i t a t i o n f or each month READ (5,FMT) P conditional p r o b a b i l i t i e s f or wet days subscript 1 for wet condition, 2 for dry condition DO 60 I = 1, 2 A 3 DO 60 J = 1, 2 60 READ (5,FMT) ( W ( I , J , K ) , K = 1, 12) C WRITE (6,20) NYEAR, T I T L E WRITE (6,21) FA, VFA, FAM, VFAM, FAN, VFAN WRITE (6,22) W REG, WAVE, WVAR, WLA, WAA, WAB, WSE WRITE (6,23) RONE, WIEN, RANVAR DO 70 I = 1, 12 70 WRITE (6,2U) I , P (I) , W (1 , 1 ,1) , W (1 , 2 , I ) , W (2 , 1 , I) , W (2 , 2 ,1) c s k i p t o a p p r o p r i a t e f i l e on l o g i c a l u n i t 7 t o a w a i t o u t p u t CALL S K I P (N F I L E , 0, 7) c u n f o r m a t t e d d a i l y d a t a o u t p u t r e q u i r e s (365*4) b y t e s LEN = 1«460 c s e t i n t e r v a l f o r p r i n t e r o u t p u t GAP = NPRINT c c i n i t i a l i z e c o n d i t i o n s f o r t e m p e r a t u r e s i m u l a t i o n c c o n v e r t v a r i a n c e s i n t o s t a n d a r d d e v i a t i o n s WVAR = SQRT (WVAR) RANVAR = SQRT (RANVAR) VFA = SQRT (VFA) VFAM = SQRT (VFAM) VFAN '= SQRT (VFAN) c a c c u m u l a t e p r o b . and f i n d l o w e s t c l a s s w i t h p o s i t i v e p r o b . DO 180 I = 1, 17 c i n d e x i n d i c a t e s w h e t h e r l o w e s t c l a s s w i t h p o s i t i v e p r o b . c has been f o u n d . i f s o , i t i s s e t t o 1 INDEX = 0 I F (PR AH ( I , 1) .EQ.0.0) GO TO 150 LOW (I) = 1 INDEX = 1 150 DO 180 J = 2, 23 K = J - 1 PRAN ( I , J ) = P R A N ( I , J ) + PRAN(I,K) I F ( P R A N ( I , J ) .EQ.0.0) GO TO 180 I F (PRAN ( I , J ) .GT. 1. 0) P R A N ( I , J ) = 1. I F (INDEX. EQ. 1) GO TO 180 LOW (I ) = J INDEX = 1 180 CONTINUE c i n i t i a l i z e tm (temp r e s i d u a l f o r markov c o m p o n e n t ) , c and wmean (mean t e m p e r a t u r e o f l a s t tamp, c y c l e ) TM = 0. WMEAN = 0. c c s a t t i m e c o u n t e r s 'myear', 'month 1 and 'mday' c i n i t i a l i z e c o n d i t i o n s f o r p r e c i p i t a t i o n s i m u l a t i o n MYEAR = 1 MONTH = 1 MDAY = 2 c c s e t c o n d i t i o n s f o r f i r s t t w o d a y s , a s s u m i n g .5 p r e c . p r o b . c a wet day h a s i n d e x 1, and a d r y day i s 2. i i i and i i A 4 c refer to the l a s t 2 days and now refers to the present day II = 2 PP(1) = 0. IF (RANUM (1) .LT.O. 5) GO TO 190 PP(1) = EXPDIS(P(1)) II = 1 190 SOW = 2 PP(2) = 0. IF (RANUM (1) .LT.O. 5) GO TO 200 NOW = 1 PP(2) = EXPDIS(P(1)) c c simulation of d a i l y temperature 200 FAVE = GAUSS(FA, VFA) FPHI = GAUSS(FAN, VFAN) FAMP = GAUSS(FAM, VFAM) c CALL TEMPER (FAVE, FAMP, FPHI, WREG, WAVE, WVAR, WLA, WAA, 1WAB, WSE, RONE, RANVAR, PRAN, LOW, TMAX, TMIN, TM, WMEAN) c c simulate daily p r e c i p i t a t i o n 210 MDAY = MDAY + 1 c reset past condition for updating III = II II = NOW Y = W (II, III ,MONTH) Z = RANUM (1) IF (Z.LE.Y) GO TO 220 NOW = 2 PP(MDAY) = 0. GO TO 250 220 NOW = 1 PP(MDAY) = EXPDIS (P(MONTH)) c check month counter 250 TIME = MDAY MONTH = TIME / 30. 4 + 1.0 IF (MONTH.LE. 12) GO TO 210 c c output daily data from l o g i c a l unit 7 c maximun and minimum data CALL WRITE (TMAX, LEN, 0, LNR, 7, &610) CALL WRITE (TMIN, LEN, 0, LNR, 7, S620) c p r e c i p i t a t i o n data CALL WRITE (PP, LEN, 0, LNR, 7, &630) END FILE 7 c c check for printer output IF (NPRINT.EQ.0) GO TO 500 c check i f i t i s the appropriate year for output PTIME = MYEAR PTIME = AMOD(PTIME, GAP) IF (PTIME.NE.0.) GO TO 500 C A 5 c output da i l y values i n graphical form c output temperature CALL GRAPH (1, TM AX, T H I N , TITLE) CALL GRAPH (2, PP, PP, TITLE) c c update year counter 500 MYEAR = MYEAR + 1 MONTH = 1 MDAY = 0 IF {MYEAR. L E.N YEAR) GO TO 200 c c put in an extra e n d - o f - f i l e mark from l o g i c a l unit 7 END FILE 7 REWIND 7 c c l i m i t s of f i l e s punched from l o g i c a l unit 7 I = NFILE + 1 J = NFILE • NYEAR WRITE (6,10) TITLE, NYEAR, I, J STOP c c error returns 610 WRITE (6,11) MYEAR STOP 620 WRITE (6,12) MYEAR STOP 630 WRITE (6,13) MYEAR STOP c EH D SUBROUTINE TEMPER (FAVE, FAMP, FPHI, WREG, WAVE, WVAR, 1 WLA, WAA, WAB, WSE, RONE, 8 ANVAR, PRAN, LOW, 1 T, TMIN, TM, WMEAN) c c simulation of d a i l y maximum and minimum temperatures c input includes: c f o u r i e r parameters - mean (fave), amplitude (famp) and c phase angle (fphi) c wave parameters - autocorrelation for mean wave tamp, c (wreg), mean S standard deviation of mean wave temp, c (wave and w a r ) , mean wave-length (wla), regression c slope (wab) & intercept (waa) , and standard error for c length-amplitude r e l a t i o n (wse) c markov component - f i r s t order markov c o e f f i c i e n t (rone) c residual white noise component - zero mean and standard c deviation (ranvar) c output includes maximum (t) 5 minimum (train) temperatures, c mean wave temp, of the l a s t period (tm), and c markov component of the l a s t day (wmean) A 6 c DIMENSION T (365) , TMIN (365) , PRAN (17, 23), L0W{17) c c i n i t i a l i z e 1 = 0 DAY = 0. GO TO 15 c generate temperature waves c generate the l a s t point on temperature wave to be averaged c with f i r s t point of next wave 10 TWLAST = WMEAN + WAMP c cycle lengths generated according to a truncated poisson c d i s t r i b u t i o n (truncation at three day level) 15 L = POISN (WLA) IF (L.LT.3) GO TO 15 WLEN = L c regression to obtain wave amplitude from wave length TEMP = GAUSS (0., WSE) WAMP = WAA + WAB*WLEN + TEMP c generate mean wave temperature c find autoregressive r e l a t i o n with mean temp, of la s t wave WMEAN = WMEAN*WREG c add residual component (according to normal distribution) TEMP = GAUSS (WAVE, WVAR) WMEAN = WMEAN + TEMP c c obtain daily temperature form a linear model: c t = t f + tw + tm + e WL = 0. c daily increment, one day i s two-pi/365 20 1 = 1 + 1 DAY = DAY + 0.017214 c fourier expansion TF = FAVE + FAMP * SIN (DAY+FPHI) c expand components of temperature waves c check i f i t i s the f i r s t day of a cycle IF (WL.GT.0.0 .OR. I.EQ. 1) GO TO 30 c average day zero of a wave with l a s t day of previous wave TW = WMEAN + WAMP TW = (T W+T WL AST) * 0. 5 GO TO 35 30 TW AVE = WL / WLEN TW = WAMP * COS(TWAVE*6.2832) + WMEAN 35 T(I) = TF + TW C c markov component TM = TM * RONE C c white noise component WIENER = GAJSS (0., RANVAR) TM = TM + WIENER T(I) = T (I) + TM c A 7 c expand daily temperature into max and min c find appropriate temperature c l a s s from -10 to 22 deg, c at 2 deg. i n t e r v a l s J = (T (I) + 10.0) * 0. 5 + 1. IF (J.LT. 1) J = 1 IF (J.GT. 17) J = 17 c find lowest class with positive probability MIN = LOW (J) Z = RANUM (1) DO 40 JJ = MIN, 23 IF (Z.LE.PRAN (J, JJ) ) GO TO 50 40 CONTINUE c obtain appropriate temp, range; c subtract one since class 1 i s 0 50 R = JJ - 1 HALF = R * 0.5 TM IN (I) = T (I) - HALF T(I) = T (I) + HALF c IF (I.EQ. 1) GO TO 80 c reset the impossible case when max. of one day i s lower c than the min. of previous day and vice versa II = I - 1 IF (TMIN (I) .LE.T (II) ) GO TO 60 c min. of day i greater than max. of day (i-1) TEMP = (TMIN (I)-T (II) ) * 0.5 TMIN (I) = TMIN (I) - TEMP T(II) = TMIN (I) GO TO 70 60 IF (TMIN (II) .LE. T (I) ) GO TO 70 c min. of day (i-1) greater than max. of day i TEMP = (TMIN (II) -T (I) ) * 0.5 TMIN (II) = TMIN (II) - TEMP T(I) = TMIN (II) c c check day counter 70 IF (I.EQ.365) GO TO 100 c chenk i f cycle expires 80 WL = WL + 1. IF (WL.GE.WLEN) GO TO 10 c i f not GO TO 20 c 100 RETURN END SUBROUTINE GRAPH (NTYPE, T, TF, TITLE) c c subporgramme to produce printer plot c input: ntype s p e c i f i e s the type of data A 8 c c c c c c c c 1 5 6 c c c 10 15 16 17 c c 20 30 1 for temperature, 2 for p r e c i p i t a t i o n , 3 for discharge. (t) i s max. temp., precip. or discharge, depending on type of output. (tf) i s minimum temp, or a dummy variable for other types, a short t i t l e (3a4) labels the base-station one year of data requires 3 consecutive pages of output DIMENSION T(365), TF(365), P(50,122), TITLE (3) DATA BLANK, POINT, STROKE, DASH / 1H , 1H*, 1H], 1H-FORMAT (1H1, 5X, 'SIMULATED TEMPERATURE PLOT',/ 2X, 1'BASED ON DATA FROM ', 3A4, 6X, 'DAYS », 13, « TO ', 13, 1// 2X, TEMPERATURE IN DEG. C'/) FDRMAT (1H1, 5X, 'SIMULATED PRECIPITATION PLOT',/ 2X, 1'BASED ON DATA FROM ', 3A4, 6X, 'DAYS ', 13, • TO ', 13, 1// 2X, * PRECIPITATION IN MM'/) FORMAT (1H1, 5X, 'SIMULATED DISCHARGE PLOT'/ 2X, 1'BASED ON DATA FROM ', 3A4, 6X, 'DAYS ', 13, • TO », 13, 1// 2X, 'DISCHARGE IN MM*/) FORMAT (1X, 13, 2H |, 122A1) FORMAT (5X, 1H|, 122A1) JL = 1 JU = 122 JJJ = 122 LOOP = 1 check for the type of data to ba plotted IF (NTYPE - 2) 15, 16, 17 WRITE (6,1) TITLE, JL, JU GO TO 20 WRITE (6,2) TITLE, JL, JU GO TO 20 WRITE (6,3) TITLE, JL, JU clear plot f i e l d DO 30 I = 1, 50 DO 30 J = 1, 122 P(I,J) = BLANK DO 80 J = 1, JJJ L = JL + J - 1 IF (NTYPE.NE.1) GO TO 50 set dash for temperature time axis P(35,J) = DASH scale max. temp, for p l o t t i n g , 1 d i v i s i o n to 1 deg. c K = 36. - T (L) IF (K.LT. 1) K = 1 IF (K. GT.50 ) k = scale minimum temp. KK = 36. - TF (L) IF (KK. LT. 1) KK = 1 IF (KK.GT.50) KK = 50 set symbol between max and min temp values 50 s i m i l a r l y k 9 DO 40 M = K, KK 40 P(M,J) = STROKE GO TO 80 c c scale p r e c i p i t a t i o n and discharge for p l o t t i n g ; c one d i v i s i o n to 3 mm 50 IF (T (L) . EQ.0.0) GO TO 75 K = 50.0 - T (L)*0.333333 IF (K. LE. 1) K = 1 DO 60 M = K, 50 60 P(M,J) = POINT GO TO 80 c dry condition, set axis 75 P(50,J) = DASH 80 CONTINUE c c plot graph DO 100 I = 1, 50 IF (I.EQ.5 . OR. I.EQ.20 '.OR. I.EQ.35 .OR. I.EQ.50)SO TO 90 WRITE (6,6) (P (I , J) , J= 1, JJJ) GO TO 100 90 IF (NTYPE. NE. 1) GO TO 92 II = 35 - I GO TO 95 92 II = (50 - I) * 3 95 WRITE (6,5) I I , (P(I,J), J = 1, JJJ) 100 CONTINUE c c increment counters JL = JU + 1 LOOP = LOOP + 1 IF (LOOP - 3) 110, 120, 130 110 JU = JL + 121 GO TO 10 120 JU = JL + 120 JJJ = 121 GO TO 10 130 RETURN END FUNCTION RANUM (IDUMMY) C c function subprogramme to generate uniformly distributed c random numbers by mixed multiplicative congruential method c r e f . : carnahan et a l i a , 'applied numerical methods', p 549 c DATA M, N, J / 1048576, 566387, 1027 / DATA XM / 1048576.0 / c N = MOD (J*N, M) A 10 RANUM = FLOAT (N) / XM RETURN END FUNCTION GAUSS (EX, STAN) c c function subprogramme to generate random numbers following c gaussian d i s t r i b u t i o n c r e f . : naylor et a l i a , •computer simulation tech.', p 95 c SUM = 0. DO 5 I = 1, 12 5 SUM = SUM + RANUM (1) GAUSS = STAN * (SUM - 6.0) + EX RETURN END FUNCTION POISN (AVE) C c function subprogramme to generate random numbers following c poisson d i s t r i b u t i o n c r e f . : naylor et a l i a , 'computer simulation tech.', p 114 POISN = 0. B = EXP (-AVE) TR = 1. 5 Z = RANUM (1) TR = TR * Z IF (TR.LT.B) GO TO 10 POISN = POISN + 1. GO TO 5 10 RETURN END FUNCTION EXPDIS (EX) C c function subprogramme to generate random numbers following c exponential d i s t r i b u t i o n c ref: naylor et a l i a , 'computer simulation tech.', p 84 c Z = RANUM (1) EXPDIS = - EX*ALOG (Z) RETURN END A 11 c programme to simulate snowpack conditions and da i l y runoff c c c inputs: c c control card: c c o l . 1-6 no. of years of record to be simulated c c o l . 7-12 no. of a l t i t u d i n a l zones (up to 100) c c o l . 13-18 i d e n t i f i c a t i o n for snowpack or runoff to c be punched from l o g i c a l unit 7: c 1 for snowpack, 2 for runoff c c o l . 19-24 no. of f i l e s to be skipped on l o g i c a l c unit 1 (input) c c o l . 25-30 no. of f i l e s to be skipped on l o g i c a l c unit 7 (output) c c o l . 3 1-36 time i n t e r v a l (in years) at which pack c information i s to be presented (set to c blank i f output i s not required) c c o l . 37-42 Say i n t e r v a l at which snowpack i s to be c plotted (ignore i f plot i s not required) c c o l . 43-48 height of base-station i n metres (f6. 1) c c o l . 49-54 height of lowest zone i n metres (f6.1) c c o l . 55-60 zonal height range i n whole number of c metres (f6.1) c c o l . 6 1-72 name of station based upon whose data c simulation i s to be performed c c data input: c max. and min. temperatures and pr e c i p i t a t i o n are read i n c as unformatted data from l o g i c a l unit 1 (these may be c outputs from temp.-precip. simulation programme) c c output: c one year of information w i l l occupy one f i l e (unformatted c data output from l o g i c a l unit 7) c f i r s t output the t o t a l no. of zones involved c then output either open and forest snowpack conditions, or c generated runoff from these sectors c snowpack output indicates the zone numbers of the lowest c zones to which the pack extends, followed by zonal pack c water equivalents (one block per month) c generated runoff output indicates the zone numbers of the c highest contributive zones, followed by runoff data (also c one block per month) c c INTEGER*2 LEN1, LEN7 DIMENSION 10(2,100), TF (2,100), TMAX(365), TMIN(365) DIMENSION P(365), P0(100), PF(100), Q0(100), QF{100), 1 PACKO(100), PACKF(100), HT(100) DIMENSION PP(100), HEAD (100), TITLE (3) , O0T(6262) c A 12 DATA BLANK, STROKE / 1H , 1H| / c 1 FORMAT (716, 3F6. 1, 3A4) 2 FORMAT (1H1, / 32X, 'PARAMETERS INPUT', // 18X, 1'SIMULATION OF SNOWPACK AMD GENERATED RUNOFF', /// 5X, 1'NO. OF YEARS TO BE SIMULATED' , 6X, I'4 , // 5X, 1 'BASED ON DATA FROM',12X,3A4,//5X,'HEIGHT OF BASE-STATION * 1, 9X, F7. 1, ' M.', /// 5X, * MID-ZONE ALTITUDES (IN M.):*, 1// 5X, 5 ('ZONE ALTITUDE ')) 3 FORMAT {/ 5 (18, 2X, F 8. 1)) 4 FORMAT (1H1,5X, 'SUCCESSFUL SIMULATION BASED ON DATA FROM ' 13A4,//4X,•OUTPUT FROM LOGICAL UNIT 7 WERE STORED IN FILES' 1, 17, • TO', 17//) 5 FORMAT (1H1,// 6X, 'SIMULATED SNOW W.E. (IN MM), YEAR', 114, / 6X, 'BASED ON DATA FROM ', 3A4,// 5X, 1'0 > 0, 1 >100, 2 >200, 3 >300, 4 >400, 5 >500' 1, ', 6 >600, 7 >700, 8 >800, 9 >900,', / 5X, 1'A >1000, B >1100, C >1200, D >1300, E >1400, «, 1»F >1500, G >1600, H >1700, I >1800, J >1900', // 2X, 1'DAYS SINCE*, 10X, 'ALTITUDINAL ZONES,', F8.0, • TO', 1F8.0, • M.«, / 3X, '1ST OCT.', 2X, 100A1) 6 FORMAT ('4X, 14, 5X, 100A1) 7 FORMAT (// 5X, 'ERROR OCCURRED AT YEAR *, 14, 1' DURING INPUT OF MAX. TEMPERATURE' //) 8 FDRMAT (// 5X, 'ERROR OCCURRED AT YEAR', 14, 1« DURING INPUT OF MIN. TEMPERATURE' //) 9 FORMAT (// 5X, 'ERROR OCCURRED AT YEAR', 14, 1' DURING INPUT OF PRECIPITATION'//) 10 FORMAT (// 5X, 'ERROR OCCURRED AT YEAR', 14, « DAY', 14, 1 ' DURING OUTPUT OF GENERATED RUN OFF * //) 11 FORMAT (// 5X, 'ERROR OCCURRED AT YEAR', 14, ' DAY', 14, 1« DURING OUTPUT OF SNOWPACK' //) C c read i n control card READ (5,1) NYEAR, NZONE, NSQ, NFILE1, NFILE7, NPRINT, 1 NPKDAY, HTBASE, HTM IN, R AN 31E, TITLE ZONE = NZONE c c skip to appropriate f i l e s for input and output IF (NFILE1.EQ.O) GO TO 20 CALL SKIP (NFILE1,0,1) 20 IF (NFILE7.EQ.O) GO TO 30 CALL SKIP (NFILE7, 0, 7) c c find mid-zone a l t i t u d e of lowest zone 30 LOW = HTMIN / RANGE HB = LOW HB = HB*RANGE +RANGE*0.5 c i n i t i a l i z e snowpack and find mid-zone a l t i t u d e s DO 200 J = 1, NZONE HT(J) = HB HB = HB + RANGE PACKO(J) = 0. A 1 3 200 PACKF (J) = 0. c s e t l o w e s t l i m i t o f s n o w - l i n e t o a z o n e a b o v e t h e b a s i n MINPKO = NZONE + 1 MINPKF = NZONE + 1 c c p l o t t i n g s nowpack c c h e c k i f p r i n t e r p l o t i s r e q u i r e d NPLOT = 1 IF (NPRINT.EQ.O) GO TO 300 c i f day i n t e r v a l f o r p a c k o u t p u t i s b l a n k , s e t t o 1 I F (NPKDAY,EQ.O) NPKDAY = 1 c i n c r e m e n t a l u n i t s f o r p l o t t i n g snowpack c o n d i t i o n s INCPL = 5 I F (NZONE.GT.20) INCPL = 3 IF (NZONE.GT.30) INCPL = 2 I F (NZONE.GT.50) INCPL = 1 MAXPLT = INCPL * NZONE c s e t h e a d i n g f o r snowpack p l o t and c l e a r f i e l d f o r p l o t t i n g DO 250 J = 1, 100 HEAD (J) = BLANK 250 P P { J ) = BLANK DO 260 J = I N C P L , MAXPLT, INCPL 260 HEAD (J) = STROKE c c p r i n t o u t i n p u t d a t a 300 WRITE (6,2) N YEA R, T I T L E , HTBASE WRITE (6,3) ( J , H T ( J ) , J = 1, N20NE) c DO 2000 MY = 1, NYEAR c c r e a d i n u n f o r m a t t e d d a t a f r o m l o g i c a l u n i t 1 CALL READ (TMAX, L E N 1 , 0, LNR, 1, S2010) CALL READ (TMIN, LEN1, 0, LNR, 1, S2020) CALL READ ( P , L E N 1 , 0, LNR, 1, S2030) c c punch t o t a l number o f z o n e s t o be computed LEN7 = 4 CALL WRITE (ZONE, LEN7, 0, LNR, 7, &2050) c i n i t i a l i z e o u t p u t i n d e x LOUT = 0 c c c h e c k f o r p l o t t i n g 400 NP = 1 IF (NPRINT.EQ.O) GO TO 450 NPLOT = MY - (MY/NPRINT)*NPRINT IF (NPLOT.GT.O) GO TO 450 WRITE (6,5) MY, T I T L E , H T ( 1 ) , HT(NZONE), HEAD c c s e t day c o u n t e r * 450 1 = 0 c s e t . m o n t h c o u n t e r MONTH = 1 c A 14 5 00 1 = 1 + 1 c c compute zonal temperatures CALL ZTEMP (TMAX (I), TMIDi (I) , HTBASE, NZONE, HT, TO, IF) c c compute zonal precipitations CALL ZPBEC (P(I), HTBASE, NZONE, HT, TO, PO, PF, MAXR) c c check for plotting IF (HPLOT.GT.0) GO TO 600 NP = I - (I/NPKDAY) *NPKDAY c c update snowpack in open and compute generated discharge 600 CALL ZSNOW (NZONE, PO, MAXR, PACKO, MINPKO, MAXRO, QO, TO, 1 NP, INCPL, PP) c update snowpack in forest and compute generated discharge CALL ZSNOW (NZONE, PF, MAXR, PA3KF, MINPKF, MAXRF, QF, TF, 1 1, INCPL, PP) c IF (NSQ.NE.2) GO TO 1100 c output zonal generated runoff (out) from l o g i c a l unit 7 c f i r s t output zonal runoff from open, then from forest c output no. of zones that generated runoff LOUT = LOUT + 1 OUT(LOUT) = MAXRO c i f there i s no generated runoff IF (MAXRO.EQ.O) GO TO 1020 DO 1010 L = 1, MAXRO LOUT = LOUT + 1 1010 OUT(LOUT) = QO(L) c output forest generated runoff 1020 LOUT = LOUT + i OUT (LOUT) = MAXRF IF (MAXRF.EQ.O) GO TO 1040 DO 1030 L = 1, MAXRF LOUT = LOUT + 1 1030 OUT(LOUT) = QF(L) 1040 MTEMP = FLOAT (I)/30. 4 IF (MTEMP.NE.MONTH) GO TO 1200 c monthly data output c length of unformatted output i s •len7' bytes LEN7 = LOUT*4 CALL WRITE (OUT, LEN 7, 0, LNR, 7, &204 0) LOUT = 0 MONTH = MONTH + 1 c 1100 IF (NSQ.NE.1) GO TO 1200 c output snowpack information (out) from l o g i c a l unit 7 c f i r s t output zonal pack condition in the open LOUT = LOUT + 1 OUT(LOUT) = MINPKO c i f there i s no snow i n the open IF (MINPKO.GT.NZONE) GO TO 1120 k 15 M = MINPKO 1110 LOUT = LOUT + 1 OUT(LOUT) = PACKO (M) M = M + 1 IF (M.LE.NZONE) GO TO 1110 c output forest snowpack 1120 LOUT = LOUT + 1 OUT(LOUT) = MINPKF c i f there i s no snow in the forest IF (MINPKF.GT.NZONE) GO TO 1150 M = MINPKF 1130 LOUT = LOUT + 1 OUT(LOUT) = PACKF (M) M = M + 1 IF (M.LE.NZONE) GO TO 1130 1150 MTEMP = FLOAT (I)/30.4 IF (MTEMP.NE.MONTH) GO TO 1200 c monthly data output c length of unformatted output i s 'len7* bytes LEN7 = LOUT*4 CALL WRITE (OUT, LEN7, 0, LNR, 7, S2050) LOUT = 0 MONTH = MONTH + 1 c c plot pack 1200 IF (NP.GT.0) GO TO 1500 IF (MINPKO.GT.NZONE) GO TO 1400 WRITE (6,6) I, PP GO TO 1500 1400 WRITE (6,6) I c c check day counter 1500 IF (I.LT.365) GO TO 500 c c put an end f i l e mark to the end of a year END FILE 7 c c read beyond endfile mark on l o g i c a l unit 1 CALL READ (TMAX, LEN1, 0, LNR, 1, 52000) c 2000 CONTINUE c REWIND 1 c put an extra endfile mark to output unit END FILE 7 REWIND 7 I = NFILE7 * 1 J = NFILE7 * NYEAR WRITE (6,4) TITLE, I, J STOP c c error returns 2010 WRITE (6,7) MY A 16 STOP 2020 WRITE STOP (6,8) MY 2030 WRITE STOP (6,9) MY 2040 WRITE STOP (6,10) MY, I 2050 WRITE STOP (6,11) MY, I c END SUBROUTINE ZTEMP (TMAX, TMIN, HTBASE, N, HT, TO, TF) c c programme for zonal temperature evaluation c input base-station max (tb(1)) and minimum (tb(2)) temp, c ht. of base-station (htbase) , mean ht. of each zone (ht) c up to (n) zones c returns max. & min. temperatures of various zones i n open: c t o ( i , j ) , i = 1,2 as max and min temp., j = 1, n zones c max. and min. temp, i n forest of various zones (tf ( i , j ) ) c DIMENSION HT(100), TO (2, 100) , TF(2,100), TLRSAX (9) , 1 TLRMIN (8) , FSLOPE (2) , FORIG (2) C c lapse rate for max and min temperatures DATA TLRMAX /-0.0069, -0.0107, -0.0088, -0.0078, -0. 0095, 1 -0.0 108, -0.0094, -0.0084, -0.0066/ DATA TLRMIN /-0.0028, -0.0047, -0.0058, -0.0072, -0.0055, 1 -0.0050, -0.0052, -0.0077/ c regression c o e f f i c i e n t s for forest-open relationships c for max. and min. temperatures DATA FSLOPE /0.9 11, 0.957/ DATA FORIG /0.40, 1.61/ c c select appropriate lapse rate c for max. temp, a l l temp, below 0 deg use rate 1, otherwise c change for every 3 deg except for over 27 deg. J = TMAX * 0.333333 + 1. IF (J. LT. 1) J = 1 IF (J.GT.9) J = 9 c for min. temp, a l l temp below -6 deg use rate 1, otherwise c change for every 3 deg except for over 18 deg. K = (TMIN+6.) * 0. 33333 + 1. IF (K. LT. 1) K = 1 IF (K.GT.8) K = 8 RATEUP = TLRMAX (3) RATEDN = TLRMIN (K) c c route through a l t i t u d i n a l zones & 17 1 = 0 20 1 = 1 + 1 c c find maximum temperature lapse rates TO (1,1) = (HT (I) - HTBASE) * HA TE UP + TM AX c find minimum temperature lapse rates TO (2,1) = (HT (I) - HTBASE) * RATEDN + TMIN c c ensure max temp i s not lower than min temp IF (TO (1,1) .GE.TO (2,1) ) GO TO 40 DIFF = ABS ((TO (1,1) - TO (2,1) ) * 0.5) TO (1,1) = TO (1,1) • DIFF TO (2,1) = TO (1 ,1) 40 DO D5 K = 1, 2 c find temp, i n forest based on temp, i n open 45 TF(K,I) = FORIG(K) + FSLOPE (K) *TO (K, I) IF (TF (1,I).GE.TF(2,I)) GOTO 50 DIFF = ABS((TF(1,I) - TF(2,I)) * 0.5) TF (1 , I) = IF (1 ,1) + DIFF TF(2,I) = TF(1,I) c c check i f a l l zones have been routed through 50 IF (I.LT.N) GO TO 20 c RETURN END SUBROUTINE ZPREC (PB, HTBASE, N, HT, TO, PO, PF, MAX) c c programme to find zonal precipitation c input p r e c i p i t a t i o n (pb) and ht (htbase) of base station c mean heights (ht) of various zones, up to (n) zones c and max. and min. temp, of various zones (to) c output zonal p r e c i p i t a t i o n in open (po) and in forest (pf) c max i s highest zone with r a i n f a l l c max i s 0 when a l l zones have snow, and below 0 i f a l l c zones are dry ( i . e . dry day) c DIMENSION HT(100), TO (2, 100), PO(100), PF(100), PSLOPE(8), 1 RB(8), TSNOW(9), TMIN (9) c c slopes of regression between alt i t u d e and pre c i p i t a t i o n DATA PSLOPE /0.00 18, 0. 025 1, 0.0377, 0. 0791, 0. 1240, 1 0,1465, 0.0691, 0.0757/ c lower bounds for p r e c i p i t a t i o n classes used in regression DATA RB/0., 10., 20., 30., 40., 60., 90., 120./ c slopes for t h r o u g h f a l l - r a i n f a l l and throughfall-snowfall DATA FRAIN, FSNOW / 0.7355, 0.0101 / c cumulative prob. of snowfall given a certain min. temp. DATA TSNOW / .032, .095, .161, .254, .368, .507, .666, A 18 1 .797, 1.0 / c min. daily temp, below which various snowfall prob. apply DATA THIN /-4., -3.5, -3., -2.5, -2., -1.5, -1., -0.5, 0./ c IF (PB.GT.0.0) GO TO 20 c set a l l values to zero i f i t i s a dry flay DO 10 I = 1, N PO(I) = 0. 10 PF{I) = 0. c max. ra i n zone i s set to negative for a dry day MAX = - 1 RETURN c c find suitable p r e c i p i t a t i o n c l a s s (j) 20 J = 1 DO 30 I = 2, 8 IF (PB. LE. RB (I) ) GO TO 40 J = I 30 CONTINUE 40 SLOPE = PSLOPE(J) c c route p r e c i p i t a t i o n i n open by a l t i t u d i n a l zones DO 50 I = 1, N PO(I) = PB + SLOPE* (HT (I)-HTBASE) c set to zero i f negative value generated IF (PO (I).LT.0.) PO(I) = 0. 50 CONTINUE c c check for ra i n or snow MAX = N c min. temp, for highest zone >0 deg., a l l zones have ra i n IF (TO (2 , N) . GT. 0.0) GO TO 260 c a l l zones have snow when lowest zone has min. temp. <-4.5 MIN = 1 IF (TO (2, 1) .LE.-4. 5) GO TO 150 c c otherwise f i n d snow-rain boundary. once a zone has snow, c a l l zones above i t w i l l assume snow condition c c find temperature at which snow occurs TEMP = RANUM (1) c find appropriate min. temperature to match the prob. class LOW = 1 100 IF (TEMP. LE.TSNOW (LOW) ) GO TO 110 LOW = LOW + 1 GO TO 100 c set 'temp* as the min. temperature below which snow occurs 110 TEMP = TMIN (LOW) c check i f the highest zone has temp, above c r i t i c a l value IF (TO (2, N) .GT.TEMP) GO TO 260 120 IF (TEMP.GE.TO(2,MIN)) GOTO 150 MIN = MIN + 1 GO TO 120 A 19 c c find p r e c i p i t a t i o n in forest for a l l snow zones 150 DO 200 K = HIS, N PF(K) = PO (K) * PO (K) *FS NOW c set to zero i f negative value generated IF (PF (K) . LI. 0 .) PF(K) = 0. c set upper bound for snow i n forest to that of the open IF (PF (K) .GT.PO (K) ) PF(K) = PO{K) 200 CONTINUE c c 'max* i s set as the highest zone with r a i n MAX = MIN - 1 c i f a l l zones are snowing, return IF (MIN. EQ. 1) RETURN c c find p r e c i p i t a t i o n i n forest for a l l r a i n zones 260 K = 0 300 K = K + 1 PF (K) = PO (K) *FRAIN c set to zero i f negative r a i n generated IF (PF (K) .LT.0.0) PF (K) = 0. c set upper bound for rain i n forest to that of the open IF (PF (K) . GT.PO (K) ) PF(K) = PO (K) c check zone counter IF (K. LT.MAX) GO TO 300 c RETURN END SUBROUTINE ZSNOW (N, P, MAXR, PACK, MINPK, MAXRUN, RUNOFF, 1 T, NP, INCPL, PP) c c subprogramme to update snowpack condition, compute runoff c and to set symbols for plot t i n g snow d i s t r i b u t i o n c c inputs c no. of zones 'n', zonal p r e c i p i t a t i o n 'p*, c maximum rain zone 'maxr*: maxr<0 when i t i s a dry day, c maxr=0 when a l l zones are snowing, and maxrX) when i t i s c raining up to zone 'maxr' c zonal pack condition 'pack', snow le v e l 'minpk', minpk>n c i f basin i s snow-free c subroutine updates 'pack* and 'minpk* c t ( i , j ) i s temp, for zone j , i=1 for max. c 'np'=0 f o r pack graphic output, with ' i n c p l ' c o l . per zone c subroutine returns pack symbols 'pp*, zonal generated c runoff 5 the highest zone that contributes runoff 'maxrun' c DIMENSION SYMBO(20), PP(100), T{2,100), P(100), PACK (100), 1 RUNOFF (100) A 2 0 c c symbols for snowpack plot DATA BLANK / 1H / DATA SYMBO / 1 HO, 1H1, 1H2, 1H3, 1H4, 1H5, 1H6, 1H7, 1H8, 1 1H9, 1 HA, 1HB, 1HC, 1HD, 1 HE, 1 HF, 1 HG, 1HH, 1HI, 1H J / c slopes and origi n s for clear-weather 5 rain melt equations DATA SLOPEC, SLOPER, ORIGC, ORIGR /. 0724 , . 278 ,-0 . 73 ,- 1. 23/ c c check highest r a i n f a l l l e v e l IF (MAXR - 0) 100, 200, 300 c c maxr negative indicates a dry day c the highest zone that generates runoff i s i n i t i a l i z e d to 0 100 MAXRUN = 0 c snowline above basin (no snow l e f t ) IF (MINPK.GT.N) RETURN c c zones below 'tninpk' w i l l contribute no runoff MIN = MINPK - 1 IF (MIN.LE.0) GO TO 110 DO 105 I = 1, MIN 105 RUNOFF (I) = 0. c c some snow i n basin, evaluate clear weather melt 110 DO 180 I = MINPK, N TM AX = T (1,1) c check i f max. temp. > 0 IF (TMAX.LE.0.0) GO TO 150 c compute melt SM = ORIGC + SLOPEC * TMAX *TMAX IF (SM.LE. 0.0) GO TO 150 IF (SM. GT. PACK (I) ) GO TO 130 c update pack and set runoff to snowmelt PACK(I) = PACK (I) - SM RUNOFF (I) = SM GO TO 14 0 c more potential melt than i s available 130 RUNOFF (I) = PACK (I) PACK (I) = 0. c update highest zone that generates runoff 140 MAXRUN = I GO TO 180 c i f no melt, set runoff to zero 150 RUNOFF (I) =0. 180 CONTINUE GO TO 500 c c case when a l l zones are snowing (maxr=0) c highest zone that generates runoff i s set to maxr 200 MAXRUN = MAXR c lower snowline to zone (maxr+1) MINPK = MAXR + 1 c add new snow to pack A 21 DO 220 I = MINPK, '8 220 PACK (I) = PACK (I) + P ( I ) c c h e c k i f any m e l t g o e s on GO TO 110 c c c a s e when some z o n e s a r e s n o w i n g , some a r e r a i n i n g c c h e c k p o s i t i o n o f r a i n b e l t w i t h r e s p e c t t o e x i s t i n g p a c k 300 I F (MINPK,LE.MAXR) GO TO 350 c r a i n i n g b e l o w e x i s t i n g s n o w l i n e , r u n o f f g e n e r a t e d DO 310 I = 1, MAXR 3 10 RUNOFF (I) = P (I) c add new snow t o z o n e s a b o v e maxr l e v e l , e t c I F (MAXB.LT.N) GO TO 200 c i f a l l z o n e s a r e r a i n i n g MAXRUN = MAXR RETURN c c r a i n i n g a bove e x i s t i n g s n o w l i n e , so s e t h i g h e s t z o n e w i t h c g e n e r a t e d r u n o f f a s 350 MAXRUN = MINPK - 1 I F (MAXRUN.EQ.O) GO TO 400 c p u t r a i n f a l l i n t o r u n o f f up t o zone (minpk-1) DO 360 1 = 1 , MAX RUN 360 RUNOFF(I) = P ( I ) c e v a l u a t e r a i n m e l t f r o m z o n e minpk to maxr 400 DO 480 I = MINPK, MAXR RAIN = P (I) c compute r a i n m e l t a s f u n c t i o n o f r a i n * m a x . t e r n p SM = ORIGR + SLOPER * R A I N * T ( 1 , I ) c c h e c k i f computed m e l t i s p o s i t i v e IF (SM.GT.O.O) GO TO 410 c i f n o t , add r a i n t o p a c k PACK (I ) = PACK (I) + RAIN RUNOFF (I) = 0. GO TO 4 80 c c c h e c k i f computed m e l t e x c e e d s r a i n f a l l 410 I F (SM - RAIN) 420, 430, 440 c r a i n f a l l e x c e e d s m e l t , some e x c e s s i v e r a i n g o e s i n t o pack 420 EXCESS = RAIN - SM PACK (I) = PACK (I) + EXCESS c m e l t w a t e r g o e s t o r u n o f f 430 RUNOFF (I) = SM GO TO 4 60 c more m e l t w a t e r t h a n r a i n f a l l 440 AVAIL = PACK (I) + RAIN c c h e c k i f a v a i l a b l e w a t e r s u p p l y meets m e l t demand IF (SM. GE. AVAIL) GO TO 450 c pack i s r e d u c e d , and r u n o f f e q u a l s m a l t PACK (I) = PACK (I) - (SM-RAIN) GO TO 430 c i n s u f f i c i e n t w a t e r a v a i l a b l e , p a c k d e p l e t e d 450 P A C K ( I ) = 0. A22 RUNOFF (I) = AVAIL c reset highest zone with generated runoff 460 MAXR UN = I 480 CONTINUE c c check snowline 500 DO 550 MINPK = 1, N IF (PACK (MINPK) .GT.0.0) GO TO 600 550 CONTINUE c i f none of the zones has snow MINPK = N + 1 RETURN c c check for p l o t t i n g 600 IF (NP.GT.0) RETURN c i n i t i a l i z e p l o t t i n g position ID = 0 c i f a l l zones have snow IF (MINPK-LE.1) GO TO 630 c set a l l snow-free zones to blank M = MINPK - 1 DO 620 I = 1, M IL = IU + 1 IU = I * INCPL DO 620 K = IL, IU 620 PP(K) = BLANK c c set symbols 630 DO 650 I = MINPK, H IL = IU + 1 IU = I * INCPL c set symbols i n the range of 1 to 2000 mm (100 mm interval) J = PACK (I) *0.01 + 1.0 IF (J.GT.20) J = 20 DO 650 K = IL, IU 650 PP(K) = SYMBO(J) c RETURN END FUNCTION RANUM (IDUMMY) C c function subprogramme to generate uniformly di s t r i b u t e d c random numbers by mixed multip l i c a t i v e congruential method c r e f . : carnahan et a l i a , 'applied numerical methods', p 549 c DATA M, N, J / 1048576, 566387, 1027 / DATA XM / 1048576. 0 / c N = MOD (J*N, M) RAN UM = FLOAT (N) / XM RETURN END 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items