g ^ J- i v ARO PRE-PROPOSAL FOR TERRESTRIAL SCIENCES RESEARCH AREA 3 - ENVIRONMENTAL SCIENCES MONITORING AND FORECASTING ATTRIBUTES OF LAND SURFACE HYDROLOGY THAT ARE RELEVANT FOR CHARACTERIZING A POTENTIAL BATTLEFIELD A. W. (Tony) England* and Roger DeRoo** *Professor of Atmospheric, Oceanic, and Space Sciences *Professor of Electrical Engineering and Computer Science **Assistant Research Scientist, Space Physics Research Laboratory University of Michigan Ph (734) 763-5534; Fx (734) 647-2106; e-mail england@umich.edu December 10, 2002 Goals and Objectives The long-term goal of this project is to develop a method to monitor and forecast the attributes of land surface hydrology that are important for characterizing a potential battlefield. Relevant attributes, for purposes of this proposal, are vertical profiles in snow and soil of volume percent liquid water and of temperature. These attributes will be derived from Soil-VegetationAtmosphere Transfer (SVAT) models forced either by observed weather data or by data from Numerical Weather Prediction (NWP). The SVAT models will be embedded in simple regional land surface hydrology models and supported by assimilation of near-daily UAV-based microwave brightness measurements of near-surface moisture or by a combination of less frequent UAV-based measurements and satellite measurements of near-surface soil moisture. Our goal is a cumulative rms accuracy for a resolution cell of -5 volume percent liquid water and -2 Kelvin, respectively. Accuracies of < 3 volume percent liquid water and < 1 Kelvin, respectively, to depths of 60 cm have been demonstrated at plot scale for grassland and agricultural terrains where weather forcing, soil texture, and vegetation canopy column density are known. Similar accuracies should be possible for a snowpack because emission and absorption processes in moist soil and wet snow are similar. The resolution cell will be -30 meters for our prototype system but could be -10 meters for an operational system. The project is divided into three phases of which this proposal addresses Phase 1. There are two parallel objectives within this phase. Because either or both may be funded, their costs are presented separately but funding both will move the project most rapidly toward the long-term goal. These objectives are: 1) To use our vegetation-over-soil SVAT/Radiobrightness model, develop a snow-over-soil SVAT/Radiobrightness model, and assimilate 1.4 GHz and/or 6.9 GHz brightness to estimate plot scale volumetric liquid water and temperature profiles in soil and snow-over-soil with rms accuracies of 5 volume percent and 2 Kelvin, respectively, and 2) To complete development of the 1.4 GHz STAR-Light imaging radiometer for use on a light aircraft to calibrate regional hydrology/SVAT/Radiobrightness models, demonstrate attainment of the project's long-term goal, and serve as a prototype for a UAV-based sensor. Background Data from the new generation of satellite radiometers will allow reliable monitoring of daily changes in the freshwater hydrologic cycle as expressed in temporal and spatial patterns of moisture stored in soil and snow. The capability will be limited in this decade to spatial scales >10 km by the low spatial resolution of satellite sensors. This resolution limitation can be 1 RL-1010 = RL-1010

reduced to -10 meters for regions of interest with near-daily over-flight of, for example, an Unpiloted Airborne Vehicle (UAV) carrying a compact 1.4 GHz imaging radiometer. The concept is summarizing in Figure 1. Effectively, this is an expansion of Numerical Weather Prediction (NWP) to include detailed forecasts of land surface hydrology. Products will be snow and soil liquid moisture and temperature profiles to depths of a meter or more. We have demonstrated the concept at plot scale for prairie grassland (Galantowicz, 1995; Galantowicz and England, 1997; Judge, 1999; Lou et al, 1999; Judge et al, 1999; Judge et al, 2001), but recovering snowpack wetness remains a research objective. Atmospheric Model Airborne 1.4 GHz Imaging Radiometer Weather & downwellingradiance Tb(observed SVAT Model | - - Assimilate Tb(observed) - Tb(model) Temperature & Moisture Profiles T(model) Radiobrightness Model Figure 1. Strategy for remotely monitoring water stored in soil and snow. The Atmospheric Model (or observed weather) provides atmospheric lower boundary weather and downwelling short- and longwavelength radiation as forcing to the Soil-Vegetation-Atmosphere Transfer (SVAT) Model. The SVAT Model provides soil, vegetation, and snow temperature and liquid moisture profiles that are used by the Radiobrightness Model to predict emitted microwave brightness, Tb. At low microwave frequencies, e.g. 1.4 GHz, the difference between observed brightness from an airborne sensor and predicted brightness represents an error in the liquid water content of near-surface soils, vegetation, or snow. This error is assimilated by the SVAT Model to improve the current estimate of near-surface moisture. Over several assimilation cycles, a calibrated SVAT model converges to an accurate moisture profile at depths below those sensed by the microwave radiometer. The project is outlined in Figure 2. Phase 1, this proposal, includes (a) use of vegetationover-soil and development of snow-over-soil SVAT/Radiobrightness models, and (b) fabrication and test of 1.4 GHz STAR-Light imaging radiometer. Phase 1 - 2 yr Develop Snow-over-Soil SVAT/R Model Fabricate and Test STAR-Light Phase 2 - 1 yr Embed SVAT/R in Regional Hydrology Model Mount STAR-Light on Light Aircraft.at - 3 - 2 y Phase 3 - 2 yr Validate Concept: Recover Moisture and Temperature Profiles Using Reaional Model and STAR-Lieht Figure 2. Project Structure. Phase 1 uses vegetation-over-soil SVAT/R model and develops snow-oversoil SVAT/R model, and fabricates and tests 1.4 GHz STAR-Light imaging radiometer. Phase 2 embeds SVAT/R models in simple regional land surface hydrology model, and mounts and tests STAR-Light on light aircraft. Phase 3 integrates model and airborne sensor to validate our concept for recovering moisture and temperature profiles in vegetation-over-soil and snow-over-soil at a spatial resolution of -10 meters. 2

(a) Vegetation-over-Soil and Snow-over-Soil SVAT/Radiobrightness Models Calibrated SVAT/Radiobrightness models have been a focus of our group's research for more than a decade. Our prairie grassland Land Surface Process/Radiobrightness (LSP/R) model exhibits excellent performance as illustrated in the White Paper submitted to ARO prior to this proposal. When forced by downwelling long- and short-wavelength radiation and observed weather, it maintains reliable temperature and moisture profiles of the upper meter of prairie soil, and accurately predicts SSM/I radiobrightness (e.g., Liou et al, 1999; Judge et al, 1999). With imperfect weather data, modeled profiles and predicted brightness will deviate from observations. Assimilation of differences between predicted and observed radiobrightness at near-daily intervals will allow modeled profiles to converge to observed profiles (e.g., Entekhabi et al, 1994; Galantowicz et al, 1999). SVAT/Radiobrightness models emerged as powerful tools for prairie and agricultural hydrology when experimental data from a series of field campaigns enabled tests of schemes for modeling energy and moisture transport in prairie soils. Notable among these field campaigns have been the First ISLSCP Field Experiment (FIFE) in 1987 (e.g., Betts et al, 1996), Washita'92 (Jackson et al, 1995), and the Southern Great Plains Hydrology Experiments - SGP'97 and SGP'99. We participated in SGP'97 with REBEX-5, our 5th Radiobrightness Energy Balance Experiment (Judge et al, 1999). Particularly relevant for this proposal, we have extended the approach to terrains beyond those of southern prairie. REBEX-1 was a 7 month experiment in grassland near Sioux Falls, South Dakota, during fall and winter of 1992-1993 (Galantowicz, 1995; Galantowicz and England, 1997). REBEX-3 was a 12 month experiment in tussock tundra near Toolik Lake, Alaska, during 1994-1995 (Kim and England, 1998; Kim, 1999, Kim and England, 2002). REBEX-4 was a summer study in northern prairie grassland near Sioux Falls during 1996 (Judge et al, 1999; Judge et al, 2001). REBEX-7 and -8 were summer studies in corn near Ann Arbor, Michigan (Hornbuckle and England, 2002 a,b). The Army's Cold Regions Research and Engineering Laboratory (CRREL) recently funded us for the Colorado Cold Lands Processes Experiment (CLPX) in winter and spring of 2003. For this project, we will use our current LSP/R model for vegetation-over-soil and develop a. snow-over-soil SVAT/R model for 1.4 and 6.9 GHz brightness. The new snow-over-soil model will include the soil model from LSP/R and a version of the CRREL SNTHERM snow model. We used SNTHERM.89 in an early version of our LSP/R model (Galantowicz, 1995) and have received a recent version of SNTHERM through the gracious assistance of Rachel Jordan at CRREL. To take advantage of our soil model's skill with both heat and moisture transport in unsaturated soil, we will modify SNTHERM to allow exchange of moisture between soil and snowpack. The vegetation-over-soil model will be calibrated with data from REBEX-7 and -8, and the snow-over-soil model will be calibrated with data from CLPX. We do not anticipate that this phase of the project will require additional field work. (b) Fabrication and Test of the 1.4 GHz STAR-Light Imaging Radiometer Only one airborne 1.4 GHz imaging radiometer is currently available for model calibration in North America - the 1-dimensional Synthetic Thinned Array Radiometer (STAR) technology instrument, ESTAR, which flies on the NASA P-3, a 4-engine turboprop aircraft. Three groups worldwide are developing airborne instruments based upon the 2-dimensional STAR technology. While none has yet actually recovering an airborne brightness image, ESA has sufficient confidence in 2-d STAR technology to base development of their SMOS satellite instrument 3

upon it. The first airborne instrument, by the Helsinki University of Technology (HUT) for ESA, is to be flown on a 2-engine turboprop aircraft. The second, by the Goddard Space Flight Center (GSFC) for NASA, is to be flown on their P-3 aircraft. The third, STAR-Light by our group with initial funding by the NSF Office of Polar Programs, is to be flown on a Super Cub class aircraft (Figure 3) - an operationally far less expensive alternative and the only alternative that offers a development path to the limited weight and power capabilities of most UAV. Figure 3. Aviat Husky Al-B. As the current equivalent to the venerable Piper Super Cub, this aircraft serves as our design baseline for STAR-Light. STAR-Light is a 1.4 GHz STAR technology imaging radiometer for use on light aircraft. It offers an alternative to the cost, complexity, and limited availability of the other instruments. It is simpler than either the ESA HUT or the NASA GSFC instruments, but is technically the most advanced. Its three innovative technologies are 1) 2-dimensional STAR architecture, 2) Direct Sampling Digital Radiometer receivers (a technology pioneered by our group), and 3) co-aligned band-definition in the digital domain (also a technology pioneered by our group). These technologies permit a compactness not realized with other 1.4 GHz imaging radiometers. The 10 STAR-Light antenna/receiver elements, shown in Figure 4, follow the "Y" pattern of the Very Large Aperture (VLA) radio telescope at Socorro, NM, and the ESA SMOS radiometer to offer maximum resolution with a minimum number of elements. STAR-Light is built in two modules: The Sensor Module (Figure 4) mounted under the belly of the aircraft, and the autonomous Control and Data Module replacing the aft seat in the two-place aircraft. Figure 4. Sensor Module. Arrow indicates flight direction. The tail of the arrow is in the slot for a glycol coolant radiator. The face of the instrument will be covered by a radome and, except for the individual antenna apertures (denoted by the 10 rings), by thermal insulation We were awarded $1.38M to develop STAR-Light by the NSF Office of Polar Programs in January, 2001. With our guidance, engineers of the Space Physics Research Laboratory (SPRL) here in the College of Engineering are developing STAR-Light. SPRL has an impressive history of developing over 100 instruments that have operated successfully on balloons, rockets, and spacecraft. Two are currently in Earth orbit and one is on the way to Saturn. The STAR-Light Preliminary Design Review (PDR), held on 3 October 2001, was attended by SPRL engineers, electrical engineering faculty, and Niels Skou, an intellectual leader of ESA's development of 4

the SMOS instrument - also an aperture synthesis radiometer. The Critical Design Review (CDR), held on 21 March 2002, was attended by SPRL engineers, electrical engineering faculty, and 4 external referees. The external referees, Alan Tanner and Mark Fischman from the NASA Jet Propulsion Laboratory, and Dave LeVine and Ed Kim from the NASA Goddard Space Flight Center, are internationally recognized radiometer engineers. All reviews endorsed the STARLight design. Except for electrical power and ground service systems, all interface specifications, schematics, parts lists, and layout CAD drawings for STAR-Light are complete. All electronic components for the STAR-Light receivers were purchased, bench tested, and thermally cycled during summer of 2002. Design of the antenna was contracted to a local electromagnetics house. As validation of the receiver design, we have used versions of the STAR-Light receiver in the 1.4 and 6.7 GHz dual-polarization radiometers of our Truck Mounted Radiometer System (TRMS3) now deployed as part of CLPX in Colorado. Design instrument performance is summarized in the following table. Angular resolution is determined in post-processing as a tradeoff between spatial resolution and radiometric noise. STAR-Light Design Performance 1 Frequency / Bandwidth 1.4135 GHz / 21 MHz NEAT / Accuracy 0.5 K / 2 K Across & Along Track Scans +/- 350 Angular Resolution 160 < 11/2 < 240 Weight 197 lbs Power 200 W STAR-Light was conceived as a "skunk works" style development much like the several 1.4, 6.7, 19, 37, and 85 GHz radiometers our group has built. Once development began, STAR's inherent complexity and our need for compactness forced us to a more systematic strategy like that used by the aerospace industry to develop satellite instruments. This "aerospace" development style greatly reduces development risk but at nearly twice the original estimate of cost. Difficulties experienced by other groups developing 2-dimensional STAR instruments support our decision to use an "aerospace" style development. SPRL performed a bottom-up cost-to-completion assessment using a multi-level Work Breakdown Structure involving each project engineer. Their no-contingency estimate, including indirect charges, is $1.258 M over 1.25 years. A 25% margin is realistic for a new development that has passed CDR. With this contingency, our estimate of cost to complete STAR-Light is $1.579 M over 2 years. Phase ] Costs Period of performance 2 Years May 1, 2003 through April 30, 2005 Cost for SVAT/R model development $ 372,309 Cost for STAR-Light radiometer development $ 1,579463 TOTAL $ 1,951,772 5

References Cited England, A.W., 1990. Radiobrightness of diurnally heated, freezing soil, IEEE Trans.Geosci. Remote Sensing, 28, pp.464-476. England, A. and J.F. Galantowicz, 1995. Observed and modeled radiobrightness of prairie grass in early fall," Proc. of IGARSS'95, IEEE GRSS, Florence, Italy, July 10-12. Entekhabi, D., H. Nakamura, and E. G. Njoku, 1994. Solving the inverse problem for soil moisture and temperature profiles by sequential assimilation of multifrequency remote sensing observations, IEEE Trans. of Geoscience and Remote Sens., 32, pp. 438-448. Fischman, M., and A.W. England, 1998. A direct-sampling receiver for synthetic thinned array radiometry, Proc. of IGARSS'98, 3, IEEE GRSS, pp. 1711-1713, July. Fischman, M., and A.W. England, 1999. Sensitivity of a 1.4 GHz direct-sampling digital radiometer, IEEE Trans. Geosci. Remote Sensing, 37, pp. 2172-2180. Fischman, M.A., A.W. England, and C.S. Ruf, 2002, How digital correlation affects the fringe washing function in L-Band aperture synthesis radiometry, IEEE Trans. Geosci. Remote Sensing, 40, pp. 671-679. Galantowicz, J.F., 1995. Microwave Radiometry of Snow-covered Grasslands for the Estimation of Land-atmosphere Energy and Moisture Fluxes, Ph.D. dissertation, University of Michigan, 221 pages. Galantowicz, J. F., and A.W. England, 1997. Seasonal snowpack radiobrightness interpretation using a SVAT-linked emission model, J. of Geophysical Research, 102(D18), pp. 21933 -21946. Galantowicz, J. F., D. Entekhabi, and E. G. Njoku, 1999. Tests of Sequential Data Assimilation for Retrieving Profile Soil Moisture and Temperature from Observed L-Band Radiobrightness, IEEE Trans. on Geoscience and Remote Sens., vol. 37, pp. 1860-1870. Hornbuckle, B., D.L. Boprie, and A.W. England, 2000. Micrometeorology theory and instrumentation for land surface process modeling, Radiation Laboratory Tech. Note, RL-968. Judge, J., 1999. Land Surface Process and Radiobrightness Modeling of the Great Plains. Ph.D. Dissertation, University of Michigan, 213 pages. Judge, J., and A.W. England, 1997. LSP/R models for bare soil and grassland in the northern Great Plains, Proc. of IGARSS'97, Singapore, Malaysia, August 3-8. Judge, J., A. England, C.L.W. Crosson, B. Hornbuckle, D. Boprie, E. Kim, and Y. Liou, 1999. A growing season Land-Surface Process/Radiobrightness model for wheat-stubble in the Southern Great Plains. IEEE Trans. Geosci. Remote Sensing, 37, 2152-2158. Judge, J., J.F. Galantowicz, and A.W. England, 2001, A comparison of ground-based and satellite-borne microwave radiometric observations in the Great Plains, IEEE Trans. Geosci. Remote Sensing. 39, pp. 1686-1696. Kim, E.J., 1999. Remote Sensing of Land Surface Conditions in Arctic Tundra Regions for Climatological Applications Using Microwave Radiometry, Ph.D. Dissertation, University of Michigan, 172 pages. Kim, E.J., and A.W. England, 1998. Land surface process modeling and passive microwave remote sensing of arctic tundra regions, Proc of IGARSS'98, Seattle, WA, July 6-10. Kim, E.J., and A.W. England, 2001, Spatial Scaling of satellite 19 and 37 GHz brightness observations of the Alaskan North Slope, submitted to J. Geophys. Res., March, 2001. Liou, Y.-A., J. F. Galantowicz, and A. W. England, 1999. A Land Surface Process/ Radiobrightness Model with Coupled Heat and Moisture Transport for Prairie Grassland, IEEE Trans. Geosci. Remote Sensing, 37, pp. 1848-1859. Liou, Y.A., E.J. Kim, and A.W. England, 1998. Radiobrightness of prairie soil and grassland during dry-down simulations, Radio Science, 33, pp. 259-265. 6

Date: Fri, 6 Dec 2002 17:50:30 -0500 From: Virginia To <vto@hpti.com> To: "Rollett, Tony" <rollett@andrew.cmu.edu>, "Huskamp, Jeffrey C." <HuskampJ@MAIL.ECU.EDU>, Kamal Sarabandi <saraband@umich.edu>, katehi@purdue.edu, "Nystrom, Nicholas A." <nystrom@psc.edu> Cc: Pet-it <pet-lt@hpti.com>, cgouin@ARL.ARMY.mil, bclarke@ARL.ARMY.mil, "Jennifer Moses (HPTilcrimmins)" <jmoses@ARL.ARMY.mil>, Charles Macon-Contractor <cmacon@ARL.ARMY.mil>, ccornwell@ARL.ARMY.mil, Mark Zottola - Contractor <mzottola@ARL.ARMY.mil>, drichie@ARL.ARMY.mil, "Sotirelis, Paul" <sotireli@asc.hpc.mil>, "Blaudeau, Jean" <blaudejp@asc.hpc.mil> Subject: February PET Review and FY04 White Papers [ Part 1.1, Text/PLAIN (charset: ISO-8859-1 "Latin 1 (Western ] [ Europe)") 235 lines. ] [ Unable to print this part. ] [ The following text is in the "iso-8859-1" character set. ] [ Your display is set for the "US-ASCII" character set. ] [ Some characters may be displayed incorrectly. ] Dear FAPOCs, Here is some additional guidance and information regarding the upcoming the FY04 White Paper/ Project Proposal submissions and the February PET Technical Review. I realize there is alot of information here, but much of it you have already seen in draft form. Please feel free to contact George, Ralph or me if you have any questions. We will be taking to you soon. ginny <?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" /> White Papers Attached are the government^Os Selection and Evaluation Plan for PET FY2004 Projects and HPTi^Os Selection/Recommendation Plan. Please note that the White Papers are due to the government on February 14, 2003 ^0 so there is a small window to incorporate new ideas received from the February Technical Review. Project proposals that cut across CTA areas, including areas like CE, ET and CDLT are encouraged. The government has agreed that in this case, the budget can be split across the CTA areas. The dominant area will have project responsibility for the project as a whole. Additionally, we highly encourage collaboration with and using CDLT

tools for seminars, meetings, training should also be encouraged, core or otherwise. Please remember that the next set of proposals is for 8 months only in the first year (I presume there is no objection to submitting multi-year proposals this year). This will allow us to synchronize the project cycle with the contract cycle. Finally, think of ways to involve the MSIs and interact with Pat on identifying MSIs who can materially contribute to their projects. We will use the same HPTi PET web site process for white paper and proposal submission. February Technical Review The Technical Review will be held at Clark Atlanta University. To maximize cross functional area interaction, we encourage your attendance for the whole agenda. All PET project PIs and on-sites are invited to attend, but not required. For the review, plan on spending roughly half your time on achievements to date, and the rest of the time laying out what you think is important to do next. Of course, given the timing, you will already know what White Papers have been submitted, so this review is a chance to promote those ideas and see what kind of reception you get. Do not forget that you will need to get information and represent the work of other PIs in your functional area. Attached please find the agenda, instructions and review presentation template. Also, please note the presentations. schedule of events below regarding the PET Technical Review: Feb. 3-7, 2003

Activity Start Date End Date Responsibility Provide Review & JTC instructions to FAPOCs 6-Dec 6-Dec To Invitations to on-sites and Project PIs 11-Dec 11-Dec To Provide presentations to HPTi (C Krahulec, C Gouin) 13-Jan 13-Jan FAPOCs First review - initial tech/formatting edit 13-Jan 17-Jan C Krahulec, C Gouin Second review - technical content 15-Jan 20-Jan Roskies, Lea, To Provide Feedback to FAPOCs 20-Jan 21-Jan Roskies, Lea, To Finalize presentations

22 —Jan 30 —Jan FAPOCS Burn CDs and hardcopies if necessary 31 —Jan 1-Feb C Gouin Arrive for review 2-Feb 3-Feb All Attend technical review 3-Feb 6-Feb All Joint Technology Council Meeting The main point of the JTC meeting is to explore cross-CTA ideas that may not have been mature enough to discuss in the open review. Hopefully, as we all will have sat through 3 days of briefings by our colleagues, some new ideas for cross-CTA collaboration will have arisen. Each FAPOC is encouraged to bring along a speculative cross-CTA idea and bring it up for discussion. Also, this will be our first opportunity to get together without having the government, or even CPOCs, present. We should use the opportunity to have some constructive discussion about how to improve the program from the FAPOC perspective. To encourage frank discussion, no material will be included in any report that will emerge from this meeting without the consent of all attending FAPOCs. If you have any agenda items for the JTC, please send them to Ralph Roskies (roskies@psu. edu)

The Joint Technology Council Meeting will be held at a conference room at the Executive Conference Center at the Hartsfield Airport (see directions below) and will commence at 1:30pm. Buffet lunch will be served at the conference center starting 12:30pm. Please do not book departing flights before February 6, 5:30pm. February Logistics at Atlanta We will book a block of rooms at a hotel in Atlanta for the HPTi team. Instructions for room reservations will be forthcoming. Directions: To the Executive Conference Center, (Airport Atrium) from Ground Level. After your arrival at Hartsfield Airport, come to the second (2nd) level. Pass through Ticketing or Baggage Claim and come to the area where the Car Rental agencies are located, (between North and South Terminals). Notice signs above for the ^OAtrium,^ (an area similar to a shopping mall). Inside the Atrium, catch the elevator to the 3rd floor. The Executive Conference Center is located beside Houlihan Os Restaurant. [ Part 2, Application/MSWORD 44KB. ] [ Unable to print this part. ] [ Part 3, Application/VND.MS-EXCEL 21KB. ] [ Unable to print this part. ] [ Part 4, Application/MSWORD 41KB. ] [ Unable to print this part. ] [ Part 5, Application/VND.MS-POWERPOINT 1MB. ] [ Unable to print this part. ] [ Part 6, Application/MSWORD 728KB. ] [ Unable to print this part. ] [ Part 7, Application/MSWORD 161KB. ] [ Unable to print this part. ]

THE UNIVERSITY OF READING Modelling the Effect of Vegetation on the Seasonal Snow Cover Submitted for the degree of Doctor of Philosophy by Melody Jane Tribbeck Environmental Systems Science Centre April 2002

Abstract A new Snow-SVAT (Snow-soil-vegetation-atmosphere transfer scheme), 'SNOWCAN', has been formulated by coupling a three-parameter optical and thermal canopy radiation model to a physically based snow energy budget model. The coupling approach has allowed feedback mechanisms to be simulated in a physically reasonable manner. Multiple scattering of radiation between the canopy and snow surface is simulated by SNOWCAN. Forest processes, including the interception of snow and reduction of wind speed by the canopy, are represented within SNOWCAN. Other processes, such as the degradation of surface albedo for a thin snowpack, are considered to be negligible for this study. Sensitivity studies have investigated the relative importance of model parameters, and uncertainties in the forcing data. SNOWCAN was tested with data from the BOREAS campaign, and simulated snow depth showed good agreement with observations. A new experiment was designed to collect data, to validate SNOWCAN, and was carried out within the Reynolds Creek Experimental Watershed, Idaho, during water year 2000-2001. Validation of SNOWCAN against sub-canopy radiation, snow depth and snowpack profile measurements has demonstrated that the model simulates observations well, generally to within the measurement error. However, following substantial melt, the simulated snowpack retained liquid water not observed in the field. Comparisons with observations have suggested that the rate of snow grain growth is underestimated by SNOWCAN. An energy balance comparison has shown that snowmelt at the Reynolds Mountain East fir site in year 2001 was dominated by longwave radiation, in contrast to the open site, which was dominated by sensible heat transfer. SNOWCAN has been used to investigate canopy spatial variability effects, and the diurnal and seasonal change in above-canopy albedo.

Acknowledgments This project was funded by the Natural Environment Research Council grant number GT04/98/TS/247, in conjunction with a CASE award from the British Antarctic Survey. Thank you to my two supervisors, Prof. Robert Gurney and Prof. Liz Morris. Their enthusiasm, knowledge and wisdom have been so inspiring and instrumental to the fantastic training that I've had during the past 3 years. Thank you to Dr. Danny Marks and all staff at the Northwest Watershed Research Center for all your hard work in the collaborative fieldwork. Thanks to Danny for processing of reflectance data and general discussions of snow model development. Thank you to Adam Winstral and Tim Link, also for snow model discussions and field data provision and proccessing. Thanks to Adam Winstral, Mark Morehead, Mark Murdoch and Tim Link for giving up valuable time in the field and to Jim Windom for your technical expertise with field equipment. Thank you to all staff and students at ESSC, particularly Ian Davenport, who is helpful and resourceful, helped in the field and in taking laboratory reference measurements. Thank you to Dave Pearson, who made the RM code available, and for useful SVAT discussions. Thanks to Dave Cobby, for picking up all the items that never reached the bin and for listening. Thank you to other snow researchers, particularly to John Pomeroy, Janet Hardy and Rachel Jordan for discussions. Thanks to Jim Foster and survivors of North Dakota '99 for fieldwork experience. i

Declaration I confirm that this is my own work and the use of all other material from other sources has been properly and fully acknowledged. Melody Tribbeck ii

Contents Acknowledgments i Declaration ii List of Figures viii List of Tables xii Notation xiii Chapter 1 Introduction 1 1.1 Project Justification................................. 1 1.2 A W alk in the Woods.............................. 4 1.3 M ethod of Approach............................... 8 1.4 Thesis Description................................ 10 Chapter 2 Snow Energy Budget Model 12 2.1 Introduction.................................... 12 2.2 Model Description................................ 14 2.2.1 M odel Setup............................... 15 2.2.2 Boundary Conditions............................ 16 2.2.3 Snow Model Assumptions............................. 22 2.2.4 Physical Principles............................ 24 2.2.5 Numerical Solution Techniques......................... 30 2.2.6 Time Step.................................. 34 iii

2.2.7 Model Outputs.............................. 35 2.3 SNOWCAN and SNTHERM Differences..................... 36 2.3.1 Model Algorithms............................. 36 2.3.2 Numerical Solutions............................. 38 2.3.3 Change in Style................................... 39 Chapter 3 Snow Model Sensitivity 41 3.1 Introduction............................. 41 3.2 Synthetic Snowpack Description.................................. 43 3.3 Meteorological Data Input........................... 44 3.4 Model Sensitivity to Parameterisation..................... 46 3.4.1 Surface Energy Balance Parameters...................... 46 3.4.2 Internal Snow Physics Parameters..................... 52 3.4.3 Site-to-site Variation.......................... 60 3.4.4 Physical Constants............................ 63 3.4.5 Model Discretisation........................... 66 3.4.6 Model Accuracy Criteria......................... 68 3.5 Sensitivity to Boundary Conditions............................... 69 3.6 Sensitivity to Initial Conditions............................... 70 3.7 Conclusions.................................... 72 Chapter 4 Vegetation Effects 75 4.1 Introduction.................................... 75 4.2 Canopy Radiation Model............................ 76 4.2.1 Introduction.............................. 76 4.2.2 Physical Principles................................ 78 4.3 Other Canopy Effects.............................. 80 4.3.1 Intercepted Snow............................... 80 4.3.2 Sub-Canopy Wind Speed................................. 84 4.3.3 Forest Debris................................ 84 4.3.4 Thin Snowpack Albedo......................... 86 4.4 Canopy Effects Not Modelled............................... 87 iv

4.5 Model Implementation................ 4.6 Summary....................... Chapter 5 Snow-SVAT Simulation I: Boreal Forest 5.1 Introduction...................... 5.2 Site Description.................... 5.3 M odel Inputs..................... 5.3.1 Model Initialisation............. 5.3.2 Canopy Parameterisation.......... 5.3.3 Forcing Data................. 5.4 Simulation of Snow Depth.............. 5.5 Conclusions...................... Chapter 6 Fieldwork 6.1 Introduction...................... 6.2 Site D etails...................... 6.3 Parameterisation Data................ 6.3.1 Canopy Parameters............. 6.3.2 Snow Parameters............... 6.3.3 Atmospheric Parameters........... 6.4 Snowpack Initialisation................ 6.5 Climate (Forcing) Data............... 6.6 Validation data.................... 6.6.1 Sub-Canopy Radiation............ 6.6.2 Snow Profiles................. 6.6.3 Snow Depth.................. 6.7 Summary..................................... 88.............. 88 92. 92. 93. 93. 93. 95. 96. 97. 101 103. 103. 104. 108. 108. 111. 114. 117. 118. 125. 125. 135. 139. 140 142. 142. 143. 143 Chapter 7 Snow-SVAT Simulation 7.1 Introduction........... 7.2 Open Site Simulation...... 7.2.1 Model Inputs...... II: ~. Reynolds Creek........................... v

7.2.2 Snow Depth Validation................................... 144 7.2.3 Snowpack Profile Validation............................ 145 7.2.4 Model Sensitivity................................... 151 7.3 Parameter Space Comparison.........................153 7.4 Fir Site Simulation................................154 7.4.1 M odel Inputs.......................... 154 7.4.2 Snow Depth Validation...............................155 7.4.3 Sub-Canopy Radiation Validation...................156 7.4.4 Sub-Canopy Wind Speed.............................165 7.4.5 Snowpack Profile Validation......................168 7.5 Canopy Parameter Sensitivity.......................... 174 7.6 Conclusions...............................175 Chapter 8 Snow-SVAT Simulation III: Diagnostic Tool 178 8.1 Introduction.............................................178 8.2 Surface Energy Balance........................... 179 8.2.1 Inter-Site Comparison.................................. 179 8.2.2 Thermal Contributions From the Sky, Vegetation and Snow........ 187 8.3 Canopy Variability................................ 187 8.3.1 LAI-Dependent Ratio of Sub-Canopy Radiation Components.... 187 8.3.2 Variation of Snow Depth.......................... 190 8.3.3 Above-Canopy Albedo........................ 191 8.4 Reversal in Site Melt Order......................... 194 8.4.1 Site Mass Differences........................... 194 8.4.2 Other Causes..........................................196 8.5 Conclusions........................................... 197 Chapter 9 Conclusions 199 9.1 Introduction................................ 199 9.2 Project Accomplishments........................................200 9.3 Consideration of Propositions....................................202 9.4 Model Limitations............................ 206 vi

9.5 Model Improvements...............................207 Appendix A Code Inconsistencies 210 vii

List of Figures 1.1 Conifer forest in RCEW, Idaho............................. 4 1.2 Snow intercepted on a conifer canopy in RCEW, Idaho, USA........ 6 1.3 Photograph of tree well............................. 7 1.4 Heterogeneity of surface and sub-surface organic matter............ 8 2.1 Simulated snowpack adaptive grid....................... 15 2.2 Comparison of SNTHERM and SNOWCAN compaction algorithm..... 37 3.1 Initialisation profiles for sensitivity analysis.................. 43 3.2 Simulated depletion of snow cover with standard model parameterisation. 45 3.3 Structure of the lower layer of the atmosphere for steady horizontal flow over a uniformly rough snow surface (from Morris (1989)).......... 50 3.4 Model liquid water fraction temperature dependence (from Jordan (1991)). 56 4.1 RM model canopy fluxes............................. 79 4.2 Snow deposition patterns beneath the trees of the taiga........... 81 4.3 SNOWCAN model structure.......................... 89 5.1 Maps showing location of BOREAS old jack pine site within the Southern study area..................................... 94 5.2 Initialisation profiles for BOREAS simulations................ 96 5.3 Simulated snow depth beneath a jack pine boreal forest canopy....... 98 5.4 Simulated snow depth beneath a jack pine boreal forest canopy [reconstructed precipitation data]........................... 99 5.5 Simulation of boreal snow depth with canopy effects removed........ 101 viii

6.1 Map showing location of Reynolds Creek Experimental Watershed within Idaho State.................................... 105 6.2 Map showing elevation of RME Study area......................106 6.3 Image of Reynolds Mountain East Watershed.....................107 6.4 Image of fir needles and branchlets used to determine canopy single-scatter albedo....................................... 110 6.5 Laboratory reflectance spectra of fir needles and branchlets............111 6.6 Photographs of snow surface during April field campaign.............113 6.7 Snow surface reflectance variation with snow depth.................114 6.8 Tracking and measurement of the sun using a Sun Photometer..........116 6.9 Ridge snow profile initialisation data...................... 117 6.10 Fir snow profile initialisation data....................... 118 6.11 2001 Meteorological data for open (ridge) site.....................120 6.12 2001 Meteorological data for RMSP..........................121 6.13 2001 Meteorological data for fir climate station...................122 6.14 Accumulation on meteorological station at snow pillow site............123 6.15 Shadowband radiometer installed at the ridge site.................124 6.16 Diffuse radiation measured by a shadowband radiometer.............125 6.17 Proportion of direct and diffuse radiation determined from a radiometer array under a stand of aspens.............................127 6.18 Upward and downward looking pyranometers and pyrgeometers at the fir sitel28 6.19 Up- and downwelling thermal and solar sub-canopy radiation at the RME fir site....................................... 129 6.20 Illustration of reflectance measurement error due to cloud cover.........130 6.21 Operator holding reference panel under spectroradiometer for reflectance measurement......................................133 6.22 Decay of snow reflectance at a fir site..........................135 6.23 Validation data from ridge snow pits..........................137 6.24 Validation data from fir snow pit: 16th March 2001..................138 6.25 Insertion of temperature array into snow pit wall..................138 6.26 Temperature profile of snowpack under fir canopy..................139 ix

6.27 Snow depth depletion at RME ridge, snow pillow and fir sites......... 140 7.1 Comparison between observed and simulated snow depth at RME ridge site 144 7.2 Comparison of snow temperature, density and grain size profiles at the RME ridge site...................................... 146 7.3 Ridge snow temperature and grain size profiles: larger initial grains....... 148 7.4 Comparison of simulated partial ice density with measured snow density.. 151 7.5 Simulated and observed snow depth beneath a fir canopy in RME....... 155 7.6 Sub-canopy downwelling solar radiation........................ 157 7.7 A comparison between simulated and measured sub-canopy solar radiation. 158 7.8 Comparison between modelled downwelling sub-canopy solar radiation with spatially averaged observations...............................160 7.9 Reduction of solar radiation beneath the canopy................ 162 7.10 A comparison between simulated and measured sub-canopy thermal radiationl63 7.11 Comparison of sub-canopy wind speed estimated by the model with observations beneath a fir canopy............................... 166 7.12 Simulation of fir depth driven by observed sub-canopy wind.......... 167 7.13 Simulation of snow depth beneath fir canopy with no reduction of wind speed168 7.14 Comparison between simulated and measured temperature, density and grain size profiles.............................. 169 7.15 Simulated temperature profile comparison with thermistor array measurements............................................. 172 7.16 Temperature profile simulation at the thermistor array location........173 8.1 Cumulative Energy Contributions at RME Ridge and Fir Sites......... 180 8.2 Relative contributions of radiative and turbulent energy contributions at RME ridge site.............................. 181 8.3 Relative contributions of radiative and turbulent energy contributions at RME fir site................................... 182 8.4 Relative contributions of radiative and turbulent energy contributions at RME ridge site: melt period days 77-85........................ 183 x

8.5 Relative contributions of radiative and turbulent energy contributions at RME fir site: melt period days 77-85.................................... 184 8.6 Relative contributions of thermal energy from the snow, canopy and sky.. 186 8.7 LAI dependent net thermal and net solar radiative energy during one melt period at a fir site (Julian days 77-84)............................. 188 8.8 Variation of thermal to solar radiative energy ratio with leaf area index. This is time and location specific................................. 189 8.9 Simulated variation of snow depth from canopy variability................ 190 8.10 Seasonal and diurnal variation in above-canopy albedo................. 192 8.11 Simulated snow depth at RME ridge site, with fir canopy radiative and turbulent effects added, compared to the snow depth measured at the open site............................................................ 195 xi

List of Tables 2.1 Wetbulb temperature dependent new snow density.............. 38 3.1 Error functions from model sensitivity to parameterisation..............47 3.2 Error functions from model sensitivity to errors in forcing data measurement. 69 3.3 Error functions from model sensitivity to errors in forcing data measurement. 71 4.1 Sublimation losses (% of above-canopy snowfall)............... 82 7.1 Sensitivity of model at RME ridge site..........................152 7.2 Correlation coefficient comparison between measured and simulated subcanopy solar radiation..................................159 7.3 Sensitivity of model to the canopy parameters at RME fir site...........175 xii

Notation Al freezing curve parameter [Unitless] AL Leaf area index [m2 m-2] B lower boundary (snow) reflectance CE bulk transfer coefficient for latent heal [Unitless] CH bulk transfer coefficient for sensible heat [Unitless] ck heat capacity of constituent k [Jkg-1K-1] CkT variation of saturation vapour pressure relative to phase k [Nm-2K-1] De effective diffusion coefficient [m2s-1] Deo standard diffusion coefficient (at 1000mb and 273.15K) [0.92 x 10-4m2s-1] E3 Third exponential integral EEO windless convection coefficient for latent exchange over snow [fHm-2K-1] EHO windless convection coefficient for sensible exchange over snow [Wm-2K-1] fl fraction of liquid water [unitless] ficap fraction of liquid water available for capillary motion [unitless] Fm, Fs error functions for sensitivity analysis frh fractional humidity [unitless] fvis proportion of visible light to near-infrared radiation [unitless] g acceleration due to gravity (9.81 ms-2) Gdry, Gwet dry and wet grain growth parameters (5 x 10-7) and (4 x 10-12) [m4kg-1] hk enthalpy of constituent k [Jkg-1] lat surface latent heat flux [Wm-2] Iprcp heat convected by precipitation [Wm-2] IswN Isw downwelling / upwelling solar radiation [Wm-2] Isen surface sensible heat flux [I7m-2] Isurface surface energy flux [Wm-2] ILWN ITW downwelling / upwelling longwave radiation [Wm-2] LW'1 LW xiii

k von Karman coefficient (0.4) ke effective thermal conductivity [WK-lrn-1] kk thermal conductivity of constituent k [WK-lm-1] Kv canopy light extinction coefficient Lli latent heat of fusion (3.335 x 106Jkg-1) Lvi latent heat of sublimation (2.838 x 106Jkg-1 @273.15K) Lvi latent heat of vapourisation (2.5045 x 106Jkg-1 @273.15K) Mk,k' mass rates of melt, evaporation and sublimation [kgm-3s-1] Patrnos atmospheric pressure [mb] PI intrinsic water pressure [mb or Nm-2] Pv,air vapour pressure of air [mb] Pvk,air vapour pressure of air at snow surface [mb] Pv0,sat saturation vapour pressure at 273.15K [mb] Q energy gain [Wm-2] Qrh relative humidity [unitless] Qtz soil quartz content R gas constant (461.296 Jkg-1K-1) R canopy reflectance coefficient Ri Richardson number [unitless] s saturation [unitless] Se effective saturation [unitless] serr effective saturation error criterion [unitless] Si irreducible saturation (minimum amount of water retained by a substance under drainage) [unitless] T nodal temperature [K] T canopy transmission coefficient Tair air temperature [K] TD depression temperature (= 273.15 - T) [K] Terr temperature accuracy criterion Uk mass flux of component k [kgm-3s-1] vk velocity of component k [ms -1] w, Wo, Wf -wind speed, wind speed at an open site, wind speed beneath a forest canopy [ms-1] z height of node [m] Z' distance to snow surface [m] zo surface roughness length [m] xiv

a Pnir 7crit fYk At Az x r, ro Kmax A A A /o P0 "E, <IH, 'M Osd Pk Tatm Ts Ok ~d Or shortwave albedo [unitless] bulk extinction coefficient for solar energy at near infrared wavelengths [m-1] bulk extinction coefficient for solar energy at visible wavelengths [m-1] critical density parameter for metamorphism compaction [kgm-3] partial density of component k [kg-n3] time step [s] layer width [m] number of days between precipitation events snow emissivity [unitless] vegetation emissivity [unitless] viscosity coefficient, viscosity coefficient at 0~C and Yw = 0 [N s m-2] hydraulic permeability [mn2] saturation permeability [7n2] fractional litter coverage [m2 m-2] radiation wavelength [nm] Daily litter rate [m2 m-2 day-1] Cosine of solar zenith angle Dynamic viscosity of water at 0~C (1.787 x 10-3Ns m-2) Metamorphism compaction parameter porosity [unitless] stability functions for transfer of water vapour, sensible heat and momentum solid porosity i.e. fractional volume of gaps in the solid (ice plus dry soil) matrix [unitless] intrinsic density of component k [kgm-3] Stefan-Boltzmann constant (5.669 x 10-8W m-2 K-4) atmospheric optical depth optical mass-dependent depth of snow fractional volume of component k [unitless] grain diameter [m] grain radius [m] XV

Subscripts The physically based snow model assumes 5 separate components: dry soil, ice, liquid water, water vapour and air. The following subscripts are used throughout the text to indicate which component the symbol represents. Where these quantities are generalised, the subscript k is used. Replacement of k with any other subscript throughout an equation gives separate component equations. a air d dry soil i ice k General indication that quantity or equation refers to one component 1 liquid water p precipitation v water vapour w total water (ice + liquid) Additional subscripts are as follows: diff diffuse solar radiation dir direct solar radiation f beneath forest canopy LW thermal radiation nir near-infrared radiation o at an open site SW solar radiation vis visible radiation xvi

Acronymns FOV GCM LAI NWRC OJP RCEW RME RMSP Snow-SVAT SSA SVAT SWE USDA Field of view General Circulation Model Leaf Area Index Northwest Watershed Research Center Old Jack Pine site Reynolds Creek Experimental Watershed Reynolds Mountain East Reynolds Mountain snow pillow site Snow-Soil-Vegetation-Atmosphere Transfer Scheme Southern Study Region Soil-Vegetation-Atmosphere Transfer Scheme Snow Water Equivalent United States Department of Agriculture xvii

Chapter 1 Introduction 1.1 Project Justification Snowmelt provides a significant proportion of the yearly water budget in areas such as Northwest America, Northern Eurasia, Austria and Switzerland. Water plays a vital role in the terrestrial biome both as a resource and as a carrier for nutrients in solution and in suspension. Prediction of the timing and magnitude of melt water runoff is vital to local communities to enable this resource to be used most efficiently, yet the effect of forest cover on snowpack ablation is poorly understood. Some researchers have found that the snow beneath a forest melts before the snow at an open site [Golding and Swanson (1986)], whereas others have found that the opposite is true [Rouse (1984), Hardy and Hansen-Bristow (1990)]. The dominant energy term that contributes to snowpack ablation beneath a forest canopy is also unclear, but generally has been found to be radiative, either longwave [Adams et al. (1996)] or shortwave [Hardy et al. (2000)]. 1

1.1. Project Justification Chapter 1. Introduction General circulation models (GCMs) are used to examine climate trends such as possible global warming. Parameterisation of snow-covered areas has a dramatic impact on model performance [Henderson-Sellers and Wilson (1983), Thomas and Rowntree (1992), Marshall et al. (1994)], particularly in forested areas [Betts (2000), Viterbo and Betts (1999)]. Many models assume a constant above-canopy albedo, but the seasonal and diurnal changes in this parameter have not been studied in detail. There is a clear need to incorporate snow-forest dynamics into GCMs. This problem is especially important, because vast areas of forested land exist, which are snow-covered for substantial amounts of time. For example, the boreal forest is the world's largest terrestrial biome, covers 15% of the Earth's land surface and is itself snow-covered for 6-8 months of the year. Simulation of the evolution of snow in forested areas such as the boreal forest would improve both local and large scale hydrological modelling, and enable improved parameterisation of forested snowpacks for climate studies. Accurate simulation of the accumulation and ablation of the snowpack beneath a forest necessarily has a number of requirements and assumptions. Ideally, the model should be capable of, or demand the following: * Examination of the optical and thermal radiation interaction between the forest cover and the snow cover. * Accurate simulation of snowpack physical properties such as temperature and density profiles, in addition to bulk properties such as snow depth. * A good physical basis for both snow and canopy models. * A minimal extra number of parameters to represent the forest cover, and these should 2

1.1. Project Justification Chapter 1. Introduction be easily obtainable in the field or determined from remote sensing. Minimising computational speed and efficiency is crucial for large-scale models such as GCMs. * Model forcing data that are common meteorological variables measured outside the canopy (either above the canopy or at a nearby open site). Field sites rarely have meteorological stations below the canopy. * Simulation of melt water runoff for use in hydrological models. No existing model has captured all of the features described above. This thesis describes a new Snow-SVAT (snow-soil-vegetation-atmosphere transfer scheme) developed within this project, by the coupling of an optical and thermal radiation canopy model with a physically-based snow model. The Snow-SVAT incorporates all of the outlined features, and has been used to examine the following propositions: 1. Snowmelt beneath the forest is dominated by thermal radiation 2. Complete ablation of the snowcover beneath a fir canopy occurs prior to complete ablation at an open site 3. Natural canopy variability leads to significant differences in snow depth, from changes in the radiative energy input 4. Changes in the snow albedo have a significant effect on the above-canopy albedo 5. Above-canopy albedo exhibits strong seasonal and diurnal variation, and cannot be represented with a single, constant value in general circulation models. 3

1.2. A Walk in the Woods Chapter 1. Introduction As will be discussed in chapter 8, propositions (1) and (2) are site-specific and may depend on canopy characteristics, location and time of year. However, a snow-covered forest is far from simple to model. In the following section, natural forest features and processes which could decrease model performance are described. Although the model is currently one-dimensional, these features are two- or three-dimensional and are discussed as such. It is acknowledged that a complete model including all these features could have an increased computation time. 1.2 A Walk in the Woods Figure 1.1: Conifer forest, Reynolds Creek Experimental Watershed, Idaho, USA A natural forest has a high degree of spatial variability. The trees are generally distributed in a non-uniform fashion, governed by biochemical sustainability. In addition the tree heights differ within a forest, due to tree age and other limiting factors. Both 4

1.2. A Walk in the Woods Chapter 1. Introduction these factors lead to a heterogeneous canopy, which in turn affects the amount of radiation reaching the snow surface. Shadows cast by the canopy onto the surface below indicate the spatial variability of radiation. Not only is tree-to-tree variability evident, but also the internal structure of individual tree canopies can be seen clearly. Bright patches on the surface occur, where radiation has an uninhibited path through the canopy, whereas shadows form when radiation is intercepted by the canopy. The shadows change over the course of the day, and so instananeous radiation variability may be much greater than the diurnal spatial variability of radiation beneath a continuous canopy. Radiation is additionally affected by the topography of the landscape. Sloped surfaces influence the radiation budget by shading and by reflection to and from adjacent slopes. The underlying ground surface also affects the accumulation distribution. Localised features in the topography result in greater accumulation in dips and basins and decreased accumulation on bumps and hills. The distribution of accumulated snow is affected by the forest. Falling snow is intercepted by the forest canopy and may remain balanced on, or bonded to, the canopy for some period of time. Figure 1.2 shows a large clump of snow intercepted on the canopy. Intercepted snow is removed from the canopy by sublimation and by slipping to the snow surface when subject to partial melt. Excess precipitation falls though the canopy once the branch loading capability has been exceeded. As snow has a much higher albedo than green biomass, intercepted snow may affect the radiation balance. Precipitation reaching the forest floor is distributed heterogeneously as a result of the interception of snow by a variable canopy. Accumulation beneath the tree canopy, 5

1.2. A Walk in the Woods Chapter 1. Introduction Figure 1.2: Snow intercepted on a conifer canopy in RCEW, Idaho, USA in general, is reduced through canopy interception losses, and the snow surface takes the form of a well around the tree trunk [Hardy et al. (1997a)]. Snow depth increases to a maximum at the canopy edge. The heterogeneous distribution of snow beneath the canopy may be enhanced by thermal emission from tree trunks, as illustrated in figure 1.3. Snow within the tree well experiences a much sharper temperature gradient than deeper snow further away from the trunk, which changes the characteristics of the snowpack considerably. Leaves or needles and other organic debris (including fallen logs) from the forest accumulate on the snow surface and may be subsequently buried within the snowpack by accumulated snow [Hardy et al. (2000)]. Figure 1.4 shows the presence of litter on the snow surface. Litter has a lower albedo than a pure snow surface, hence absorbs more solar radiation and subsequently emits more thermal radiation. The presence of organic litter may substantially alter the radiation balance at the snow surface, and within the 6

1.2. A Walk in the Woods Chapter 1. Introduction Figure 1.3: Photograph of tree well around aspen in Reynolds Creek Experimental Watershed snowpack. Understory (or sub-canopy) vegetation also affects the energy balance. Protruding vegetation (as shown in figure 1.4), has a lower albedo than snow and similarly to surface debris, absorbs more solar radiation and emits more thermal radiation. Submerged and partially submerged vegetation enable pockets of air to form within the snow cover, allowing greater transport of heat and mass thereby reducing stratification and local heterogeneities in snow characteristics. Problems associated with understory vegetation can be greatly reduced for validation of the model by appropriate choice of field site. Operational forecasting of melt events in a catchment with a high degree of sub-canopy biomass, however, requires closer examination of this matter. Finally, one further complication arises from the wildlife prevalent and specific to the area. Examples of these include beaver dams in Vermont, which change biomass in a 7

1.3. Method of Approach Chapter 1. Introduction Figure 1.4: Heterogeneity of surface and sub-surface organic matter. local area, affecting both radiation and snow accumulation from drift. Another example is vole tunnelling within the snowpack. These tunnels provide preferential pathways for melt water during the ablation phase. For an unsaturated snowpack, melt water flows around the tunnel walls, but water will flow through the tunnel under saturated conditions, and melt water is removed from the snowpack in a shorter time [Jordan (1995)]. 1.3 Method of Approach The need for a model that was capable of accurate simulation of the snow cover beneath a forest canopy was discussed in section 1.1. Within the model, radiation interaction between the canopy and snow surface should be well represented, as should the internal snow processes of the snowpack. Two models were chosen to form the Snow-SVAT basis: (i) a multi-layer phys 8

1.3. Method of Approach Chapter 1. Introduction ically based snow model capable of representation of snowpack processes such as metamorphism and compaction and (ii) an optical and thermal radiation canopy model. The canopy radiation model was coupled to the snow model in such as way that the simulated snowpack properties affect, and are affected by, the sub-canopy radiation balance. The forest effect on turbulent transfer of energy to the snow surface was simulated through a reduced wind speed algorithm. The sublimation loss of snow intercepted in the canopy was assumed to be a simple function of precipitation. The Snow-SVAT was initially tested through simulation of snow depth during accumulation and ablation periods beneath a boreal pine forest during water year 1993 -1994.1 Then snow cover evolution beneath a fir forest in Idaho, USA2 was simulated for water year 2000-2001 and the model was validated against measurements of snow depth, sub-canopy up- and downwelling optical and thermal radiation and profiles of density, temperature and grain size within the snowpack. Spatial variability of the forested snowpack was examined by sensitivity analysis of the Snow-SVAT to the canopy parameterisation. Comparisons were made between model simulations at an open, non-forested site and a forested site. The Snow-SVAT was used to examine the propositions stated in section 1.1. A summary of the contents of each chapter of this thesis and how they contribute to this method of approach follows in the next section. lData were collected as part of the BOREAS campaign, as described in chapter 5. 2Reynolds Mountain East Experiment 2000-2001 was designed specifically to collect validation data, and is described in chapter 6. 9

1.4. Thesis Description Chapter 1. Introduction 1.4 Thesis Description The need for a new Snow-SVAT has been introduced in Chapter 1. Natural forest features that complicate the energy and mass balance of the snowcover have been described and the method of approach of this project has been outlined. A snow model has been chosen to be incorporated within the Snow-SVAT, and this model is described in Chapter 2. The description includes the model methodology and boundary conditions. Chapter 3 describes an investigation of the sensitivity of this snow model to its parameterisation and to errors in the initialisation and forcing data, using a homogeneous snowpack under idealised meteorological conditions. A number of parameters and initial conditions have been selected for further study with real data, because the model sensitivity was high, because the uncertainty was high, or because the model was expected to be more sensitive under real conditions. The vegetation radiation model chosen to couple to the snow model is described in Chapter 4. The coupling method of approach, and the advantages over existing pre-processing approaches is also discussed. Assumptions made in the formation of the Snow-SVAT are justified, which include the reduction of wind speed beneath the canopy, sublimation of intercepted snow, and reduction of snow surface albedo with litter and for a thin snowpack. Simulations of the snow depth beneath a boreal coniferous forest canopy are presented in chapter 5. Data taken during the BOREAS campaign have been used to test the Snow-SVAT, but the lack of validation data has meant that a new experiment was required to collect data for model validation. Chapter 6 describes an experiment that was carried out in water year 2000-2001 10

1.4. Thesis Description Chapter 1. Introduction within Reynolds Creek Experimental Watershed, Idaho, USA. Data were collected to initialise and drive the Snow-SVAT. Validation data collected includes sub-canopy solar and thermal radiation, snowpack profiles of temperature, grain size and density, and automatic measurements of snow depth. Problems and difficulties experienced in the field are also discussed in Chapter 6. Data collected during the Reynolds Mountain East field experiment have been used to simulate the accumulation and ablation of snowcover beneath a fir canopy. The snow model parameterisation was derived from a validated open site simulation, therefore no model calibration was required at the fir site. Model simulations generally agreed well with observations, and are demonstrated in Chapter 7. Because the model has been proven to represent physical processes will beneath a fir canopy at Reynolds Mountain East, Chapter 8 examines physical processes outside the limits of the RME data. The effect of canopy density on the radiative energy balance and subsequently snow depth is addressed, and also the effect of canopy density and snow albedo degradation on the diurnal and seasonal variation of above-canopy albedo is discussed. The thesis propositions are addressed in Chapter 9, and a summary of this thesis, and directions for future work are outlined. 11

Chapter 2 Snow Energy Budget Model 2.1 Introduction The need for physically based snow energy budget models in forested areas was described in section 1.1, and focussed on hydrological management and climate modelling. Physicallybased snow models also have a role in avalanche prediction. Adequate representation of the physics that cause snowpack instabilities reduce the need for in situ measurements, which are often difficult or dangerous to take. Detailed avalanche model requirements are less important here, but the snow model required for this project should use physical principles to represent energy and mass transfer within and to the snowpack, and show good agreement with field observations. A snowpack is extremely dynamic, undergoing compaction, densification, phase change and crystal growth. The seasonal snowpack is stratified, with layers of snow crystals corresponding to different precipitation events, which are modified by post-depositional energy and mass transfer. Each layer may have a different density, temperature, average 12

2.1. Introduction Chapter 2. Snow Energy Budget Model crystal size and other physical characteristics that are initially governed by the meteorological conditions as the layer is deposited. Energy and mass transfer between each of the layers and between the atmosphere and snow surface result in the evolution of the snowpack with varying meteorological conditions. The snow albedo describes the proportion of solar radiation reflected from the snow surface. Reflection from the snow surface is non-isotropic, exhibiting strong scattering in the forward direction. The proportion of solar radiation reflected through any specific angle can be fully described as a function of the incident radiation angle by the bi-directional reflectance distribution function (BRDF). Absorption of radiation within the snowpack is grain size and frequency dependent. Longwave radiation is attenuated within the top few millimeters of the snowpack whereas shortwave radiation penetrates up to -20cm SWE [Wiscombe and Warren (1980)]. The model chosen should capture the above processes to a good degree of accuracy. Degree-days models [summarised by Melloh (1999)] do not represent the physical processes within the snowpack and require calibration for the specific location simulated. Single- and dual-layer models [Kondo and Yamazaki (1990), Tarboton and Luce (1996), Marks and Winstral (2001)] were considered to be too simplistic for this project, since the stratigraphy of the snowpack is to be simulated. Two suitable multi-layer models that simulate physical processes within the pack existed at the project outset: CROCUS [Brun et al. (1989)] and SNTHERM [Jordan (1991)].1 CROCUS does not consider water vapour phase change, diffusive transport or turbulent surface mass loss, nor does it model the soil. SNTHERM89v4 was chosen in 1Other suitable physically-based models have been developed over the past few years, and are described in the Snow Model Intercomparison Project (SNOWMIP) on the web page http://www.crm.rneteo.fr/snowmip/model-review.htm [13

2.2. Model Description Chapter 2. Snow Energy Budget Model this project since it modelled processes not incorporated within CROCUS. In addition, SNTHERM is well-known, is freely available and good support is provided for its users. In this project, SNTHERM was re-coded from Fortran 77 into Fortran 90. This was undertaken in order to make use of improved language features, to remove code that is extraneous to this project (and therefore to increase model efficiency) and to aid model comprehension.2 The physical principles behind SNTHERM were retained, but changes were made to specific algorithms. The differences between the resulting model, "SNOWCAN", and SNTHERM are highlighted at the end of this chapter (section 2.3). These differences include changes to the compaction algorithms and new algorithms to determine fresh snow properties and the solar zenith angle, different numerical solution algorithms for the solution of tridiagonal and cubic equations, and a change in style enabled through the use of Fortran 90. SNOWCAN is described in the following section. 2.2 Model Description SNOWCAN is a one-dimensional, discrete, multiple layer model. The model is driven by commonly measured meteorological variables and uses conservation of mass, energy and momentum to simulate the change in physical state of the snowpack over time. The initialisation of the snowpack and the layer system are described in section 2.2.1, boundary conditions are discussed in section 2.2.2, model assumptions are presented in section 2.2.3 and the physical principles of the model are stated in section 2.2.4. The solution methods 2A new version of the model, SNTHERM99 became available during the course of this project, but time constraints limited the model update in this project. SNTHERM99 contains improved representation of turbulent transfer mechanisms, and additionally simulates the effect of litter and a thin snowpack on the snow albedo and the effect of wind on snow densification and mass transport. The improvements to SNTHERM89 incorporated in SNTHERM99 are described by Jordan et al. (1999) 14

2.2. Model Description Chapter 2. Snow Energy Budget Model and time step procedure are described in sections 2.2.5 and 2.2.6 respectively, and model outputs are shown in section 2.2.7. The notation used in this thesis is listed on page xiii. 2.2.1 Model Setup The simulated snowpack and underlying soil are split into discrete, horizontal layers,3 as shown in figure 2.1. The snowpack layers typically correspond to layers observed in the field, as described in section 2.1, so that the stratigraphy of the snowpack is simulated. X.XX..... 3.. X.. * rT r r lll.,l..,?,.....,............_....... Figure 2.1: Simulated snowpack adaptive grid SNOWCAN is run when observed snow is already in situ, therefore snowpack layers are defined at the start of the model, and layers are added by the model as precipitation accumulates on the surface. Initialisation of the snowpack requires the width, temperature, density and approximate grain diameter for each layer. Near the surface, changes in snowpack properties are generally greater and faster. Model accuracy and stability are increased by smoother transitions, achieved by an increase in the number of surface layers., with a decrease in layer width closer to the surface. The model could be initialised prior to the first snowfall, which would require initialisation of soil temperature and water content profiles. However, the accuracy of 3SNTHERM terminology denotes layers as either snow or soil type. Here, layers are used to describe grid elements within each material type. 15

2.2. Model Description Chapter 2. Snow Energy Budget Model the model would be much more dependent on the accuracy of the model precipitation algorithms. The simulation of an in situ snowpack is considered to be a more accurate and stable approach. The simulated snowpack is based on a non-Eulerian grid, which deforms with the compacting snowpack, as shown in figure 2.1. In this way, layers largely retain the original sample of snow throughout the model run, with the exception of liquid and vapour transport to and from the layer, as described in section 2.2.4. Height is defined as positive upwards from the snow-soil interface, hence fluxes within the snowpack (and soil) are also defined as positive upwards. This is in constrast to the surface fluxes, which are defined as positive downwards by meteorological convention. 2.2.2 Boundary Conditions Once the snowpack has been initialised, meteorological data are required to drive the model. These forcing data are: * Incident solar radiation * Incident thermal radiation * Air temperature * Wind speed * Relative humidity * Precipitation The need for these forcing data stems from Neumann boundary conditions that constrain the physical equations described in section 2.2.4. The boundary conditions 16

2.2. Model Description Chapter 2. Snow Energy Budget Model chosen for SNOWCAN, following the work of Jordan (1991) are described below. Lower boundary mass transfer Liquid water and vapour transport within the soil is not modelled, therefore the lower boundary conditions for mass transfer are set at the snow-soil interface, where excess mass is drained artifically. This not only affects the rate of heat transport within the soil and between the soil and snowpack, but also degrades the model as a hydrologic runoff tool. Dirichlet lower boundary conditions are specified as: Uz<o = 0 (2.1) Uv Iz< = 0 (2.2) Upper boundary mass transfer Mass transfer to the snowpack surface results from precipitation (liquid and solid water), and sublimation / condensation processes result in mass loss to, or gain from, the atmosphere. Precipitation mass flux, Up, is defined as: Up -— 7p SWEp (2.3) where the density of precipitation, -p, is determined from the wetbulb temperature, as described in section 2.3. Precipitation snow water equivalents, SWEp (m s-1), are converted from the precipitation forcing data (SWE m hr-'). Latent exchange of energy is described by equation 2.17 on page 22. The mass transfer at the surface by turbulent processes (k = i, 1) is: EEO + EE w) Ulat = -EO (PEv,air - frh Pvk,sat) (2.4) Lvk 17

2.2. Model Description Chapter 2. Snow Energy Budget Model Lower boundary energy transfer The heat flux from the soil constitutes the energy transfer from the lower boundary condition. A Dirichlet condition maintains constant heat flux by specification of a constant temperature: 9T T =0 (2.5) at soil base The constant temperature of the lowest soil node is set by the initialisation of the snowpack (see section 2.2.1). Upper boundary energy transfer The energy transfer at the snow surface is a dynamic interaction between the snowpack and the atmosphere. A Neumann boundary condition is used to specify heat transfer at the upper boundary: ke =Ts Issurface (2.6) ez z= surface The surface energy balance is defined as: surface -= ISW + ISW + LW +- I+W + Isen + Ilat + Iprcp (2.7) By meteorological convention, fluxes are designated as positive downwards (into the pack). 1. Solar radiation. Incident shortwave radiation (300-3000nm) is input directly from meteorological measurements. The proportion of incident solar radiation that is reflected is given by the snow albedo, as. As described in section 2.1, snow reflectance is non-Lambertian 18

2.2. Model Description Chapter 2. Snow Energy Budget Model and is also frequency dependent. A model that simulated the bi-directional reflectance distribution function (BRDF) of solar radiation from a snowpack was previously coupled to SNTHERM by Glendinning (1997). However, that version of SNTHERM was not used here, since validation data are difficult to collect, particularly beneath forest canopies, as described in section 6.6.1. The simpler 2-stream albedo algorithm(2.12), where snow grains are assumed to be spherical, was used instead [Hardy et al. (1997a), adapted from Marks (1988), which was based on Marshall and Warren (1987)]. Direct and diffuse albedo for clear days is given separately for visible and nearinfrared (the notation used in this thesis is given on page xiii): avis, diff = 1 - 2 / (2.8) anir, diff = 0.85447 x e-21.23 (2.9) Cvis, dir vis, diff + 1-575 V ("- - - 0) (2.10) Canir, dir = nir, diff + (2 Gr + 0.12) (- o) (2.11) The shortwave albedo on a clear day is therefore given as: is SW, diff + a dir SW, dir O~s - Ivis, ciff -— iI- + Q-vis, dir i --- f vis SW_ SW d + aOnir, diff SW -diff+ ~nir, dir SW i (1 - fvis) (2.12) sw sw / where the fraction of visible to near-infrared radiation is given by fvis = 0.436 for clear days, and the proportion of direct to diffuse light is either input from shadowband radiometer measurements, or estimated from the atmospheric optical depth (1sw, dir = Isw x e-7atm/l'o). It should be noted that fvis may be different 19

2.2. Model Description Chapter 2. Snow Energy Budget Model beneath a canopy because of preferential wavelength absorption by the canopy, which is not considered explicitly in the canopy model described in chapter 4. The solar radiation that is not reflected penetrates the snowpack. Solar radiation absorption within the snowpack is governed by the snow grain size, snow density and radiation frequency. Absorption of near-infrared radiation occurs near the surface of the snowpack, assumed to be within the top 2mm in this model. At visible wavelengths, the radiation penetrates much deeper within the pack, so is described by physics within the model. The description is nevertheless included here within the boundary conditions section as it is analogous to the near-infrared absorption. At visible wavelengths, the absorption coefficient is estimated as: 0.003795 w (2.13) Pvis V~ (2.13) 01 At near-infrared wavelengths, the absorption coefficient, /3nir, is set by the user, and all near-infrared radiation absorption occurs near the surface of the snowpack (assumed to be within the top 2mm). t3nir effectively allows the user to determine the proportion of solar radiation that is absorbed at the surface, which was investigated by Brandt and Warren (1993) and Koh and Jordan (1995). An absorption coefficient of /3nir = 400 m-1 corresponds to ~ 55% of net solar energy absorption occuring at the surface of the snowpack. The energy gain from solar radiation absorption for the surface layer is: Qsw, surface = Is (1- a) (1 - - vi- e- ' -0002 nir) (2.14) The energy gain from solar radiation absorption for layers within the snowpack is: Qs W -= I 1 (1 -e/-visAz)j (2.15) 20

2.2. Model Description Chapter 2. Snow Energy Budget Model where Az is the layer width and j refers to the layer number, so I$w l represents the solar radiation penetration from the layer above. All radiation is assumed to be absorbed within a maximum mass of snow. For snowpacks with less mass than this maximum, all remaining radiation that has penetrated the snowpack is absorbed at the soil surface. 2. Thermal radiation. Meteorological measurements of downwelling thermal radiation incident on the snow surface, Ijw, are used to drive the model. Most of the incident thermal radiation is absorbed at the surface, given by the snow emissivity (E5), and the remaining small proportion is reflected. Thermal radiation emitted from the snow surface depends on the surface temperature, and is estimated from Kirchoff's Law for emission of thermal energy (EaT4). The total upwelling thermal radiation is therefore: I -W =-( ) IV -e - w T (2.16) Snow has similar thermal characteristics to a black body, therefore an emissivity of 0.99 was chosen in this model to parameterise the thermal radiation component of the upper boundary condition for the energy balance equation. 3. Turbulent transfer. Turbulent transfer of heat at the snow surface occurs through convection sensible heat and sublimation / evaporation (latent heat) processes. The turbulent fluxes of 21

2.2. Model Description Chapter 2. Snow Energy Budget Model latent and sensible heats are given by Andreas and Murphy (1986) as Ilhat = (EEO +- 1 CE W)(Pvair - frhP,sat) (2.17) EE Isens = (EHO + PairCairCH W)(Tair-Tn) (2.18) EH where the windless exchange coefficients for latent and sensible heat, EEO and EHO are input by the user. The bulk transfer coefficients, CE and CH are estimated from the roughness length (input by user) and measurements of wind speed, air temperature, relative humidity and the heights at which these measurements are taken. A full explanation of the turbulent transfer algorithms used in SNOWCAN is included in Jordan (1992). 4. Precipitation heat. When snow starts to fall, a new layer is created within the model, since the precipiation may have totally different physical properties to the snow beneath. Only rain, which is treated in SNOWCAN by analogy to liquid water flow from a layer above, increases the energy in the existing top node by: Iprcp = -c1 Up Tp (2.19) This term is negative, since fluxes within the pack are defined as positive upwards, in constrast to surface energy fluxes. 2.2.3 Snow Model Assumptions Assumptions have been made regarding the geometry and physics of the snowpack, as described below. Truncational and numerical assumptions are also listed in this section. 22

2.2. Model Description Chapter 2. Snow Energy Budget Model 1. Snow grains are assumed to be spherical. This affects solar radiation scattering and absorption both within the snowpack, and at the surface, as described in the upper boundary condition for energy in section 2.2.2. 2. Solar energy incident on the snow cover is split into 2 broad-bands: Visible (300 -700nm) and near-infrared radiation (700-3000nm). 3. Ice crusts due to melting and subsequent refreezing and liquid water ponding on the top of impermeable ice layers are not taken into consideration in this model. This is a particularly important effect for snowpacks overlying a sloped topography, where substantial amounts of liquid water runoff may occur over the top of the impermeable layer. However, sloped surfaces are not simulated in SNOWCAN. 4. Soil is assumed to be immobile and incompressible. Vapour and liquid water flow are disallowed in soil in this model. Liquid water and water vapour are artificially drained at the snow-soil interface, as described in the lower boundary conditions for mass (section 2.2.2). 5. Flows are constrained to the vertical direction in this 1-dimensional model. Water flow is assumed to be gravitational, since the capillary pressure effect of the gaseous constituents is not taken into consideration for this model. 6. Horizontal homogeneity is assumed, so that bulk densities and source terms are constant over the element volume, AV. The formation of vertical flow pipes within the pack, which increases the rate of liquid water flow to the base of the pack, is not considered in this model. 23

2.2. Model Description Chapter 2. Snow Energy Budget Model 7. Liquid water is assumed to have the same velocity as the compacting ice matrix. This results in the conservation of the partial masses of ice and liquid water during compaction: (7kAz)t= ( kAZ)t+At. 8. The densities of water vapour and dry air remain invariant during matrix deformation, resulting in the expulsion of a portion of gas from the contracting volume. Temperature effects on gaseous density are also not considered. 9. The air is assumed to be at atmospheric pressure throughout the pore volume. This necessarily implies that all pores are interconnecting. 10. All constituents within a layer are assumed to be in thermal equilibrium. 11. Water vapour density reaches equilibrium relative to ice when the fractional volume of liquid water is less than or equal to an arbitrary value of 2% [Jordan (1991)] and with liquid water for a fractional volume greater than 2%. 12. Gaseous components constitute less than 1% of the total mass. Therefore their contributions to the total bulk density, total specific heat and total thermal conductivity are neglected in this model. 13. For model discretisation, it is assumed that diffusive/conductive processes are centralised difference and convective processes are stepwise. 2.2.4 Physical Principles There are three fundamental conditions that must be satisfied in the physically based snow model: conservation of mass, conservation of momentum and conservation of energy. These conservation equations are given below. The boundary conditions required to form a 24

2.2. Model Description Chapter 2. Snow Energy Budget Model closed set of equations were stated in section 2.2.2 and the method of solution is described in section 2.2.5. As stated in section 2.1, this snow model is based on SNTHERM. Jordan (1991) and Jordan et al. (1999) contain a more rigorous explanation of the snow physics used in this model. Conservation of Mass A fundamental assumption made by mixture theory is that the conservation of mass can be applied simultaneously to the total medium and to the separate constituents independently. The equations for mass balance can be written: Total Medium: t JptdV Z Uk.dS (2.20) at Jk Water Constituents ( k, k' = i, 1, v) Ot PkdV = -j Uk dS + Mk,k,dV (2.21) at Jk' Dry Air: 9 | pdV = - Ua dS (2.22) The phase change term is the last term in equation 2.21. Mk,k, denote the mass rates of melt, evaporation and sublimation. Note that Mk,k' = -Mk',k and that Mk,k = 0 Water vapour flux, U, is determined from the diffusion equation: OT Uv = De, snow CkT oT~ (2.23) The liquid water flux, U1 is determined from conservation of momentum. 25

2.2. Model Description Chapter 2. Snow Energy Budget Model Conservation of Momentum Under the assumption that the snowpack is isotropic and the water is a Newtonian fluid undergoing negligible divergence, the conservation equation for the flow of water through the snow pack is:: - ilvldV+ lvl ~ dS - 0 Pli - dS- J — ig+ I dV ot Nv s s v inertial term convection pressure gravity viscous stress /v(Miv + i -Mv V V )dV phase change (2.24) The inertial, convective, phase change and capillary terms are negligible in comparison to gravitational and viscous stress [Jordan (1991)], hence equation 2.24 simplifies to give the expression for the liquid water mass flux, equation 2.25. Ui -= -p g (2.25) which is the gravitational form of the empirical Darcy Law. Liquid water flux can also be represented in terms of effective saturation (se = 1-) using the Brooks and Corey (1964) relation for hydraulic permeability (equation 2.26). Ni = lmaxSE (2.26) where e is a quantity which depends on the pore size distribution and is suggested by Colbeck and Anderson (1982) to be approximately equal to 3, from their experimental data. hmax is the saturation permeability and is approximated from the grain size and density 26

2.2. Model Description Chapter 2. Snow Energy Budget Model by Shimizu (1970) as Kmax 0.077 (id e(-00078t ) (2.27) The effective saturation is the independent variable, which is determined from the solution to the mass balance, as described in section 2.2.5. Conservation of Energy Equation 2.28 describes the conservation of energy for each layer, where the external energy flux is defined as positive downwards. a J pthtdV = - Ukhk dS + ktVT.dS Ot kv"JsJ mass flux conduction + I - dS (2.28) Js external flux The external energy flux for the top layer was determined from the boundary conditions (I = I0op), as described in section 2.2.2. For internal layers, only solar radiation contributes to the radiative term, since all other external energy is absorbed at, or close to, the surface. The absorption and scattering of solar radiation within the snowpack was also included in section 2.2.2 The independent variable in the energy balance equation is the temperature, and solution of the energy conservation equations is discussed in section 2.2.5. Compaction and Metamorphism of the Snowpack The snowpack is subject to the temperature and pressure gradients across it. With time, the snow crystals change and undergo metamorphism. Newly fallen snow crystals typi27

2.2. Model Description Chapter 2. Snow Energy Budget Model cally have sharp edges. Since the vapour pressure over edges is much higher than over flat surfaces, the sharp edges start to sublime away. Over time, crystals subject to low temperature gradients will tend to form spherical grains. Compaction due to the weight of overlying snow occurs. Fresh dendritic snow has a greater strength than other crystal forms, as the sharp grain shapes lock in to each other, which reduces the compaction rate. As low temperature gradient or equitemperature metamorphism tends the snow crystals towards spherical shapes, this strength is lost, and the snow grains can geometrically fit closer together. Anderson (1976) proposed an empirical compaction rate due to this process, which has been modified in SNOWCAN (equation 2.29) to improve the simulated snow depth agreement with field observations. A comparison of the compaction algorithms used in SNOWCAN with those used in SNTHERM is presented in section 2.3. 1 O -0.04(273.15-T) (2.29) AZ Ot metamorphism where ~ = 2.0835 x 10-5 for the simulations presented in this thesis, and may vary for different sites. This value is greater than the SNTHERM parameterisation of ~ = 2.778 x 10-6 [Jordan (1991)], which is discussed in section 2.3. If the bulk ice density is above a critical value4, ycrit, the compaction rate is slower and equation 2.29 is multiplied by a factor e-0~046(7i-7ycrit). For wet snow, (-yi > 1%) the compaction is estimated to occur twice as fast and equation 2.29 is modified by a factor of 2. After the new snow has settled, the compaction process is due mainly to the pressure of the overlying snow. Anderson (1976) suggested that this much slower process 4This is set by the user, but is suggested by Jordan (1996) as 100 Kg m-3 28

2.2. Model Description Chapter 2. Snow Energy Budget Model is described by: IOAz -P (2.30) A\z at overburden t7 where 7 is the viscosity coefficient and is defined in the snow model (and the value of 770 is input by the user) as: 71 = 770 x e0 08(273.15-T) x e0.023p, (2.31) Snow crystals which are subject to high temperature gradients tend to form a different structure called depth hoar. The high temperature gradient causes sublimation from the top of the crystal which then refreezes on the bottom of the crystal above, forming a regular structure. In this way, crystals grow in size and the snow can be less dense as a result of the hollow structure of depth hoar crystals. The snow grain size increase with temperature gradient is modelled in SNOWCAN by equation 2.32. Sublimation is a surface-area and temperature dependent process. Smaller crystals and crystals with a high surface structure (e.g. dendritic) have a higher surface area to volume ratio than large, rounded crystals, and therefore greater sublimation per volume. Vapour diffusion from areas of high vapour concentration and subsequent condensation means that large crystals will grow larger at the expense of smaller crystals. This will have two effects: firstly the average particle size will increase and secondly (due to the upward vapour flux) the distribution of larger grain sizes will move upwards. The proposed description of this mechanism (Jordan (1991)) is: Od= Gdy for dry snow (2.32) at Od ad = e(0 + 0.05) for wet snow (01 < 9%) (2.33) atd = (0.14) for wet snow (0 > 9%) (2.34) at~ - we (0.14) f or wet snow (0, > 9%) (2.34) 29

2.2. Model Description Chapter 2. Snow Energy Budget Model where Od is the mean grain diameter and U, is the average mass vapour flux. The inverse diameter dependence models the reduction in grain growth with increased grain diameter. Equation 2.33 models a linear increase in grain growth rate with liquid water content. Another type of metamorphism which may occur is due to the melting of the snow during the day and the refreezing of the water during the night. This leads to large clusters of strongly bonded round crystals. As the liquid water content increases later on in the melt season the bonds become weaker. However, the effect of melt-freeze events on snow grain sizes are not simulated in this snow model. 2.2.5 Numerical Solution Techniques Simulation of a snowpack occurs over a period specified by the meteorological forcing data, which are typically provided at hourly intervals. SNOWCAN has an adaptive time step to optimise both accuracy and computational efficiency, as discussed in section 2.2.6. Meteorological data (downwelling solar and thermal radiation, air temperature, wind speed and relative humidity) are interpolated between data points to give a more continuous data record. Precipitation is not interpolated since measurements are hourly averages, therefore precipitation is assumed to be constant over the hour. The model solves the mass balance equation first to determine the layer effective saturations and then solves the thermal balance equation to determine the layer temperatures. The accuracy of the solutions are checked, and if they do not pass a user-determined threshold then the time step is reduced, and the mass and energy balances are recalculated at the smaller time step. If the solutions are accurate enough then the mass and energy balances are calculated at the next time step. 30

2.2. Model Description Chapter 2. Snow Energy Budget Model Mass Balance Solution The change in layer mass and density is determined by a number of stages: 1. Precipitation is added to the snowpack. A new layer is created if snow has just started to fall. If precipitation is continuing, snow is added to the top layer at stage 8, until a maximum layer width is reached. New layers are created to accommodate further snowfall. Rain is added to the existing top node as flow from a theoretical layer above. 2. Densification by compaction, as described in section 2.2.4, is determined by a change in nodal width. Each layer retains its sample of snow (constant mass) during compaction, and the new densities are calculated from the change in layer width. The rate of compaction is calculated here, but actual changes are made to the node in stage (8). 3. Solar absorption within the layers and the surface energy fluxes are estimated, as described in the upper boundary condition for energy (section 2.2.2). This energy governs the amount of melt within a time step. 4. Layers that have sufficient liquid water for flow to begin, or where liquid water is already flowing are identified from the effective saturation. A residual amount of liquid water is retained within the ice matrix, bonded to the surface of ice grains, and at ice grain boundaries, by tension (temperature dependent). Here, the amount of residual water is assumed to be constant (set by the model user). Once the amount of liquid water in the layer exceeds the residual requirement (Se > si), liquid water flow begins. 31

2.2. Model Description Chapter 2. Snow Energy Budget Model 5. Liquid water fluxes are determined from equations 2.21 (liquid water constituent) and 2.25, where the effective saturation is the independent, unknown variable. Each region of flow contains one or more layers, with a corresponding equation for each layer that depends on its effective saturation and the effective saturation of the node above (flow across the layer upper boundary). The tridiagonal matrix solution method can be used to solve a set of equations of the form: AA XXJ + A Xi + A -1-B = (2.35) where X is the independent variable. The effective saturation is linearised to reduce the flow equations to the form of equation 2.35. Layers at the edge of the flow boundaries are subject to rapid changes in flow. Here, the effective saturation is estimated from the solution to the cubic equation, Se,est (assuming the liquid water flow from the node above has not changed from the previous time step). The true effective saturation is linearised as: S (Se,est) Se (2.36) Away from the flow boundaries, changes in liquid water flow are much slower. The effective saturation is linearised using a Taylor series expansion as: Se -2 ( -At)3 3( -At)2 S (2.37) The layer effective saturations are subsequently calculated using the tridiagonal matrix method of Press et al. (1992), and the liquid water fluxes are determined. 6. Vapour fluxes are determined from the diffusion equation (2.23). 32

2.2. Model Description Chapter 2. Snow Energy Budget Model 7. The increase in the average grain size by metamorphism is estimated for each layer from the algorithms in section 2.2.4 (equations 2.32 to 2.34). 8. Finally, mass changes are made to the layers. Precipitation is added to the top node, if required. The mass of water is recalculated as a result of liquid water and vapour flow. The nodal width and density are determined after compaction processes. All nodal properties (e.g. porosity, liquid water content) are recalculated to reflect changes to the mass. Energy Balance Solution The nodal temperatures are calculated by the following procedure: 1. Nodal thermal properties are calculated from the new density. Effective thermal conductivity is estimated as [Jordan (1991)]: ke = ka + (7.75 x 10-5 yw + 1.105 x 10-6 ) (ki - ka) (2.38) The effective thermal conductivity includes heat transfer by vapour transport, but does not take into account the proportion of liquid water in the layer. Glendinning (1997) suggested that this algorithm (equation 2.38) overestimated the effective thermal conductivity. 2. Liquid water content of a node is normally calculated by this model as a function of temperature. However, changes in the liquid water content of snow are rapid at temperatures close to the melting point. A relatively small uncertainty in the layer temperature would result in a large uncertainty in the liquid water content. For nodes close to the melting point, the energy balance equation (2.28) is redefined in 33

2.2. Model Description Chapter 2. Snow Energy Budget Model terms of the melt variable to increase model accuracy. A full description of the 'melt zone' is given by Jordan (1991). 3. The energy balance equation 2.28 has only one non-linear term: the thermal radiation emitted from the snow surface. Although a Taylor's expansion should be used to linearise the T4 component, this was not included in SNOWCAN because of the form of the thermal component of the canopy model, which was coupled to this snow model in chapter 4. The time steps were assumed to be small enough to minimise this error, and the thermal emission of radiation was estimated from the snow temperature at the previous time step. This simplification both degrades the model and decreases model efficiency. 4. The set of energy balance equations are solved using the tridiagonal matrix method. The temperature of nodes close to melting (stage 2) are retrieved using a NewtonRaphson iterative scheme. 5. Surface fluxes and snow heat fluxes are determined from the new layer temperatures. 6. Liquid water content is estimated from the nodal temperatures and other bulk properties are recalculated. 2.2.6 Time Step The need for an appropriate model time step was introduced within the description of numerical methods (section 2.2.5). SNOWCAN has an adaptive time step to optimise model stability, accuracy and efficiency, following the work of Jordan (1991). A CrankNicolson scheme is used for discretization of time variation. The size of the next time step 34

2. 2. Model Description Chapter 2. Snow Energy Budget Model is governed by the accuracy of the solutions at the current time step. If the accuracy is exceeds a threshold defined by the user, the time step is increased. If the solutions do not reach the required accuracy, the time step is reduced and the calculations at that time step are repeated (time is not advanced to the next step until the desired accuracy has been reached). 2.2.7 Model Outputs To summarise the model capabilities, the model simulates the following properties, which can be validated against surface, above surface and within snowpack measurements: * Surface fluxes: - Reflected solar radiation (or surface albedo) - Upwelling thermal radiation - Sensible heat - Latent heat * Snow profiles: - Temperature - Density (or integrated SWE) - Grain size - Liquid water content - Snow layer widths 35

2.3. SNOWCAN and SNTHERM Differences Chapter 2. Snow Energy Budget Model 2.3 SNOWCAN and SNTHERM Differences This section highlights areas where the snow model used in this project differs in style from SNTHERM as a result of the change in computer language. Replaced algorithms and numerical solution methods are described. Appendix A describes inconsistencies in the SNTHERM code. The complete code of SNOWCAN is included on a CD, which accompanies this thesis. 2.3.1 Model Algorithms Compaction Section 2.2.4 described the compaction algorithms, which differ from the SNTHERM parameterisation. The compaction due to melt implemented in SNTHERM is neglected in SNOWCAN, because a 1i factor makes the model unstable at low bulk ice densities. The change in the metamorphism compaction, and the parameterisation of both compaction processes implemented in the simulations presented in this thesis probably compensate for the compaction clue to melt, providing the time step is small enough. As described on page 28, the rate of settling compaction is greater than that used in SNTHERM, but in contrast, the rate of compaction by overburden pressure is decreased in SNOWCAN through an increase in the viscosity coefficient, 770o. In the simulations presented in this thesis, o0 - 9 x 107 N s m-2, compared to literature values of 1.7 x 105 < 70O < 1.9 x 106 N s m-2 [Kojima (1955)], 710 = 3.6 x 106 N s m-2 [Anderson (1976)] and 07o = 1.0 x 106 N s m-2 [Brun et al. (1989)]. However, the compaction parameterisation could be a source of error within the model. A comparison between SNTHERM and SNOWCAN compaction algorithms is 36

2.3. SNOWCAN and SNTHERM Differences Chapter 2. Snow Energy Budget Model presented in figure 2.2. The initial compaction curve is sharper, but is shallower towards the end, which offers an improved representation of field observations of snow depth in chapters 5 and 7. 0.2................................................................................................................................................. 018 0.16 5 0.14 - \ SNTHERM (100 kg m"-3) \ SNOWCAN (lO kg mJ-3) | \ --— SNTHERM (220 kg m^-3) 0.12. 0.2- \SNOWCAN (220 kg m^-3) 0.1 0.08 0.06 -.... e. e... e e - I - X.... i. i. i I i. i. I. i. I 0 10 20 30 40 50 60 70 80 90 100 110 120 Time (days) Figure 2.2: Comparison of SNTHERM and SNOWCAN compaction algorithm Precipitation In SNOWCAN, precipitation density and liquid water content are determined from the wetbulb temperature, rather than the air temperature. A summary of the algorithm [from D. Marks, pers. comm., but calculated from dewpoint temperature] is presented in table 2.1 below. 37

2.3. SNOWCAN and SNTHERM Differences Chapter 2. Snow Energy Budget Model Wetbulb temperature New snow density Liquid water content (K) (kg m-3) (%) < 268.15 75 0 268.15 - 270.15 100 0 270.15 - 271.65 150 0 271.65 - 272.65 175 0 272.65 - 273.15 200 25 273.15 - 273.65 250 75 > 273.65 RAIN 100 Table 2.1: Wetbulb temperature dependent new snow density Solar Zenith Angle The solar zenith angle influences snow albedo, ratio of direct to diffuse solar radiation and the transmission of direct solar radiation through the canopy. A new algorithm, which determines solar zenith angle from the field site latitude and time of day and year has been introduced in SNOWCAN [Jacobson (1999)]. 2.3.2 Numerical Solutions In SNOWCAN, a new algorithm has been used for the solution of the tridiagonal matrix [Press et al. (1992)] to improve the accuracy of calculated water saturation and temperature profiles. Determination of the estimated liquid water saturation at the wetting boundaries has also been improved by a different cubic solution algorithm [Press et al. (1990)]. This required re-definition of variable 'ap' in function FSSEST, used in the estimation of liquid water saturation at the wetting boundaries (equation 105b, [Jordan (1991)]) as: ap = -(cl * dzdt * sso - unold - upo/2d0 + bmelts) /dum (2.39) so that 'ap' and 'p' could be used as the coefficients in the cubic solution method. 38

2.3. SNOWCAN and SNTHERM Differences Chapter 2. Snow Energy Budget Model 2.3.3 Change in Style Fortran 90 (F90) allows definition of TYPEs, which are used extensively in this model. A variable is described by its TYPE e.g. INTEGER, REAL or CHARACTER. Definition of a new TYPE allows the variable to be composed of, say 3 REAL numbers and 2 INTEGERS. In this way, a NODE type was defined, that carried all variables related to a particular node, and an array of nodes was built. This contrasts with the SNTHERM approach, where variables are retained in either node or layer5 arrays. In a similar manner, TYPEs were defined for the meteorological data6, external model parameters (user-defined)7, flags and calculated fluxes. Several advantages are gained by use of TYPE definitions: 1. Very few variables are passed between subroutines, since only the TYPE arrays are required, and this minimises the risk of misnamed variables with variable transfer. 2. Old nodal values are easily assigned 3. Nodes are easily deleted, combined and subdivided However, one disadvantage occurs when nodes are created (by precipitation). Since layer properties are stored for each node, the new node requires these to be assigned. This problem can be overcome by initially creating a node identical to the one below it. Great care must be taken to update all nodal specific properties to ensure only layer remnants of the node below remain. 5Layers in SNTHERM are defined by the material it contains e.g. snow or soil. The layer is then divided into nodes. In this model, only the individual nodes are defined, so that any variables common to the 'layer' are stored in the nodal array. Although this duplicates information, the reduction in structure complexity is beneficial. 6All data were read in at the start of the model. 7In contrast, all internal model parameters are grouped within one module, and all are labelled with capital letters to aid identification as internal parameters. 39

2.3. SNOWCANV and SNTHERM Differences Chapter 2. Snow Energy Budget Model This change in program style is particularly important, not only for the data handling advantages highlighted above, but also because it changes the way the model works. 40

Chapter 3 Snow Model Sensitivity 3.1 Introduction A description of the physically-based snow model, SNOWCAN, was given in chapter 2. The model is constrained by a number of parameters, boundary conditions and initial conditions. In this chapter, the effect of uncertainty in the parameter space on the rate of snowmelt has been investigated, and the most important parameters have been identified. The importance of sensitivity studies for physically-based snow models was highlighted by Tuteja and Cunnane (1997), who examined the sensitivity of the UCGVDSM model. Model accuracy is also affected by errors in the measurement of data used to initialise and drive the model. The sensitivity of the model to errors in the forcing data and initialisation profiles are described in sections 3.5 and 3.6 respectively. The model was initialised with an isothermal, homogenous snowpack (described in section 3.2), so that the model sensitivity was not affected by initial snowpack stratification. Uniform meteorological data (described in section 3.3) were used to drive the 41

3. 1. Introduction Chapter 3. Snow Model Sensitivity model, to maintain a linear rate of change within the snowpack (as far as possible). Two error functions were defined in order to quantify the model sensitivity to a set of parameters, Pi, from model outputs, Vp, such as snow depth and surface energy balance components. In this chapter, the error functions are defined relative to the model output, Vsd, where the model has the standard parameterisation given in table 3.1, and is initialised and driven by data presented in sections 3.2 and 3.3. Vp also applies to uncertainties in input data. In chapter 7, the error functions are defined relative to observations, Vobs, rather than a simulation under standard model parameterisation (Vsd). The error functions are given below: Fm (Pi) Pi - Vsdj (3.1) =1(Vpi-Vd)i F (Pi) ( E - Vd) (3.2) i=i (Vsd)i Equation 3.1 gives the magnitude of the sensitivity and equation 3.2 gives the general sign of the sensitivity. A negative value of Fs given by equation 3.2 indicates snow depth is less than that of the standard simulation - either because the compaction is increased, or because the change in parameter produces faster snow melt. It is expected that for the most part of this sensitivity study, changes in model parameters will act only to increase or decrease the rate of melt, rather than produce a mixed response. Under these circumstances, the two error functions Fm and Fs will be equal in magnitude. For the simulation of a natural snowpack (chapters 5 and 7), the response of the model to changes in the parameter space may not be linear. The error function Fm is regarded as a more accurate indicator of the model error, and this function is used in chapters 5 and 7 to quantify the model accuracy. 42

3.2. Synthetic Snowpack Description -Chapter 3. Snow Model Sensitivity 3.2 Synthetic Snowpack Description The snow cover changes more rapidly during the ablation period than the accumulation period, therefore this sensitivity analysis focussed on a warm snowpack that melted during the simulations. An initial snowpack temperature of 272K was chosen, so that a small period of warming was studied, followed by the active melt period. The lowest soil layer was initialised at 273.15K, since the Dirichlet boundary condition constrains the temperature of the lowest layer to be constant, and a temperature lower than the melting point would cool the snowpack. The partial water density of all layers was initialised at 250 kg m-3, and the average snow grain diameters were initialised at 1mm. The input files used in this simulation are included on the CD that accompanies this thesis. i:;;itial tnltempeat p rd al l; -; - deniritieInai rneal plriled -.. *.-.- 02 -.-~ * —,=^ — i ' - o ano S *;;P a 0 o.................................... 0.05 a -J 0.05 * 0.0'50r 'i^:^ * X 'if i^ * as^ ^ ~ f -:.-.-. -'.:-A.-.-. -*- i if- i'iaii| < j; 0.05-w L i iL1 0 l: -00 5- W -- 271.6 272 2725 27 248 24- 2-50 251 2,52 0. 0.9 1 1.1 Figure 3.1: Initialisation profiles for sensitivity analysis Figure 3.1 shows the initial snowpack profiles used in this simulation. The total snow depth at the start of the simulation was 0.2m (the approximate average depth of the 43

3.3. Meteorological Data Input Chapter 3. Snow Model Sensitivity real snowpacks studied in this thesis). 3.3 Meteorological Data Input Three different meteorological scenarios were used to drive the model. As explained in section 3.1, the meteorological data used to drive the model were constant, so that the snowpack melted at an approximately constant rate, and the melt period was examined as changes are much more rapid during the ablation period. For most simulations, the normal meteorological input data (referred to as 'N' in table 3.1, which summarises the sensitivity analysis results) were as follows: * Downwelling longwave radiation (350 W m-2) * Air temperature (280K) * Windspeed (1 m s-1) * Relative humidity (50%). A low windspeed was used in order to minimise turbulent transfer (which is expected beneath a canopy). Solar radiation was set to zero to avoid diurnal cycle effects, particularly on the snow albedo. Precipitation was set to zero to allow the snowpack to melt without mass influx.1 Parameters that affected the solar radiation input to the snowpack were tested under similar meteorological conditions to those stated above, with an additional solar radiation input of 200 W m-2. This is not strictly a constant input, since the model 'Although mass can be added to the snowpack by condensation from the atmosphere, this influence is much smaller than precipitation accumulation. 44

3.3. Meteorological Data Input Chapter 3. Snow Model Sensitivity assumes zero solar radiation when the sun is below the horizon. This meteorological scenario is denoted as E in the sensitivity analysis results. The final set of meteorological conditions were used to test the model sensitivity during precipitation events. This scenario, denoted by *, is similar to scenario X, but with a constant air temperature of 274K, to ensure that precipitation was in solid form. A constant precipitation rate of 0.0001 m hr-1 water equivalent was input to the model, except during the model spin-up time; no precipitation was input during the first two timesteps. 0,.2 0,2.,,_ 0.15 N 0.1 L' 0.05 0 ' " ".....'..., "...."" 'I" I"","'" l'" I"lll"ll I' ","," 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 Time (days) Figure 3.2: Simulated depletion of snow cover with standard model parameterisation The melt profiles of the snowpack described in section 3.2, under the meteorological conditions stated above (N, I and *), are shown in figure 3.2. The set of parameters used in this simulation are similar to the set of parameters that have been used in the 45

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity simulations of chapter 7 and in the model included on the CD that accompanies this thesis. The melt profiles shown in figure 3.2 are therefore regarded as the 'standard', to which all other simulations are compared using the sensitivity functions (equations 3.1 and 3.2). 3.4 Model Sensitivity to Parameterisation Model parameters were varied by ~10% relative to the standard parameterisation used in chapter 7, and the error functions defined in equations 3.1 and 3.2 were calculated from the snow depth depletion profiles. Table 3.1 summarises the error functions for each parameter in descending order of magnitude. In the following sections, the reasons for the model sensitivity to the parameterisation are investigated. Parameters that affect the surface energy balance are examined in section 3.4.1 and those that affect the internal snow physics are studied in section 3.4.2. The sensitivity to differences between study sites is approached in section 3.4.3 and model errors arising from uncertainty in physical constants are investigated in section 3.4.4. Sections 3.4.5 and 3.4.6 examine the model sensitivity to the discretisation scheme and user-defined accuracy thresholds respectively. 3.4.1 Surface Energy Balance Parameters Thermal Radiation Thermal radiation, solar radiation and turbulent transfer of energy at the surface were all affected by parameters in the sensitivity study. The largest variation in snow depth, shown in Table 3.1 occurred through sensitivity to the snow emissivity, where a -10% variation from the standard value (0.99) resulted in an error function (Fm) of up to 0.15. Decreasing 46

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Parameter [Psd] Met. F+ Fm_ F+ FCs [0.99]t N 0.15 -0.15 Es [0.99]t ' 0.12 0.12 Ycrit [150] * 0.095 0.085 -0.095 0.085 Latitude [53.91634] ' 0.090 0.091 0.090 -0.091 EHO [2.0] N 0.080 0.063 -0.080 0.063 amax [0.8] ' 0.00070 0.048 0.00067 -0.048 Ycrit [150] N 0.040 0.030 -0.040 0.030 EEO [2.0] N 0.0060 0.013 0.0060 -0.013 Pvo,sat [6.112]t N 0.0067 0.012 0.0066 -0.012 Al,s [100]t N 0.0083 0.0034 0.0082 -0.0034 AZmax [0.016667]t N 0.0067 0.00070 -0.0067 -0.00012 Terr [0.05] N 0.0062 0.0029 -0.0062 0.0028 Gwet [4 x 10-12]t N 0.00054 0.0044 -0.00050 -0.0044 Tatm [0.07] ' 0.0044 0.0041 -0.0044 0.0041 k [0.4]t N 0.0044 0.00041 -0.0044 -0.00031 ki [2.29]t N 0.00093 0.0040 0.00092 -0.0039 AZzin [0.002]t N 0.0018 0.0040 0.0018 -0.0040 Gwet [4 x 10-12]t ' 0.0023 0.0038 0.0023 0.0038 Qtz [0.4] N 0.0034 0.0012 -0.0032 0.0010 Gdry [5 x 10-7]t ' 0.0023 0.00018 0.0023 0.000083 3rnir 0.0013 0.0020 -0.0013 0.0020 rs [25]t 'P 0.0017 0.00075 0.0017 -0.00070 ka [0.023]t N 0.00022 0.0013 -0.00015 0.0013 zo,s [0.005] N 0.0013 0.00074 -0.0013 0.00072 r7o [9.00 x 107] * 0.00081 0.0012 -0.00080 0.0012 Patm [945.12] N 0.0012 0.00089 -0.0012 0.00084 si [0.04] N 0.00074 0.0011 -0.00062 0.0011 Pi [0.001787]t N 0.00026 0.0011 -0.000084 0.0010 CDN/CHN [1.0] N 0.00091 0.00057 0.00090 -0.00050 Ca [l000]t N 0.00088 0.00023 -0.00083 -0.000096 Ci,o [2117]t N 0.00000078 0.00070 -0.00000078 0.00052 Deo [0.00009]t N 0.00052 0.00038 0.00050 -0.00029 7r0 [9.00 x 107] N 0.0000089 0.0000096 0.0000079 -0.0000096 CDN/CEN [0.7] N 0.0000018 0.0000018 0.0000018 -0.0000018 Gdry [5 x 10 —7] N 0.00000037 0.00000037 -0.00000035 0.00000035 Azincr [0.04]t * 0 0 0 0 serr [0.001] N 0 0 0 0 Table 3.1: Error functions from model sensitivity to parameterisation, where superscript (+) and (-) denote sensitivity to a + and - 10% change in parameter. t denotes that the parameters are internal. The meterological conditions applied to the model are described in section 3.3. 47

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity the snow emissivity resulted in an increase in the time taken to melt the pack, since less incident thermal radiation was absorbed at the surface. This effect was partially offset by a decrease in emission of thermal energy from the snowpack, but the model is still very sensitive to this parameter. Fortunately, the snow has similar thermal characteristics to a blackbody, and thus the emissivity is estimated to vary from 0.99 by only <3%. In addition, the thermal energy contribution to the surface energy balance in this study was high, therefore the model is expected to be much less sensitive to the snow emissivity under observed meteorological conditions. Solar Radiation Solar radiation absorption within the snowpack is governed by the snow albedo, which is snow grain size and solar zenith angle dependent, as described in section 2.2.2. The model allows an upper limit on the snow albedo, am(tx, which is defined by the user. An upper limit on the snow albedo increases flexibility for the user at present, but will be removed once the simulation of snow albedo by SNOWCAN has been validated. Chapter 6 discusses surface reflectance measurements that will be used in the future to validate the snow albedo algorithms incorporated in SNOWCAN. The model appears to be sensitive to the snow albedo upper limit. This is not surprising for a -10% change in snow albedo, since this increases the energy input to the snowcover, and the rate of snowmelt is faster. However, the sensitivity to a +10% increase in snow albedo is unexpected, since the albedo is limited by grain size and solar zenith angle, and under the meteorological conditions applied in this case, no fresh snow was added to the snowpack. 48

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity The model sensitivity to the albedo limit lies in the timing of the combination of layers within the model. The model determines the albedo of the top layer only, and the albedo of the underlying layers are not calculated until they are exposed to the surface. In addition, the snow albedo can only be determined during the day, when the solar zenith angle is positive. When layers are combined, the adaptive time step is very small. The snow albedo is determined at the end of the short time step, and the error caused by the momentarily high albedo is small. However, if the layers are combined during the night, the time step at the start of the day may be much longer, and then a high albedo may have a larger impact on the model performance. Under these simulated meteorologial conditions, solar radiation is set at 200 W m-2 as soon as the sun rises above the horizon, and so the model is sensitive to the albedo limit. Under observed meteorological conditions, though, the solar flux is very small as the sun rises, and so the model sensitivity to the albedo limit may be expected to be less during simulations driven by observed meteorological data. Turbulent Transfer of Energy Turbulent transfer of energy at the surface uses five parameters that are defined by the user: windless exchange coefficients for sensible and latent heat, roughness length and ratios of sensible and latent heat transfer coefficients to the eddy transfer coefficient. Zero wind transfer is governed by the windless exchange coefficients2, EEO and EHO, as described in equations 2.17 and 2.18. Table 3.1 shows that under the meteorological conditions applied, the model sensitivity to a ~10% change in the windless exchange 2The current recommendation for the windless exchange coefficients used in SNTHERM is EHO = 1 and EEO = 0 [Jordan et al. (1999)]. 49

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity coefficients is among the highest for the parameters studied. In addition, Jordan (1991) suggests a 25% error on these empirically derived coefficients. At higher wind speeds, the windless exchange coefficients become less important than the windy exchange coefficients, EH and EE, but beneath a forest, where wind speeds are low, the windless exchange coefficients are expected to be important parameters. Turbulent transfer during windy conditions depends not only on the wind speed, but also on the atmospheric stability. Objects protruding from the mean surface level provide frictional drag on the wind. The effect of surface features on the atmospheric stability is demonstrated in figure 3.3 [from Morris (1989)]. A laminar sub-layer may form above a smooth surface, whereas a rough surface extends into the turbulent flow layer. Figure 3.3: Structure of the lower layer of the atmosphere for steady horizontal flow over a uniformly rough snow surface (from Morris (1989)) The stability of the atmosphere is determined from the wind speed in terms of 50

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity the bulk Richardson number3, which is defined as: 2 g ZT (Tair - Tsurf) (33) Ri r1_ -27 -r - (3'3) W (Tair + Tsurf) A neutral atmosphere is characterised by Ri = 0, a stable atmosphere, where warm air lies above cold air, is characterised by 0 < Ri < 0.2, and an unstable atmosphere, where cold air lies above warm air, is determined as Ri > 0. Turbulent transfer of energy during windy conditions is determined from the user-defined parameters roughness length and ratios of sensible and latent heat transfer coefficients to the eddy transfer coefficient. A neutral atmosphere is initially assumed and correction factors are subsequently applied for non-neutral atmospheres. The stability functions for neutral and stable atmospheres are assumed to be equal to 1, whereas for an unstable atmosphere, an iterative technique is required to estimate the stability functions from the Richardson number and estimated bulk transfer coefficients, which is fully documented in Jordan (1992). The bulk transfer coefficients used to determine the surface turbulent transfer of energy, described in chapter 2, are determined as: CD = CDN( 1 I CDN ) CE = CEN CD 1-[ C keJ (3.4) C DN L V(7DN^ where the coefficients at neutral stability are defined as: k2 CDN = [n (Z/ZO)]2 3The Richardson number was determined from the absolute measurement height for the sensitivity study. This was changed to the relative height for the experimental results. 51

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity kx/CDN CHN k CD In (Z/ZT) CEN = C ) (3.5) in (Z/ZQ) (3.6) Table 3.1 indicates a low sensitivity of the model to the roughness length. However, a low wind speed was chosen for the simulation, which minimises turbulent transfer. It is therefore expected that the model sensitivity to the roughness length would be increased for simulations driven by observed meteorological data. Table 3.1 also shows that the model is more sensitive to the roughness length than to the neutral atmosphere ratios of bulk transfer coefficients of sensible and latent heat to the drag coefficient. It is interesting to note that increasing the roughness length increases the rate of snow melt through increased turbulence, whereas an increase in the neutral bulk transfer ratios has the opposite effect. 3.4.2 Internal Snow Physics Parameters The simulated energy input to the snowpack is determined by the boundary conditions applied to the model, but the transfer of mass and energy within the snowpack is governed by the internal snow physics. In this section, the empirical parameterisation of the internal snow physics is investigated, and the impact of parameter uncertainty on compaction, snow grain growth, liquid water content and transfer, vapour diffusion and radiation transfer is discussed. 52

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Snow Compaction As described in section 2.2.4, two parameters are used to describe the compaction of the snowpack. The first is a critical density, 7Ycrit, and beyond this density limit, the compaction that occurs as a result of destructive metamorphism is reduced with increasing snow density. The second parameter, 710, limits the rate of compaction from the pressure of the snow above. A full sensitivity analysis of the model to these parameters was not carried out, because of the computational expense, but should be carried out in the future as part of an investigation into a linked compaction and grain growth model, which is discussed in chapter 9. The sensitivity during the melt period is examined here, but the model is expected to be more sensitive during the accumulation period. Despite the limitation of the sensitivity study, the critical density for compaction appears to be one of the most important parameters governing model performance. This is true for both cases studied: constant precipitation and precipitation-free meteorological conditions. An increase in the density limit produces a corresponding increase in the rate of densification. This not only affects the snow depth directly, but also governs the rate of heat transfer within the pack, which is density-dependent. The model uses the following algorithm to determine the thermal conductivity of the snow: ks = ka + (7.75 x 10-5yw + 1.105 x 10-6y) (ki - ka) (3.7) although the model combines heat transport through conduction and vapour diffusion, and both processes are represented with a single effective thermal conductivity: ke = ks + LviDeCkT (3.8) The effect of increasing the density threshold for compaction not only increases 53

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity the density of the snowpack, but also the rate of heat transfer within the pack. The effect of increasing the thermal conductivity of the snowpack is examined in section 3.4.4. Increasing the viscosity coefficient has conflicting effects on the rate of snowmelt for the different meteorological conditions studied, as demonstrated in table 3.1. For the non-precipitation case, an increase in the viscosity coefficient decreases the compaction rate through the 1 / 770o dependence. Therefore, the depth of the snowpack is increased relative to the standard parameterisation and the rate of snowmelt is slower. The model is more sensitive to this parameter as mass is added to the pack, and the pressure is increased on the snow below, but an increase in 770o produces an apparent increase in the rate of snow depth depletion. This counter-intuitive model behaviour results from the density dependence of r7, which increases exponentially with increasing density, as described by equation 2.31. The decrease in densification with an increase in 70O described above conflicts with the increase in densification for lower density snow. Whilst this effect is minimised for a single sample of snow (because the density would increase to a point where the 1 / 70 effect dominates at higher densities), the model behaviour here arises because fresh, low density snow is continually added to the snowpack. The representation of the compaction processes in the snow model is crucial to model performance. Here, the empirical algorithms used require two user-defined parameters to calibrate the compaction routine. The snow model would benefit greatly from a more detailed physically-based compaction approach, which incorporates the metamorphism of the snow crystals explicitly with the pressure of the overlying snowpack. 54

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Snow Grain Growth The mechanism of grain growth within a snowpack was described in section 2.2.4, and is represented empirically within the model. Two parameters describe the rate of grain growth: one for dry snow and one for wet snow. Under the conditions applied, the model is most sensitive to grain growth parameter Gwet, because the snowpack is melting for most of the simulation. The model shows some sensitivity to the grain growth parameter Gdry because the snowpack is not initially melting, and some grain growth occurs under low liquid water content whilst the snowpack is warming. As described in section 2, the snow grain size affects the reflection and absorption of radiation and the rate of gravitational liquid water flow. Increasing the rate of grain growth under conditions of no solar radiation increases the rate of melt. Liquid water flow is enhanced by grain growth, and heat transfer within the snowpack is enhanced by the increased mass transfer. However, removal of liquid water from a layer also affects the layer density and thermal conductivity. Enhanced grain growth with the addition of solar radiation acts to reduce the rate of melt. Less solar radiation is reflected from larger grain sizes, therefore more energy is input into the snow cover. The bulk extinction coefficient at visible wavelengths is reduced with increased grain size (equation 2.13), therefore a smaller proportion of energy is absorbed in each layer. Whilst the overall energy input at the surface is increased, more radiation may be absorbed by the soil, rather than within the snowpack, which in turn reduces the rate of melt. The model sensitivity is expected to differ for an optically semi-infinite snowpack. 55

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Liquid Water Content The liquid water content of the snow is determined from its temperature, using a semiempirical function (given below for snow): 1_ 1 7w 1 + A (273.15 -T)2 (3.9) The form of this function is shown in figure 3.4, which was taken from Jordan (1991). Figure 3.4: Model liquid water fraction temperature dependence (from Jordan (1991)) The slope of the curve in figure 3.4 is sharpened by increasing parameter Al in equation 3.9. This has the effect that the snowpack persists for longer under melting conditions, as the water remains frozen at higher temperatures. Table 3.1 indicates that the freezing curve parameter, Al, is more important than many user-defined parameters. The standard value was arbitrarily chosen [Jordan (1991)], and therefore a reasonable degree of uncertainty is associated with this parameter (> 10%). However, the model sensitivity to the liquid water content parameterisation will be greatly reduced outside the melting period. 56

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Liquid Water Transport The irreducible saturation describes the amount of liquid water retained by a substance after drainage. Using this parameter, the effective saturation, which is the independent variable in the mass balance solution, is defined as: s- Si Se = s i (3.10) 1 - Si An increase in the irreducible saturation means that the snowpack retains more liquid water. The reduction of heat transfer from the snow by liquid water flow results in a faster rate of snowmelt. Sensitivity of the model to a ~10% change in this userdefined parameter is relatively low under the conditions applied, but the uncertainty in this parameter could be greater than 100% [Kattelmann (1986), Colbeck (1974), Jordan (1991)]. The gravitational drainage of water is also dependent on the viscosity of water (equation 2.25), although the sensitivity analysis showed that the model was less responsive to changes in this parameter than to changes in the irreducible saturation. An increase in the water viscosity reduces the rate of liquid water flow, which increases the rate of melt. The expected parameter variation is low (< 2%), but non-negligible since the viscosity of water is pressure-dependent. The model is expected to be less sensitive to this parameter over a seasonal simulation. Vapour Diffusion Heat and mass transport by vapour diffusion was introduced in chapter 2. The mechanism of vapour transport is assumed to occur by sublimation from one side of a crystal, vapour diffusion across the air void in the the direction of the downwards temperature 57

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity and concentration gradient and condensation onto the next crystal. This process occurs for all crystals, and whilst the process is a diffusion process, the same molecules are not transported through the pack, hence this process is represented with an effective diffusion coefficient. The effective diffusion coefficient in equation 2.23 is estimated from measurements at 1000mb pressure and 0~C, and from pressure and temperature relations as: (1000 T \6 'De ( )273.15) (3.11) P' — a 2 ( T"3-15 The model was found to be relatively insensitive to the diffusion coefficient De,0. Dorsey (1940) reported a diffusion coefficient of water vapour in air of 0.000022 m-2 at 0~C and atmospheric pressure, and a temperature dependence of T1'75 from the compilation of Boynton and Brattain (1929). The effective diffusion coefficient is expected to be higher than the value of Boynton and Brattain (1929), since the diffusion path length is effectively much shorter for the crystal-to-crystal vapour transport. Radiation Transfer Radiation extinction within the snowpack, as described in section 2.2.2, is governed by the bulk extinction coefficients for visible and near-infrared radiation, which are calculated from the snowpack properties according to equation 2.13 (f/3is) and are input by the user (/3nir). Definition of the bulk extinction coefficient by the user allows the division of the solar radiation into a proportion of energy that is absorbed at the surface, and the remainder that penetrates into the pack and is absorbed lower down in the pack. This approach was taken by Koh and Jordan (1995), which was based on the spectral model of Brandt and Warren (1993). 58

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity An increase in the bulk extinction coefficient for near-infrared radiation increases the amount of solar radiation absorbed at the surface, but also decreases the energy absorption lower in the snowpack. Radiation that penetrates the snowpack to the soil is absorbed at the soil surface, and hence reduces the energy available for snowmelt. The penetration of radiation is discussed further in this section. In addition to this effect, heat is more easily transported downwards through the snow during the melting period, because conduction, diffusion and liquid water flow all transport heat downwards, whereas only conduction and diffusion mechanisms carry heat upwards (for a non-isothermal snowpack). This will be discussed in greater depth in section 3.4.4. Therefore, an increased proportion of energy absorbed at the surface enhances the rate ol melt, relative to the case where the same energy is absorbed lower in the snowpack, or in the soil. Absorption of radiation that penetrates the snowpack is grain size and density dependent. Although the transmission of radiation through the snowpack is generally modelled as an exponential function, this model assumes a maximum mass-dependent depth, at which all remaining radiation is absorbed. The model was found to be less sensitive to this parameter than to the bulk extinction coefficient for near-infrared radiation. An increase in the maximum radiation penetration depth was found to decrease the rate of snow depth depletion. This effect could be due to a similar mechanism to the bulk extinction parameter effect, as the increase in parameter lowers the depth at which a proportion of energy is absorbed, and the energy transport through the snowpack is affected. As the depth is decreased to below the penetration depth, the radiation is absorbed by the soil, hence the model is not sensitive to this parameter at low snow depths. 59

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity The sensitivity of the model to the maximum radiation penetration depth should be examined carefully in the context of a seasonal snowpack. Exponential extinction of radiation within the snowpack implies that some radiation should always reach the underlying soil. If this is the case, then measurements of surface reflection should be affected by the depth of snow. However, measurements taken during the course of this study, as explained in chapter 6, indicate that this is not the case, and in fact, the surface reflectance is only significantly by the depth of snow at very low snow depths (<-5cm). 3.4.3 Site-to-site Variation Site Location The model sensitivity to site location was found to be relatively high under the conditions studied. Latitude is used by the model to calculate the solar zenith angle. The increase in latitude has two immediate effects: the day is shorter, so less solar energy is input to the pack, and secondarily, the albedo is decreased, because of the increased solar zenith angle. Although these factors are contradictory, the decrease in energy input from decreased length of day overwhelms the increase in energy input from decreased albedo, and the snowpack melt rate is reduced. However, this effect is exacerbated by the meteorological conditions chosen. The model is expected to be much less sensitive under diurnal cycle conditions, and in addition, the uncertainty associated with measurement of site latitude is generally negligible. 60

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Underlying Soil In the model, two soil types are parameterised: clay and sand. A complete sensitivity analysis of all soil parameters was considered to be unnecessary for this study, but the sensitivity of the model to soil type was investigated, and is described below. The model sensitivity to the lowest soil layer temperature was also investigated, since it affects the heat transfer to the snowpack, and is described within the sensitivity to initial conditions section (3.6). A change in soil type from clay to sand resulted in sensitivity functions Fm = 0.033 and Fs == -0.033. The model is reasonably sensitive to soil type, and the depletion of the snowcover over a sandy soil was quicker than the depletion of the snowpack over clay soil. Since vapour diffusion is disallowed for soil, and the liquid water is artificially drained at the snow-soil interface, the main effect on soil is through the heat transfer process. The thermal conductivity for sandy soil is higher than for clay soil, therefore the heat transfer to the warming snowpack is more rapid for the sandy soil. The model appears to show similar sensitivity to the soil quartz content. This parameter is not usually measured in the field, and it is important to gauge the accuracy required for this parameter, which is suggested to range between 0 and 0.45 by Jordan (1996). An increase in the quartz content of the soil increases the soil thermal conductivity, which increases the rate of snow melt under these initial conditions. 61

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Meteorological Conditions The atmospheric optical depth is used by the model to estimate the clear sky direct / diffuse split of solar radiation, if measurements are unavailable, according to: II - i e-Tatm// o SW,dir - SW 'SW,diff - iSW (l- )e (3.12) The model shows relatively high sensitivity to the optical depth, and an increase in the optical depth was found to decrease snow depth relative to the standard simulation. The proportion of direct to diffuse radiation affects the snow albedo. The albedo for direct solar radiation is greater than the albedo for diffuse radiation, as demonstrated by equations 2.8 to 2.11. An increase in the atmospheric depth increases the proportion of diffuse radiation, which lowers the snow albedo and enhances the rate of melt of the snowcover. The atmospheric optical depth is not simple to measure at a field site, because remote sensing equipment and data processing are required. Measurements cannot be made during cloudy conditions, which affect the direct and diffuse proportions, and the atmospheric optical depth can vary over time, particularly if the field site is close to high pollution areas. This sensitivity study highlights the benefits of diffuse solar radiation measurements as model forcing data. Determination of the model sensitivity to the atmospheric pressure is particularly important as this parameter is assumed to be constant during a model simulation. The atmospheric pressure is used in the calculation of vapour diffusion (the air within the snow is assumed to be at atmospheric pressure); air density, which affects the turbulent transfer (62

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity of sensible heat; and the temperature of precipitation. An increase in the average atmospheric pressure decreases the diffusive heat transport (equation 3.11), but increases the sensible heat transfer to the surface (equation 2.18). The additional energy input from turbulent transfer with increased atmospheric pressure enhances the rate of melt of the snow, which was observed in this sensitivity study. The atmospheric pressure can vary by 10%, but the model was found to be less sensitive to this parameter than to many internal parameters for the conditions studied. However, a low wind speed was used in this simulation, and the model sensitivity is expected to increase with greater rates of turbulent transfer. 3.4.4 Physical Constants Saturation Vapour Pressure Water vapour pressure affects the latent heat transfer at the surface and diffusion processes within the snowpack, for example, a low density of water molecules relative to the maximum (saturation vapour pressure) enhances net sublimation. The upper limit of vapour pressure is the saturation vapour pressure, which is temperature dependent. Vapour pressure is often defined in terms of the saturation vapour pressure at a fixed temperature, which is 0~C here, as: Pv = frh X f{T} Pvosat (3.13) Pvk,sat The model demonstrated high sensitivity to the saturation vapour pressure, and an increase in this parameter was found to decrease the rate of snow melt. From equation 2.17, it can be seen that an increase in the saturation vapour pressure decreases the rate of latent heat transfer to the surface. Diffusion within the snowpack (equation 2.23) is 63

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity increased as the slope of the saturation vapour density vs. temperature curve is increased: CkT Pvk sat l osat (Lvk - ) e-( RT )e(TO) (3.14) =T RT2 - RT 2 but the model sensitivity resulted from the energy decrease at the surface. Fortunately, the saturation vapour pressure at 0~C is well known, but the algorithm used to determine the saturation vapour pressure at a different temperature (from Buck (1981)) has not been tested in this study. Von Karman Coefficient The model shows moderate sensitivity to the von Karman constant, which is crucial in the determination of turbulent transfer. A theoretical increase in the von Karman constant increases the drag coefficient (equation 3.6). This also increases the turbulent transfer of latent and sensible heat, since the ratios of sensible and latent heat to the drag coefficient at neutral stability are constant in the model. The rate of snow melt in the model increases with the additional energy input to the snow. Values for the von Karman constant vary in the literature by approximately 5% from the value used here [Jacobson (1999) suggests a value between 0.35 and 0.43], and the sensitivity of the model to the uncertainty in this parameter is expected to increase with higher wind speeds. However, the model error that results from the uncertainty in the von Karman constant is insignificant in comparison to the uncertainty in the algorithms used to represent the turbulent transfer. 64

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Thermal Conductivity A decrease in the thermal conductivity of ice reduces the rate at which energy is transported away from the surface, hence melt occurs earlier at the surface, but the snow pack as a whole is cooler. However, liquid water flow from the surface through the pack warms the snow quickly, and the rate of melt for the snowcover is increased. The thermal conductivity of ice is temperature dependent, and increases by -12% over a - 20~C change. This is not important for a melting snowpack, since the temperature is constant and an appropriate value is used for the thermal conductivity of ice. However, for a cold snowpack, the thermal conductivity may be different, although the model is expected to be less sensitive to the thermal conductivity, since the height at which the snow starts to melt is not relevant. A much smaller increase in the thermal conductivity of the snow cover, such as a 10% increase in the thermal conductivity of air, has the opposite effect to larger increases in the thermal conductivity. Heat is distributed more evenly through the pack, and the total snowpack melts more quickly. The thermal conductivity of air is pressure, and therefore altitude dependent, and the resultant uncertainty could be as high as 15%. However, a larger source of error could lie in the liquid water content of the snowpack, since the thermal conductivity of liquid water (0.561 W m-1 K-1 at 0~C) is not explicitly considered in the snow thermal conductivity algorithm. 65

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity Heat Capacity The model is relatively insensitive to 10% changes in the heat capacity of both ice and air. The heat capacity of ice at 0~C is used to determine the temperature limits of the meltzone, and hence affects the temperature solution accuracy. The heat capacity of ice is temperature-dependent, and is estimated for a snow layer from the following algorithm: ci = 7.8 T - 13.3 (3.15) The transfer of sensible heat at the surface is dependent on the heat capacity of air, but this parameter is approximately constant for the range of expected temperatures. 3.4.5 Model Discretisation In this section, the model sensitivity to the layer size limits was investigated. The model should theoretically become more accurate as the layers reduce in size, and representation of snowpack profiles are smoother in the model. Maximum Layer Width A limit is placed on the maximum size of the top two nodes to improve model accuracy, since snowpack changes are most rapid at the surface as a result of the boundary conditions. In this sensitivity study, the width limit for the top two nodes was changed by ~10%, for both layers at the same time. The standard value for the top layer is given in table 3.1, and the standard value for the layer beneath the top is 0.033m. This width choice was based on SNTHERM numerical studies and was found to be the most efficient width (R. Jordan, pers. comm.). 66

3.4. Model Sensitivity to ParameterisationC Chapter 3. Snow Model Sensitivity A change in the top layer width has two immediate effects: the energy absorbed in a thin layer produces a higher rate of change of layer temperature than for a thicker node, and the change in temperature profile affects heat transfer through the snowpack; and the timing of layer division and hence solution accuracy is affected. It is not easy to separate out these effects to predict or explain the model response to a change in the maximum layer width. However, the upper limit on layer width is a balance between thinner layers, which give a smoother response to the surface energy input, and thicker layers, which decrease computational expense. The surface layer width cannot be very thin, since the assumption of energy absorption at the surface would cause model instability. Minimum Layer Width The model is less sensitive to the minimum layer width under the conditions studied than to the maximum layer width. The effect of nodal width on heat transport was discussed above, but the minimum layer width governs the timing of layer combination rather than division. Precipitation Layer Width The model showed no sensitivity to the maximum width for a fresh layer of precipitation for the conditions studied. This resulted from insufficient precipitation to overcome the mass loss from melt, and therefore fill the layer to the maximum width. The sensitivity of the model to the precipitation layer width under true accumulation conditions is unknown, and cannot be easily tested by this sensitivity approach. 67

3.4. Model Sensitivity to Parameterisation Chapter 3. Snow Model Sensitivity 3.4.6 Model Accuracy Criteria At the end of a time step cycle, the accuracy of the mass balance solution and energy balance solution is tested. If the solution accuracies are within the threshold levels set by the user, the solutions are accepted and the time step is advanced. If either solution accuracy does not satisfy the accuracy criterion, the time step is shortened and the calculation for that time step is repeated. Under the conditions tested, the temperature accuracy threshold was exceeded, but the saturation threshold was not. An increase in the temperature error limit resulted in an increase in the snow melt rate, whereas a decrease in the temperature error limit had the opposite effect. This points to an overestimation in the modelled snow temperatures, and suggests that the model may benefit from a reduced accuracy threshold, despite the increase in computational expense. Linearisation of the temperature dependent thermal emission of radiation may aid the model performance in this area. The time step is increased after a number of successive accurate time step calculations, which is set by the user. The model sensitivity to this parameter was tested, by +1, from the standard number (2 time steps). An increase in the number of time step was found to increase the rate of melt, as Fm = 0.0068 and F8 = -0.0068, whereas a decrease in the number of time steps had the opposite effect, with Fm = 0.0046 and F8 = 0.0046. If the model increases the time step after one good calculation, but subsequently is not accurate enough, the time step is lowered to below that of the first time step, and the solution is more accurate. Whilst the time step may be increased more frequently than for a greater number of time steps criteria, the average time step may be lower, as the model repeats more time step calculations. A lower tinie step, and more accurate calcula 68

3.5. Sensitivity to Boundary Conditions Chapter 3. Snow Model Sensitivity tions reduce the overestimation of temperature by the model. However, the calculation is generally robust with the time step. An overview of the order of calculation is presented in figure 4.3 in chapter 4, which also includes the vegetation components. 3.5 Sensitivity to Boundary Conditions The model sensitivity to errors in the measurement of forcing data were investigated under the meteorological conditions 1, described in section 3.3. The error functions (equations 3.1 and 3.2) are summarised in table 3.2 below. The uncertainty in thermal and solar radiation can be high, as demonstrated by Field et al. (1992), who compared different 4-component radiometers from several manufacturers at the same site and time, and found measurement differences of up to 80 W m-2, and systematic-root-mean-square instrumental differences of 15 to 35 W m-2. Here the error estimate of 20 W m-2 is more conservative, since Eppley radiometers were used in the collection of fieldwork data described in chapter 6, which are considered to be among the highest quality radiometers available. Forcing data [error] Met. F+ F | F+ FT I I 1t [20 W m-2] 4> 0.38 0.72 -0.38 0.72 Tair [1 K] 4 0.11 0.088 -0.11 0.088 Jrh [10%] ( 0.079 0.059 -0.079 0.059 I1 [20 W m-2] ~ 0.071 0.053 -0.071 0.053 w [0.5 m s-1] ( 0.014 0.0083 -0.014 0.0083 Table 3.2: Error functions from model sensitivity to errors in forcing data measurement. Several points are of interest to note. An increase in all forcing data examined increased the energy input to the snowpack, and therefore the rate of melt. Turbulent 69

3.6. Sensitivity to Initial Conditions Chapter 3. Snow Model Sensitivity transfer of latent and sensible heat are opposite in magnitude, but sensible heat dominates here, because of the values used for the user-defined coefficients CDN/CHN and CDN/CEN. An error in the thermal energy input to the snowpack has the greatest effect, for the meteorological conditions studied. This effect is greater than the same error applied to the solar radiation. This is because the model sets the solar radiation to be zero when the sun goes below the horizon. Therefore the measurement error for solar radiation is applied for less than half of the day. The handling of radiation within the snowpack will also differentiate the model sensitivity to the measurement error of thermal and solar radiation. An error in thermal radiation measurement primarily affects the energy balance at the surface, whereas the same error in solar radiation immediately affects the energy absorbed within the snowpack. An increase in the fractional humidity decreases the amount of evaporation / sublimation from the surface. However, since this process cools the snowpack, an increase in the fractional humidity therefore suppresses energy -transfer from the pack, and the rate of melt increases. A comparison of tables 3.1 and 3.2 indicates that the model is working reasonably, since the model is generally much more sensitive to uncertainties in the forcing data than to uncertainties in the parameter space. 3.6 Sensitivity to Initial Conditions The sensitivity of the model to measurement errors in the initial snow profiles of temperature, density and grain size was tested to determine the effect of uncertainty in the initial conditions. The error functions for the snow depth depletion are given in table 3.3. 70

3.6. Sensitivity to Initial Conditions Chapter 3. Snow Model Sensitivity Initial profile [error] Met. F+ F F+ F-?yw [50 kg m-3] [ 0.44 0.43 0.44 -0.43 T1 [1K] 4 0.29 0.095 -0.29 0.095 Ts [1K] 4D 0.085 0.014 -0.085 0.014 (d [0.5 mm] 3 ( 0.077 0.081 -0.077 0.081 -(w,soil [50 kg m-3] I 0.010 0.018 [ 0.0095 -0.015 Table 3.3: Error functions from model sensitivity to errors in forcing data measurement. Table 3.3 shows that of the initial snowpack profiles studied, the model is by far the most sensitive to the initial density profile. As discussed in section 3.4.4, an increase in the density of the snowpack increases the snowpack effective thermal conductivity, which reduces the rate of melt. An underestimation in the initial snowpack density profile may be compensated in the model through the refinement of the compaction algorithm. On the other hand, an overestimation of the initial snow density may not be corrected easily, since densification is largely irreversible. The model is much more sensitive to the temperature of the lowest soil node than to the temperature profile of the remaining soil and snow layers. The Dirichlet lower boundary condition for energy, as described in section 2.2.2, assumes that the temperature of the lowest soil layer is constant. Therefore a measurement, or averaging error, in the temperature of this layer is propagated throughout the model simulation. The temperature profile of a snowpack generally follows the temperature of the air, with a lag time related to the thermal conductivity of the snow and temperature height within the snowpack. The initial temperature profile is therefore lost relatively quickly - within a few hours - and the model sensitivity to the initial temperature profile is low. Under these initial and meteorological conditions, however, the model sensitivity is higher, since the time taken to melt the snowpack also includes the time taken for the snowpack 71

3.7. Conclusions Chapter 3. Snow Model Sensitivity to ripen. The model sensitivity to the initial grain size is important because, similarly to the density profile, an overestimation in the initial snow grain size cannot be undone, whereas an underestimation can be compensated through increased grain growth. Here, the model is relatively insensitive to the initial grain size profile, but as the snow grain size governs solar energy absorption within the snow, the model is expected to be more sensitive to the initial grain size profile for a longer time series simulation. 3.7 Conclusions In this chapter, the model sensitivity to parameter space, uncertainties in forcing data and the initial conditions was examined in order to highlight the most important areas where uncertainties in parameters and input data affect the model performance. An ablation period was examined in this sensitivity study, because the snowpack changes the most rapidly during this period. The model was generally less sensitive to uncertainties in the parameter space than to measurement errors in the forcing data. A number of parameters have been selected for a sensitivity study at an open site, using observations to initialise, drive and validate the model. Some parameters with high model sensitivity during this study, as shown in table 3.1, were not selected for further investigation because the model was expected to be less sensitive over a seasonal simulation, or because the parameter uncertainty is much less than 10%. Further sensitivity studies are reported in chapter 7, where the model sensitivity to the following parameters and initial conditions has been investigated for the accumula 72

3.7. Conclusions Chapter 3. Snow Model Sensitivity tion and ablation of a real snowpack, driven by measured data: * Windless exchange coefficients EHO and EEO. * Compaction critical density Ycrit. * Compaction viscosity coefficient 710. * Bulk extinction coefficient for near-infrared radiation 3nir* Soil quartz content. * Atmospheric optical depth, Tatm. * Initial snow density profile. * Initial lowest soil layer temperature (Dirichlet boundary condition). As described in chapter 1, the wind speed beneath a forest canopy is greatly reduced relative to the wind speed at an open site. This sensitivity study highlighted the windless exchange coefficients as important parameters that affect model perfomance during the ablation period, under the low wind speeds studied. The uncertainties associated with the windless exchange parameters are also high. The densification algorithms affect the heat and mass transport within the snowpack, and therefore the parameterisation of the snowpack compaction will be examined further. The model sensitivity to the initial snowpack density profile is high, and the measurement error warrants further sensitivity investigations (chapter 7). Unknown parameters such as the bulk extinction coefficient for near-infrared radiation, atmospheric optical depth and particularly the soil quartz content all have a high degree of uncertainty. The model sensitivity to these parameters over a seasonal 73

3. 7. Conclusions Chapter 3. Snow Model Sensitivity simulation (chapter 7) will indicate the suitability of parameter choice, and the accuracy required for these parameters. The lower boundary energy condition requires a constant soil temperature, which is given by the initial condition data. Investigation of the model sensitivity to this condition is highly beneficial. If an average temperature over the season is required, then the operational use of this model is limited. However, if the accuracy required for this temperature is low, or the temperature can be estimated from typical early or late season temperatures, then the application of this model can be more widespread. The model tendency to overestimate the snow temperature was highlighted in section 3.4.6. This indicates a possible benefit in the reduction of the temperature accuracy threshold in the model, despite the increase in computational expense. This should be looked at in greater detail once the linearisation of thermal radiation has been addressed, which was discussed in chapter 2. This chapter examined the sensitivity of the model to its parameters under one set of idealised meteorological conditions, and therefore offers a snapshot of the model reponse to changes in the parameter space. A much more comprehensive sensitivity study is planned for the future, and directions for this future work are given in chapter 9. The sensitivity study presented in this chapter indicated the suitability of SNOWCAN to be incorporated within a Snow-SVAT. The next chapter describes the vegetation radiation model that has been coupled to this snow model, the method of approach and other forest effects that are represented within the model, or are considered to be negligible for this study. 74

Chapter 4 Vegetation Effects 4.1 Introduction As described in section 1.2, the forest cover has a number of effects on the snow cover. These include radiation shielding and emittance, reduction of precipitation reaching the snow surface by intercepted snow sublimation and reduction of snow albedo by leaf litter within the snowpack. A canopy optical and thermal radiation model was chosen and coupled to the physically based snow model that was described in chapter 2, and other forest processes were examined, to form the Snow-SVAT (snow-soil-vegetation-atmosphere transfer) scheme, 'SNOWCAN'. 75

4.2. Canopy Radiation Model Chapter 4. Vegetation Effects 4.2 Canopy Radiation Model 4.2.1 Introduction The attenuation of radiation by a canopy has been previously modelled with different degrees of complexity, ranging from a simple Beer's Law approach, which ignores scattering, to geometrical optical techniques. Longwave and solar radiation are often treated independently. Beer's Law assumes exponential decay of solar radiation. Wigmosta et al. (1994) adopt this approach, using a constant attenuation coefficient. In reality, attenuation of radiation varies with path length through the canopy and is therefore a strong function of solar zenith angle. Some models (Yamazaki and Kondo (1992), Hellstrom (1999)) have addressed this issue by using a solar zenith angle dependent attenuation coefficient. However, over-simplified treatment of radiation, such as the use of Beer's Law, is unjustified where there are more physically reasonable models available. Geometrical models [Li and Strahler (1986)] approximate tree elements as simple geometrical structures (cones and ellipsoids). The bidirectional reflectance is calculated under all illumination and viewing angles by assuming that radiation is not transmitted through the canopy crowns. This approach has subsequently been improved in hybrid models, where radiative transfer theory is used to model the transmission of radiation within canopy crowns [Li et al. (1995)]. However, there are a number of problems with geometrical-optical techniques: * Initialisation of the model is complex as the distribution of tree stems must be known * Tree parameters required are difficult and laborious to measure in the field (a large 76

4.2. Canopy Radiation Model Chapter 4. Vegetation Effects number of samples must be taken to adequately represent the canopy) * The technique is computationally expensive A vegetation model should include a physically based optical and thermal radiation model, that is also reasonably simple, thus minimizing the number of parameters required to initialise and drive the model. The radiation model "RM" [Pearson et al. (1999)] approximates the canopy as a series of discrete, randomly orientated scatterers. Radiation incident on the canopy is attenuated by striking the canopy elements. The radiation is either absorbed or scattered iostropically with equal fluxes above and below the canopy. By assuming a proportion of the incident radiation is intercepted (dependent on the canopy structure), the physical depth of the canopy can be ignored. The proportion of intercepted radiation which is scattered is governed by the leaf area index and by the single-scattering albedo at shortwave wavelengths and the emissivity for thermal radiation. Thermal radiation is treated in the same way as shortwave radiation. The model allows infinite scattering between the canopy understory and the ground beneath in order to retain conservation of energy. It is driven by incident solar and thermal radiation. Diffuse radiation from the sky is assumed to be isotropic and direct radiation is approximated as a delta-function beam in the direction of the sun. Thermal radiation is emitted from the atmosphere, canopy and snow surface. RM was chosen to form the canopy radiation component of the SNOWSVAT because of its physical basis and elegant mathematics. This model has many advantages over other models: * RM models both optical and thermal radiation 77

4.2. Canopy Radiation Model Chapter 4. Vegetation Effects * Initialisation of the model is simple * The canopy parameters required by RM are simple to measure and can be determined in the field or readily obtained using remote techniques * Computational efficiency is relatively high This model has already been coupled to a soil model as part of a SVAT (soil vegetation atmosphere transfer scheme) and was validated with a high degree of success [Pearson et al. (1999)]. 4.2.2 Physical Principles The vegetation model RM assumes the canopy to be composed of a series of spherically distributed, randomly orientated discrete scattering components. Radiative transfer theory is used to define transmission and reflectance coefficients, which enable the total upwelling and downwelling fluxes to be determined above and below the canopy. This is shown schematically in figure 4.1, where Io is the above-canopy downwelling flux, T is the transmission coefficient, R is the canopy reflectance coefficient and B is the lower boundary reflectance - the snow albedo. In order to satisfy conservation of energy, multiple reflections of radiation between the snow surface and canopy underside are accounted for. The total sub-canopy downwelling and upwelling fluxes are given by equations 4.1 and 4.2. Io =- TIo+RB T Io+(RB)2 TIo-+ =[1- RB] T Io (4.1) ItTo = [1 -RB]-1 B TIo (4.2) 78

4.2. Canopy Radiation Model Chapter 4. Vegetation Effects CANOPY `'I"~t t K1t v t:1 v (R71 Ut 1i Figure 4.1: RM model canopy fluxes The reflection and transmission coefficients R and T are defined in equations 4.3 and 4.4 for shortwave (Rsw, Tsw) and thermal (RLWV, TLW) radiation. A full derivation can be found in Pearson et al. (1999). The coefficients are included here for reference. Solar radiation: RSw = - a [1 - 2E3(Tv) TSw = 2E3(T) + Rsw (4.3) Thermal radiation: RLW = -(1- E) [1- 2E3(7T)] TLW = 2E3(T) + RLW (4.4) where a, is the vegetation albedo, E3(.) is the third exponential integral (E3(x) = e —Xt f~ F dt, x > 0), ev is the vegetation emissivity and Tv is the vegetation optical depth. Tv is determined from the leaf area index AL (equation 4.5) for the canopy approximation 79

4.3. Other Canopy Effects Chapter 4. Vegetation Effects as spherically distributed, randomly orientated discrete scatterers. A full justification for this relationship is detailed in Pearson et al. (1999). T AL (4.5) 2 4.3 Other Canopy Effects In addition to the scattering and absorption of radiation by the canopy, other forest processes affect the mass and energy balance of the snowpack. The representation of these processes by SNOWCAN is described below. 4.3.1 Intercepted Snow Interception of snow by the canopy significantly affects the distribution and mass of snow accumulated on the ground, as shown in figure 4.2 [taken from Sturm (1992)]. Therefore, interception processes are immensely important for the mass and energy balance of a winter forest. The mechanisms by which snow is intercepted on, and subsequently lost from the canopy, were introduced in section 1.2, and are described in greater detail below. Canopy Loading The loading mechanism by which snow accumulates on a canopy is described by Schmidt (1991). At the beginning of accumulation the ice crystals generally fall though the canopy. A very thin layer accumulates along the sterns of each branch. As snow collects at the base of the needles, the interception efficiency increases. Bridges form between the needles and subsequently between branches and larger scale structures, until the maximum loading capability is achieved. 80

4.3. Other Canopy Effects Chapter 4. Vegetation Effects Figure 4.2: [From Sturm (1992)] Snow deposition patterns beneath the trees of the taiga. (a) leaning birch or aspen with enlarged snow cone; (b) vertical aspen or birch with snow cone; (c) white spruce with large tree well; (d) black spruce, covered with snow, forming chimney; (e) black spruce with small tree well Over 50% of cumulative snowfall can be stored as interception within boreal, montane and subalpine forests in mid-winter [Pomeroy and Gray (1995)]. Measurements of cumulative intercepted snow in a boreal pine canopy ranged between 30 and 44% throughout the winter period in two separate years [Pomeroy et al. (1998)]. This was consistent with measurements of boreal pine throughfall [Pomeroy and Dion (1996)], which was found to be 68% of the above-canopy snowfall. Intercepted snow may remain on the canopy for up to a month in cold environments [Hedstrom and Pomeroy (1997)], and leaves the canopy by unloading, or by sublimation. 81

4.3. Other Canopy Effects Chapter 4. Vegetation Effects Canopy Unloading Snow intercepted on a canopy may be removed by changes in wind, or air temperature [Sturm (1992)]. Radiative heating may cause partial melt, which causes water to drip from the canopy and wind disturbance may allow the snow mass to slip to the surface [Hellstr6m (1999)]. Most of the intercepted snow in forests of the taiga is removed from the canopy by unloading mechanisms [Sturm (1992)], but in boreal, montane and subalpine forests, sublimation losses from snow intercepted on the canopy are high [Pomeroy and Gray (1995)]. Mass Loss by Sublimation Table 4.1 summarises field measurements of sublimation losses from intercepted snow, as a percentage of the total above-canopy snowfall, or snowfall at an open site. Source Loss (%) Canopy type Pomeroy and Schmidt (1993) >30 boreal Pomeroy and Gray (1995) 25-45 boreal, montane and subalpine Harding and Pomeroy (1996) 30 boreal Hardy et al. (1997b) 35 boreal black spruce Pomeroy et al. (1998) 30-32 boreal mature pine Pomeroy et al. (1998) 38-45 boreal black spruce Table 4.1: Sublimation losses (% of above-canopy snowfall) No physically-based model exists that determines the mass and energy balance of the intercepted snow. However, the intercepted mass and its subsequent sublimation have been simulated using a process-based model [Pomeroy et al. (1998)], which showed good agreement with observations. In that model, canopy unloading was parameterised with a constant unloading coefficient that was determined from field measurements in the 82

4.3. Other Canopy Effects Chapter 4. Vegetation Effects boreal forest [Hedstrom and Pomeroy (1998)]. The process-based algorithms developed by Pomeroy et al. (1998) require additional parameters to those used within SNOWCAN, and the extra parameters cannot be determined using remote sensing methods. An interception model was not included within SNOWCAN for these reasons. Based on field measurements detailed in the literature (see table 4.1), SNOWCAN assumes 30% of precipitation is intercepted by the canopy, and that all snow intercepted on the canopy is lost to the atmosphere by sublimation. Albedo of Canopy with Intercepted Snow Snow has a much higher albedo than the canopy, therefore the presence of intercepted snow in the canopy affects the radiation balance of the canopy and underlying snowpack. Measurements of above-canopy albedo showed little change when snow was present on the canopy [Leonard and Eschner (1968), Harding and Porneroy (1996)]. Harding and Pomeroy (1996) concluded that solar radiation was trapped beneath the canopy. Snow intercepted on the canopy is assumed here to have two effects on the radiation balance: the penetration of solar radiation through the canopy is decreased (leaf area index is effectively increased) and that less radiation is absorbed by the canopy (canopy albedo is increased). These assumptions have opposite effects on the radiation balance. Given the measurements made in the field and the assumptions made above, the effect of intercepted snow on the canopy radiation balance is assumed to be negligible, and is not modelled by SNOWCAN. 83

4.3. Other Canopy Effects Chapter 4. Vegetation Effects 4.3.2 Sub-Canopy Wind Speed Turbulent transfer of energy to the snowpack surface is dependent on the wind speed, which is greatly reduced beneath the canopy. Since SNOWCAN is forced with meteorological data measured outside of the canopy, the shielding effect of the canopy must be considered. Link and Marks (1999) suggested an algorithm to estimate the sub-canopy wind speed from the above-canopy wind speed as: Wf= 5 (4.6) where wf is the sub-canopy wind speed and Wo is the open wind speed. The minimum sub-canopy wind speed is limited to 0.2 m s-1 in order to maintain model stability. This algorithm was implemented in SNOWCAN to estimate the sub-canopy wind speed from the wind speed measured at an open site, so that turbulent energy exchange beneath the forest was not overestimated. No attempt has been made to relate the stability of the atmosphere to the presence of the canopy in SNOWCAN, except through the reduction in wind speed. The algorithms used are the same as those over an open snowpack. Turbulent transfer over open areas is now well understood [Jordan et al. (1999)], but much more research is required beneath forest canopies. Windpumping through the forest can transfer a large amount of energy, and this process is not simulated in SNOWCAN. 4.3.3 Forest Debris Snow beneath a vegetation canopy becomes increasingly dirty as debris falls from the canopy and accumulates on the snow surface. This, in turn, affects the snow albedo and 84

4.3S. Other Canopy Effects Chapter 4. Vegetation Effects resulting energy balance. The amount of surface litter varies with time and precipitation, as recent snowfall covers the surface litter and hence the albedo increases to that of fresh snow. As the snow melts, litter buried under later precipitation re-emerges with the diminishing snowpack. The effect of litter on the snow albedo must be taken into account in order to represent the forest-snow energy balance in a physically realistic manner. The method for estimating the distribution of litter within the snowpack, and its effect on the snow surface albedo are descibed below. Fractional Coverage of Litter Hardy et al. (2000) modelled the surface and within-pack fractional coverage of litter and its effect on the surface albedo. The algorithm (equation 4.7) presented by Hardy et al. (2000) to determine the percent coverage of surface and within snowpack litter has been implemented in SNOWCAN. A = 1 -- (1 - A)x (4.7) where A is the fractional litter coverage (rn2m-2), A is the daily litter rate (m2m-2day-) and X is the number of days since the last precipitation event (for surface litter) or number of days since the first precipitation event (for total snowpack litter). Effect of Litter on Snow Albedo Forest debris has a lower albedo than the albedo of pure snow, therefore the distribution of litter within the snowpack lowers the surface albedo. Here, a simple algorithm is proposed, which lowers the surface albedo as a function of surface litter coverage: as = A av + (1 - A) a, (4.8) 85

4.3. Other Canopy Effects Chapter 4. Vegetation Effects where as ' is the degraded surface albedo. This does not take into account the distribution of litter within the snowpack, multiple scattering of radiation within the snowpack, the change in spectral reflectance and therefore broadband albedo, or the aging, and therefore browning of the litter, which changes the albedo from that of the canopy, av. 4.3.4 Thin Snowpack Albedo Solar radiation at visible wavelengths penetrates the snowpack to a depth of up to ~20cm SWE [Wiscombe and Warren (1980)], although this is grain size and wavelength dependent. For a thin snowpack containing less than 20cm SWE, solar radiation may penetrate to the soil beneath the snowpack. A smaller proportion of radiation is scattered from the soil surface than would be scattered by snow, hence, the overall albedo for a thin snowpack may be lowered through absorption of radiation within the soil. Snow albedo measurements of thin snowpacks above a black surface by Giddings and LaChapelle (1961), showed no change in albedo at wavelength A = 590nm, for snow depths greater than 2cm (not SWE). Warren and Wiscombe (1980) suggested that the asymptotic value for albedo was lowered by impurities within the pack, and Wiscombe and Warren (1980) calculated the semi-infinite depth for clean snow as 4 times that measured by Giddings and LaChapelle (1961). Spectral albedo measurements of finite snowpacks show characteristic flattening of the albedo at visible wavelengths [Wiscombe and Warren (1980)]. Measurements made in the field of the spectral reflectance of shallow snowpacks, as described in section 6.3.2 estimated the semi-infinite depth to be ~5cm. From the experimental work of Giddings 86

4.4 Canopy Effects Not Modelled Chapter 4. Vegetation Effects and LaChapelle (1961) and fieldwork that has been described in chapter 6, SNOWCAN assumes a negligible change in surface albedo of a thin snowpack as a result of radiation absorption by the underlying soil. 4.4 Canopy Effects Not Modelled In this study, several canopy effects were not modelled. The canopy absorbs solar and thermal radiation and can therefore have a different temperature to that of the air. An energy balance for the canopy could be used to estimate the temperature of the canopy. The alternative used in this study was to input the measured canopy temperature where available, or to assume the canopy was at the same temperature as the air temperature. The thermal radiation from vegetation stems was not simulated. However, measurement of leaf area index typically includes canopy branch and stem structure. The tree trunk is effectively incorporated in the canopy parameterisation, from which the emission of thermal radiation from the tree elements is calculated. Therefore it is considered unnecessary to explicitly model the thermal emission from the canopy trunk. Vegetation and forest floor debris such as fallen branches, buried within the snowpack were not considered. These would affect the radiation balance (through increased absorption) and would have cavity effects, impacting transport processes. The forest floor at the main study site described in chapter 6 and modelled in chapter 7, was free from vegetation and large-scale debris, therefore this effect is assumed to be negligible at this point, but would need to be considered for a stand-scale simulation or greater. 87

4. 5. Model Implementation Chapter 4. Vegetation Effects 4.5 Model Implementation The method used to couple the physically based snow model and vegetation canopy model is shown in figure 4.3. After interpolation of met data for a particular time step, the vegetation model is used to determine sub-canopy optical and thermal radiation fluxes. The mass balance of the snow cover is computed first, followed by the energy balance. Finally the snow albedo is computed as a result of the change in snowpack properties. The snow albedo is a very important feedback parameter as it both governs and is governed by the energy balance. Therefore the process of coupling the two models is important. This differs from the pre-processed forcing data approach taken by many models to adjust the radiation balance for the presence of the canopy, which cannot simulate this feedback mechanism. 4.6 Summary Here, a summary of the forest processes modelled by SNOWCAN is presented, and the assumptions made are reiterated. Canopy processes modelled: 1. Optical and thermal radiation. The sub-canopy down- and up-welling optical and thermal radiation fluxes are modelled by the coupling of the vegetation model RM with the physically based snow model SNOWCAN. Multiple scattering beneath the canopy is modelled. 2. Radiative heating of canopy. An energy balance of the canopy is not included 88

4.6. Summary time rtep.a d repe at Chapter 4. Vegetation Effects Initialisation FiMet Datak Figure 4.3: SNOWCAN model structure 89

4.6. Summary Chapter 4. Vegetation Effects at present, but the canopy temperature may be included as forcing data, where available, or assumed to equal the air temperature otherwise. 3. Interception loss. 30% of the precipitation measured above the canopy (or at an open site) was assumed to be intercepted by the canopy. All intercepted snow was assumed to be lost by sublimation. 4. Wind speed. Shielding by the forest was assumed to reduce the wind speed beneath the canopy by 80%. 5. Forest debris. Litter was assumed to fall at a constant rate, and the effect of litter on surface albedo was represented by a simple, linear algorithm1. Canopy processes assumed to be rnegligible, or not considered: 1. Intercepted snow albedo. Canopy albedo is assumed to be unchanged as a result of the snow intercepted on it. 2. Atmospheric stability. Algorithms to determine turbulent transfer over open snowpacks were also applied to forested snowpacks. 3. Thin snowpack albedo. The albedo change as the snowpack thins was assumed to be negligible. 4. Thermal emission from stems. This effect was assumed to be incorporated by canopy parameterisation and the thermal emission model for the general canopy. Although this process is included in SNOWCAN, it has not been implemented for any of the simulations presented in this thesis, because insufficient data were available, or because the amount of litter deposition was low. 90

4. 6. Summary Chapter 4. Vegetation Effects 5. Large contaminants. Vegetation and large-scale debris within the snowpack were assumed to be absent. Three extra parameters are required to represent the forest cover in SNOWCAN: Leaf area index, canopy albedo and canopy emissivity. Canopy temperature forcing data are desirable, but are estimated from the above-canopy air temperature if the data are unavailable. 91

Chapter 5 Snow-SVAT Simulation I: Boreal Forest 5.1 Introduction SNOWCAN was tested with data taken during the BOREAS campaign, which was an international, multi-disciplinary experiment designed to investigate processes affecting global change on a wide range of scales [Sellers et al. (1997)]. The experiment was conducted over 5 years within the Canadian boreal forest, which is snow covered for 6-8 months of the year. Data collected during BOREAS are now freely available via the internet, and can be obtained from the Oak Ridge National Laboratory website (http:\\www-eosdis.ornl.gov). Data collected by several groups of investigators were used for the canopy parameterisation, snowpack initialisation and to drive SNOWCAN, as described in section 5.3. The snow cover beneath an old jack pine canopy site, described in section 5.2, was simulated for the period 1 Jan - Apr 20 1994, and the simulated snow depth was compared 92

5.2. Site Description Chapter 5. Snow-SVAT Simulation I: Boreal Forest with snow depth measured in a small canopy gap. The comparison between simulated and observed snow depths is presented in section 5.4. 5.2 Site Description A number of sites with different canopy characteristics (species, LAI, height) were studied during BOREAS. The old jack pine site (OJP) within the Southern Study Area (SSA) was selected to test SNOWCAN during the winter period 1993-1994. This data set was chosen for two reasons: the range of data collected during this period, and the previous research that has already focussed on this year and site [Hardy et al. (1997a)]. The location of several BOREAS experimental study sites, including the old jack pine site (OJP) studied in this chapter, are shown in figure 5.1. 5.3 Model Inputs 5.3.1 Model Initialisation Detailed snow pit measurements were unavailable for model initialisation, therefore the snow temperature, density and grain size profiles were estimated from available data and typical observations. The snow depth at the start of the simulation was determined from the automatic depth gauge, and the snowpack was arbitrarily split into 6 layers, with thinner layers towards the snow surface. A logarithmic density profile was applied to the snowpack, to give an integrated density of 220 kg m-3, which was consistent with the measured integrated density of 231 kg m-3, taken during the standard snow course on 14th January 1994 by BOREAS group HYD04. Within the snowpack, the temperature 93

*0 0 LA. 0 C 0 E CD 0 C CO 0 a. C.), kiz -------- - t t.% 1 1, I; iii -A. P,. 4 t III * It IA i ~0 ci) 0 +-o 0 4c a) -oI 0 Cr1 0n 2t 0 z a) 4 0 4 0 a) C) 0 a) a) 0 -SC C) -o 0 S a) -e C 0 V Cl) - Ci) Q <*5 S 4a) A 4-D CZ,-.;;........... ---- ----- - 1-1-1..........-...... 11...

5.3. Model Inputs Chapter 5. Snow-SVAT Simulation I: Boreal Forest profile was assumed to be linear between the soil surface and air temperature. Temperatures within the soil were measured by BOREAS group HYD01 at 10cm, 20cm and 50cm from the surface. These measurements were used to define the soil layer boundaries, and the temperatures of the layer midpoints were determined by linear interpolation from the point measurements. The lowest soil layer temperature was set to 274K to satisfy the Dirichlet boundary condition. HYD01 determined the soil to be a clay type, and the partial water density of the soil was arbitrarily initialised to be 100 kg m-3. This value is between the residual water content of 3% and the saturated water content of 40%, which were both measured by HYD01. The sensitivity study presented in chapter 3 showed that the model was less sensitive to the soil water content than to uncertainties in other initial conditions, forcing data and many model parameters. Snow grain size information was unavailable, therefore commonly observed snow grain sizes of 1mm were used to initialise the snowpack. No data were available to describe the deposition of litter on the snow surface, therefore this process was not modelled during the BOREAS simulation. The soil and snowpack initialisation profiles are summarised in figure 5.2, and the initialisation file is included on the CD that accompanies this thesis. 5.3.2 Canopy Parameterisation SNOWCAN requires three parameters to represent the canopy: leaf area index, emissivity and albedo. Leaf area index was determined as 2.5 by Chen et al. (1997), at the end of the 1993 growing season from measurements taken with a Decagon Ceptometer1 (see 1http://www. deca.gon.com/ 95

5.3. Model Inputs Chapter 5. Snow-SVAT Simulation I: Boreal Forest Figure 5.2: Initialisation profiles for BOREAS simulations section 6.3.1 for further information about this instrument) and a LI-COR LAI-2000 2. Canopy single-scatter albedo was assumed to be similar to that of a fir canopy within the Reynolds Creek Experimental Watershed, ID. The single-scatter albedo (Cv = 0.25) was determined in the laboratory and has been described in greater detail in chapter 6. Vegetation thermal characteristics are generally similar to a blackbody, therefore a canopy emissivity of 0.977, which was given by Monteith and Unsworth (1990) for a poplar canopy., was used in this simulation. 5.3.3 Forcing Data Meteorological variables (above canopy temperature, relative humidity, wind speed and downwelling shortwave and longwave radiation) were measured by instruments mounted above the canopy on a tower. Precipitation was measured with a Belfort weighing pre2http://env.licor.com/products/LAI2000/2000.htm 96

5.4. Simulation of Snow Depth Chapter 5. Snow-SVAT Simulation I: Boreal Forest cipitation gauge mounted with an Alter-shield. This instrument was situated beneath the canopy 42.7m east-southeast of the tower. Snow depth was measured with an ultrasonic depth gauge situated in a small canopy gap 37m west-northwest of the tower. Cloud cover was assumed to be negligible throughout the duration of the simulation. Missing data were linearly interpolated for time spans of one hour or less. However, data were missing from 30 Mar 23:30 to 01 Apr 00:15. For this period, data were estimated from data collected at the old aspen site (approximately 21 km SW of OJP site). The meteorological variables at the two sites exhibited a high degree of correlation, typically with an r2 value of 0.9 (for the variables used). Measured precipitation data are unavailable from 31 Mar (Julian day 90) to the end of the season, therefore precipitation water equivalent was also estimated from data collected at the old aspen site. 5.4 Simulation of Snow Depth Snow depth in the canopy gap was simulated as demonstrated in figure 5.3, with an r2 correlation coefficient of 0.91. Figure 5.3 exhibits large increases in simulated snow depth not observed in the field, particularly on days 37, 49, 52, 57 and 79. These differences result from delayed throughfall of intercepted snow. Since the snow depth is measured beneath a small gap in the canopy, the measured snow depth increases concurrently with precipitation events. However, the snow depth sensor is situated beneath the canopy where snow may remain intercepted on the canopy until it is lost by sublimation or falls to the ground either by partial melt causing 'snow slip' [Hellstrom (1999)] or by the intercepted mass exceeding 97

5.4. Simulation of Snow Depth Chapter 5. Snow-SVAT Simulation I: Boreal Forest 0.4 0.35 0.3 E 0.25 r" s 0.2 I 0.15 0.1 0.05 0................................................................... + Measured I- Simulated 0 10 20 30 40 50 60 70 80 90 100 110 Julian Day Figure 5.3: Simulated snow depth beneath a jack pine boreal forest canopy the canopy storage capacity. Increases in the snow depth directly beneath the canopy lag behind snow depth increments in the snow depth beneath the canopy gap, and may be smaller in magnitude due to canopy losses. In order to illustrate the time lag experienced by the precipitation gauge, the large precipitation events observed on Julian days 37, 49, 52 and 79 were assumed to occur on clays 31, 48, 50 and 78 respectively. Figure 5.4 demonstrates the simulated snow depth, driven by reconstructed precipitation data. The mass added to the snowpack remained the same, as only the timing of precipitation events was altered. The large precipitation event simulated on day 57 is particularly interesting as the measured snow depth profile suggests a large precipitation event on day 50, followed by several days of either light precipitation or interception fall through the canopy. As yet, no physically-based model exists that represents gradual throughfall over a period of 98

5.4. Simulation of Snow Depth Chapter 5. Snow-SVAT Simulation I: Boreal Forest time for the input data, so this event was not altered in any way. By taking into account the mass input timing differences caused by snow interception in the vicinity of the precipitation sensor, figure 5.4 shows that the model predicted and measured snow depths exhibit good agreement, with a regression coefficient of 0.94. The regression coefficient may not be the most appropriate statistical test of agreement, because the results are not independent, but serves as a quick guide to how well the model is working. 0.45 W 0.4-..Measured E0.3 0.25 - 0.1 0.05 -0 10 20 30 40 50 60 70 80 90 100 110 Jutlian Day Figure 5.4: Simulated snow depth beneath a jack pine boreal forest canopy [reconstructed precipitation data] Differences between the measured and modelled snow depth could arise from a number of causes. The mass input to the snowpack beneath a forest is affected by the interception of snow on the canopy. Here, canopy throughfall was used to drive the model, rather than above-canopy precipitation (or precipitation measured at an open site), 99

5.4. Simulation of Snow Depth Chapter 5. Snow-SVAT Simulation I: Boreal Forest therefore the uncertainty in the mass input is decreased. However, the spatial variation of interception loss is unknown, and the mass input to the surface in the vicinity of the precipitation gauge may be different to the mass input in the gap between individual tree canopies, where the snow depth was measured in this experiment. The leaf area determined by the BOREAS team is representative for the site, but is probably less in a canopy gap, which would affect the radiative energy balance at the snow surface. The canopy single-scatter albedo and emissivity were both estimated from different canopy type measurements, and may be different for the boreal Jack Pine canopy. A sensitivity study of SNOWCAN to the canopy parameters is presented in chapter 7, where the effect of these parameters on the radiation balance is studied in greater detail. The final melt profile of the jack pine snowpack, days 98-106 in figure 5.4, is represented well by the model. This indicates that the assumptions regarding the albedo decrease from litter distributed throughout the snowpack, or radiation penetration to the underlying surface, compaction processes and the lower Dirichlet boundary condition, are probably reasonable assumptions. The benefits of the approach of this SVAT model can clearly be seen in figure 5.5, where the same input data (including throughfall precipitation data) were used to drive the model, but with all other canopy effects (i.e. radiative and wind speed reduction) removed from the model. This is for illustration purposes and is not necessarily physically realistic. The regression coefficient between measured and modelled snow depth is r2 = 0.29. 100

5.5. Conclusions Chapter 5. Snow-SVAT Simulation I: Boreal Forest 0.45 0.4 0.35 E 0.3 ~- 0.25 0~ c 0.2 015 0.15 0.1 0.05 0 0 10 20 30 40 50 60 Julian Day 70 80 90 100 110 Figure 5.5: Simulation of boreal snow depth with canopy effects removed 5.5 Conclusions This chapter describes how the physically-based Snow-SVAT was used to simulate the evolution of the snowpack beneath a boreal jack pine forest, using data taken during the BOREAS experiment. Simulated snow depth fitted observed snow depth, with a correlation coefficient, r2 = 0.94. However, validation of any Snow-SVAT against measurements of snow depth is difficult, since both the energy balance and the snow compaction, which is represented empirically in the model, have a large impact on the snow depth profile. Several parameters that affect the energy balance have a large uncertainty associated with them, such as the leaf area index and lower boundary soil temperature. The canopy measurements were not taken at the same location as the point measurement of validation, and neither were the mass forcing data measurements. 101

5.5. Conclusions Chapter 5. Snow-SVAT Simulation I: Boreal Forest More data are required to test SNOWCAN fully. Ideally, measurements of subcanopy radiation and snowpack physical properties are required to examine all components of this model. Appropriate validation data were not collected as part of the BOREAS campaign, therefore a new experiment was designed to collect a complete data set, which could be used to validate any Snow-SVAT. This experiment was carried out at a fir site within the Owyhee Mountains, Idaho, USA during water year 2000-2001, and is described in chapter 6. 102

Chapter 6 Fieldwork 6.1 Introduction A new experiment was designed to collect essential data to validate SNOWCAN, which were not available from previous experiments. Fieldwork took place during January-May 2001 within the Reynolds Creek catchment in Idaho, USA. This experimental watershed is owned by the US Department of Agriculture (USDA) and is managed and maintained by the Northwest Watershed Research Center (NWRC) in Boise, Idaho. This site was selected after an extensive search of possible sites because of its comprehensive existing background observations, because of the presence of both vegetation-free and vegetated sites, and because of the interest from the scientists at NWRC. Details of the site are given in section 6.2. Two main periods were studied, during which continuous measurements of validation data were taken. These were 15-22 March 2001 and 13-22 April 2001. Several separate periods of 1-3 days duration were also used to collect necessary data throughout 103

6.2. Site Details Chapter 6. Fieldwork the season. Two sites were studied: an open site and a fir site. An open site simulation was carried out in chapter 7 to allow model calibration of compaction algorithms, and to determine the accuracy of the snow model decoupled from the canopy model. The open site parameterisation of the snow cover was then used for the simulation of the snow cover beneath a fir site. The methods used to determine atmospheric, canopy and snow parameters are described in section 6.3. The data used to initialise the snowpack are described in section 6.4, and the meteorological data used to drive the model are shown in section 6.5. The collection of validation data is discussed in section 6.6. The data described in this chapter are used to simulate the seasonal snow cover at an open site, and beneath a fir canopy, and the simulation results are presented in chapter 7. 6.2 Site Details The fieldwork site is located within the Reynolds Creek Experimental Watershed (RCEW), Idaho, USA. Maps showing the location of Idaho are shown in figure 6.1. Reynolds Creek Experimental Watershed has been instrumented since 1964. Details concerning the history of measurements are available from Hanson et al. (2000). The area chosen for the fieldwork is the Reynolds Mountain East (RME) subcatchment. Reynolds Mountain East occupies 0.36 km2 and ranges in elevation from 2024m to 2140m, as shown in fig 6.2. Two climate stations have been operational within RME since 1968. The first, an open site known as the "ridge site" or "site 176", is situated on the edge of the catchment. This site is exposed, and high wind speeds are characteristic at this site. The second site is situated approximately centrally within the catchment, in 104

6.2. Site Details Chapter 6. Fieldwork Figure 6.1: Map showing location of Reynolds Creek Experimental Watershed within Idaho State an open area adjacent to a stand of aspen trees. This site is known as the "snow pillow site", because several different snow pillows have been installed at this site. The snow pillow site was also used for an intercomparison of snow pillows during the 2000-2001 winter season [Johnson and Shaefer (2001)]. A third station was installed within the fir stand prior to the experiment, which is now known as the "fir site". The ridge and snow pillow climate stations provide the forcing data used to drive the model, which are shown in section 6.5. A panoramic image of the Reynolds Mountain East catchment is shown in figure 6.3. 105

6.2. Site Details Chapter 6. Fieldwork AI. [J.. v. I L:. 5 i ^ ' a^ II I~lL f.7 Il '' s, 7<t2: 7 a9~ t.. v'W * * E,P ~: 4 St L:DI t Q i/ 4 IL:' VAT I I...,.:: 4.:. 4., i. Il4alml. t. E ii NII Figure 6.2: Map showing elevation of RME Study area Figure 6.2: Map showing elevation of RME Study area 106

6.2. Site Details Chapter 6. Fieldwork 4)Z AW.0! t, ~.:$ Is 5 *5 / / tf 1^ Figure 6.3: Image of Reynolds Mountain East Watershed 107

6.3. Parameterisation Data Chapter 6. Fieldwork 6.3 Parameterisation Data In this section, measurements taken to determine atmospheric, canopy and snowpack parameters are described. The measurement methods are discussed and refinements are suggested for the future. 6.3.1 Canopy Parameters As described in chapter 4, only three parameters are required to represent the forest canopy in SNOWCAN. These are leaf area index (LAI), emissivity and canopy singlescatter albedo. Leaf area index was determined above the fir site snow depth sensor, and around the fir stand using a Decagon SF-80 Sunfleck Ceptometer. The ceptometer has 80 light sensors at 1cm intervals along a probe. The probe is placed horizontally beneath a canopy and measures photosynthetically active radiation (PAR) i.e. radiation in the wavelength band 400-700nm, averaged over all sensors. In addition, the ceptometer determines the ratio of unshaded to shaded sensors, which gives the canopy gap fraction. Further details about this instrument are available from Decagon Devices (1987). The experiment followed the method of Pierce and Running (1988), adapted as described below. Canopy transmittance was measured on 3 separate clear days between noon and 2pm during the winter season. A reference measurement of PAR (Qo) was made at the start and end of each sampling period at a nearby open location. Measurements were made at 15~ increments around a full circle centred around the snow depth sensor and subsequently averaged over the circle. Leaf area index was derived from the canopy transmittance [Pierce 108

6.3. Parameterisation Data Chapters 6. Fieldwork and Running (1988)] as: A In (Iswf/Ijs, o) (6.1) AL - K where the canopy transmittance Iswf / Isw,o is averaged over the circle and the light extinction coefficient Kv is given as 0.52 from Pierce and Running (1988) for a coniferous canopy. Difficulties encountered with the measurement of leaf area index included the lack of clear weather to take measurements, which was also a problem for measurements of atmospheric optical depth (section 6.3.3) and surface reflectance (section 6.6.1). This problem was compounded by data loss during bumpy SnowCat journeys. Five measurements of LAI (these ranged from 3.9 to 5.6) were saved within 5m of the location of the fir climate station. At the point location of simulation, a barbed wire enclosure protected the climate station and prevented access. The LAI required at this point was determined by taking measurements at 0.5m increments around three sides of the perimeter of the enclosure (the fourth side was blocked by a tree trunk). The tip of the ceptometer was directed towards the centre of the enclosure at all times. Using this adapted method, the LAI was determined to be 4.0 at this point location. An additional problem experienced when taking ceptometer measurements was that the screws of the Ceptometer were magnetic, which hindered measurement of the circular increments with a compass. Care was taken by the ceptometer operator to ensure that personal shadows did not fall on the ceptometer light sensors, and the bullseye spirit level was used to ensure that the ceptometer was kept as level as possible. The canopy single-scatter albedo was determined from laboratory measurements of the spectral reflectance of needles and branchlets taken from the fir stand at RME. 109

6.3. Parameterisation Data Chapter 6. Fieldwork The fir needles and branchlets, shown in figure 6.4, were placed on a sheet of black paper to eliminate reflectance from the underlying surface. The spectral reflectance of the fir elements under laboratory illumination, which mimics the solar spectrum, is demonstrated in figure 6.5. The broadband single-scatter albedo of the fir canopy has been determined as 0.25 by D. Marks, by a solar spectrum weighted average of the spectral reflectance from 350nm to 2500nm [ D. Marks, pers. comm.]. ii/i,,, Figure 6.4: Image of fir needles and branchlets used to determine canopy single-scatter albedo As in the BOREAS simulation (chapter 5), the canopy emissivity was assumed to be 97.7%, which was the value given by Monteith and Unsworth (1990) for a poplar canopy. Although a preferable approach may be to determine the canopy emissivity in the field, measurements of canopy thermal emission and canopy temperature are not easy to take, and have large errors associated with them. The impact of this parameter assumption on the model performance is examined in chapter 7. 110

6.3. Parameterisation Data Chapter 6. Fieldwork 4. 3' -- Solar Constant -- Lab Reference illumination '" \....Target 01 1 E 3........ Target 012 XC Target 013 25 Target 014 Target 015 0 x Target 016 t — 0 500 1000 1500 2000 2500 3000 Wavelength (nm) Figure 6.5: Laboratory reflectance spectra of fir needles and branchlets (from D. Marks, pers. comm.) 6.3.2 Snow Parameters The rate of litter deposition on the snow surface was determined following the method of Hardy et al. (2000). Although this may be regarded as a canopy parameter, it is included here with a discussion of its effect on the snow albedo. Daily digital photographs were taken of the snow surface in the forest plot during periods of intense measurement; some examples of these are shown in fig 6.6. Percentage of litter coverage was determined from these images. Photographs of the snow surface were converted to binary images using image processing software. The proportion of black to white pixels were counted and the resulting fractional coverage of litter determined. 111

6.3. Parameterisation Data Chapter 6. Fieldwork Details about the modelling of surface and within pack litter are given in section 4.3.3. Figure 6.6 illustrates the change in snow surface over five days during the ablation. In particular, fig 6.6(e) and fig. 6.6(f) show the massive changes that can occur during one day. Changes during ablation are much larger than changes during the accumulation. In addition to new litter falling on the surface from the canopy, old litter becomes re-exposed to the surface as the snowpack melts and sublimates. However, figure 6.6(f) shows the maximum observed surface coverage of litter throughout the 2000-2001 season. From equation 4.8, a litter coverage of 7% lowers the snowpack albedo by <5%. Thus, the deposition of litter is assumed to have a negligible effect on the snow albedo throughout the winter season, and litter accumulation throughout the snowpack is not modelled in simulations at the RME fir site for the water year 2000-2001. The decrease in surface albedo for a thin snowpack was examined in the field, using the portable field spectrometer described in section 6.6.1. The spectral surface reflectance for different snowpack thicknesses are shown in figure 6.7. Although an albedo decrease was observed for a snowpack thickness 8cm, the snow surface at the measurement location was noted to be dirty, and that the ground could not be seen. The albedo decrease for a dirty 8cm snowpack is of similar magnitude to a snowpack of thickness 2cm, but is likely to be less for a clean surface. Based on these observations, the decrease in surface albedo with a decline in snow depth was assumed to be negligible in the simulations presented in chapter 5 and 7. 112

6.3. Parameterisation Data Chapter 6. Fieldwork (a) 13th Apr. 0% litter (b) 15th Apr. 0.01% litter (c) 16th Apr. 0.3% litter (d) 17th Apr. 1% litter (e) 10:00 18th Apr. 1% litter (f):L7:30 18th Apr. 7% Litter Figure 6.6: Photographs of snow surface during April field campaign. Litter percent coverage (m2rr-2 x 100%) is given for each photograph. 113

6.3. Parameterisation Data Chapter 6. Fieldwork mid-dirt( 56 cm depth | (u I ^........;. 0,8 0,2 0 500 000 1500 2000 25WX ____ Wavelength (rnm 40 t......... 0.8 0,2 0 t0 m0.6 I ^n ^ S end of drift; 24 cm depth *' 500 1000 1500 2000 2500 Wavetegth inm) I m N of dift end 0.8 I 0.6 OA 0.2 0 2 im N of drift e:t 20 cm depth 0 500 100lo0 1500 2000 2W50 Wavelength (nmo S edge of drift; 8 cm depth i0*")' M Ai Q -4- - 0 500 1000 1500 2000 250 Wavelength (nr) S edge of drift: 2 cm depth 0 Ig 0., 0 0 0 -: sM A I C 1000 1500 2000 25 Wave4ength (nml ' r 0.8 I 0.6 0 l.46 0.2 0o 0 500 000 1500 2 2500 2 Wavelength (nm) Figure 6.7: Snow surface reflectance variation with snow depth (A. Winstral, pers. comm.) 6.3.3 Atmospheric Parameters Spectral measurements were made of shortwave radiation from the sky using a CIMEL CE 318-3 Sun Photometer to determine the atmospheric optical depth, which is used to determine the proportion of direct and diffuse solar radiation when shadowband radiometer measurements are unavailable (see section 6.5). Details of this instrument are available from the NERC EPFS (Natural Environment Research Council Equipment Pool for Field 114

6.3. Parameterisation Data Chapter 6. Fieldwork Spectroscopy) website: http://www.soton.ac.uk/-epfs/. As discussed below, useful data from this instrument have not become available during the course of this project. The discussion of this instrument and its use within this experiment is discussed here to provide insight into the experiment philosophy and to indicate future work. The Cirnel CE 318-2 is a portable instrument that automatically computes the position of the sun and tracks its movement, and measures sun and sky luminances in 8 filters over visible to near infrared wavelengths. The CE 318-2 is fitted with filters at 440 nm, 670 nm, 870 nm, and 1020 nm for measuring atmospheric aerosol optical thickness, and a filter at 936 nm for measuring atmospheric water vapour. The CE 318-2 is additionally fitted with 3 polarised filters at 870 nm. The sun photometer was installed at the ridge site and was operated during days of little or -no cloud cover, as shown in figure 6.8 within the periods of continuous measurement. The sun and sky luminances were measured automatically in the principal plane (principal plane scenario) and azimuthally in the plane of the sun (almucantar scenario). The measurement scenarios were run approximately every 2 hours, because the instrument required the scenarios to be started manually. However, a significant problem was encountered during sun photometer operation. The sun photometer detects the brightest point in the sky in the approximate direction of the sun, calculated from latitude, longitude and time of day. Under conditions of light cloud cover or greater, the sun photometer often detects the brightest spot as a gap in the cloud, which it assumes to be the sun, and all measurements are then meaningless. This process could not be prevented manually. Hence, the predominantly cloudy conditions throughout the season limited the use of this instrument. 115

6.3. Parameterisation Data Chapter 6. Fieldwork Figure 6.8: Tracking and measurement of the sun using a Sun Photometer In addition to the operational problems of the sun photometer, the software required to convert the intrument outputted voltages to irradiances and to determine the optical depth (and other useful atmospheric parameters) is not readily available. For the purposes of this study and the simulations presented in chapter 7, the atmospheric optical depth under cloud-free conditions was assumed to be constant, and was estimated as 0.07, which is typical for a clean atmosphere such as at RME. This parameter is only required when shadowband data are unavailable, but may be a source of error in the model. It is hoped that the atmospheric parameters can be determined from the sun photometer measurements in the near future. 116

6.4. Snowpack Initialisation Chapter 6. Fieldwork 6.4 Snowpack Initialisation Snow pits were dug at the fir site on 16th February 2001 and at the ridge site on 31st January 2001 for model initialisation. Snow density samples were taken every 10cm vertically, using a density sampler of known volume and weight, and were weighed with an electronic balance. Temperature profiles were measured with a long digital temperature probe, which was inserted into the snow at 10cm vertical intervals. After profile measurements were taken at the ridge site, the snow pit wall was shaved back and the measurements were repeated to gauge local scale spatial variability, and therefore temperature and density measurement uncertainties. At the fir site, additional measurements of the thickness of different snow layers were made. Crystal sizes and shapes were examined with a pocket microscope. The snowpack profiles measured at the ridge and fir sites are shown in figures 6.9 and 6.10 respectively. Ridge Temperature Profite Ridge Density Profile 30/01/01 30/01/01................................................... 7 -.... -.. 1 60 - 60 -I50 50 40 - E 40 30 3 030 - -.-Pitl 20 20 10 10 - 30o -10 -8 -6 -4 -2 0 200 300 400 Temperature (deg C) Density (kg m. -3) Figure 6.9: Ridge snow profile initialisation data 117

6.5. Climate (Forcing) Data Chapter 6. Fieldwork Snow pit profiles at the ridge site are shown in figure 6.9. The time difference between measurements was less than 30 minutes, and the spatial difference was approximately 50cm. From the two sets of measurements taken in the snow pit, the uncertainty in density measurements was estimated to be ~ 50 kg m-3, which is approximately the same as the measurement instrument error, and the uncertainty in temperature measurements was estimated to be ~ 1K. The temperature, density and grain shape profiles at the fir site are presented in figure 6.10. The snow stratigraphy has been plotted using the SnowPit98 macro program developed by Shultz and Albert (1998). Fir snow pit 16/02/01 4 5................................................... 40 -35 30 25 ~ 20 15 / 10 —. 5. 0 -I ------ 0 200 400 Density (kg m'-3) Fir snow pit 16/02/01 Fir snow pit 16/02/01 45 - 40 - 35 - -30 c = 25 8 20 -15 I 10, 0 El 0 Q Q 47 Julian Day................................... Figure 6.10: Fir snow profile initialisation data 6.5 Climate (Forcing) Data Meteorological climate data, taken automatically by the climate stations installed at the ridge and snow pillow sites are presented in this section. The data are hourly averages 118

6.5. Climate (Forcing) Data Chapter 6. Fieldwork and are transferred by telemetry to the USDA-NWRC daily. Meteorological data taken by the fir climate station1 are included within this section for comparison with the ridge and snow pillow climate data. Although the fir meteorological data may be used to drive some Snow-SVAT models, SNOWCAN uses above-canopy or open site data to drive the model. Sub-canopy measurements, such as those taken by the fir climate station, are used to validate SNOWCAN in chapter 7. The climate stations at the ridge, snow pillow and fir sites measured downwelling solar radiation, relative humidity, wind speed and air temperature. Wind direction was measured at the ridge and snow pillow sites, and thermal radiation was measured at the snow pillow station. Precipitation was measured using shielded and unshielded Nipher gauges at the ridge and snow pillow stations. Snow precipitation is notoriously difficult to measure, and errors occur mainly because of gauge undercatch [Pomeroy and Goodison (1997)]. The true value of precipitation SWE was estimated from the shielded and unshielded gauge measurements using the unpublished method of Hanson [A. Winstral, pers. comm.]. Battery problems were experienced at the fir climate station, which were fixed on Julian day 60. Data acquired prior to day 60 should therefore be regarded with suspicion. Data from the fir climate station were only used to compare the simulated sub-canopy wind speed with measurements made at this climate station (section 7.4.4). Climate data for the ridge, snow pillow and fir sites are shown in figures 6.11, 6.12 and 6.13 respectively. Correlation of air temperature, relative humidity and solar radiation between the ridge and snow pillow sites were good, with r2 correlation coefficients ranging 1The fir meteorological data are not transferred by telemetry - the data storage module is collected on a two week basis. 119

6.5. Climate (Forcing) Data Chapter 6. Fieldwork 2m Air I 3 0 0............................................................................................................................................ 290 20 210 270 4 260 250 "; 25... i,,,,i.;.,,,. i,.. ~.,,.,..;.......;...,. i... I....;... Temperature,,:W.:.\.\,.......:...: '''.1''''1 '''1 '''1''I''I''' l I''I''I'''''"' " ' " ' I" "'''' 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian Day 120 - 6 100 -E BO - - S60 40 20 -It 2m Relative Humidity, r1>............................................................................... 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian Day 2m Wind Speed E (11 -a Tur o I 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Solar Radiation 1000 E oo 600 c 400, 200 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day... re ip tati.................................................................................................................................................. 7 -.............. 6 -L-; 5 -l I i L. LI L. L., i 'i * i.J.. 1 1 1... I.L,...I. 0 -- '... ",....,...,.....,J.,,.- L...;....... 5., i, M... d...... L..... 1 I - I.... I I. I 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 Julian day 75 80 85 90 95 100 105 110 115 120 125 Figure 6.11: 2001 Meteorological data for open (ridge) site 120

6.5. Climate (Forcing) Data Chapter 6. Fieldwork --------— — 2m Air Temperature 30. 290 280 260 250 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian Day 2m Relative Humidity 1 2 0.......................................................... 100 E 80 40 a 20 - / 20 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day 2m Wind Speed E 10 -5 0 c 0 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Solar Radiation 1 2 0 0............................................................................................................................................................................................................................................................................................................ 1000 E 200 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Thermal Radiation Precipitation 300 3: 25- 0 200 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Figur e 6.12: 2001 Meteorological dation E3 E 2 o 5 10 1 20 25 30 35 40 45 50 5,., 60 65 70 75 80 85 90 95 100 105 110 115 120 125 JuLian day Figure 6.12: 2001 Meteorological data for RMSP 121

6.5. Climate (Forcing) Data Chapter 6. Fieldwork 2m Air Temperature 300 -................................................................................................................................................................................................... 290 280 270 260 250.. 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day 2m Relative Humidity 120 100..l E 80 20 - I 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day 2m Wind Speed 15...-................................ 15 -0,.-,h.' A, |, 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Sub-canopy Solar Radiation..........1200.............. 1000 -E J 500 X 400 200 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian dayt Figure 6.13: 2001 Meteorological data for fir climate station between 0.88 and 0.98. A comparison between solar radiation measured at the ridge and snow pillow sites in figures 6.11 and 6.12 indicates periods where the snow pillow pyranometer may be snow-covered. Icing problems were often observed at the snow pillow site, and one such event is shown in figure 6.14. At the ridge site, higher wind speeds decreased the precipitation accumulation on the meteorological instruments. Solar radiation data from the snow pillow site were considered to be much less reliable than data from the ridge 122

6.5. Climate (Forcing) Data Chapter 6. Fieldwork site, and only solar radiation measured at the ridge site was used to drive SNOWCAN (chapter 7). Figure 6.14: Accumulation on meteorological station at snow pillow site However, although the pyrgeometer at the snow pillow site is subject to the same accumulation problems as the pyranometer, the data do not show this effect clearly because the diurnal cycle is less pronounced. Since the pyranometer at the snow pillow site is the only pyranometer installed within RME, comparisons cannot be made with other thermal radiation data. The longwave radiation data used to drive the model is used with caution, given the icing problems observed at the snow pillow site. Diffuse radiation was measured at the ridge site with a shadowband radiometer, shown in figure 6.15. This instrument requires daily manual adjustment of the metal frame, 123

6.5. Climate (Forcing) Data Chapter 6. Fieldwork to follow the solar declination angle. Data from the shadowband radiometer during periods where no operator was present, were therefore discarded. In addition, on 15th March 2001, the shadowband radiometer was found to be incorrectly aligned (the shadowband shield faced magnetic north, rather than true north), and all data prior to this date were also discarded. Diffuse radiation measured over the season at the ridge site is shown in figure 6.16, with periods of bad data shown in dark. Although it is theoretically possible to retrieve some of the data, time constraints and the lack of adjustment when no operator was present means that it is not feasible to determine the proportion of diffuse radiation for the periods of bad data. Figure 6.15: Shadowband radiometer installed at the ridge site 124

6.6. Validation data Chapter 6. Fieldwork 74 a81 102 112 7 0 0..........................Shaidowband Operator Oper tor 600. 400500 not present i not pr sent E 500 iz c 400 -~ii 3 00 -rrii 1 1 1 ~ti 30 200 1 30 40 50 60 70 80 90 100 110 120 Julian day Figure 6.16: Diffuse radiation measured by a shadowband radiometer 6.6 Validation data In this section, data taken to validate the snow model are described. Measurements were taken of sub-canopy solar and thermal radiation, snow profiles and snow depth to test both canopy radiation and snow components of the Snow-SVAT model. These validation data are presented below. 6.6.1 Sub-Canopy Radiation An array of radiometers was installed within the fir site during periods of intense measurement. The array consisted of 10 Eppley pyranometers and 3 Eppley pyrgeometers. The radiometers were placed randomly around the snow surface to measure sub-canopy downwelling radiation, and were moved and relevelled on each sunny day. A typical2 diurnal cycle of solar radiation measurements from all pyranometers is demonstrated in figure 6.17. The upper measurement envelope represents direct radiation from the sky 2These radiometer array measurements were made under a stand of aspens during the RME 2001 Experiment, with a different radiometer array. Section 6.7 highlights additional validation data collected that are not described or used elsewhere in this thesis. 125

6.6. Validation data Chapter 6. Fieldwork that has penetrated the canopy (a sunlit radiometer), whereas the lower measurement envelope represents the proportion of diffuse radiation that has penetrated the canopy (a shaded radiometer) [D. Marks, pers. comrn.]. The physical meaning of these envelopes is illustrated by the inclusion of the shadowband measurement, and the ridge site solar radiation measurement for that particular d(lay, which was sunny. Prior to the second period of intense measurement, one pyranometer and one pyrgeometer were inverted at the fir site to measure sub-canopy upwelling solar and thermal radiation. The radiometers were mounted approximately Im above the snow surface and because of the hemispherical nature of these instruments, surrounding objects (including vegetation) in the field of view and shadows of the instruments on the snow surface all contribute to radiation measurement errors. A second pyranometer and pyrgeometer were mounted on the same apparatus, but were upward-looking and hence measured downwelling radiation. The remaining radiometers were scattered around the snow surface, as before. The arrangement of mounted radiometers is shown in figure 6.18, and the down- and upwelling solar and thermal radiation measurements from these radiometers are presented in figure 6.19. Downwelling thermal radiation in figure 6.19 shows that the trees are heating up during days 106-108. In addition to the upwelling solar radiation measured by an inverted Eppley pyranometer at the fir site, measurements of the spectral reflectance of the snow surface (both in the open and under forest canopies) were made throughout the experimental period under a range of meteorological conditions using a GER3700 handheld spectroradiometer. The spectroradiometer measures radiance in the range 350-2500nm in 650 narrow spectral bands. Details about this instrument can be found on the GER website 126

6.6. Validation data Chapter 6. Fieldwork (a)l Sunny day - approximately 38% of total solar radiation is diffuse 1000 900 - -. - Radiometer array, 800 exposed radiometer E 700 D. 600. * Shadowband 600 radiometer 5050 400 - 300 1 0r 200 "I 0 0 5 10 Hour 15 20 (b) Cloudy day- approximately 82% of total solar radiation is diffuse Figure 6.17: Proportion of direct and diffuse radiation determined from a radiometer array under a stand of aspens 127

6.6. Validation data Chapter 6. Fieldwork Figure 6.18: Upward and downward looking pyranometers and pyrgeometers at the fir site (http://www.ger.com/ground.html). The surface reflectance was determined from the measurement of radiance of a BaSO4 standard reflectance panel and the subsequent measurement of the radiance of the snow surface. The surface reflectance is calculated as the ratio of these measurements at each wavelength. There are, however, a number of problems associated with this measurement technique, which are detailed below. A General problems 128

6.6. Validation data Chapter 6. Fieldwork 500 450 400 - < 350 E 3 300 - x | 250 - l4 al.- 200 -.-5 - 15Q100 - 50 - 0 - - Downwelling solar - Upwelling solar - Downwelling thermal - Upwelling thermal - %, I: j f. i I 'ir' * I A...:'Q,, i 4...,,, on\-.t^' W 1,N 1,85 -I l i l h 1 I! %. I I 1 I I I 1 1 ~ 9 2 94 96 98 100 102 104 Julian Day 106 108 110 112 Figure 6.19: site Up- and downwelling thermal and solar sub-canopy radiation at the RME fir A-I Cloud cover. Reference panel and target surface measurements are not simultaneous, each having a measurement time of approximately 10 seconds. This can result in significant errors on cloudy days, where lighting conditions change rapidly through the movement of clouds. These changes may not be perceptible to the eye. Figure 6.20 illustrates the magnitude of the cloud effect. Measurements were taken on a diffuse day using a standard panel for both the reference and the target, hence an expected result would be perfect i.e. 100% reflectance. However, even though cloud changes were barely observed with the naked eye, a change in reflectance of up to 25% was measured using the spectroradiometer. 129

6.6. Validation data Chapter 6. Fieldwork ----------- ----------- ------ -- -----— ~ --------- ------ --— ---- -- ---- - — -- — --- - --- Effect of cloud cover on reflectance measuremerts 130 E' 120......... 100 - I.~ 90 8o 70 60 0 500 1000 1500 2000 2500 3000 Wavelength (nm).........._-...........-......._ _ -..........::..:.......:w:.. -.........-_..:.........:.....-.......::...... --- —::...... ----:...............:_._ _.......:.......-.....:::_. _-. ---- -w::..__._......:.::_-:_. ----:-:-.:_.:.............-.....:: —.:.._............:.-......: -_- _-....... — Figure 6.20: Illustration of reflectance measurement error due to cloud cover This error becomes even more important over highly reflective surfaces such as snow, where relative changes in reflectance are small but have a large impact on the surface energy balance. The cloud cover error illustrated in figure 6.20 is relatively small in the visible region of the spectrum, where snow reflectance is less variable. This effect is greater towards the infra-red region of the spectrum, where snow reflectance is more spectrally dependent. A-II Instrument Field of View. The reflectance panel used is 30x30 cm2 in area - this limits the height at which reflectance measurements can be made. This is probably the largest panel that can be used in the field as it is heavy to hold and transport, is expensive and must be kept clean and damage-free for optimum reflectance. 130

6.6. Validation data Chapter 6. Fieldwork The measurement height governs the instrument field of view (FOV). Two different optics were used in this experiment - the standard optic (3~ FOV) and fibre optic (23~ FOV) attachments. These are not easily interchangeable in the field and hence were used on separate days. By simple geometry, to achieve a FOV < 30cm diameter, measurement distances of < 35cm and 2.85m are required using the fibre and standard optics respectively. These are the maximum heights that can be used to measure the radiance of the reference panel. The minimum measurement distance is limited by the instrument shadow on the reflectance panel, which depends on the solar geometry relative to the instrument set-up. This effect is easily observed on sunny days, but extra care must be taken on days where the proportion of direct radiation is small, but still significant, as the operator may be unaware that the equipment (particularly the fibre optic) is casting a shadow on the reference panel. The most appropriate optic to use is dependent on several factors. * Surface heterogeneity * Type of measurement * Environmental * Convenience. A-III Solar Geometry. Snow is a non-Lambertian reflector, exhibiting strong forward scattering [Glendinning (1997)]. Therefore it is important when comparing reflectance measurements to ensure that differences caused by sun movement are minimised. This can be done by taking geometry consistent measurements 131

6.6. Validation data Chapter 6. Fieldwork e.g. at the same time of day, in constant diffuse conditions (problems are caused by cloud cover) or by using a wider FOV optic (which would measure radiance from the snow surface at a range of angles). Where the bidirectional reflectance distribution function (BRDF) is to be represented, a narrower optic will allow more accurate measurements. A-IV Surface Heterogeneity. Surfaces with small scale topography, such as a windblown snow surface may vary in reflectance, depending on the local slope of the surface measured. This effect may be minimised by appropriate measurement techniques such as choice of optic to average small scale variations or to eliminate larger scale variations. The consideration of heterogeneity scale is also relevant to surface contamination, as described below. A-V Anthropogenic Change in Lighting Conditions. Perhaps the most surprising error is the effect that human movement can have on the radiance measurements. This is most relevant to measurements over snow, where the surface reflectance is high (errors have a more significant effect) and the target surface cannot be disturbed. The reference panel is slotted under the optic at an appropriate distance and must be held level for at least 10 seconds, as shown in figure 6.21. However, this is very uncomfortable for the reference panel holder. As a result, it is a natural tendency for the operator to move away from the equipment as the reference panel is moved away. It was continually observed that the change in lighting conditions caused by the movement of the operator resulted in target radiance measurement greater than the reference panel i.e. 'reflectance' of 132

6.6. Validation data Chapter 6. Fieldwork Figure 6.21: Operator holding reference panel under spectroradiometer for reflectance measurement >100%, even though no shadow was cast on the reflectance panel. This error was subsequently reduced by maintaining operator position but removing the reference panel from the FOV for the target measurement by sliding the panel as horizontally as possible. B Canopy specific problems B-I Canopy Structure. The radiation incident on the snow surface can be extremely heterogenous, because the structure of the canopy governs the amount of direct radiation that penetrates the canopy to the surface below. The radiation heterogeneity is indicated by the shadows cast on the snow surface, which can also change with wind-induced canopy movement. The spatial variation in downwelling radiation causes problems when 133

6.6. Validation data Chapter 6. Fieldwork attempting to determine the surface reflectance. If the canopy is relatively dense, then the reflectance of a shadow may be measured, or the reflectance of a light patch may be determined if there is a sufficient gap in the canopy. However, neither of these approaches will give a surface reflectance that is representative of the general forest surface. The reflectance of a surface that has a mixture of sunflecks and shadows is extremely difficult to measure. There is a finite height difference between the reference panel and snow surface, which cannot be disturbed. Therefore the shadows that fall on the reference panel are different to those that fall on the surface, and the surface reflectance measurements may be dominated by differences in incident radiation, rather than differences in surface reflectance. B-II Surface Contamination Surface contamination, such as leaf litter, absorbs more radiation than the snow cover, and lowers the surface albedo. Contaminants are distributed heterogeneously, and the measurement of surface albedo may vary greatly. This is even more important for a smaller field of view. Despite the numerous difficulties in the measurement of snow surface reflectance, the spectral reflectance of the snow cover at the ridge and fir sites was measured throughout the season. Two decay periods were captured, where measurements were made daily after a fresh snowfall event. One such period is shown in figure 6.22. Additional measurements were made using a similar instrument, the ASD FieldSpec Pro. This instrument only measured radiation from 250nm-1100nm, but showed much less sensitivity to the problems highlighted above. The data from this instrument therefore allows diagnosis of spectral anomalies shown by the GER intrument. 134

6.6. Validation data Chapter 6. Fieldwork 0.9.. 0.6, 0.5 < 4 - Apr-14 -Apr-15 0.2 Apr-17 0 500 1000 1500 2000 2500 3000 Wavelength (nm) Figure 6.22: Decay of snow reflectance at a fir site measured by a GER handheld spectroradiometer. The plots are averages of several measurements, each of which require assessment for measurement problems. However, each reflectance spectrum in each daily set of spectra at each site requires careful consideration of the possible errors. For validation of the Snow-SVAT against surface albedo, the reflectance spectra need to be averaged to determine the broadband albedo for comparison with the simulated albedo. The averaging is not linear, since the solar energy spectrum is wavelength dependent. This work is currently being undertaken by Dr. Danny Marks at USDA-NWRC, but the broadband albedo change over the fieldwork season is not yet available for model validation. 6.6.2 Snow Profiles As described in section 6.4, snow profile measurements from snow pits were taken at the ridge and fir sites for model initialisation. Further snow pits were dug at the open and fir 135

6.6. Validation data Chapter 6. Fieldwork sites in order to determine density, temperature and grain size profiles for model validation. Figure 6.23 shows temperature, density and grain shape profiles taken throughout the accumulation and ablation season at the ridge site. Figure 6.24 shows temperature, density and grain shape profiles taken on 16th March at the fir site. A snow pit was not dug at the fir site during April because the snowpack was too shallow. Continuous temperature profiles at the fir site were detemined in addition to the snow pit temperature profile described above. Thermistors were mounted every 10cm on a wooden frame, and the fir thermistor array was installed during the accumulation period by gently sliding the array into a vertical snow pit wall, as shown in figure 6.25. Great care was taken during this procedure to prevent the formation of air gaps around the thermistors and to disturb the snowpack as little as possible. Data from the fir temperature array are shown in figure 6.26. The thermistors were blue in colour, and were coated in a clear resin, both of which affected the radiation balance in the surrounding snow. The tips of the thermistor array, except for the thermistors, were covered in white tape to minimise radiative heating of the thermistor array. Fig. 6.26 shows several typical features. The surface variation in temperature is damped and lagged by the snow, which depends on the snow thermal conductivity and therefore density. Another feature observed is the ripening of the snowpack. On day 65, the snowpack becomes isothermal. Although the temperatures shown are not at 0~C, they are within the 0.5~C error associated with the snow temperature array set up. For model validation, all measured temperatures from the thermistor array were increased by a constant 0.35K to correct for the array temperature offset. 136

6.6. Validation data Chapter 6. Fieldwork Ridge snow pit 17/02/01 80......................... 70 60 c50 E I40 \ 30 - 20 10 250 300 350 400 450 Density (kg m' -3) Ridge sr 80 70 I 60 u50-.r tL 40 - 3 30 20 10 o now pit 17/02/01 a *I Julia Day -S 1 (a) 17th Feb Ridge snow pit 16103/01 | 9 0............................. 80 70 60 ' 50 30 20 10 0 200 400 600 | Density (kg rn'-3) Ridge snow pit 16/03/01 90 80 7(0 6(0 c 501 3( 2(O 10 21 -3 - 20t __ _ _ * w * TII 3 JuliP Day -3 -2 -1 0 Temperature (deg C) (b) 16th Mar Ridge snow pit 17/04/01............................... 6..0 l < 50. 40 30 20 loI -0.6 -0.4 -0.2 0 Temperature (deg C) Ridge snow pit 17/04/01 60... 60( I I 50 40 -. 30 20 10 - I Ridge snow pit 17/04/01 O I...........-............ 0 200 400 600 Density (kg m'-3) (c) 17th Apr Figure 6.23: Validation data from ridge snow pits 137

6.6. Validation data Fir sniow pit 16/03/01 F ir sw pit16/031 Fir 60 60., __.~, I. --- —---- ( 60 50- 50 5 40 40 40 E./ E 30 30 30 20 20 I I 20 10 m10 o.10 0 ----—.ue...... 1 0... 0 200 400 i -10 -5 0 Density (kg m'-3) Temperature (deg C Figure 6.24: Validation data from fir snow pit: Chapter 6. Fieldwork snow pit 16/03/01 * 0? Q... — - ^........ — -.. 75 Julian Day.. th M arch 2001.................... 16th March 2001 Figure 6.25: Insertion of temperature array into snow pit wall 138

6.6. Validation data Chapter 6. Fieldwork Temperature array data 10 ------- 10cm 20cm ' '"30cm... 40cm (in aW) -15, " not al 8i 8 l 8 8 ~ coDay Figure 6.26: Temperature profile of snowpack under fir canopy Subsequent melting-out of the thermistors can be seen on days 78, 80 and 82 for probes 30cm, 20cm and 10cm from the ground respectively. On these days, the thermistors became exposed to the air and temperatures followed the 40cm probe, which was exposed to the air throughout the season. 6.6.3 Snow Depth Figure 6.27 shows the snow depth at the ridge, fir and snow pillow sites. The snow depth was measured using an ultrasonic depth gauge at each of the three climate stations, and the data were transferred back to USDA-NWRC by the same method as the climate data described in section 6.5. 139

6. 7. Summary Chapter 6. Fieldwork Ridge Site -1 E 1 'O 042 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 B0 65 90 95 100 105 110 115 120 125 Julian day Snow Pillow Site 0 o.8 o 0.2 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Fir Site 0.8 0.6 - 0 5 10 15 20 25 30 35 40 45 50 55 60 55 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Figure 6.27: Snow depth depletion at RME ridge, snow pillow and fir sites 6.7 Summary In this chapter, measurements taken as part of the Reynolds Mountain East Experiment 2000-2001 were described. More data were taken over the whole RME basin than has been described here, including two catchment depth surveys. Validation measurements, similar to those made at the fir site, were also made within a stand of aspen trees. Further details about this experiment can be found on the experiment website: http://www.nercessc.ac.uk/-mjt/rme. Data from this experiment have been used to initialise and drive SNOWCAN. The seasonal snowcover at the ridge and fir sites were simulated for the water year 2000 - 140

6.7. Summary Chapter 6. Fieldwork 2001, and the model was validated against the measurements described above. These model simulations are presented in chapter 7. 141

Chapter 7 Snow-SVAT Simulation II: Reynolds Creek 7.1 Introduction The data collected as part of the BOREAS experiment were insufficient to validate SNOWCAN, as described in chapter 5. A new experiment was designed and carried out within the Reynolds Creek Experimental Watershed, Idaho, to allow a full test of SNOWCAN to be carried out. The experiment, and data collected were described in chapter 6. In this chapter, the accumulation and ablation of the snowcover at two sites within the Reynolds Mountain East catchment is simulated. Model validation at an open site (section 7.2) is used to determine the best parameter set, and the parameters are compared to those used in the BOREAS simulation (section 7.3). The accumulation and ablation of the snowpack at the fir site is simulated (section 7.4), using the open site parameterisation of the snowcover. SNOWCAN is validated against measurements of sub 142

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek canopy radiation, snow depth and snowpack profiles of temperature, density and grain size. The model sensitivity to the canopy parameters is also investigated (section 7.5). 7.2 Open Site Simulation This section describes the simulation of the seasonal snow cover at the RME ridge site. The inputs used to initialise and drive the model are described in section 7.2.1. The model is validated against measured snow depth and temperature and density profiles in sections 7.2.2 and 7.2.3. Chapter 3 describes the sensitivity of the snow model to its parameterisation, and several parameters are highlighted, for which further investigation of model sensitivity with non-idealised input data could be beneficial. The model sensitivity to these parameters is studied in section 7.2.4, using the RME ridge site input data. 7.2.1 Model Inputs Model Initialisation The model was initialised for the open site from snow pit measurements made at 15:30 on 30th Jan 2001, as shown in figure 6.9 (section 6.4). No measurements of grain size or layer information were taken from this snow pit, therefore each snow layer was arbitrarily initialised with 10cm width, except the two upper layers, which were initialised with 2cm and 1cm widths. The average grain size of the lower half of the snowpack was assumed to be 1mm and the average grain size of the upper half of the snowpack was assumed to be 0.5mm. The forcing data used to drive SNOWCAN were described in section 6.5. Thermal radiation was not measured at the ridge site, therefore measurements of longwave radiation 143

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek at the snow pillow site were used to drive the model. The input files used in this simulation are included on the CD that accompanies this thesis. 7.2.2 Snow Depth Validation 1...................................... 0.9- - 0.8 e 0.7 - o 0 n + Measured | * \ -~0.6 -0.4 - c: 0.3 -0.2 - Simulated 0. 1 30 40 50 60 70 80 90 100 110 120 Julian day Figure 7.1: Comparison between observed and simulated snow depth at RME ridge site Snow depth simulated at the RME ridge site is shown in figure 7.1. This graph shows that there are obviously problems with the precipitation forcing data. Precpitation is difficult to measure, particularly in a windy environment such as at this site. Nonetheless, measured and simulated snow depths show good agreement, with a regression coefficient of 0.94. The final melt period (days 112-120 in figure 7.1) is modelled less well than the 144

7. 2. -Open Site Simulation Chapter 7. Snow-SVAT Simulation I: Reynolds Creek rest of the ablation period. Without examination of the snow physics within the pack, it is difficult to determine whether the difference in snow depth occurs in the mass balance e.g. compaction or water flow, or in the energy balance, such as a decrease in snow albedo for a thin snowpack, or in the lower boundary Dirichlet condition. The internal snow physics can be investigated through the comparison of simulated snow temperature, density and grain size profiles with the field measurements described in section 6.6.2. 7.2.3 Snowpack Profile Validation A comparison of the simulated snow temperature, density and grain size profiles with field measurements is demonstrated in figure 7.2. The agreement between simulated and measured profiles is generally excellent, and to within the measurement error. This gives an indication that the model simulates the seasonal snowcover well, both in accumulation and ablation periods. A number of features shown in figure 7.2 are extremely interesting. The temperature profile during the accumulation period, illustrated on day 48, is remarkably accurate. Differences in the surface temperature are discussed later in this section. However, in the ablation period, the simulated and measured temperature profiles show less agreement. The lower half of the simulated snowpack retained liquid water, which was not observed in the field. An examination of the grain size profiles in figure 7.2 indicates an underestimation in the average grain diameter for the lower half of the snowpack on day 75. Liquid water flow is dependent on the grain size. If the grain size is underestimated by a factor of 2, which is conceivable from figure 7.2, the liquid water mass influx to the layer below is 145

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek iI I I,,.T. 1. t T i I t i e p,c - r, C - r '.I uCJ] jL4i;~| I - - T X,,, I r,,, 8 '4 r' tN I I,-; 'l, e- -i o i I - 41 IC, J I 0 -,.. ) jr., _(O II) I1|RfH I fP I I,-.I M bn P 2 E 0rn". - I.,;. i J I: ' i I ' ' ' - ''..... I t....) I'*"-.,:,,,-J;,-~. iL,._ T i r_l ('1) 14,~! H C;.4b rn in rt r4 o -13 0 0 0 (LLI) 145~i.H M,2 IT a, E. r,} (ll ] -l~l lfH I (N m J~: p - W i v ll t 'J l d d d o o.1s:I.. I,. A / sf' ) 1) () - I | I I i II. (u mkisi i............................... t- 0 0 JJ, - l. i- i. u.2 I, -h ( ) 14, L" I o Figure 7.2: Comparison of snow temperature, density and grain size profiles at the RME ridge site 146

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek reduced by a factor of 4, since U1 oc ~d (equations 2.25 to 2.27). The rate of liquid water removal from the snowpack is determined by the wavefront or seepage velocity, and is dependent on the saturation and thermal state of the snow in addition to the local water flux, as residual liquid water and heat deficits must be satisfied prior to liquid flow. The underestimation in grain size could, in part, arise from the initial conditions, since the grain size profile was not determined from the first ridge snow pit. The grain size profile on day 48 indicates that the grain size in the lower half of the snowpack is already underestimated in the accumulation period. To investigate whether the grain size underestimation results from the initial conditions or from the model, the simulation was repeated, but the snow grain diameters in the lower 30cm of the snowpack were initialised at 1.5mm, snow grain sizes in the 20cm above were initialised at 1mm, and the top 13cm of snow was initialised with an average 0.5mm snow grain diameter. Figure 7.3 shows the temperature and grain size profiles for days 48 and 75. In the simulation presented in figure 7.3, the average grain size diameter is larger in the lower half of the snowpack, than was measured in the field on day 48. However, by day 75, lower within the snowpack, the average grain diameter simulated by the model was still less than the average grain diameter measured in the field. The differences between simulated and measured temperature profiles were not resolved by this change in initial conditions. Figures 7.2 and 7.3 indicate that the model underestimates the rate of grain growth. The snow crystals at this depth were observed to be faceted, with slight rounding, and therefore indicate that the grain growth from temperature gradient metamorphism is not modelled accurately under these conditions. 147

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek. Julian day 48 ---— I8Measured 0.8 E.......... - Simulated 0.6 I- 0.4.......0 A | 0.2 266 267 268 269 270 271 272 273 274 Temperature (K) 0.8 - E 0.6 - 0.4 -0.2 - 0 Julian day 48......................................................................................................................................................................... I Measured.............. —.~............... IE............... Simulated,.............,....................:"................................................ m.. _:!- -9 - - -*...... -,j 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 Average grain diameter (m) I I, t - Julian day 75 Julian day 75 1 - 0.8 - E 0.6 -~ 0.4 - 0.2 -. a. Measured.................................... -; -e- Simulated..................K.......4................... 0 I 11. I..............I.....I 266 267 268 269 270 271 272 273 274 Temperature (K) i, i - I. I 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 Average grain diameter (m) Figure 7.3: Ridge snow temperature and grain size profiles: larger initial grains The removal of liquid water from the snowpack is enhanced by the formation of flow fingers, which are not modelled by SNOWCAN. However, flow fingers were not observed within the snow pit. Preferential water flow could occur around (unsaturated) or through (saturated) air pockets within the snowpack [Jordan (1995)]. The ridge site is heterogeneous., and has a large proportion of sage brush coverage, which is buried within the snowpack for most of the season. The locations of the snow pits were chosen to avoid sage brush within the snow, but local vegetation may affect the 3-D flow of liquid water. Another source of error in the model could be in the parameterisation of the irreducible saturation, which is dependent on the physical state of the snow, but is assumed to be constant in the model. 148

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek On day 107, both simulated and measured snow temperatures on day 107 reveal an isothermal, ripe pack, which matched visual observations of high liquid water content. The measured temperature profile suggests a thermometer offset of-0.5~C. This effect was not due to coupling of measurement with air temperature, because the air temperature was measured to be +8.2~, using the same temperature probe at this time. The temperature measurements on day 48 and 75 potentially have the same offset, with measured temperatures up to 0.5~C less than the actual temperatures. An underestimation in both simulated density and snowpack height, such as on day 48 in figure 7.2, indicates a difference in snowpack mass. This problem is greatest during the accumulation period, which suggests that the error is more likely to result from the spatial variability of the snowpack between the depth sensor (where the point model is focussed, as is discussed below) and the snow pit locations (where the model is validated), rather than from errors in the model mass balance. However, the model performance is degraded by uncertainties in the upper mass boundary forcing data, which are difficult to measure. The difference in location between the point simulation and the measurements used to initialise and validate the model causes other problems. For example, the difference in snowpack heights affect the temperature gradient within the snowpack. Day 48 in figure 7.2 shows minor surface differences in temperature. Snow profiles are sometimes compared from the surface downwards, since changes are most rapid at the surface. The advantage of this method is that model agreement with surface measurements may be improved, but the snow-soil interface may be at inconsistent depths. Here, the profiles are plotted from the snow-soil interface upwards, so that the height measurements are true 149

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek to the simulation and observations. This requires no visual processing on behalf of the reader and avoids height confusion. The simulated density profiles show good agreement with measured density profiles, to within measurement error, except for the top 20cm of the pack on days 48 and 107. Site 176 (ridge site) is an exposed site, with relatively large wind speeds, and features characterised by windblown or wind-scoured snow are often evident at this site. The model SNOWCAN does not simulate compaction due to wind effects, therefore an underestimation in snow density is expected at the surface of a simulated snowpack under high wind speed conditions. However, SNOWCAN simulated accurate densities lower in the pack, which should be borne in mind as the underestimation in surface densification may in part be compensated by an overestimation in either the pressure densification or destructive metamorphism. The snowpit on day 75 immediately followed a fresh precipitation event. The simulated density of new snow is computed from an algorithm using the wetbulb temperature. As compaction is a time-dependent process, it is unsurprising that the snow density profile is accurate, even at the surface, as the wind compaction is negligible. This demonstrates the effectiveness of the algorithm that computes the new snow density. The density profile on day 75 supports the argument of suppressed simulated water flow relative to that observed in the field. To illustrate this, the simulated partial ice density was compared to the measured snow density in figure 7.4 (the snow density would be equal to the partial ice density if all liquid water was removed from the snow sample). Figure 7.4 shows that the simulated density profile shows better agreement with field observations if the liquid water is assumed to flow more rapidly. 150

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek Julian day 75 1 0.8 I 0.6 _r to..3 0.4 I O.m 0.2 0 0 100 200 300 400 Density (kg m' -3) 500 600 Figure 7.4: Comparison of simulated partial ice density with measured snow density The validation of the snow model in this section indicates that the model is working well, particularly given the wide diversity in meteorological conditions experienced during the season simulated. However, as described in section 3, the model performance depends on its parameterisation. The model sensitivity to several important parameters was investigated further in the following section. 7.2.4 Model Sensitivity Chapter 3 examined the model sensitivity to the parameters used to constrain the model, under constant melting conditions. Several parameters were highlighted in section 3.7 for further investigation, because the model showed relatively high sensitivity, the parameter uncertainty was large or the parameter was expected to have a greater effect under realistic meteorological conditions. Table 7.1 shows the sensitivity functions (equations 3.1 and 3.2) to a ~ 10% change in parameter value, based on the comparison of simulated snow depth with observed 151

7.2. Open Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek snow depth. The sensitivity of the model to a ~ 50 kg m-3 change in the initial snow density was investigated, to determine the senitivity of the model to the uncertainty in the density measurement. The sensitivity of the model to the Dirichlet constant temperature lower boundary condition was examined for a ~ 1K change in temperature. The regression coefficient for the snow depths is included in table 7.1. Parameter [std] F+ F~ F Fs (r2)+ (r2)EHlO [2.0] 200 210 110 121 0.9304 0.9352 EESO [2.0] 210 210 140 130 0.9367 0.9335 Ycrit [150] 190 230 28 190 0.9323 0.9387 r70 [9 x107] 210 210 130 120 0.9356 0.9321 /3nir [500] 210 210 140 110 0.9340 0.9293 Qtz [0.4] 210 210 140 140 0.9359 0.9387 Tatm [0.07] 200 210 120 130 0.9304 0.9345 Initial 7y profile 360 240 350 -120 0.9472 0.9168 T1 [273.6K] 250 230 -100 180 0.8830 0.9402 Table 7.1: Sensitivity of model at RME ridge site. For the simulation presented in section 7.2.2, r2 -= 0.9402, Fm = 220 and Fs -= 160. The model is relatively insensitive to the parameters examined in table 7.1. The sensitivity functions are dominated by the last melt period in figure 7.1, which began on day 112, because the difference between simulated and measured snow depth is greatest during this period. A better agreement in the snow depths from day 112 until the complete ablation of the snowcover lowers the sensitivity functions, but may not reflect a general improvement in the simulation of snow depth by the model. The regression coefficient is regarded as a better test for the sensitivity of the snow model for the entire seasonal simulation, but should be used with caution, since snow depth data are not independent, as discussed in chapter 5. Table 7.1 shows that for most parameters studied, a change in parameter de 152

7. 3. Parameter Space Comparison Chapter 7. Snow-SVAT Simulation II: Reynolds Creek creased the model accuracy. This indicates that a good set of parameters was chosen for this simulation. C)ne interesting result is that out of the parameters and initial conditions studied here, the simulation is only improved by a 50 kg m-3 increase in the initial snow density profile (the error on the snow density measurement). The same parameters used to simulate the accumulation and ablation of snow at the ridge site will also be applied to the model for the simulation of the evolution of snow cover beneath the forest canopy at the RME fir site in section 7.4. 7.3 Parameter Space Comparison Most of the model parameters were kept the same for simulations of the seasonal snowcover at both the BOREAS jack pine site and the RME ridge site. The main difference in parameter space affects the model representation of snow compaction. The parameter is ycrit A higher value was used for the RME ridge site, which increased the rate of compaction relative to the BOREAS simulation. The boreal forest is an extremely cold environment during the winter months. By contrast, the RME field site is further south and is much warmer. Several melt periods were experienced at this site in water year 2000-2001. Compaction is a temperature dependent process, and whilst this is reflected in the model, a higher threshold snow density in a warmer climate enhances the temperature dependence. In addition, the model does not explicitly model the change in compaction rate with liquid water flow, so the higher parameter value used at the RME site reflects the observed increase in compaction at this site, especially from liquid water flow. The model would benefit greatly from a more vigorous treatment of snow compaction, to avoid necessary site-to-site calibration. 153

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek 7.4 Fir Site Simulation The seasonal snow cover was simulated beneath a fir canopy at the RME fir site during water year 2000-2001. In this section, the model is validated against measurements of sub-canopy snow depth, sub-canopy up- and downwelling solar and thermal radiation and measurements of the snow temperature, density and grain size profiles. The results indicate how well this Snow-SVAT represents the radiative interaction between the canopy and the snowpack, and the other physical processes occurring within the snowpack. 7.4.1 Model Inputs The model parameters, with the exception of the canopy parameterisation, were determined at the ridge site in section 7.2. The parameters were taken directly from the ridge site study and were applied here to the fir site. Parameters required to represent the canopy were determined experimentally in the field or in the laboratory, or taken from literature, as described in section 6.3.1. The model was initialised from measurements of temperature, density and grain size profiles taken from a snow pit on Julian day 47, as described in section 6.4. A finer grid was used at the fir site than at the open site to optiinise model accuracy. Model forcing data (downwelling thermal radiation, air temperature, wind speed, relative humidity and precipitation) were measured by the climate station at the snow pillow site, and these data are presented in section 6.5. The snow pillow site experienced problems with equipment icing, which is evident from the solar radiation data presented in figure 6.12. For this reason, the downwelling solar radiation data used to drive the model are taken from climate station measurements at the ridge site. All model input files are included on the 154

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek CD that accompanies this thesis. 7.4.2 Snow Depth Validation The accumulation and ablation of snow beneath the fir canopy at the RME fir site was simulated during water year 2000-2001. Simulated and observed snow depths are shown in figure 7.5. 0.7 0.6 0.5 0.4 C~0.3 C: 0.2 0.1 0 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Figure 7.5: Simulated and observed snow depth beneath a fir canopy in RME Figure 7.5 shows that the model simulates the snow depth during both accumulation and ablation periods extremely well. The excellent agreement between simulated and measured snow depths is shown by the regression coefficient, r2 = 0.99. In particular, the assumption that 30% of snow is intercepted by the canopy, and that all intercepted snow is sublimated, appears to be reasonable from figure 7.5. However, to demonstrate that the model captures the physical processes accu 155

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek rately, the model must be tested more rigorously. One of the main aims of this project, as discussed in chapter 1, is to model accurately the interaction of radiation between the canopy and snow beneath. The radiation component of the Snow-SVAT was tested in the following section (7.4.3). 7.4.3 Sub-Canopy Radiation Validation Solar Radiation A comparison between simulated sub-canopy down- and upwelling solar radiation with the shortwave radiation flux measured by the mounted upward and downward looking pyranometers is shown in figure 7.6 for the entire period of sub-canopy radiation measurement. A shorter period of five days was chosen to focus more closely on the model validation of the solar radiation component, as shown in figure 7.7. This period of five days consisted of two sunny days and three diffuse days, and are the only five days where a complete set of diurnal downwelling sub-canopy solar radiation measurements were taken by an array of pyranometers beneath the fir canopy, as described in section 6.6.1. Figures 7.6 and 7.7 compare simulated sub-canopy radiation with radiation measured by the single upward and downward mounted radiometers. A representative measurement of solar radiation beneath a canopy is not easy to acquire. Shadows from the canopy move across the pyranometer with the change in sun angle, and thus the measured solar radiation is highly variable. Figure 7.7 shows a measured radiation minimum around the middle of the day on sunny days for both upwelling and downwelling radiation. This was caused by the shadow of a tree trunk. A statistical test of agreement between the modelled solar sub-canopy radiation and observed 156

01 0 0 zj350 - E= 300 -X 250 -m 200 - MV 150 - ~t100 - 0 * I I I I I% I i 71 U I I I II I I iki p 7II I LIAJ - -in- -ffieaured 1-2 IL-& K WI m - I -,T-.,, i I i. - i i -.T-., -T 92 93 94 95 96 97 98 99 'OD 1D' 102 103 1]4 105 106 107 108 109 110 111 112 113 114 JL~iafl cay 3W 2 ~-W ---asured 2C0 150 76 ico 0 92 93 094 95 96 97 98 99 100 101 1C2 103 104 iC5 106 107 108 109 110 111 '12 113 1P4 Jula day m I C-i Cb CD 0 Cl) 0 CD 0. C) CD

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek 250 --- E 200 x I 50 10 -107 108 109,i 110 Julian day 111 112 (a) Sub-canopy downwelling solar radiation 2 5 0 -...................................................................................-....................................................................................................................................................................... l 200 -i- -Measured E I --- Simulated x 150 I 100 0 107 108 109 110 111 112 Jutian day (b) Sub-canopy upwelling solar radiation Figure 7.7: A comparison between simulated and measured sub-canopy solar radiation 158

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek solar sub-canopy radiation over the total period of measurement reveals a correlation coefficient of r2 - 0.46 for the downwelling radiation and r2 = 0.63 for the upwelling radiation, which is surprisingly low, given the apparent visual agreement in figure 7.6. Also, the downwelling radiation correlation is lower than the upwelling correlation. This is rather unexpected, because the model simulates multiple reflections between canopy and snowcover, and the model error in downwelling solar radiation is propagated to the upwelling radiation (as is indicated in table 7.2). The low correlation for downwelling radiation suggests periodic snow settlement on the pyranometer, which would not occur for the inverted pyranometer that measured upwelling radiation. The correlation between measured and simulated sub-canopy downwelling solar radiation was improved for a five day period, when the radiometer was known to be snowfree. The radiometer array was operated during this five day period (Julian day 107-111 inclusive). Out of the five days, two were sunny (107, 108) and three were cloudy (109 -111). The measured and simulated sub-canopy radiation was compared and the correlation coefficients for the total, sunny and cloudy periods are given for the downwelling and upwelling radiometers and for the radiometer array average in table 7.2. Met. Conditions (Days) Downwelling Upwelling Downwelling (array average) Mixture (days 107-111) 0.70 0.60 0.82 Sun (days 107-108) 0.60 0.52 0.77 Cloud (days 109-111) 0.93 0.88 0.90 Table 7.2: Correlation coefficient comparison between measured and simulated sub-canopy solar radiation for Julian days 107-111 inclusive. As described above, the correlation coefficient for the upward looking radiometer is improved for a period when the radiometer is known to be snow-free. The correlation 159

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek coefficient for the downward looking radiometer is relatively unchanged. In addition, the correlation coefficients for the upwelling radiation are less than the downwelling radiation, because of the propagation of model errors with the reflection of radiation. Table 7.2 shows that the model simulates the observed solar radiation extremely well during diffuse periods, but less well in periods of direct sunlight, as the measurement shows high variability over time. The agreement is improved if a spatial averaged measurement is used as the basis for comparison. Measurement variability is decreased, since at any given time, a proportion of radiometers will be in direct sunlight, and some will be shaded. Figure 7.8 shows the array averaged sub-canopy downwelling solar radiation measurement with the simulated radiation. The radiation minimum observed with a single radiometer, as shown in figure 7.7, is removed in a spatially averaged measurement. 250 - -- - Average Measurement E 200^ -s —Simulated E 200 -!; 150 CZ 6 50 107 108 109 110 111 112 Julian day Figure 7.8: Comparison between modelled downwelling sub-canopy solar radiation with spatially averaged observations 160

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek However, a model comparison to spatial averaged measurements gives a lower correlation coefficient on diffuse days. Under cloudy conditions, the measurement is much less variable, and one radiometer gives a representative measurement at that point location. The radiation beneath a canopy is dependent on the canopy characteristics. An averaged measurement is representative of the canopy as a whole, but not necessarily at the point location. For these reasons, a single radiometer measurement is more appropriate on diffuse days (if the radiometer is placed close to the point location) and an array averaged measurement is more appropriate on sunny days. In spite of the measurement difficulties experienced in the field, the model simulates the observed sub-canopy radiation well. Figure 7.9 shows the downwelling measured and simulated sub-canopy solar radiation in comparison to the open site incident solar radiation used to drive the model. The shielding effect of the canopy is clearly seen in figure 7.9, where the model replicates an observed factor of 4-5 reduction of solar radiation beneath the canopy relative to the open site. Thermal Radiation The thermal radiation component of SNOWCAN is demonstrated in figure 7.10, which shows both simulated and observed sub-canopy longwave radiation. The downwelling thermal radiation measured at the snow pillow site, which is used to drive the model, is also shown in figure 7.10. Figure 7.10a clearly shows the importance of thermal radiation treatment within Snow-SVATs, since downwelling thermal radiation beneath the canopy may be over 100 W m-2 higher than that in the open. Simulated sub-canopy thermal radiation appears to follow the measured radiation 161

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek -1000 -0 E 900 -x 700 -> 600 -co t 500 - 0 - 400 - - 300 - I 200 - 0 -* 10 -iiri i.. tIt I r X a $: I - X R X -. - '^ y v l A& a p ji I' I I I )7 ~l I Imm...... I -— mm F —. F... - 1 -.... -.11.-... 108 109 '110 Jutian day 111 112 Figure 7.9: Reduction of solar radiation beneath the canopy closely, which is quite different from the forcing data. However, the simulated radiation generally exceeds the measured radiation by approximately 10-20 W m-2, which is within the measurement error. Pyranometer and pyrgeometer design introduces an error called the 'Dome Effect', where the radiometer dome alters the radiation balance between the target and the sensor. This leads to an underestimation in longwave radiation by up to 10 W m-2 and up to 25 W rm-2 underestimation in solar radiation (Ji and Tsay (2000)). This could be one reason for the discrepancy between measured and simulated thermal radiation. An extra bias could be introduced from non-level radiometers. In addition, the radiometer has a hemispherical field of view and therefore all objects in the hemisphere have an effect. This includes tree trunks and people and specifically for the downward 162

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek (a) Sub-canopy downwelling thermal radiation 4300 N 3150 a) 250 o 2'00 - - -*' -Measured I150 - Simulated 92 94 96 98 100 102 104 106 108 110 112 114 Julian day (b) Sub-canopy upwelling thermal radiation Figure 7.10: A comparison between simulated and measured sub-canopy thermal radiation 163

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek looking radiometer, the supporting apparatus, wires and the case containing the dataloggerl. Aside from the offset between measured and simulated sub-canopy thermal radiation, the model simulates the sub-canopy thermal radiation balance extremely well. Modelled and observed radiation exhibit a correlation coefficient of 0.96 for the downwelling radiation and 0.84 for the upwelling radiation. The agreement is less good for the upwelling radiation, in part because the measurements exceed the simulation on days 107 and 108, in contrast to the general offset. Although the upwelling thermal radiation exceeds the emission of radiation by a black body at 273.15K by up to 3.1 W m-2, this occurs when the incident thermal radiation is also high. The snow is assumed to have an emissivity of 0.99 in this model, which would account for a 3.5 W m-2 increase in upwelling radiation (in addition to the ~315 W m-2 emitted from the melting snow surface), from reflection of downwelling thermal radiation of magnitude 350 W m-2. The periods where the measured upwelling thermal radiation exceeds the simulated radiation on days 107 and 108 therefore could result from reflection of downwelling thermal radiation rather than from emission of foreign objects in the field of view, in which case the snow emissivity may be set slightly too high in this model. On the other hand, the radiative heating and emission of surfaces other than snow, such as forest debris, which is known to be in the field of view and to have a lower emissivity than the snowcover, have not been studied in detail as part of this project, and warrants further investigation in the future. 1This was buried in the snow, but may have been exposed during melt 164

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek Another interesting feature shown in figure 7.10 is the decreased variability in the forcing data between days 97 and 104. Comparison of the climate data at the snow pillow site, as shown in figure 6.12 with the climate data measured at the ridge site (figure 6.11) indicates that the pyranometer was probably covered in snow, in which case, the pyrgeometer was almost certainly also covered in snow. No thermal radiation data were available from other sites to provide forcing data during periods when this pyrgeometer was snow-covered. However, this error is relatively small, since the thermal radiation component from the sky is much less than the emission of thermal radiation from the canopy. The relative components of the energy balance is examined in greater detail in chapter 8, but summarised here, the simulated sub-canopy radiation shows extremely good agreement with the measured sub-canopy radiation, despite the errors in forcing data, because the downwelling thermal radiation is dominated by thermal emission from the vegetation. 7.4.4 Sub-Canopy Wind Speed The radiation component of the model has been tested fully and was found to simulate subcanopy observations with good accuracy. The canopy also affects the turbulent transfer of energy at the snow surface, which is suppressed by reduced wind speed beneath the canopy. An algorithm was used to estimate the sub-canopy wind speed from the work of Link and Marks (1999), which was described in section 4.3.2. Figure 7.11 shows a comparison of the open site wind speed, which is used to drive the model, the wind speed estimated by the algorithm and hence used to determine turbulent transfer of energy, and the true sub-canopy wind speed, which was measured by the fir climate station. 165

7.4. Fir Site;5imulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek. 4 - Link E Marks wind speed - Open site wind speed -3 i -- 3 Sub-canopy wind speed CD 2 rr F. 0" / L i 107 108 109 110 111 112 107 108 109 110 111 112 Julian day Figure 7.11: Comparison of sub-canopy wind speed estimated by the model with observations beneath a fir canopy Figure 7.11 presents a five day comparison of measured and modelled sub-canopy wind speed. The comparison for the entire period may be viewed from the CD that accompanies this thesis. The agreement between measured and simulated sub-canopy wind speed is not; high (r2 = 0.47), but is only slightly less than the agreement between wind speed observed at an open site and that observed beneath the canopy (r2 = 0.50). The maximum wind speed observed beneath the fir canopy throughout the season was less than 3 m s-l, therefore the turbulent energy transfer is expected to be quite low. Unfortunately no eddy correlation measurements were made during the RME 2000-2001 field experiment2, so the turbulent transfer component of this model cannot be tested fully at this stage. However, the error in simulated wind speed is expected to have a low impact on the model performance because the turbulent transfer of energy is low. This 2Although the loan of equipment was offered, no skilled operator was available to run the equipment 166

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek supposition was tested by a simulation of the fir site snow cover, driven by observed subcanopy wind speed. Snow depth simulated over the season is compared to the observed snow depth in figure 7.12. 0.7 0.6 -;. Measured 0.5 -— Simulated 0.4 - 0 -0.3 0.2 0.1 10 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Figure 7.12: Simulation of fir depth driven by observed sub-canopy wind Measured snow depth, and snow depth simulated with observed sub-canopy wind speed show excellent agreement (r2 = 0.99). A comparison of figures 7.5 and 7.12 demonstrates that the wind speed errors incurred through the use of the wind speed algorithm (equation 4.6), as opposed to forcing the model with observed sub-canopy wind speed, have little impact on the simulation of the snow cover. The reduction of wind speed beneath the forest canopy is an important component of the model, and the model performance is affected if this process is not taken into account. Figure 7.13 shows the simulated snow depth, when the wind speed was not reduced beneath the canopy. The agreement between simulated and measured snow depth 167

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek is less good (r2 into account. 0.98) if the reduction of wind speed beneath the fir canopy is not taken 0.7 0.6 0.5 E^ Q- 0.4 a) ~0.3 CD _0 ^:0. 3 0 c LN 0.2 0.1 0 4.5 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 Julian day Figure 7.13: Simulation of snow depth beneath fir canopy with no reduction of wind speed 7.4.5 Snowpack Profile Validation Sections 7.4.3 and 7.4.4 demonstrated that the model simulation of energy input to the snow surface is reasonable. In this section, the snow component of the model is tested by comparisons of snowpack profiles of temperature, density and grain size. Unfortunately, only one set of snow pit validation profile measurements were made at the fir site, because of the unusual shallowness of the snowpack in water year 2000-2001. The profiles simulated by the model are compared to the snow pit measurements in figure 7.14. Figure 7.14 shows that the model simulated the snow density profile to within the measurement error. The snow pit used to initialise the model was located at a similar 168

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek Figure 7.14: Comparison between simulated and measured temperature, density and grain size profiles 169

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek distance to the climate station, but diametrically opposite to the validation snow pit. The initialisation data required the addition of 10cm to the base of the snow pit so that the initial snow depth was consistent with the snow depth measured by the snow depth sensor, therefore the validation measurements may be offset from the simulated densities by a similar amount. However, the differences are not significant given the measurement error, and the general simulation of density profile appears to be reasonable from figure 7.14. The temperature profile simulated by the model at the fir site on day 75, as shown in figure 7.14, is similar to the temperature profile simulated at the open site, shown in figure 7.2. The simulated snowpack appears to retain liquid water, which was not observed in the field. Whilst insufficient data are available to determine whether, as at the open site, the grain size is underestimated by the model, observations in the field indicate large clusters of snow grains. Clusters of snow crystals are not modelled by SNOWCAN, so their effect on liquid water flow has not been considered here. The effect of melt-freeze cycles on snow crystal growth is also not explicitly simulated by the model. However, several melt-refreeze events were observed at the RME field sites during water year 2000-2001. Both of these physical processes affect the flow of liquid water in the natural snowpack, and would benefit from further investigation in the future. Evidence of water flow was observed in the lower 3-10cm of the fir snowpack on day 75. This suggests that drainage pipes may have been formed. SNOWCAN does not simulate liquid drainage via flow pipes, which could also explain why the modelled snowpack retains liquid water not observed in the natural snowpack. As discussed in section 7.2.3, the model assumption of constant irreducible saturation and choice of pa 170

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek rameter may be inadequate and could result in the retention of liquid water in the model, if the parameter is assumed to be too large. The simulated temperature profiles at the open and fir sites on day 75 have highlighted that improvements need to be made to the model in two areas to simulate the observed drainage of water: accurate simulation of grain size and clustering, and the formation of flow pipes in the snowpack. However, the simulation of snowpack temperature at the ridge site prior to liquid water flow is good, as demonstrated in figure 7.2. At the fir site, the simulation of the snowpack temperature profile throughout the season was tested with data from the thermistor array. A comparison between the simulated and measured temperatures is presented in figure 7.15. The simulated snow temperatures follow the trend of measured temperatures until the snowpack ripens, although the simulated temperatures are generally warmer than the measured temperatures. The difference between the continuous measured and simulated snow temperatures, as shown in figure 7.15 mainly results from the difference in location between the point simulation and snow pit measurement site. The distance between the two locations is less than 2m, but the snow depths at initialisation differed by 10cm. The snowpack at the snow temperature array location was simulated (i.e. the 10cm additional snow added to the base of the initial snow profile was removed), and the simulated temperatures at the snow-soil interface, 10cm, 20cm and 30cm from the snowpack base are demonstrated in figure 7.16. From figure 7.16, the model spin-up time appears to be of the order 3-4 days. It is more likely that the difference between the simulated and measured temperatures results from errors in the initial density profile rather than the Dirichlet lower boundary temper 171

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek lI I p. r-. g I" R F (a,.L. 1 -E % 0 r = I* I I I '- | I r r-? *,2?, r.. r.i r.i r \. r-.i ()')1 1 1, @IB~ U, iI C) LL ) '.11 1 c ip Lu! 1),U CL t "I I le I'll I r It, I T -0 ~o;I I.I, L. -1... I l. (i) wn~i~jdwuQ I, -I 'It N C' CI il e i LI. %7j 13 1 If r-. qIt.= p I. I II M CL/ t r (. r ~j N J. 2.. r-.~ u1 ~J '~ rLr Irl I I LU I'.6 1C \ I 1.J Figure 7.15: Simulated temperature profile comparison with thermistor array measurements. Measurements are given for the snow-soil interface, 10cm, 20cm and 30cm from the base of the pack. The shaded light blue represents the thermistor measurement error (0.5K). 172

7.4. Fir Site Simulation Chapter 7. Snow-SVAT Simulation II: Reynolds Creek.16J m.M4 ~21 8 rn -~r., I m 1 U-,J r-. '-4~ I',. LM W L44 II (,I);.mri4iv.i;.irJLLi;.ij i -IL LL Lu 4~4 -t- I':' I I I. L71 '.P. CU 4-, -n LL74 L44 LAl r — 41\ 3 '14 4-1. \t l l I I \ I0:Il l,f c 4 (U 14 2 - I I Di) iijiJwi 8 uri r — uLL r, LIw LL) C LU-.4i4.14 "In.1I 1i 2i rif 'V I VI (0 14a if I 44. L i LL LLCi - u~j u.I '.1 I d; %1 I I I '1' I,3 C..I r -. -IL. %Z (. J c "I Figure 7.16: Temperature profile simulation at the thermistor array location 173

7.5. Canopy Parameter Sensitivity Chapter 7. Snow-SVAT Simulation II: Reynolds Creek ature, because the difference appears immediately and is consistent at all measurement heights. Between the end of model spin-up time and the start of significant phase change to liquid (around day 64), the model simulates the observed temperatures well, predominantly within the measurement error of the thermistor array. The regression coefficients between days 51 and 64 are as follows: * Snow-soil interface: 0.93 * 10cm: 0.97 * 20cm: 0.95 * 30cm: 0.90 However, as discussed earlier, the snowpack does not simulate the observed rate of drainage of water, because the representation of grain geometry in the model does not reflect field observations adequately, because liquid water drains rapidly via flow pipes in the field, which are not simulated by the model and because the assumption of constant irreducible saturation may not be valid. The modelled temperatures remain at or just below the freezing point once melt has occurred, because of the liquid water retention by the model. 7.5 Canopy Parameter Sensitivity The model sensitivity to a ~10% change in canopy parameters was tested by comparison of simulated and observed snow depths at the RME fir site. The correlation coefficient and sensitivity functions for the change in parameterisation are presented in table 7.3. 174

7.6. Conclusions Chapter 7. Snow-SVAT Simulation II: Reynolds Creek Parameter [std] F+ F I F+ F- (r2) (r2)LAI [4.0] 92 140 54 130 0.9860 0.9868 av [0.25] 100 150 -27 140 0.9839 0.9866 Ev [0.977] 500 500 0.9648 Table 7.3: Sensitivity of model to the canopy parameters at RME fir site. For the simulation presented in section 7.4, r2 = 0.9857, Fm = 98 and Fs = 56 Table 7.3 shows that the model is relatively insensitive to the canopy leaf area index and to the canopy albedo. The model is more sensitive to a 10% decrease in canopy emissivity, but the uncertainty in this parameter is only - 2%. The canopy radiation component of the model, described in chapter 4, determines the sub-canopy radiation balance from the third exponential integral of the leaf area index. Given the shape of this function, it is likely that the model is more sensitive to the leaf area index at lower values of this parameter. 7.6 Conclusions This chapter described how SNOWCAN was used to simulate the snow cover at a fir site during both accumulation and ablation periods. Data collected as part of the Reynolds Mountain East Field Experiment 2000-2001 were used to initialise, drive and validate the Snow-SVAT. The model parameterisation was determined largely through the validation of the snow component mode at an open site, and additional canopy parameters were determined experimentally and from literature. Validation of the model at an open site presented in section 7.2 showed that the snow component model simulates the snow cover well, both in accumulation and ablation periods, despite the difficulty in the measurement of precipitation. The regression 175

7.6. Conclusions Chapter 7. Snow-SVAT Simulation II: Reynolds Creek coefficient between simulated and observed snow depths is 0.94. A comparison between simulated and measured density profiles indicates that the empirical compaction functions generally represent the densification of snow cover to within the measurement error, although the observed surface compaction was modelled less well, because the increased densification at high wind speeds is not considered by the model. The simulated snow temperature profile prior to a significant period of melt has excellent agreement with the temperature profile observed from a snow pit. However, after a period of snow melt (e.g. Julian day 75), the model retains liquid water, which was not observed in the field. A comparison between simulated and observed grain size profiles indicates that the model underestimates snow grain size by as much as a factor of 2, and thus the rate of liquid water flow is underestimated by up to a factor of 4 (dependent on the residual water and thermal deficits of the snowpack). At the fir site, the model simulation of snow depth improved from the open site simulation, and the regression coefficient between simulated and observed snow depth is 0.99. The model simulated the sub-canopy radiation balance well, and compared to the upwelling and downwelling solar and thermal radiation with high r2 coefficients, despite measurement errors caused by canopy shadows. Simulated temperature profiles at the fir site, generally agreed with continuous measurements from a thermistor array to within the measurement error, between the end of model spin-up time (4 days) and the first appearance of significant liquid water. Similarly to the open site, the model retained liquid water. Field observations suggested that rapid liquid water drainage in the field through flow pipes may contribute significantly 176

7. 6. Conclusions Chapter 7. Snow-SVAT Simulation I: Reynolds Creek to the discrepancy between observed and simulated temperature profiles. The Snow-SVAT developed as part of this project has been shown to simulate the seasonal snow cover well, over a range of meteorological conditions. The accurate representation of physical processes by SNOWCAN means that this model can now be used to address questions related to the heterogeneity of the field site, and to the energetic contributions to snow melt. The use of SNOWCAN as a diagnostic tool is investigated in chapter 8. 177

Chapter 8 Snow-SVAT Simulation III: Diagnostic Tool 8.1 Introduction In this chapter, the validated Snow-SVAT is used to address general questions regarding the effect of the forest canopy on the snowcover beneath it. Data from both the BOREAS campaign and from the Reynolds Mountain East Experiment 2000-2001 will be used to expand on the studies carried out in chapters 5 and 7, and to examine the spatial variability in snow depth that results from spatial variability in the canopy density. Specifically, the following questions will be addressed: 1. Does the snow beneath the forest melt faster than the snow at an open site? 2. How do the relative contributions of energy terms compare between the forest and open site? 178

8.2. Surface IEnergy Balance Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool 3. Could a variation in leaf area index change the dominant radiative energy term? 4. Under what circumstances would the opposite scenario to question 1 occur? This chapter cannot answer these questions definitively, as that would require substantial amounts of validation data from numerous sites. However, the data already available and the Snow-SVAT model can illustrate a line of approach, and indicate possible conclusions in the future. 8.2 Surface Energy Balance Chapter 1 discussed observations that demonstrated that the snow beneath a forest melted before the snow at an open site, or vice versa. In this section, SNOWCAN is used to look at the particular case of Reynolds Mountain East, where the snow beneath the forest does indeed melt before the snow at an open site, but conditions that could cause the snow at an open site to melt before the snow at a forest site are discussed in section 8.4. There are two main causes of more rapid melt at the forest site: the relative energy balances of open and forest sites, and mass difference by interception processes. A comparison between the energy balance at the open site and beneath the fir canopy is given in section 8.2.1, and the role of snow interception by the canopy in the site melt order is discussed in section 8.4.1. 8.2.1 Inter-Site Comparison Figure 8.1 shows the model cumulative radiative and turbulent energy contributions for the Reynolds Mountain East ridge and fir sites. Several points of interest are demonstrated by figure 8.1, 8.2 and 8.3. The melt period from day 77-85 is additionally demonstrated 179

8.2. Surface Energy Balance Chapter 8. Snow-SVAT Simulation IIl: Diagnostic Tool 5.0E+08 4.0E+08 3.0E+08 E 2.0E+08 x 1.0E+08 U., O.OE+00 I > -1.OE+08 7 -2.0E+08 E U -3.0E+08 -4.0E+08 -5.OE+08 45 YZ7.>- Net Solar Radiation Net Thermal Radiation Sensible Heat... ---..i.. Latent Heat —... Total Radiative and Turbulent;....x 0 5 60 6 70 7 80 8 0 9 00 0 0 1 20 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 Julian day................ -- -................ - -....... (a) Open Site 7.0E+07 -.. 6.0E+07 -I 5.0E+07 -E 4.0E+07 L _ 3.0E+07 LL.i 1, 2.OE+07 I, 1.OE+07.-. X O.OE+00.... 3 | -1.0E+07 -2.OE+07 -t -3.OE+07 - -. ---. Net Solar Radiation - Net Thermal Radiation -- Sensible Heat.-b —....Latent Heat - Total Radiative and Turbulent..................................................................... 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115; A-,., 120 (b) Fir Site Figure 8.1: Cumulative Energy Contributions at RME Ridge and Fir Sites. Open site ordinate scale is an order of magnitude greater than for the fir site. Melt periods can be seen from snow depth measured and simulated at these site in figures 7.1 and 7.5. 180

8.2. Surface Energy Balance Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool Figure 8.2: Relative contributions of radiative and turbulent energy contributions at RME ridge site 181

8.2. Surface Energy Balance Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool Figure 8.3: Relative contributions of radiative and turbulent energy contributions at RME fir site 182

8.2. Surface Energy Balance Chapter 8. Snow-SVAT Simulation II: Diagnostic Tool Total Radiative and Turbulent 750.................................................................................................................................................................................... 550 -" I 350 -150 -50 L I450 77 78 79 80 81 82 83 84 85[ Julian day Net Solar Radiation...................................................................................................................................................................................................................................................................... I N 300 -.E 200 100 -100.,oo I| -200 77 78 79 80 81 82 83 84 85 Julian day Net Thermal Radiation N 0 30 i E 2. - i 1 7( 7 - 2100 -l1t(1 77 78 79 80 81 82 83 84 85 Julian day Sensible Heat 1E 2(00 l -oo.................................................................................................................................................................. --—................. 20X....................... —..................................................... 77 78 79 80 81.2 83 84 85 Julian day Latent Heat 100 - 100 2: - oo.................................................................................................................................................................................................................................................................... 7;' 78 79 80 81 82 83 84 85 Julian day Figure 8.4: Relative contributions of radiative and turbulent energy contributions at RME ridge site: melt period days 77-85 183

8.2. Surface Energy Balance Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool JierSoRadiaton -25 77 78 79 80 81 Julian day 82 83 84 WL-t Tarsa-.l D.A.I#-,r 75 50 -25 N E i - o -.... -.. I.......-......:........................................................... -25 - -- - --............... 77 78 79 80 81 Julian day UL_....82 83 84 85 82 83 84 85 Sensible Heat 75 50 (x^, 25 E I o -r' -25 '....I....................................t... 77 78 79 80 81 Julian day 82 83 84 85 Latent Heat 75...... 50 E 0 0 -25 --- 77 -- --... - --, - - - - "' — " — - ---— j - '' - -- - -, LI, - i, -'~la 78 79 80 81 Jutian day 82 83 84 85 Figure 8.5: Relative contributions of radiative and turbulent energy contributions at RME fir site: melt period days 77-85 184

8.2. Surface Energy Balance Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool in figures 8.4 and 8.5. Firstly, the magnitudes of energy terms are much greater at the open site, and approximately three times as much energy is required to melt the open site snowpack. Secondly, the higher wind speeds measured at the ridge site resulted in a higher rate of energy transfer from convection than from solar heating. In addition, the heat lost by thermal emission from the snowpack exceeded the heat gain from solar heating. The periods of snowmelt at the open site, days 77-85, 106-109 and 112-120 from figure 7.1 on page 144 were largely driven by sensible heat transfer. This contrasts with snow melt at the fir site, where the snow melt periods on days 64-67, 77-84, 105-109 and 112-116, as shown in figure 7.5 on page 155 were mostly driven by thermal heat transfer, with a significant contribution from sensible heat. Chapter 3 discussed the model uncertainties, which were relatively high for all aspects of the turbulent transfer of energy. These uncertainties are more important for the open site simulation, where the snow melt is driven by the sensible heat transfer. This is also limited by the formulation of turbulent transfer, and is probably the cause of the less accurate simulation of snow depth at an open site (r2 = 0.94) compared to the simulation of snow depth at the fir site (r2 = 0.99), as reported in chapter 7. This accuracy in simulated snow depth was achieved at the fir site, despite icing problems with the pyrgeometer that provided the thermal radiation forcing data for the model. Given that at the fir site the thermal radiation dominated the provision of energy for snow melt, it is remarkable that the errors in forcing data did not appear to significantly affect the model performance. This is explained in the next section. 185

8.2. Surface Energy Balance Chapter 8. Snow-SVAT Simulation II: Diagnostic Tool Ri' 400 E o - Total -50 150 Atmospheric iC 50 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 12C Julian day I 4 55 0 - -- —............. -.................................................................................................................................................................................................................................................... 400 E 350:300 on 250 3c200 a-.T otal 150 ~-"'~Atmospheric (reflected from snow surface) c ---- Vegetation (reflected from snow surface) Y100 now -..........______............ Sno 50 0 -45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120i Julian day Figure 8.6: Relative contributions of thermal energy from the snow, canopy and sky 186

8. 3. Canopy Variability Chapter 8. Snow-SVAT Simulation II: Diagnostic Tool 8.2.2 Thermal Contributions From the Sky, Vegetation and Snow The relative thermal contributions from the sky, vegetation and snow are presented in figure 8.6. It can be seen from figure 8.6 that the downwelling sub-canopy thermal radiation is dominated by the thermal emission from the canopy, whereas the upwelling sub-canopy thermal radiation is dominated by thermal emission from the snow surface. The atmospheric thermal radiation contribution to the thermal radiation balance beneath the fir canopy at this site is small, and for this reason, equipment icing errors in the measured thermal radiation used to drive this model had little impact on the model performance. 8.3 Canopy Variability In section 8.2, the thermal radiation balance was shown to be much more important than the solar radiation component beneath the RME canopy. However, as described in chapter 1, many studies have suggested that canopy shielding from solar radiation is the dominant effect of a canopy on the snow beneath it. In this section, the canopy leaf area index dependence of the relative contributions of solar and thermal radiation to snow ablation was investigated, to determine at what value of LAI the solar radiation input becomes greater than the thermal radiation input. 8.3.1 LAI-Dependent Ratio of Sub-Canopy Radiation Components The net solar and thermal radiative flux at the snow surface is shown in figure 8.7 for the second melt period at the fir site (Julian days 77-84). Figure 8.7 also demonstrates the radiative energy contributions for a more sparse canopy (LAI = 1.0). As can be seen from figure 8.7, the radiative energy input is dominated by shortwave radiation for a thin 187

8.3. — Canopy Variability Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool LAI = 1.0 200 - - - Net Solar - - - - Net Thermal 2 150 E B x 100 Ua 5 z / I 1\ I I ~ I I I \_ _ivx_ I -- - - - -- --, -- - - p s wwwr —j — 7tr -- I Ir"ll- -I 7 --— I --- —81 W i4 V 4.......... J ~....... ' ' i-_ -50 Julian Day LAI = 4.0 200 -. ' 150 X 100 - * 50. n:.- 04 4J z 7 z 7 -- Net Solar - - - - Net Thermal I IJAJJ.k 46 ~ -.! I I. I!,. 7 ------- 82 83 78 79 80 81 J84 -50 -L Julian Day Figure 8.7: LAI dependent net thermal and net solar radiative energy during one melt period at a fir site (Julian days 77-84). 188

8.3. Canopy Variability Chapter 8. Snow-SVAT Simulation IIl: Diagnostic Tool canopy, and by longwave radiation for a dense canopy. One point of interest is that there appears to be less lag between net solar and net thermal at higher LAI. This arises partly from the decreased atmospheric contribution to the thermal energy beneath a more dense canopy, and partly from the model assumption that the canopy temperature is equal to the air temperature. The LAI dependence of the ratio of total thermal energy input during this melt period to the total solar energy input is shown in figure 8.8. The assumptions of the model break down for a discontinuous canopy, and therefore should not be used for canopies with LAI <1.0. The ratio of thermal and solar radiative energy for a canopy with LAI = 0.5 is only shown in figure 8.8 to illustrate the general trend. 3 2.5 0 0 2 1.5 4- 1 0LAI -0.5 -1 -1.5 LA}................ 3 3.5 4 4.5 - ---— --- —— — — --- —------------ " "' " Figure 8.8: Variation of thermal to solar radiative energy ratio with leaf area index. This is time and location specific. Figure 8.8 shows that under these meteorological conditions, the daytime thermal 189

8.3. Canopy Variability Chapter 8. Snow-SVAT Simulation II: Diagnostic Tool heat exceeds nocturnal radiative heat loss for canopies with leaf area index greater than 1.6 for this melt period. For canopies with LAI greater than 2.7, the thermal energy input to the snowcover is greater than the solar energy contribution. 8.3.2 Variation of Snow Depth 2 o...................................................................................................................... -...... 4 -— Maximum Difference 5..... Average Difference 15j E | Minimum Difference C 10 I:.; 5 -................................................. ~....... -!0. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 LAI Figure 8.9: Simulated variation of snow depth from canopy variability The effect of canopy variability on the radiative energy balance was examined in section 8.3.1; the variation in the surface energy balance gives rise to variability in snow depth. Figure 8.9 shows the maximum, minimum and average snow depth difference from that measured at the RME fir site for the season simulated in chapter 7. The variability in snow depth. shown in figure 8.9 highlights variability from the canopy radiative effects only, and the natural snow depth variability is greater, as discussed in chapter 1, because of ground topography and the interception of snow by the canopy (which must depend on canopy structure), among other forest effects. As stated previously, the model is not valid 190

8. 3. Canopy Var~ability Chapter 8. Snow-SVAT Simulation II: Diagnostic Tool for LAI < 1.0, but the snow depth difference is shown in figure 8.9 for a canopy with LAI = 0.5, for illustrative purposes. At the RME fir simulation point location, the LAI was measured in the field to be 4.0, and figure 8.9 indicates that this is the best model parameter value. However, the differences between simulated and observed snow depths are anomalous for LAI = 4.5. The simulated snow depth is greater that simulated for a canopy of LAI = 4.0, despite an increase in energy to the snowpack surface. The cause of this anomaly is not immediately apparent, but may be affected by the reduction of solar absorption within the pack with greater canopy shielding. This argument is, however, counter-intuitive to the simulations presented in chapter 3. Above a canopy LAI of 4.5, the thermal energy input reaches an asymptote in the model, since the thermal radiation almost entirely originates from canopy emission. 8.3.3 Above-Canopy Albedo In addition to the investigation into the variability in snow depth, SNOWCAN can also be used to examine the variation in the albedo measured above the canopy. The above-canopy albedo is an extremely important variable in climate models, as it is used to estimate the solar energy input to the land surface. Above the Old Jack Pine site in the Canadian boreal forest (LAI = 2.5), which has been studied in chapter 5, BOREAS investigators measured the above-canopy albedo to be 0.15 [Betts and Ball (1997)]. Although Betts and Ball (1997) looked at seasonal variation in above-canopy albedo, they did not address the diurnal variation. This is important, as an average value may not be representative, because the incident solar radiation flux is not constant, the direct beam irradiance angle 191

8.3. Canopy Variability Chapter 8. Snow-SVAT Simulation IIl: Diagnostic Tool varies during the day, and canopy structure is important. In this section, the simulated albedo' above the RME fir canopy during the RME 2000-2001 field experiment is examined to investigate the magnitude of diurnal, seasonal and canopy variability effects. LAI = 1.0 0.35 - o.0 0.3 ~ 0.25 U & 0.2 < 0.15 bnucnh~~ t t - I t; 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 Julian day LAI = 4.0 0.15 0.148 o 0.146 s0 0.144 >,0.142 O 0.14 5 0.138 - ' 0.136 0 0.134 0.132 n 1 - l R ' L... l.. >.,, i.. III i -. I - 1- --- 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 Julian Day Figure 8.10: Seasonal and diurnal variation in above-canopy albedo Figure 8.10 shows that the variability in above-canopy albedo increases with 1The above-canopy albedo includes radiation that has been reflected from the snow surface and has penetrated back through the canopy 192

8.3. Canopy Variability Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool decreasing canopy cover. This is because the albedo of the underlying snowcover has a greater influence on the above-canopy albedo. At low canopy density, the above-canopy albedo is a maximum around the middle of the day, and decreases at low solar angles because of strong forward scattering at the snow surface. At higher canopy density, the above-canopy albedo is a minimum around midday and increases at lower sun angles, because the path length through the canopy is greater, and more radiation is intercepted by the canopy and is reflected from the canopy elements. The above-canopy albedo clearly decays between precipitation events, and is low during melt periods. For a canopy with a leaf area index of 3.0, the above-canopy midday albedo is a minimum during melt periods but a maximum otherwise. This points to a threshold canopy LAI, above which the above-canopy albedo is dominated by canopy effects rather than the snow albedo. For a canopy with a leaf area index of 2.5, such as the boreal jack pine forest that has been studied in chapter 5, the 'average' albedo (constant weighting for all times of day) simulated with the RME data is 0.160, and ranges between 0.150 and 0.171. The average simulated albedo is slightly higher than that measured by the BOREAS investigators. In addition, the average above-canopy albedo does not give a true estimate of the amount of solar energy input to the surface, because the radiation input is greater in the middle of the day. A better estimate of the average above-canopy albedo could be obtained by a diurnal cycle weighted average. For low canopy densities, the weighted average would raise the average albedo, but the average albedo would be lowered for high density canopies. It is apparent from this study that the above-canopy albedo varies with canopy density and time between precipitation events. SNOWCAN could be used to develop a 193

8.4. Reversal in Site Melt Order Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool model to describe the above-canopy albedo, for use in general circulation models. However, at this time, the lack of validation data would limit confidence in the accuracy of such a model. The use of a fixed daily albedo in climate or ecological models is clearly an unacceptable assumption, even though it is done regularly. 8.4 Reversal in Site Melt Order Investigations into the comparative energy balance at an open site and beneath a fir canopy at RME revealed that the cumulative energy input to the snowcover is much greater at the open site than at the fir site, yet the snow beneath the fir canopy melts before the snow at the open site. SNOWCAN has been used to probe whether the order of site melt arises from interception losses, as described in section 8.4.1. 8.4.1 Site Mass Differences In order to investigate the causes of complete fir site ablation prior to open site ablation, the RME ridge site simulation was repeated, with the addition of a canopy of LAI = 4.0. The wind speed beneath the canopy was also reduced, but all precipitation was assumed to fall through the canopy, so that the mass input was the same as for the open site simulation presented in chapter 7. The simulated snow depth at the open site but with radiative and turbulent canopy effects of the fir canopy included is shown in figure 8.11. Figure 8.11 demonstrates that the snow still melts beneath the canopy prior to the open site, which indicates that the difference does not arise from the mass input difference. However, the wind speed at the ridge site is higher than at the snow pillow site, which was used as forcing data for the fir site in chapter 7. The turbulent energy input in this 194

8..Reversal in, Site M~elt Order Chapter 8. Snow-SVAT Simulation III: Diagnostic Tool 1 0.9 0.8 0.7 E 0,6,r 0.5 i 0.4 C v 0.3 0.2 0.1 0 30 40 50 60 70 80 Julian day 90 100 110 120 Figure 8.11: Simulated snow depth at RME ridge site, with fir canopy radiative and turbulent effects added, compared to the snow depth measured at the open site. simulation is higher than for the fir site simulation in chapter 7, and the turbulent and radiative energy input to the snowpack shown in figure 8.11 is approximately 12% higher than the turbulent and radiative energy input at the open site. The total amount of turbulent and radiative energy required to melt the open site snow pack is 1/3 higher than would be required to melt the same snowpack beneath a canopy (interception effects excluded). As all other inputs and conditions are equal, the difference must arise from the handling of mass and energy within the snowpack and from the ground heat flux. Chapter 3 highlighted the effect of heat transport and radiation absorption within the pack and at the soil surface. Here beneath the forest, a greater proportion of energy is absorbed at the surface, since less solar radiation is incident on the snow surface, and subsequently less solar radiation is absorbed within the snowpack or at the soil surface. As described in chapter 3, this has two effects: (i) for thin snowpacks, solar 195

8.4. Reversal in Si'te Melt Order Chapter 8. Snow-SVAT Simulation Ill: Diagnostic Tool energy that may have been absorbed by the soil is reduced, and a greater proportion of energy is absorbed within the snowpack, and (ii) heat is more easily transported downwards through the snowpack during the melt period, by gravitational drainage of liquid water in addition to diffusion and conduction, than upwards by diffusion and conduction alone. Greater absorption of energy at the surface of the fir snowpack could explain why less energy is required to melt the same mass of snow. A similar investigation has been carried out in chapter 5, where the canopy radiative and turbulent effects were removed. A comparison of figures 5.4 and 5.5 suggests that the snow beneath the canopy melts before that at an open site, and that the differences do not arise from mass input differences. 8.4.2 Other Causes The question still remains as to the processes that could cause the snow at an open site to melt faster than the snow beneath a canopy. One answer could lie in the meteorological conditions. If an even greater proportion of incident radiation arose from the shortwave rather than longwave, such as under cold conditions at a site located further south (where the days are longer), then the shielding effect of the canopy may more than compensate for the energy gain from increased canopy thermal emission with greater canopy density. However, increased incident solar radiation will raise the temperature of the canopy and hence increase thermal emission of the canopy, therefore increased incident solar to thermal incident radiation is unlikely to change the site melt order. A similar argument can be applied to the atmospheric optical depth. These effects cannot be investigated until the canopy temperature is modelled in SNOWCAN. 196

8. 5. Conclusions Chapter 8. Snow-SVAT Simulation Ill: Diagnostic Tool Snow deposition mechanisms may affect the order of site ablation. At a forest site, if less precipitation is lost by sublimation from intercepted snow, then the forest snowpack will persist for longer. This could occur in warmer or windier environments, where the canopy interception is removed more quickly, or for thinner canopies, which may intercept less snow. SNOWCAN would benefit greatly from the inclusion of a physically based canopy interception model, since the assumptions made by this thesis may not be applicable in other environments. One physical process not considered by this thesis is blowing snow. At sites where this process is significant, the snow mass at an open site could become less than that beneath a canopy, even if the initial precipitation is the same. It would be extremely interesting to address these questions and assumptions further, through SNOWCAN simulation and validation from experiments similar to that conducted at Reynolds Mountain East, but at a location where the snow at an open site melts prior to the snow beneath a canopy. 8.5 Conclusions This chapter extended the use of the model from simulation at a specific site, such as Reynolds Mountain East in chapter 7, to hypothetical sites with different canopy characteristics. Comparison of the energy balance at the RME ridge and fir sites has shown that the ridge site snowmelt is driven by sensible heat transfer, whereas the fir site snowmelt is driven largely by thermal emission from the canopy, although the contribution from sensible heat transfer is non-trivial, despite low wind speeds beneath the canopy. The distinction between surface energy absorption and solar absorption within the snowpack 197

8.5. Conclusions Chapter 8. Snow-SVAT Simulation II1: Diagnostic Tool has been shown to be important. This method of approach has allowed questions to be addressed regarding the relative melt rates at open and forested sites. SNOWCAN has also been used to investigate the effects of diurnal and seasonal variations and canopy variability on the above-canopy albedo, which is extremely important for general circulation models. Although these discussions are purely theoretical at this stage, this chapter has shown how SNOWCAN, or a similar type of model, could be used as a tool to probe the balance of the physical processes operating in the environment, and to test larger scale model simplifications. However, before these questions can be answered, more data are required to validate SNOWCAN under a number of different canopy and meteorological conditions. Improvements should also be made to SNOWCAN to improve the representation of some of the physical processes in the model. Chapter 9 highlights directions for future work, discusses the model limitations and examines the model as a platform for improved physicsbased modelling. 198

Chapter 9 Conclusions 9.1 Introduction Depletion of snowcover brings about crucial changes in the land surface, through flux of water into the soil and across the land surface, and by an increase in solar energy input by a decrease in the surface albedo. This in turn leads to an increase in biological activity and an increased CO2 and water exchange with the atmosphere. Accurate modelling of snowmelt has many implications, from local water management to climate modelling and prediction. Snowcover simulations have previously been shown to represent observed snowpacks over open areas, but relatively little research has been carried out in forested environments. Large areas of forested snowpacks, such as the boreal forest, have great impact on the global energy balance. This thesis has examined the processes affecting the snow energy balance under a forest. In particular, the project focussed on the interaction of radiation between the canopy and snowcover. In this final chapter, the project achievements are summarised in section 9.2. 199

9.2. Project Accomplishments Chapter 9. Conclusions Chapter 1 stated several propositions, and these are considered in section 9.3. The limitations of the model and further work to improve the representation of physical processes by the model are discussed in sections 9.4 and 9.5 respectively. Further use of the model as a diagnostic tool and to improve existing models is approached in section 9.5. 9.2 Project Accomplishments Eight major achievements were accomplished in this project: 1. A snow model, SNTHERM, was recoded, underwent a change of style and new algorithms were included, to form the new snow model, SNOWCAN. 2. The sensitivity of SNOWCAN to the model parameters was investigated, under constant, high energy meteorological conditions. Several parameters were highlighted as particularly important because the model showed high sensitivity, because the uncertainty in the parameter was high or because the model was expected to be more sensitive under more general meteorological conditions. These parameters were the windless exchange coefficients, two compaction coefficients, the bulk extinction coefficient for near-infrared radiation, the soil quartz content and the atmospheric optical depth. In addition, the model was sensitive to the initial snow density profile and the lowest soil layer temperature (Dirichlet boundary condition). 3. The physically based snow model was coupled to a physically-based optical and thermal canopy radiation model to form the basis of a Snow-SVAT. Other canopy effects were also incorporated into SNOWCAN. 4. SNOWCAN was tested with data taken during the BOREAS campaign. The snow 200

9.2. Project Accomplishments Chapter 9. Conclusions depth simulated beneath a boreal jack pine forest has shown good agreement with observations. However, insufficient data were available to validate the model. 5. A new experiment was designed and carried out at the Reynolds Creek Experimental Watershed, Idaho, USA during water year 2000-2001. A complete set of data was taken to initialise, drive and validate SNOWCAN. 6. Data from the RME experiment were used to simulate the evolution of the seasonal snowcover at an open site. Validation data at the ridge site were used to determine the most appropriate compaction parameters and to test the accuracy of the standalone snow model. Validation of the model against measurements of grain size profiles highlighted an underestimation of grain growth by the model. Consequently, liquid water was retained by the snow in the simulation, but was not observed in the field. 7. The accumulation and ablation of the snowpack beneath the RME fir canopy has been simulated from data taken during the RME 2000-2001 field experiment. The model has been tested and validated with measurements of sub-canopy up- and downwelling solar and thermal radiation, snowpack profiles of temperature, density and grain size, and automatic measurements of snow depth. Simulations showed good agreement with observations, although the results again highlighted retention of liquid water in the simulated snowpack, that was not observed in the field. At the fir site, however, observational evidence of water flow suggested an increase in liquid water flow, through the formation of flow pipes. 8. Validation of SNOWCAN has shown that the representation of physical processes by 201

9.3. Consideration of Propositions Chapter 9. Conclusions the model is reasonable. SNOWCAN therefore has been used to investigate physical conditions outside the limits of the RME data. A comparison of the energy balances beneath the RME fir canopy and at the ridge site has shown that the fir snow melt was dominated by thermal emission from the canopy, and that complete ablation occurred prior to that at the ridge site, where the snow melt was dominated by sensible heat transfer. As a result, the simulation of snow depth was less accurate at the ridge site, because the turbulent transfer of sensible heat is less well understood than thermal processes. SNOWCAN has also been used to investigate the effect of canopy variability on snow depth, and the variability of above-canopy albedo. 9.3 Consideration of Propositions A number of propositions were stated in chapter 1, and have been approached throughout this thesis. These propositions are discussed below. Snowmelt beneath the forest is dominated by thermal radiation A simulation of the accumulation and ablation of the snowcover beneath a fir forest within Reynolds Creek Experimental Watershed showed that the dominant energy input, which caused the snow to melt, was the downwelling thermal radiation. An additional study concluded that the atmospheric thermal contribution was small, and that the thermal emission from the vegetation provided most of this energy. This was supported by measurements of the sub-canopy up- and downwelling thermal radiation, which were closely followed by the simulated radiation. The simulation also revealed that the turbu 202

9.3. Consideration of Propositions Chapter 9. Conclusions lent transfer of sensible heat provided a significant energy input, despite low wind speeds observed and simulated beneath the canopy. However, an investigation into the change in radiative energy input with canopy leaf area index, which focussed on one melt period, suggested that the solar energy input exceeded the thermal contribution below a leaf area index of 2.7. In addition, this study has shown that the nocturnal radiative cooling of the snowpack exceeds the daytime thermal heat gain for canopies with LAI < 1.6, for the melt period studied. Further investigation would be beneficial to determine whether these LAI thresholds can also be applied to different sites and under different meteorological conditions. Hence, whilst this proposition is true for the site and meteorological conditions studied, it may not be true under different conditions. This proposition is probably not true for thinner canopies, where either the solar radiative energy or turbulent sensible heat exchange may dominate. Complete ablation of the snowcover beneath a fir canopy occurs prior to complete ablation at an open site Observations made as part of the RME field experiment 2000-2001, and the validated simulations of the accumulation and ablation of the seasonal snowcovers at a ridge and a fir site have shown that this proposition is true under the conditions studied. It has also been shown that the total radiative and turbulent energy input to the fir snowpack is less than the total radiative and turbulent energy to the ridge snowpack. Whilst the difference partly arises from mass loss at the fir site as a result of sublimation of intercepted snow and from the difference in ground heat flux, the energy absorption 203

9..3. - Consideration of Propositions Chapter 9. Conclusions within the pack also has an impact. At the fir site, a greater proportion of energy is absorbed at the surface, and less energy input is required to melt the snowpack. However, the snow at an open site, could melt before the snow beneath a canopy because of topographical effects, and differences in mass input, where blowing snow may reduce the mass at an open site relative to that beneath the forest. The question as to whether the snow beneath a given forest will melt before that at an open site is open to debate, and would benefit greatly from further investigation. Natural canopy variability leads to significant differences in snow depth, from changes in the radiative energy input Figure 8.9 showed that the snow depth may vary by up to about 20cm, from radiative energy c hanges with canopy variability. However, the variation in leaf area index that causes this magnitude of change is much greater than can be expected in the field, and would demand a different type of model to represent it, as described in section 9.4. The relatively high variability in snow depth observed in the field is more likely to result from topographic and interception processes than from radiative differences. This postulate is rejected, because the observed canopy variability does not account for the observed variability in snow depth. |Changes in the snow albedo have a significant effect on the above-canopy albedo The decay of snow albedo has an obvious effect on the above-canopy albedo, which also decays between precipitation events. For a canopy with high LAI, such as 204

9.3. Consideration of Propositions Chapter 9. Conclusions at the RME fir site, however, the above-canopy albedo only changes by <0.01. For this canopy, the albedo decay has a greater effect than the diurnal variation, which is dominated by canopy radiation interception. As the LAI decreases, the scattering of radiation by the snow becomes more important, and dominates around LAI of 3.0. The above-canopy albedo of a canopy with LAI < 2.0 is more sensitive to diurnal variation than to the decay of snow albedo. Changes in the snow albedo do have a significant effect on the above-canopy albedo. The decay in snow albedo affects the above-canopy albedo, even for a relatively dense canopy. Thinner canopies are affected by the diurnal changes in snow albedo, and the effects may be larger than the albedo decrease with snow aging. Above-canopy albedo exhibits strong seasonal and diurnal variation, and cannot be represented with a single, constant value in general circulation models As described above, the diurnal variation in above-canopy albedo is small for dense canopies, but increases with decreasing canopy density. However, seasonal variation in above-canopy albedo arises from the decay of snow albedo between precipitation events. At high canopy density, a constant above-canopy albedo could be used in general circulation models, if an uncertainty of 0.01 is acceptable. A much better approach would be to account for both diurnal and seasonal variation in above-canopy albedo, and a physically based Snow-SVAT, such as SNOWCAN, could be used to develop an above-canopy albedo model. However, the effect of intercepted snow on the above-canopy albedo has not been considered here, and would need to be studied in greater detail. 205

9.4. Model Limitations Chapter 9. Conclusions 9.4 Model Limitations Despite the excellent performance of SNOWCAN, there are a few situations where use of the model may not be appropriate. The assumptions made about the canopy structure mean that the model cannot be used for a discontinuous canopy (LAI < 1), where cavity effects may dominate. This may also be a problem where the canopy exhibits a high degree of spatial variability. The postulates discussed in section 9.3 cannot be treated as hypotheses, because of the lack of validation data. The ability of the model to describe the physics under very different canopy conditions has not been proven. Although the model simulated snow depth well beneath a boreal canopy (with LAI = 2.5), snow depth is governed by snow compaction, which is represented empirically, and has been calibrated at particular sites. It would be beneficial to validate the model beneath a canopy at the limit of model assumptions i.e. LAI around or just above 1.0, by a similar approach to that used for the RME fir canopy. At present, this model has not been tested for use as an operational tool. The Dirichlet lower boundary energy condition used by this model, and the initialisation conditions make this model inappropriate for operational use, although data assimilation techniques could be used in the future to move this SNOWCAN towards operational status. Before any of these limitations are addressed, however, a number of improvements could be made to the model in its current form. These improvements are highlighted in the following section. 206

9.5. Model Improvements Chapter 9. Conclusions 9.5 Model Improvements Although the model can be taken in many directions, several areas where the model would benefit from improved representation of the physics have been discovered during the course of this project. Three of these relate to the snow model alone, and four are concerned with the Snow-SVAT as a whole. Validation of the snow model against field measurements of temperature and grain size has indicated that the model underestimates the rate of snow grain growth. This error not only affects the absorption of solar radiation within the pack, but also the rate of liquid water drainage. Comparison with temperature profiles has shown further evidence of this, since the snowpack retained liquid water not observed in the field. Liquid water drainage and heat transport within the snowpack are governed by the snow densification, which is represented empirically within the model, and is calibrated between sites. Snow densification partly arises from growth of snow grains, but these processes are not linked within SNOWCAN. This model would benefit greatly from a new approach to the physics of snow densification. At the fir site, the model again showed that the simulated snowpack retained liquid water, which was not observed in the field. However, field observations suggested that the formation of flow pipes enhanced drainage in the natural snowpack. SNOWCAN does not model rapid liquid water flow via drainage pipes, and more research is required to represent this important mechanism for water removal from the snowpack. Snow depth simulated beneath a fir canopy in chapter 7 has shown better agreement with observations than the snow depth simulated at the nearby ridge site. An energy balance investigation has shown that snow melt at the ridge site was dominated 207

9.5. Model Improvements Chapter 9. Conclusions by turbulent transfer of sensible heat, in contrast to the fir site, which was dominated by thermal radiation. The turbulent transfer mechanism is less well understood than other mechanisms of heat transfer, particularly over a smooth surface, such as snow. Further investigation into this energy transfer mechanism beneath the canopy would be beneficial, where the turbulent transfer mechanism has not been treated rigorously in this model. SNOWCAN would benefit from a canopy temperature model linked to the radiative and turbulent energy absorption within the canopy. As stated earlier, however, the canopy effect on atmospheric stability and turbulent transfer over snow needs to be studied further. The coupling of an additional canopy interception model would reduce uncertainty in the mass input to the forest surface. These two improvements probably need to be studied together, since the canopy temperature would affect the amount of intercepted snow, and vice versa. The additional representation of CO2 and methane transport in this Snow-SVAT would be beneficial for climate modelling. One other effect that should be incorporated within SNOWCAN is the effect of forest debris on the snow albedo. This process was considered to be negligible at the Reynolds Mountain East site during water year 2000-2001, because photographs of the surface taken towards the end of the ablation period indicated that litter deposition was relatively low. This assumption may not applicable for other sites, where this model could be used as a diagnostic tool. Inclusion of these improvements will lead to a good physical representational of all the components of a Snow-SVAT. These extensions would also probably reduce the number of parameters, as the processes in the snowpack are related and the extended model would take advantage of this. After testing with a variety of field data, SNOWCAN can then be 208

9.5. Model Improvements Chapter 9. Conclusions used to produce simplified, yet physically-based models for GCMs and operational snow modelling. This in turn will be a great advance over existing empirical representations and allow predictions, for instance, of possible regional and global changes, with far greater confidence. 209

Appendix A Code Inconsistencies In this section, minor code changes made either directly or indirectly (through the change in program style) to the SNTHERM code to resolve inconsistencies are described below. The relevant subroutine or function is given in square brackets. * Specific gravity g, defined as L was changed to 1.0905d0 [FILTRATE] Pi * A safety condition recalculates effective saturation for a thin surface node. Determination of liquid water flow (U = fn(Se)) was added following this safety condition. [FILTRATE] * Saturation vapour pressure was reduced from 6.1368d0 to 6.112d0 [INIT] * Variable almb was renamed albm [GETMET] * A safety condition ensures that the effective saturation does not exceed total saturation. The loop index for the effective saturation was changed from i to j (ss(i) to ss(j)) [FILTRATE] 210

Appendix A. Code Inconsistencies * Initial grain size, if not included in the input files, is estimated by function FD. This function is now called after calculation of the bulk properties in subroutine DENSITY. [SNTHERM] * Grain diameter and snow depth are returned to previous values if the solution has not converged or if a node has melted too quickly, and the time step calculation is to be repeated. [RESET] * Where subroutine DENSITY is called to reset all nodal bulk properties once the time step calculation is to be repeated, the loop index for porosity within the subroutine variable list is changed from m to i. [RESET] * Extra safety conditions have been input to set negative measured solar radiation to zero, prior to calculation of snow albedo [GETMET] 211

Bibliography Adams, R. S., Spittlehouse, D. L., Winkler, R. D., 1996. The effect of a canopy on the snowmelt energy balance. In: 64th Annual Meeting Western Snow Conference. Bend, Oregon, pp. 171-174. Anderson, E. A.., 1976. A point mass and energy balance model of a snow cover. Tech. Rep. CRREL Report 83-12, USA Cold Regions Research and Engineering Laboratory. Andreas, E. L., Murphy, B., 1986. Bulk transfer coefficients for heat and momentum over leads and polynyas. Journal of Physical Oceanography 16 (11), 1875-1883. Betts, A. K., Ball, J. H., 1997. Albedo over the boreal forest. Journal of Geophysical Research 104 (D24), 28901-28909. Betts, R. A., 2000. Offset of the potential carbon sink from boreal forestation by decreases in surface albedo. Nature 408, 187-190. Boynton, W. P., Brattain, W. H., 1929. Int. Crit. Tables 5, 62-65. Brandt, R. E., Warren, S. G., 1993. Solar-heating rates and temperature profiles in Antarctic snow and ice. Journal of Glaciology 39 (131), 99-110. 212

Bibliography Bibliography Brooks, R. H., Corey, A. T., 1964. Hydraulic properties of porous media. Hydrology Papers Colorado State University, Fort Collins, CO. Brun, E., Martin, E., Simon, V., Gendre, C., Coleau, C., 1989. An energy and mass model of snow cover suitable for operational avalanche forecasting. Journal of Glaciology 35, 333-342. Buck, A., 1981. New equations for computing vapor pressure and the enhancement factor. Journal of Applied Meteorology 20 (12), 1527-1532. Chen, J., Rich, P. W., Gower, S. T., Norman, J. M., Plummer, S., 1997. Leaf area index of boreal forests: theory, techniques and measurements. Journal of Geophysical Research 102 (D24), 29,429-29,444. Colbeck, S., Anderson, E. A., 1982. The permeability of a melting snow cover. Water Resources Research 18 (4), 904-908. Colbeck, S. C., 1974. The capillary effects on water percolation in homogeneous snow. Journal of Glaciology 13 (67), 85-97. Decagon Devices, 1987. Sunfleck ceptometer user manual. Tech. rep., Decagon Devices, Pullman, Washington, USA. Dorsey, N. E., 1940. Properties of ordinary water-substance. Reinhold Publishing Corporation, Ch. 13, pp. 72-73. Field, R. T., Fritschen, L. J., Kanemasu, E. T., Smith. E. A., Stewart, J. B., Verma, S. B., Kustas, W. P., 1992. Calibration, comparison, and correction of net radiation instruments used during FIFE. Journal of Geophysical Research 97 (D17), 18,681-18,695. 213

Bibliography Bibliography Giddings, J. C., LaChapelle, E. R., 1961. Diffusion theory applied to radiant energy distribution and albedo of snow. Journal of Geophysical Research 66, 181-189. Glendinning, J. H. G., 1997. The modelling of radiative transfer in snow at visible and near infrared wavelengths. Ph.D. thesis, The University of Reading. Golding, D. L., Swanson, R. H., 1986. Snow distribution patterns in clearings and adjacent forest. Water Resources Research 22 (13), 1931-1940. Hanson, C. L., Marks, D., Vactor, S. S. V., 2000. Climate monitoring at the Reynolds Creek Experimental Watershed, Idaho, USA. Tech. Rep. ARS Techincal Bulletin NWRC-2000 -6, USDA - Agricultural Research Service, Northwest Watershed Research Center, 800 Park Blvd, Suite 105, Boise, ID 83712-7716. Harding, R. J., Pomeroy, J. W., 1996. Energy balance of the winter boreal landscape. Journal of Climate 9 (11), 2778-2787. Hardy, J. P., Davis, R. E., Jordan, R., Li, X., Woodcock, C., Ni, W., McKenzie, J. C., 1997a. Snow ablation modeling at the stand scale in a boreal jack pine forest. Journal of Geophysical Research 104 (D24), 29397-29405. Hardy, J. P., Davis, R. E., Jordan, R., Ni, W., Woodcock, C., 1997b. Snow ablation modeling in conifer and deciduous stands of the boreal forest. In: 65th Western Snow Conference. Banff, Canada, pp. 114-124. Hardy, J. P., Hansen-Bristow, K. J., 1990. Temporal accumulation and ablation patterns of the seasonal snowpack in forests representing varying stages of growth. In: 58th Western Snow Conference. Sacramento, California, pp. 23-34. 214

Bibliography Bibliography Hardy, J. P., Melloh, R., Robinson, P., Jordan, R., 2000. Incorporating effects of forest litter in a snow process model. Hydrological Processes 14, 3227-3237. Hedstrom, N. R., Pomeroy, J. W., 1997. Accumulation of intercepted snow in the boreal forest: measurements and modelling. In: 56th Western Snow Conference. Banff, Canada, pp. 130-141. Hedstrom, N. R., Pomeroy, J. W., 1998. Accumulation of intercepted snow in the boreal forest: measurements and modelling. Hydrological Processes 12, 1611-1625. Hellstrom, R. A., 1999. Representation of forest in a physically based snowmelt model, phase I. In: 56th Annual Meeting Eastern Snow Conference. Fredericton, N.B., Canada, pp. 215-231. Henderson-Sellers, A., Wilson, M. F., 1983. Surface albedo data for climatic modeling. Reviews of Geophysics and Space Physics 21 (8), 1743-1778. Jacobson, M. Z.. 1999. Fundamentals of Atmospheric Modeling. Cambridge University Press, Ch. 10, pp. 280-282. Ji, Q., Tsay, S., 2000. On the dome effect of Eppley pyrgeometers and pyranometers. Geophysical Research Letters 27 (7), 971-974. Johnson, J. B., Shaefer, G. L., 2001. The relationship between real-time snow water equivalent pressure sensor measurements and snow/soil mechanical, thermal and hydrologic processes. In: EOS Transactions AGU,82(47), Fall Meeting Supplement, Abstract IP51A-0727. San Francisco, CA, USA. 215

Bibliography Bibliography Jordan, R., 1991.. A one-dimensional temperature model for a snow cover: Technical documentation for SNTHERM. Tech. Rep. CRREL Special Report 91-16, USA Cold Regions Research and Engineering Laboratory. Jordan, R., 1992. Estimating turbulent transfer functions for use in energy balance models. Tech. Rep. Internal Report 1107, USA Cold Regions Research and Engineering Laboratory. Jordan, R., 1995. Effects of capillary discontinuities on water flow and water retention in layered snow covers. Defence Science Journal 45 (2), 79-91. Jordan, R., 1996. User's guide for USA CRREL ONE-DIMENSIONAL SNOW TEMPERATURE MODEL (SNTHERM.89). Tech. rep., USA Cold Regions Research and Engineering Laboratory. Jordan, R. E., Andreas, E. L., Makshtas, A. P., 1999. Heat budget of snow-covered sea ice at North Pole 4. Journal of Geophysical Research 104 (C4), 7785-7806. Kattelmann, R., 1986. Measurements of snow water retention. In: Proceedings of the Cold Regions Hydrology Symposium, American Water Resources Association. pp. 377-386. Koh, G., Jordan, R., 1995. Sub-surface melting in a seasonal snow cover. Journal of Glaciology 41 (139), 474-482. Kojima, K., 1955. Viscous compression of natural snow-layer, I. Low Temperature Science Ser. A (14). Kondo, J., Yamazaki, T., 1990. A prediction model for snowmelt, snow surface tempera 216

Bibliography Bibliography ture and freezing depth using a heat balance method. Journal of Applied Meteorology 29, 375-384. Leonard, R. E., Eschner, A. R., 1968. Albedo of intercepted snow. Water Resources Research 4 (5), 931-935. Li, X., Strahler, A. H., 1986. Geometrical-optical bidirectional reflectance modeling of a conifer forest canopy. IEEE Transactions on Geoscience and Remote Sensing GE-24, 906-919. Li, X., Strahler, A. H., Woodcock, C. R., 1995. A hybrid geometric optical-radiative transfer approach for modeling albedo and directional reflectance of discontinuous canopies. IEEE Transactions on Geoscience and Remote Sensing 33 (2), 466-480. Link, T. E., Marks, D., 1999. Point simulation of seasonal snow cover dynamics beneath boreal forest canopies. Journal of Geophysical Research 104 (D22), 27,841-27,857. Marks, D., 1988. Climate, energy exchange, and snowmelt in Emerald Lake Watershed, Sierra Nevada. Ph.D. thesis, University of California, Santa Barbara. Marks, D., Winstral, A., 2001. Comparison of snow deposition, the snowcover energy balance, and snowmelt at two sites in a semi-arid mountain basin. Journal of Hydrometeorology 2 (3), 213-227. Marshall, S., Roads, J. O., Glatzmaier, G., 1994. Snow hydrology in a general circulation model. Journal of Climate 7, 1251-1269. Marshall, S. E., Warren, S. G., 1987. Parameterization of snow albedo for climate models. 217

Bibliography Bibliography In: Goodison, B. E., Barry, R. G., Dozier, J. (Eds.), Large scale effects of seasonal snow cover. Vancouver, Canada, pp. 43-50. Melloh, R. A., 1999. A synopsis and comparison of selected snowmelt algorithms. Tech. Rep. CRREL Report 99-8, USA Cold Regions Research and Engineering Laboratory. Monteith, J. L., Unsworth, M., 1990. Principles of Environmental Physics, 2nd Edition. Arnold, Ch. 6, p. 81. Morris, E. M., 1989. Turbulent transfer over snow and ice. Journal of Hydrology 105, 205-223. Pearson, D., I)aamen, C. C., Gurney, R. J., Simmonds, L. P., 1999. Combined modelling of shortwave and thermal radiation for one-dimensional SVATs. Hydrology and Earth System Sciences 3 (1), 15-30. Pierce, L. L., Running, S. W., 1988. Rapid estimation of coniferous forest leaf area index using a portable integrating radiometer. Ecology 69 (6), 1762-1767. Pomeroy, J. NV., Dion, K., 1996. Winter radiation extinction and reflection in a boreal forest pine canopy. In: 53th Annual Meeting Eastern Snow Conference. Williamsburg, Virginia, pp. 105-118. Pomeroy, J. W., Goodison, B. E., 1997. The Surface Climates of Canada. McGill-Queen's University Press, Ch. 4: Winter and Snow, pp. 68-100. Pomeroy, J. W., Gray, D. M., 1995. Snowcover Accumulation, Relocation and Management. Tech. Rep. NHRI Science Report No. 7, NHRI, Environment Canada, Saskatoon. 218

Bibliography Bibliography Pomeroy, J. WV., Parviainen, J., Hedstrom, N., Gray., D. M., 1998. Coupled modelling of forest snow interception and sublimation. Hydrological Processes 12, 2317-2337. Pomeroy, J. W., Schmidt, R. A., 1993. The use of fractal geometry in modelling intercepted snow accumulation and sublimation. In: 50th Eastern Snow Conference. Quebec City, Canada, pp. 1-10. Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling, W. T., 1990. Numerical recipies in C. Cambridge University Press, Ch. 5, p. 157. Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., 1992. Numerical recipies in Fortran. Cambridge University Press, Ch. 2, pp. 42-43. Rouse, W. R., 1984. Microclimate at arctic tree line 1. Radiation balance of tundra and forest. Water Resources Research 20 (1), 57-66. Schmidt, R. A., 1991. Sublimation of snow intercepted by an artificial conifer. Agricultural and Forest Meteorology 54, 1-27. Sellers, P. J., Hall, F. G., Kelly, R. D., Black, A., Baldocchi, D., Berry, J., Ryan, M., Ranson, K. J., Crill, P. M., Lettenmaier, D. P., Margolis, H., Cihlar, J., Newcomer, J., Fitzjarrald, D., Jarvis, P. G., Gower, S. T., Halliwell, D., Williams, D., Goodison, B., Wickland, D. E., Guertin, F. E., 1997. BOREAS in 1997: Experiment overview, scientific results, and future directions. Journal of Geophysical Research 104 (D24), 28731-28769. Shimizu, H., 1970. Air permeability of deposited snow. Tech. Rep. Contribution No.1053, Institute of Low Temperature Science, Sapporo, Japan, english Translation. 219

Bibliography Bibliography Shultz, E. F., Albert, M. R., 1998. An automated procedure for plotting snow stratigraphy. In: 55th Annual Meeting Eastern Snow Conference. Jackson, New Hampshire, pp. 147 -151. Sturm, M., 1992. Snow distribution and heat flow in the Taiga. Arctic and Alpine Research 24 (2), 145 —152. Tarboton, D. G., Luce, C. H., 1996. Utah energy balance snow accumulation and melt model (UEB). computer model technical description and users guide. Tech. rep., Utah Water Research Laboratory, Utah State University and USDA Forest Service Intermountain Research Station. Thomas, G., Rowntree, P. R., 1992. The boreal forests and climate. Quarterly Journal of the Royal Meteorological Society 118, 469-497. Tuteja, N. K., Cunnane, C., 1997. Modelling coupled transport of mass and energy into the snowpack - model development, validation and sensitivity analysis. Journal of Hydrology 195, 232-255. Viterbo, P., Betts, A. K., 1999. The forecast impact of changes to the albedo of the boreal forests in the presence of snow. Journal of Geophysical Research 104, 27,803-27,810. Warren, S. G., Wiscombe, W. J., 1980. A model for the spectral albedo of snow. II. Snow containing atmospheric aerosols. Journal of Atmospheric Science 37, 2734-2745. Wigmosta, M. S., Vail, L. W., Lettenmaier, D. P., 1994. A distributed hydrologyvegetation model for complex terrain. Water Resources Research 30 (6), 1665-1679. 220

Wiscombe, W. J., Warren, S. G., 1980. A model for the spectral albedo of snow. I. pure snow. Journal of the Atmospheric Sciences 37 (12), 2712-2733. Yamazaki, T., Kondo, J., 1992. The snowmelt and heat balance in snow-covered forested areas. Journal of Applied Meteorology 31 (11), 1322-1327. 221