Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Impact of an artifical circulation device on the heat budget of an ice-covered mid-latitude lake Rogers, Christopher K. 1993

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

Item Metadata


831-ubc_1993_spring_rogers_christopher.pdf [ 6.97MB ]
JSON: 831-1.0050455.json
JSON-LD: 831-1.0050455-ld.json
RDF/XML (Pretty): 831-1.0050455-rdf.xml
RDF/JSON: 831-1.0050455-rdf.json
Turtle: 831-1.0050455-turtle.txt
N-Triples: 831-1.0050455-rdf-ntriples.txt
Original Record: 831-1.0050455-source.json
Full Text

Full Text

IMPACT OF AN ARTIFICIAL CIRCULATION DEVICEON THEHEAT BUDGET OF AN ICE-COVERED MID-LATITUDE LAKEbyChristopher Kavanagh RogersB. A. Sc. (hons.), University of WaterlooA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESDepartment of Civil EngineeringWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1992© Christopher Kavanagh Rogers, 1992In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)Department of OWL_^--ArarAfEereptiaThe University of British ColumbiaVancouver, CanadaDate DE-6 (2/88)ABSTRACTTwo lakes located in the Southern Interior Plateau of British Columbia were selected fora field investigation in order to assess the impact of artificial circulation on the heatbudget of a small mid-latitude ice-covered lake. The heat transfer algorithm developedby Patterson & Hamblin (1988) for high latitude lakes was extended to include theimpacts of snowmelt due to rain, sediment heat transfer, snow-ice formation, and day today variations in snow density, snow conductivity and albedo. Since the lakesconsidered are nearly isothermal in winter, the new model ignores the internalhydrodynamics of the problem. This model was tested on Harmon Lake, which wasselected as a control for the study. The model was further modified to include the effectsof artificial circulation at nearby Menzies Lake. These effects include polynyadevelopment, and a substantial reduction in average temperature. Heat losses due to freeconvective evaporation and direct snowfall on open water were added to the set ofstandard aerodynamic formulae used to determine the net meteorological heat flux acrossa water surface. Turbulent heat transfer from the circulated water to the ice cover wasestimated based on an empirical surface velocity relationship derived from fieldmeasurements. The size of the polynya is estimated by means of a simple heat balancewhich also involves the surface velocity function.The Harmon Lake predictions agree well with the field data. All discrepancies could beaccounted for by parameter uncertainties and expected observation error. It was foundthat sediment heat transfer may be important in early winter in preventing a net loss ofheat from the lake water. Significant heat gains in the latter part of winter, however, areattributed to the penetration of solar radiation. Once calibrated, the Menzies Lakepredictions are also good. It was found that, over the period in which lake temperaturedropped, the average heat loss due to turbulent heat transfer between water and ice wasiithree times that across the polynya surface. The former heat flux continued to increase asthe lake warmed up again, while the latter fell, on average, over a short period untilincreased solar heating resulted in a reversal in the direction of heat transfer across thepolynya. Discrepancies between early winter ice thickness predictions and observations,could not be accounted for. It is suspected, however, that these discrepancies are a resultof the impact which the heat flux across the polynya may have on the heat flux throughthe ice and snow cover. In this thesis, it is assumed that these two fluxes are independentof one another.iiiTABLE OF CONTENTSABSTRACT^ iiLIST OF FIGURES^ viiLIST OF TABLES^ xLIST OF SYMBOLS xiACKNOWLEDGMENT^ xivINTRODUCTION^ 11.1 Winterkill  11.1.1 The Oxygen Depletion Problem^ 11.1.2 The Artificial Circulation Solution  11.2 Project Scope^ 41.2.1 Importance of the Heat Budget^ 41.2.2 Winter Stratification^ 61.2.3 The Natural State 6BACKGROUND PHYSICS & RELEVANT LITERATURE REVIEW^ 82.1 Existing Heat Budget Models^ 82.1.1 A Brief History of Lake Simulation 82.1.2 The Development of DYRESMI^ 92.1.3 DYRESMI Formulation  112.2 Heat transfer from Bottom Sediments 172.2.1 Review of Observed Rates of Heat Transfer^ 172.2.2 Methods of Estimating Sediment Heat Transfer 182.2.3 Sediment properties^  192.3 Effects of Snow Cover 202.3.1 Snow Properties 202.3.2 Estimation of Snow Albedo^ 222.3.3 Snow-ice Formation 262.3.4 Snowmelt Due to Rain Events 262.4 Effect of Artificial Circulation^ 282.4.1 Review of Systems Used Worldwide^ 282.4.2 Heat Transfer Between Water and Air 31iv2.4.3 Heat Transfer Between Water and Ice^ 332.4.4 Flow Field^ 362.5 Previous Work and the MLI Model^ 38FIELD INVESTIGATION^ 403.1 Characteristics of Study Lakes^ 403.1.1 Morphometry 403.1.2 Annual Temperature and Oxygen Budgets^ 443.2 Data Collection Procedures and Observations 493.2.1 Data Requirements^ 493.2.2 Instrumentation 513.2.3 Field Procedure 553.2.4 Field Trip Observations 58NUMERICAL MODEL^ 724.1 MLI Formulation 734.2 Bulk Aerodynamic Formulae^ 754.3 MLI: Lake in the Natural State 794.3.1 MLI Main Routine 794.3.2 MLI Subroutines^ 824.3.3 Sediment Heat Transfer Estimates^ 944.4 MLI-C: Artificially Circulated Lake 974.4.1 MLI-C Subroutines^ 994.5 Data Preparation^ 1064.6 Summary of Parameters 107RESULTS AND DISCUSSION^ 1095.1 Lake in Natural State 1095.1.1 Observation Errors^ 1115.1.2 Analysis of MLI Predictions^  1125.1.3 Components of the Heat Budget 1175.1.4 Impact of MLI Parameters on Predictions^ 1235.2 Artificially Circulated Lake^ 1365.2.1 Analysis of MLI-C Predictions^ 1365.2.2 Components of the Heat Budget 1425.2.3 Impact of Varying MLI-C Calibration Parameters^ 1485.2.4 Improving Ice Thickness Predictions 154CONCLUSIONS AND RECOMMENDATIONS^ 1566.1 General Remarks^ 1566.2 MLI^ 1576.2.1 Summary of Results^ 1576.2.2 Recommendations for Future Research^ 1606.3 MLI-C^ 1626.3.1 Summary of MLI-C Results^ 1626.3.2 Recommendations for Artificial Circulator Design^ 1646.3.3 Recommendations for Future Research^ 165REFERENCES^ 167APPENDIX ATHEORETICAL SOLAR RADIATION UNDER CLEAR SKIES^ 172viLIST OF FIGURES1.1 The Air-o-lator®^ 21.1a Polynya Generated by Air-o-lator® at Menzies Lake^ 31.2 Oxygen and Temperature Isotherms: Menzies Lake, 1991 52.1 Heat Fluxes Across an Ice and Snow Cover^ 133.1 Field Study Location^ 413.2 Bathymetry of Menzies and Harmon Lakes^ 423.2b Menzies and Harmon Lakes Hypsographs 433.3 Isotherms at Menzies and Harmon Lakes, June, 1991 - April, 1992^ 453.4 Oxygen Isopleths at Menzies and Harmon Lakes, June, 1991 - April, 1992 463.5 Total Dissolved Solids Near Lake Surface and Sediments^ 483.6 Menzies Lake Weather Station^ 523.7 Isotherms at Menzies and Harmon Lakes, December 13th - March 10th^ 603.8 Air Temperature and Solar Radiation^ 613.9 Ice and Snow Thicknesses^ 633.10 Isotherms Generated from Datalogger Data at Menzies Lake^ 663.11 Maximum and Minimum Likely Heat Content at Menzies Lake^ 673.12 Radial Jet Velocity Profile^ 714.1 MLI Flow Chart^ 804.2 SNOWDENS Flow Chart^ 834.3 SNOWICE Flow Chart 864.4 SOLAR Flow Chart^ 884.5 ALBEDO Flow Chart 894.6 MELT Flow Chart^ 914.7 ICEWATER Flow Chart 92vii4.8^Heat Fluxes Across Ice-Covered and Ice-Free Zones of Lake^ 984.9^NEWAREA Flow Chart^ 1025.1^Ice and Snow Thickness and Whole Lake Temperature at Harmon Lake:A Comparison of Results Using MLI and Patterson & Hamblin's (1988)model^ 1105.2^MU: Surface Heat Budget Components and Heat Fluxes to and from LakeWater 1165.3^a) MLI & DYRESMI: Comparison of Heat Budget Componentsb) MLI: Heat Fluxes Through the Cover; Air and Surface Temperatures ^ 1215.4^Effect of Heat Associated with the Formation of Snow-Ice on Ice andSnow and Whole Lake Temperature^ 1245.5^Effect of Increased Solar Attenuation Through Snow-Ice on Ice and SnowThickness and Whole Lake Temperature 1265.6^Effect of Increased Maximum Snow Density on Ice and Snow Thicknessand Whole Lake Temperature^  1285.7^Effect of Increased Snow Conductivity on Ice and Snow Thickness andWhole Lake Temperature 1305.8^Effect of Decreased Sediment Conductivity on Ice and Snow Thicknessand Whole Lake Temperature^ 1315.9^Effect of Reduced Albedo Decay Rate on Ice and Snow Thickness andWhole Lake Temperature  1335.10^Effect of Reduced Thermal Gradient at Ice-Water Interface on Ice andSnow Thickness and Whole Lake Temperature^ 1355.11^Snow and Ice Thickness and Whole Lake Temperature Predictions:A Comparison of MLI and MLI-C Output for Menzies Lake^ 1375.12^Predicted and Observed Polynya Radius^ 1435.13 MLI-C: Surface Heat Budget Components and Heat Fluxes to and fromLake Water at Menzies Lake^ 1445.14^MLI-C: Surface Heat Budget Components Across Polynya at MenziesLake^ 1475.15^Effect of Increased Radius of Turbulent Heat Transfer on Polynya Radiusand Whole Lake Temperature^ 1495.16^Effect of Increased Minimum Possible Ice Thickness on Polynya Radiusand Whole Lake Temperature  151viii5.17^Effect of Decreased Turbulent Heat Transfer Coefficient on PolynyaRadius and Whole Lake Temperature^ 153ixLIST OF TABLES2.1 Albedo Under Various Surface Conditions ^ 232.2 Polynya Areas in Finnish Harbours 293.1 Morphometry of Menzies and Harmon Lakes^ 443.2 Field Trip Schedule^ 563.3 Observed Average Daily Cloud Cover^ 684.1 Snow-Ice and Pure Ice Albedo 874.2 Summary of Parameters^ 108xLIST OF SYMBOLSThe following is a list of the most common symbols in this thesis.A^ spectral fraction^ areaC^ fraction of sky which is cloud-coveredCpi^heat capacity of iceCpw^heat capacity of waterCt^ turbulent heat transfer coefficientd^number of days since snowfallH^net meteorological fluxh^ thickness of mediumhen^ latest snow-ice accumulationhsm^maximum possible snow thickness10^non-reflected solar radiation at surfacelw^ solar radiation reaching waterK^conductivityL^ latent heat of fusionP^ rainfallPatm^ atmospheric pressureof^ heat flux in ice at ice-water interfacego^ heat flux at the surfacegsed^ sediment heat fluxqt^ turbulent heat transfer from water to iceqw^ heat conduction from water to iceQ^volumetric flow rateQc^ sensible heat fluxQe^ evaporative heat fluxQfc^ heat flux due to free convectionQ0^ incoming solar radiationQr^ heat due to rainfallQsi^heat due to snow-ice formationQsp^heat flux due to snowfall on waterQt^ total turbulent heat transfer to icerp^polynya radiusroc,^ rp + firRli^ incoming long-wave radiationR10^outgoing long-wave radiationsvp0^ saturated vapour pressure at surface temperaturesvpd^vapour pressure of airS^ Stefan NumberSW^observed total daily short-wave radiationSw^wind sheltering coefficientSWCS ^ total daily short-wave radiation on clear dayT^ temperatureU radial jet velocityUw^wind velocity^ lake volumez^vertical coordinate (depth)zmax^maximum depthi^mean depthSubscripts1 ^visible spectrum of radiation2^ infrared spectrum of radiatione snow-icef^ freezingi ^ icem^meltingo surfacep^polynyape^polynya edges^ snowsed^ sedimentw^waterGreek Symbolsa^ albedo0^ solar angle at solar noonAr^ radial distance over which Qt is active£^ emissivityX^ solar attenuation coefficientp^density131 ^ snow density 24 hrs after snowfallPm^maximum snow densityPn^ density of new snowa^ Stefan-Boltzmann ConstantACKNOWLEDGMENTA great deal a support was provided for this project in terms of financing, technicalconsultation and field assistance. Funding was provided by the Natural Sciences andEngineering Research Council (N.S.E.R.C), the Science Council of British Columbia,and the University of British Columbia. Equipment purchases and field expenses werepaid for by N.S.E.R.C. and the B.C. Ministry of the Environment. The author would liketo thank Dr. P.F. Hamblin for his patient explanations regarding the complexities of lakemodeling, and Dr. G.A. Lawrence for lending a sympathetic ear in times of frustration.The suggestions and comments provided by these two individuals were crucial to thesuccessful completion of this thesis. Dr. Craig Stevens and Dr. Noboru Yonomitsu arealso recognized for their help in this regard.The field study component of this project could not have been carried out without theguidance of Mr. K.I. Ashley who has been an excellent tutor in the area of fieldlimnology. The students and staff of the U.B.C. Department of Civil Engineering (aswell as many friends) who volunteered their time as field assistants are too numerous tolist here, but their help has been greatly appreciated. In particular, the author wishes tothank Mr. B. Laval for the many long days he has dedicated to the acquisition of fielddata.xivChapter 1INTRODUCTION1.1 WINTERKILL1.1.1 The Oxygen Depletion ProblemMany of the small kettle lakes of the B.C. Interior are naturally eutrophic, and thereforeexperience high oxygen depletion rates due to the decay of organic material. Over thewinter season, replenishment of oxygen is prevented by ice and snow cover. The ice andsnow constitute a physical barrier to the diffusion of oxygen into the water and can alsoinhibit photosynthetic production of oxygen due to light limitation. As a result, thedissolved oxygen content in these lakes may drop below levels which can support fishlife. Deleterious effects on aquatic life may occur below about 1 0C (pers. comm.,Ashley, 1991). The B.C. Ministry of the Environment (Fisheries Branch) has addressedthis problem, commonly termed winterkill, by implementing an extensive winter aerationprogramme throughout the Southern B.C. Interior.1.1.2 The Artificial Circulation SolutionThe aeration scheme simply involves artificially circulating the water in order to maintaina polynya (ice-free patch) over part or all of the winter season. This is often done using acommercial unit called the Air-o-lator®, which is shown in Figure 1.1. The polynyagenerated by this device at Menzies Lake near Merritt, B.C. is shown if Figure 1.1a. Thesize of the aeration equipment must be sufficient to ensure a large enough ice-free zone toallow adequate diffusion and turbulent entrainment of oxygen into the lake, but small1N2"2/31. Diffuser2. Float3. Propeller4. Draught Tube5. Motor6. Inflow7. Radial jetenough to avoid cooling the water to temperatures which threaten the resident aquaticlife. The hardware generally consists of a compressor or axial flow pump fixed inside ofa vertical draught tube, which is attached to a floatation device. The assembly isanchored as close to the deepest part of the lake as practically possible. A single 0.75 kWpump will produce an ice-free zone of 30-50 m in diameter, which is sufficient to preventwinterkill in lakes up to approximately 10 ha in surface area (Ashley, pers. comm.,1991).Water is drawn up through the draught tube, which extends 0 to 3 metres into the lake,and is projected into the air. The fountain-like jet plunges back into the water andspreads radially. The area around the device remains ice-free out to a critical distancewhere the artificially-induced turbulence and thermal energy is insufficient to preventfreezing.2Figure 1.1 The Air-o-lator®Figure 1.1a Polynya Generated by Air-o-lator® at Menzies Lake(photograph taken from north end of lake)31.1.3 Circulator EffectivenessField measurements at Menzies Lake near Merritt, B.C. in mid-February, 1991 indicatedvery effective reaeration due to artificial circulation. The weak inverse thermalstratification normally observed under the ice was mixed down to a depth of about 12metres. The depth of the pycnocline is controlled by a 3 m draught tube and the lakemorphometry as a result of its shallowness. As shown in Figure 1.2, no significanthorizontal variations of oxygen or temperature were observed other than those caused byheat and mass transfer across the sediment-water interface (see Chapter 3). Thisindicates that a one-dimensional representation of the effect of the circulator isappropriate provided a sufficient rate of pumping is maintained. The top 12 metrescontained 7 to 8 mg/1 of dissolved oxygen and was isothermal at about 1°C. Below thethermocline, the water was anoxic (below 0.5 mg/1) and near the temperature ofmaximum density. In addition to aeration, the circulator has significant cooling effect onthe lake. When ice-covered, a mid-latitude lake in its natural state will generally be closeto 4°C over the entire water column below a relatively thin boundary layer which bringsthe water to 0°C at the ice-water interface.1.2 PROJECT SCOPE1.2.1 Importance of the Heat BudgetOxygen concentrations depend on water temperature and ice cover characteristics inseveral ways. These dependencies include rates of decay of organic material, oxygensaturation concentrations, inhibition of photosynthesis through light attenuation, andduration of ice-cover. The first step in developing a comprehensive oxygen budget41.8Oxygen IsoplethsMenzies Lake, February, 19910.0^►^I^I^I^I^I8.0circulator location127^184 241^298^3553.67.010.812.6^/40.2/ lake bottom►^►^I14.4/ 116.2-44^7013355,lake bottom13^70^127^184^241^298IsothermsMenzies Lake, February, 1991I^I^I^I^I-44Distance from Circulator Location (m)5Figure 1.2. Oxygen and Temperature Isotherms: Menzies Lake, February, 1991model for the artificially circulated lake is therefore a model which simulates the heatbudget. The current investigation is limited to this first step.1.2.2 Winter StratificationAn exhaustive treatment of the heat budget would require a numerical description of theinternal hydrodynamics induced by the artificial circulation. It is suspected that selectivewithdrawal theory (see §3.2.4) may govern the fate of the winter pycnocline, but thecontinual operation, induced lake circulation, and radial jet behaviour make the problemunique. In addition, with the formation of a circulator-induced winter thermocline, thehypolimnetic water quickly becomes anoxic (this has been observed throughout thecurrent study). The anoxia leads to the release of metals which are sensitive to oxidation-reduction potential (see Mortimer, 1941; 1942; Stauffer, 1987). The build-up of solidsnear the lake bottom may actually be greater than if the lake were not artificiallycirculated because diffusion is restricted by the density gradient associated with thethermocline above. The resulting increase in solids concentration serves to stabilize thecirculator-induced stratification.The circulation equipment, however, is generally sized in B.C. winterkill applications tocause the lake to overturn continually over the winter period, without disturbing the lakesediments. Furthermore, the winter hypolimnion at Menzies Lake is of almost negligiblevolume, compared with the mixed region above (this fact is illustrated by themorphometric data given in Chapter 3). In order to simplify the heat budget problem,complete mixing will be assumed, and the internal hydrodynamics ignored.1.2.3 The Natural StateBefore examining the effects of artificial circulation, a model must first be developed topredict the thermal behaviour of the ice-covered lake in its natural state. In spite of the6existence of many hydrodynamic lake models, few have tackled the ice-cover problem.Those that have are highly simplified and/or do not deal with the particularitiesassociated with mid-latitude lakes. For this study, an existing ice and snow-cover modelwas selected and improved to account for these particularities. The product is called MLI(Mixed Lake with Ice cover). It is applied here to Harmon Lake which is situated in thesame geographical region as Menzies Lake (see §3.1). The improved version was thenmodified to account for artificial circulation. This modified version is called MLI-C, andis applied to Menzies Lake.7Chapter 2BACKGROUND PHYSICS & RELEVANTLITERATURE REVIEW2.1 EXISTING HEAT BUDGET MODELS2.1.1 A Brief History of Lake SimulationA number of water quality simulation models have been developed to address the need tomanage the quality of water stored in reservoirs. Most are fashioned after the WRE(Water Resources Engineers) model originally developed by Chen and Orlob (1975).Thermal stratification is of primary importance to water quality and therefore much efforthas been focused on its prediction. As temperature variations in small to medium sizedlakes and reservoirs is for most purposes, only significant in the vertical direction, theone-dimensional horizontal slab concept first advanced by Raphael (1962) is basic tomost of these models. The heat budget of most lakes with high retention times arebasically a function of meteorological forcing, and for simplicity, bulk aerodynamicformulae, such as those given by the Tennessee Valley Authority (TVA, 1972) aregenerally applied.The incorporation of ice-cover routines into existing models, however, is only in theearly stages of development. Recently, two independent research teams have make thisattempt using DYRESM (see Imberger & Patterson, 1981), a commercially availablestratification model which has been thoroughly tested and proven for temperate climates.Both models, however, were developed for lakes at high latitude, and have beeninsufficiently tested (Patterson & Hamblin, 1988; Gosink, 1987). They are also similar in8that they involve a solution to the heat conduction equation which includes a depthdependent heat source to account for the penetration of solar radiation through the ice andsnow. The Patterson & Hamblin (1988) model is superior, however, to the Gosink(1987) model in at least two ways. In the latter model, the heat flux through the entire iceand snow cover is assumed constant and equal to the balance of meteorological forcing.Secondly, the short wave radiation is assumed to consist of only one spectral bandcharacterized by a single attenuation coefficient for each medium through which itpasses. In the former model, however, the heat flux through the cover is a function ofdepth, and two spectral bands are used. Consequently, the Patterson & Hamblin (1988)model (DYRESMI) was selected and modified for application to this investigation. Theformulation will be described in detail following a brief overview of the development oftheir model.2.1.2 The Development of DYRESMIThe definitive thermodynamic model of surface ice formation and ablation, according thePatterson & Hamblin (1988), is that of Maykut and Untersteiner (1971). It is a 1-Dunsteady heat transfer model through a two component ice and snow cover, developedfor sea ice, but not linked to oceanic mixing. As such, the heat flux across the ice-waterinterface is left unknown and must be specified by the user in order to arrive at a solution.As an alternative to the same model, Semter (1976) showed that a quasi-steady statesolution to the heat transfer equations is reasonable if the ice is thin. The approachimplies that the temperature distribution in the ice-snow cover has reached a steady statebefore a time step has elapsed. It has been shown that a quasi-steady state solution willbe within 5% of a transient solution if the Stefan Number, S, is less than 0.1 (see Hill &Kucera, 1983):S— CPi• ATL (2.1)9where Cpi = specific heat of iceAT = characteristic temperature difference across ice thicknessL = Latent heat of fusionFor the above condition to hold, AT must be less than about 15 0C. This restriction ismore likely to be satisfied at mid-latitude lakes than at those at high latitudes. The Stefancondition may be violated, however, provided that the ice (or snow) layer is thin. If it isassumed that the temperature gradient through the is ice cover is small compared withthat through the snow cover, then the time scale for heat conduction through the snowmust be short compared to the time scale of changes in meteorological forcing. Thislimiting time scale is of order h s2hcs , where h s = snow thickness, and Ks = thermaldiffusivity of snow (see Patterson & Hamblin, 1988). Given a time step of one day (as isthe case for DYRESM and DYRESMI), the snow thickness must be less than about 20cm for the steady-state solution to be reasonable. Since the maximum snow thicknessover lake ice is limited due to the low buoyancy of ice (see §2.3.3), it is relativelyuncommon to observe thicknesses in excess of 20 cm, especially during mild winters.For a snow-free ice surface, the equivalent limiting ice thickness is about 30 cm.Patterson and Hamblin (1988) developed two main improvements to the Maykut &Untersteiner (1971) model. First, a thermodynamic link with DYRESM was established,negating the need to specify the ice-water heat flux. Secondly, the effect of partial icecover was incorporated into the model.The thermodynamic link with DYRESM is interactive; DYRESM determines the heatflux at the ice-water interface, thereby affecting the ice cover, which in turn provides aboundary condition for the mixing model. Partial ice cover is dealt with by assuming thatice will not persist below a minimum thickness, but will be transported by wind stressand accumulate elsewhere at this specified thickness. This thickness has been observedto be about 10 cm for medium-sized lakes (Hamblin, pers. comm., 1992). During periods1011of partial cover the heat transfer will be non-uniform over the lake surface. Horizontaltemperature gradients are instantly relaxed, however, so that the model remains one-dimensional.2.1.3 DYRESMI FormulationProvided that freezing conditions persist over the winter period, both a layer of ice and alayer of snow will usually be observed on a lake. The distribution of temperature in theice and snow is governed by the heat conduction equation which includes the absorptionof solar radiation as a distributed heat source. The conductivity and attenuationcharacteristics of ice and snow are significantly different, and this must be reflected in thegoverning equations. A two term absorption law was adopted to account for thedependence of attenuation on wavelength. The absorption of radiation is assumed tofollow Beer's Law, which for one spectral band is given as follows:Iz = 10 e-Xz^ (2.2)where, Io = incident (non-reflected) radiation= ( 1 - a) • Q0 , where a = albedo; Q0 = total incoming radiationIz = radiation remaining at depth z below the surface= attenuation coefficientThe attenuation coefficient is dependent on the radiation wavelength and the transparencyof the medium through which the radiation passes and the radiation wavelength. Theappropriate form for radiation characterized by two spectral bands is:Iz = 10 (A 1 e-X1z + A2 e -X2z)^ (2.3)where, A = fraction of total radiation in a given spectral band (A 1 + A2 = 1)subscripts 1 and 2 refer to the two spectral bandsPatterson & Hamblin (1988) set Al = 0.7, the fraction associated with visible radiation,and A2 = 0.3, that associated with infrared radiation, in accordance with Kirk (1983) andBoer (1980). It has been found that an increase in the number of spectral bands beyondtwo does not result in significant change in the distribution of heat below the ice andsnow cover (Patterson & Hamblin, 1988).The steady-state one-dimensional heat flow equation for a material with a uniform heatsource is:a2TK ^ + 4 = 0az2where, T = material temperatureK = material conductivity= heat generated per unit volume(2.4)For an ice and snow cover with a depth-dependent solar heat source (see Figure 2.1), theequivalent equations are:a2TKs az2 s + Al ks1 10 exP[-Xs1 (hi + hs - z)] + A2 Xs2 lo exP[-X52 (hi + hs-z)] = 0,^hi + hs^z^hi21 TiAlK.—az2 + ^Xii Io exP[?s1 hs - Xil (hi - z)] + A2 2 i2 Io exP[ -Xs2 hs - Xi2 (hi - z)] = 0,^z^0 (2.5)where, X = attenuation coefficienth = material thicknesss = snowi^= ice12Qo Qelo 1Qc RhH 1air4^snowhsI^icehiofwaterqwFigure 2.1. Heat Fluxes Across an Ice and Snow CoverThe appropriate boundary conditions for the above equations are:Ti = Tf^Z = 0Ti = Ts.m._ i^= ts.aTi .„, s ars } z . h iaz^azTs = To^z= hi+ hs^(2.6)where, To = the surface temperature of the top component of the coverTf = the freezing temperatureFor the given equations to hold, the temperature distributions must stabilize on a timescale which is considerably shorter than that of variations in To and the snow and icethicknesses. This condition is satisfied provided that the Stefan Number is less than 0.1as previously described. The thermodynamic balances at the surface and ice-waterinterface can be treated independently, and an analytic solution is possible. As the13component temperatures are not required, the equations need only be solved for q 0, theconduction of heat at the air-snow interface:aTsK (2.7go = - saz z = hi + hs^)The solution can be written as:CsKi+hiK1t^sksii  1 -exp(-2t.sihs) +^KiXiiexp(-Xsihs)[1-exPeXiihi)l}KiKs (clo-Io) = Tf - To - A go KA 21T {1 -eXP('IS2hS) + eXP(-XS2h0[ 1-eXP(-2q2h1- rt0^ (2.8)^KsXs2^KiXi2The solution reduces to the correct form in the absence of snow. It should be noted thatthere is a typographical error in this solution as given in Patterson & Hamblin (1988) (thecorrect form is in the DYRESMI code). The two terms representing the attenuation ofshort wave and infrared energy within the cover should result in an increase in q o for anincrease in any given attenuation coefficient. For this condition, the sign preceding thesetwo terms must be negative, not positive as shown in Patterson & Hamblin (1988).The components of the ice and snow cover heat budget are as shown in Figure 2.1. Thethermodynamic balance at the surface depends on the meteorological forcing and on q 0 ,the flux of heat in the top component at the interface with the air. The surfacetemperature, To, will adjust itself so that a heat flux balance is achieved. T o, however isbounded by Tm , the melting temperature. This provides the condition for surfacemelting. The surface condition may be expressed as:go(To) + H(To)^= 0 ,dh= - PLa 'To < TmTo = Tm^(2.9)14where, H(To) = net incoming meteorological flux= Rli - R10(To) + Qc(To) + Qe(To)R11 = incoming long-wave radiationR10 = outgoing long-wave radiationQc = sensible heat transfer between surface & atmosphereQe = evaporative heat fluxThe above approach differs from that of Maykut and Untersteiner (1971) in that the solarradiation term enters the surface balance implicitly through the evaluation of the heat fluxin the upper component, q0 . These authors specify a fraction of short-wave energyassumed to be available at the surface, which would produce a less accurate thermalprofile in the cover. The model does not allow for ice growth at the surface. Snowthickening can only occur due to snowfall.A criterion, based on a simple hydrostatic force balance, is given for the formation ofsnow-ice. If more snow falls than what the ice can support, the ice will bend and crack,resulting in flooding of the snow cover. In sufficiently cold conditions, this water willfreeze, creating snow-ice. This phenomenon is discussed further in §2.3.3.Ablation and accretion of ice may occur at the ice-water interface only. The flux of heatat this point, qf, depends on Tf, the freezing temperature, and the surface conditions. It isobtained from the solution of Equation 2.5:of = -Kiaz= qo - Ai lo ( 1 - exp-(Xsihs+ Xiihi)) - A210( 1 - exp-Otsihs+ Xiihi)) (2.10)Independently, the heat flux from the water to the ice, q w , depends only on the conditionsin the water column. Any imbalance between of and qw results in the freezing or meltingof ice:15aTiz=0•dhidtof - qw = Pi Li (2.11)16The flux of heat from the water to the ice may be due only to conduction if the water canbe considered stagnant. In this case:, dTw Iqw = —'w dz I z_o (2.12)where, w = waterThe solution properties associated with the given formulation are described in Patterson& Hamblin (1988).The above equations have been incorporated into DYRESMI to account for the transferof heat across the ice and snow covered lake surface. Further investigation was requiredin order to tailor the given formulation to the small mid-latitude lakes with which thisproject is concerned. Areas of particular concern include the transfer of heat stored in thelake sediments over the summer, the effect of rainfall on the ice and snow-cover, and theformation of snow-ice due to excessive snow-cover. In addition, the great variability ofalbedo observed at mid-latitudes (see Strickland, 1982), required the incorporation of anappropriate albedo model, rather than the use of a constant value, as assumed in versionof DYRESMI reported in Patterson & Hamblin (1988). (Another version of DYRESMInow exists which employs the empirical time dependent albedo relationships derived bythe U.S. Army Corps of Engineers (Hamblin, pers. comm., 1992). These relationshipsare given in §2.3.2.)2.2 HEAT TRANSFER FROM BOTTOM SEDIMENTS2.2.1 Review of Observed Rates of Heat TransferHeat transfer from the lake bottom to the overlying water has been identified as beingimportant in the circulation of small lakes over the winter period (see Birge et al., 1927;Mortimer & Mackereth, 1958; others). Considerable amounts of heat are stored belowlakes during the warm season (Ashton, 1986). Thandertz (Ashton, 1986) computed heatfluxes at Lake Velen (590N, mean depth, --z- = 6.5 m) in Sweden from measurements oflake water warming during periods of ice cover. Average heat fluxes of 1.0, 1.9 and 1.8W/m2 were computed over 3 consecutive years of observations. Maximum heat fluxes inthe order of 4 W/m2 were calculated for the month of January based on measuredsediment properties and thermal profiles. This maximum flux decreased in a linearfashion throughout the winter period. The classic work of Birge et al. (1927) regardingthe sediment heat budget of Lake Mendota in Wisconsin, U.S.A. (43 0N, i = 12.4 m)produced results consistent with the Lake Velen findings. Sediment properties andthermal profiles were also measured at Mendota, and the investigation revealed aconsiderable range in heat transfer rates. The average winter period (December 15th toApril 1St) heat transfer rate of 2.9 W/m 2 , is considerably higher than the average ratescalculated for Lake Velen due to higher summer temperatures. Mortimer & Mackereth,(1958) however, used the thermal profiles given in Birge et al. (1927) to estimate anaverage value of 1.6 W/m2 for the actual ice covered period which began on January 1st.Welch and Bergmann (1985), for Methane (z max = 14 m) and No-Fish (zmax = 13 m)Lakes (Canadian NWT, 630N) estimated a heat flux of 1.5 W/m2 based on averagetemperatures at and 2 cm below the sediment surface, assuming a thermal moleculardiffusivity of 1.2 x 10 -3 cm2/s. The range in fluxes found in the literature is consistent17with the theoretical model by O'Neill and Ashton (1981) which computes heat transferrates based on monthly average air temperatures.Estimated rates of heat transfer in the Western Basin of Lake Tornetrask in SwedishLapland (680N, zmax = 169 m) were much lower, in the range of 0.2 to 0.4 W/m2(Mortimer & Mackereth, 1958). The lower rate is expected given the great depth of thelake, and the relatively cool summer temperatures (the maximum summer temperature ata depth of 1 m was 11.8 0C on average between 1921 and 1929, compared with about260C at Lake Mendota).2.2.2 Methods of Estimating Sediment Heat TransferQuantitative sediment properties and thermal profiles are not available for the lakesinvestigated in this study. In order to estimate sediment heat transfer, use of publisheddata and empirical equations are required.In general, the shallower the lake, the more heat is stored in the sediments over thewarming season. Morphometry, altitude and latitude are also factors which influenceheat storage. As a rough estimate, however, Falkenmark developed an empiricalrelationship which predicts the ratio of sediment to lake heat budget based on the meandepth only (see Ashton's [1986] Figure 4.8). This relationship could be used to estimatesediment heat transfer rates.A second method of estimating sediment heat transfer rates is available. Thermal profilescharacteristic of lake sediments change markedly over the year in response to changinglake temperature and as a function of lake depth. Temperatures below about 2 m underthe sediment-water interface remain, however, more or less constant over time (see, forexample, Ashton, 1986). The temperature at 2 m is often close to the average annualtemperature of the lake water above. As an approximation of the thermal gradient at the18sediment-water interface, it is reasonable to employ the following formula (Hamblin,pers. comm., 1992):dT Ty - Tw dz — zsed^ (2.13)where, Ty^= the average annual whole lake temperatureTw^= the current water temperaturezsed^= distance over which sediment temperature (linearly) increases fromT to TY^w= 2 mThe above formulation implies the use of the conduction equation to estimate sedimentheat transfer. This necessitates an estimation of sediment conductivity based on dataavailable in the literature. In most studies reported in the literature, only observations ofthermal diffusivity have been made. In each of these cases, to determine conductivity, itis necessary to estimate the density and heat capacity of the sediments.2.2.3 Sediment propertiesBirge et al. (1927) made sediment temperature profile measurements which wereconsistent with the assumptions that heat was transported by molecular diffusion, andthat the thermal diffusivity, lc, of the sediments equaled that of still water. Investigationsby McGaw (O'Neill & Ashton, 1981) and Likens & Johnson (1969) are in agreementwith these observations. Hutchinson (1957), however, reported values of 32.5 x 10-8 m2/sfor Lake Mendota, or almost 2.5 times that of pure water (Kw = 13.6 x 10-8 m2/s at 50C).He also cited Neumann (1953), who determined a value of 40 x 10 -8 m2/s for sedimentswith some sand content.Urban & Diment (1988) have reported conductivities in the range of 1 to 3 times that ofwater at Clear Lake California, the higher values being associated with the coarser-grained sediments. Likens & Johnson (1969) found that the gelatinous ooze consisting of1920partly decomposed organic matter in the sediments at two Wisconsin lakes had the sameconductivity as water. These latter two results suggest that the density and heat capacityof many water saturated sediments are comparable to water.2.3 EFFECTS OF SNOW COVER2.3.1 Snow PropertiesThe physical and optical properties of snow (density, conductivity, albedo, attenuation)are extremely variable especially at mid-latitudes. Because of the great importance ofthese properties on the heat transfer to the water, average values used as universalparameters (as is the case in the DYRESMI model) are not generally appropriate for thelakes in question. Here attention will be given primarily to density and albedo. A simplerelationship is available which relates conductivity to density, and attenuation will beassumed constant owing to a lack of information.DensityThe density of snow varies widely depending on the following processes:a) Heat transfer (convection, condensation, radiation, conduction)b) The weight of overlying snowc) Wind drifting and compactiond) The temperature and variation of water content within the snow packe)^The percolation of melt waterSnow density may vary from as little as 10 kg/m 3 for dry falling snow, to as high as 380kg/m3 for wind-toughened, well compacted snow. The presence of liquid water mayincrease the density of melting snow to as high as 500 kg/m3 before it drains away(McKay, 1968). Patterson & Hamblin (1988) used a constant snow density of 330 kg/m 3for DYRESMI. The density of the snow varies directly with depth, but layering leavessignificant variations. Initially there is rapid settling following a snowfall. The highestrates of densification occur within the first few hours after deposition. Snow which hasaccumulated at an average of 46 kg/m3 has been observed to increase to 176 kg/m3 after24 hours of drifting. (McKay, 1968). Changes in form and displacement of particles areresponsible for the settlement (USACE, 1956). Much attention has been given to thedensity of newly fallen snow, because of its high variability. A simple air temperaturedependent relationship has been derived from observations in the Rocky Mountains inColorado (Grant & Rhea, 1973):Ps = (0.0785 + 0.00219 Tail- + 0.00023 Taut) Pw^ (2.14)In the range of 0 to -15 0C, this equation, on average, predicts a density of about 80kg/m3 .ConductivityThe conductivity of snow increases substantially with density (Ashton, 1986):Ks = 0.021 + 4.2 x 10 -3 p s + 2.2 x 10-9 psi^(2.15)Furthermore, conduction of heat from the underlying ice to the air will also increase asdensification proceeds because of the associated reduction in thickness of the snow layer.AttenuationFew data are available which quantify the attenuation of radiation through snow at mid-latitudes. Patterson & Hamblin (1988) use coefficients of 6 and 20 m -1 for visible andinfrared radiation respectively. For the shortest wavelengths, however, Grenfell &Maykett (1977), found coefficients of about 18 m -1 for compact dry snow in the ArcticBasin.21Da— 100_ 10 (1.05 - 0.07d)Act — 100_ 10 (0.78-0.69d)2.3.2 Estimation of Snow AlbedoIn DYRESMI, the snow albedo is set to 0.85 when the air temperature is below zero, and0.6 when it is above zero. As described earlier, however, rapid and extreme weathervariations should result in a wide range of snow albedos at mid-latitudes. A variablealbedo should therefore be incorporated into the heat transfer model. Typical ranges andfactors which determine albedo will now be described.Albedo is commonly taken to decrease exponentially from about 0.8 for fresh snow, toabout 0.4 for melting, late season snow (Gray, 1970). Grenfell, Perovish & Ogren (1981)report a significantly higher range, with values decreasing form 0.93 to 0.63. Morerecently, Henderson-Sellers (1984) tabulated a wide range of range of albedos for variouslake surface conditions (Table 2.1).A simple empirical relationship developed by the U.S. Army Corps of Engineers (seePetzold, 1977) accounts for the reduction in albedo due to compaction, fine flakedestruction, accumulation of dust and debris, and, for melting conditions, themetamorphosis from powder to granular snow. Two equations, called the USACE decayfunctions, were derived; one for an accumulating snowpack, and the second for a meltingsnowpack. Both express the change in albedo only in terms of the number of days, d,since the last snowfall:(accumulating)(melting)^ (2.16)22With the exception of accumulation of impurities, the processes listed above all tend toTable 2.1 Albedo Under Various Surface Conditions (Henderson-Sellers, 1984)Surface Condition I^Albedo1. Freshly formed, snow-free black ice 0.022. Snow-free snow-ice 0.4-0.53. Freshly fallen snow 0.8-0.954. Old snow of moist wet snow 0.7-0.85. Wet snow 0.6-0.76. Slushy, gray melted snow layer 0.2-0.67. Completely wet top surface (ponding) 0.058. Water 0.06increase the snow grain size. The following is a calibrated equation for albedo, a, interms of grain size (Bohren & Beschta, 1979):a = 1 - 5.964ji • dg^(2.17)where, ji = the wavelength dependent absorption coefficient for pure homogeneousicedg = the average snow grain diameterRobinson & Kukla (1984) monitored the albedo of a melting snowpack for over onemonth following an early February snow storm in the north-eastern U.S. For openfarmland they found that the albedo dropped at an increasing rate from a maximum of0.8. The opposite trend, observed for both urban and shrubland areas in the same study,is implicit in the USACE decay functions. Gray & Male (1981) explain that, for shallowsnow cover, the albedo drops at an increasing rate due to the increasing influence of therelatively low albedo of the underlying ground or ice. The USACE functions were23derived from deep mountain snowpack data, and are not influenced by the opticalproperties of other materials below the snow. According to Choudhury (1982), thesnowpack can be treated as semi-infinite when it is greater than only 10 cm in depth. Thislimiting depth, however, increases with snow grain size, and may be up to 50 cm for oldmelting snow (Wiscomb & Warren, 1980). This latter observation appears to beconsistent with the results of Robinson & Kukla (1984).Petzold (1977) found that USACE functions consistently under-predicted the actualalbedo for the 1969 summer field season on Meighen Ice Cap, N.W.T. He derived thefollowing quadratic correction function:3 86 + 0 380 d + 0.123 d 2 Aa — (2.18)100The above correction indicates that, as implied by Robinson & Kukla (1984), and Gray &Male (1981), a simple decay function may not be appropriate, at least for the data set towhich the correction was applied.In order to reduce the scatter produced when the USACE functions were applied to thedata set used by Petzold (1977), further corrective functions were developed in order toaccount for cloud conditions and solar angle. It is commonly reported that albedoincreases with increasing cloud cover, opacity and with decreasing cloud height. Albedodecreases, however, as the solar angle increases up to about 40 degrees, beyond which, itremains fairly constant. Strickland (1982) found that 14% of albedo variability in sunnynon-melting conditions could be explained by solar angle. Diurnal or seasonal trendscaused by changing solar angle, may, however, be obscured by variable cloud cover andaging snow (Bergen et al, 1983). The effect of solar elevation may be accounted forusing Petzold's (1979) Figure 2a, to which the following expression has been fitted:24ea — -1 .9 + 24.8•0/ 15 . 5100 (2.19)where, [3 = solar noon anglePetzold (1979) also performed a regression analysis, using data from four polar stations,to derive a simple relationship expressing the change from clear sky albedo to a valueunder a given cloud cover:Da — 0.449 + 0.0097• (10•C) 3100 (2.20)where, C = fraction of sky obscured by cloudIn general, he found that the greatest errors in the complete albedo model occurred ondays with fresh snow. He suggested that rain mixed with the falling snow may have, onsome occasions, resulted in lower than predicted albedos.When applying the USACE functions it is generally assumed that the maximum possiblealbedo is 0.84. Petzold (1977) however, suggests that a new snowfall albedo between0.84 and 0.89 is reasonable, and furthermore, this value is likely to vary between 0.90and 0.95 at low sun angles. Strickland (1982) has observed natural snow albedos up to0.95 at solar noon near Peterborough, Ontario over the 1980-81 snow season. Bergen etal. (1983) have also documented late winter albedos in excess of 0.96 on a level clearingin the White Mountains of Arizona, at an elevation of 3000 m. Over a one week period,the daily average albedo never dropped below 0.92. Snow purity also has a profoundinfluence on albedo. Fresh, pure snow may reflect more than 95% of incident radiation(see Wiscomb & Warren, 1980).Strickland (1982) observed a dramatic drop in albedo following mid-winter rainfall. Theeffect of rainfall is not considered in either the USACE functions or Petzold's (1977)computer model. Strickland (1982) found the accuracy of USACE functions to be adisappointing 48% for her data set. She implies that they may not be appropriate for25latitudes in the range of 40 to 50 degrees north, due to rapid and extreme weathervariations. Better agreement with actual albedo may be expected at the lakes consideredin the present study, however, given that the USACE equations were developed usingdata from high elevation snowpacks.2.3.3 Snow-ice FormationIf the weight of the added snow is too great, the underlying ice will crack and floodingwill occur. This results in the formation of snow-ice when the flood water freezes. Asimple force balance, based on buoyancy theory, is used to determine the maximumthickness of snow that the ice will support (Patterson & Hamblin, 1988):hSMhi•w - pi)Ps(2.21)In the DYRESMI formulation, the excess snowfall is reduced in thickness, according tothe ratio of snow to ice density, and added to the ice layer.The water which seeps through cracks in the ice caused by the weight of the overlyingsnow and into the snowpack constitutes a significant source of heat in the snow layer.Sensible heat is given up to the snow as the flood water temperature is cooled to thefreezing mark. Under most weather conditions following snowfall, the water will thenfreeze, and snow-ice is formed. This latter process releases latent heat to the surroundingice and snow.2.3.4 Snowmelt Due to Rain EventsSnowmelt due to rain may be an important factor in the depletion of a snowpack. Thispossibility was ignored in the DYRESMI formulation because it is a rare event at highlatitude lakes (Patterson & Hamblin, 1988). When rain falls on snow, the sensible heatassociated with the rain water is given up to the snow. If the heat transferred is more than26sufficient to raise the snow temperature to 0°C, then the excess heat results in snowmelt.The heat released to the snow is given as follows:Qrs = Cpw Pw Fair - Tm) • P^ (2.22)where, Qrs^= sensible heat released to snowP^= total rainfallCpw = heat capacity of waterIt is often assumed that when rain does occur, the snowpack is isothermal at 0°C and allof the heat associated with the rain is used to melt snow (USACE, 1956; Harr, 1981).Harr (1981) uses this assumption in illustrating that rain is only a significant contributorto snowmelt in forested watersheds when total daily rainfall is in the order of severalcentimetres. If a sub-freezing snow pack were assumed, melt due to rain would be shownto be even less significant. However, in this latter case, the rain water would freezeinside the snowpack, thereby releasing latent heat.The relative importance of latent heat is illustrated by the following example: 46 mm ofrain at 0°C will increase a 1.8 m snowpack (p s = 400 kg/m3) from a mean temperature of-10°C to 0°C. If, however, the same quantity of rain were at 8°C and fell on a snowpackat 0°C, less than 0.5 cm of meltwater would be produced. Although some snow is melteddirectly by the rain, the associated high humidity and air temperature results in heattransfer to the snow dominated by condensation of water vapour. When precipitation isless than 13 cm/day (less than 2 cm of rain falls on average between December andMarch in the B.C. Southern Interior Plateau), condensation is 4 times more importantthan the sensible heat released from rainfall (Harr, 1981). Condensation will occur inwarm, humid conditions when the vapour pressure of the air above the snow exceeds thesaturation vapour pressure at the temperature of the snow surface. Condensation is27handled by the bulk aerodynamic equations adopted for use in the DYRESMI heatbudget (see Chapter 4).2.4 EFFECT OF ARTIFICIAL CIRCULATION2.4.1 Review of Systems Used WorldwideArtificial circulation has long been used to prevent or remove ice cover in ports andmarinas to facilitate their use year-round. These systems, in general, differ from thatused at Menzies Lake in three ways:a) the circulated water may be salineb) the water is only a fraction of a degree above freezingc) the system is open: there is effectively an infinite body of water on which to draw.Very little quantitative information, however, is available regarding their behaviour orimpact.Eranti et al. (1983) describe the use of flow developers in Finnish ports where theaverage harbour temperature is 0.1 0C above freezing. It was found that the thermalreserve was sufficient for the prevention of freezing when surface currents are greaterthan 0.6 m/s. This is consistent with the findings of investigations of ice cover on rivers(Stigebrandt, 1978). The polynyas generated, however, varied significantly in sizedepending on the local water temperature and ice conditions. Furthermore, the surveysuggested that the presence of a warm water discharge may be more important inincreasing the polynya size than the flow developer engine power. This latter point isillustrated in Table 2.2.28Table 2 2 Polynya Areas in Finnish Harbours (Eranti et al., 1983)Engine Power (kW) I^Ice-Free Area (m2)11 17011 20015 5030 30030 100030 30001f warm water discharge (Q=0.13 m3/s, T= 70C)Predictive modelsThe only known model which predicts the size of an artificially-generated polynya is thatdeveloped by Ashton (1979). The model predicts the flow induced by a point sourcebubbler system and the associated heat transfer to the ice cover. An earlier model(Ashton, 1974) deals with the problem of a line source bubbler. In the new model, air isdischarged from a point at some depth below the water surface. Water is entrained as theplume rises to the ice-water interface at which point it is redirected and spreads radially.The temperature of the plume is of primary importance and is calculated by integratingthe product of temperature and the rate of increase of flow over the vertical distance, z,through which the plume rises:HwTwH - Tm — ,-,^ of(Tw(z) - Tm) dQ(z) dzY wH dzwhere, TwH = temperature of the plume at height H above the bubblerQwH = flow at height HTm^= melting temperature(2.23)29A coefficient is developed for the turbulent transfer of heat from the water to the iceusing the results of Gardon & Akfirat (1966), who investigated the heat transferassociated with an axisymmetric air jet impinging on a flat plate. This coefficient, h, isgiven as a function of radial distance, r, from the point of impingement:h 2.08 Kw UcH0.55 r -0.45^ 24)—^ (2.r0.55where, UCH = The centerline water velocity at the point of impingementv^= kinematic viscosityThe variation with r is an empirical result from the work of Gardon & Akfirat (1966), andis not based on boundary layer heat transfer.The actual melting of the ice cover is handled by employing a much cruder model thanthat of Patterson & Hamblin (1988). Ashton's (1979) model was applied to a mid-western location in the U.S. for a water temperature of 0.2°C. It is unclear whether thistemperature is the mean ambient temperature, the mean plume temperature, or an initialcondition. It is also unclear if buoyancy was of significance in the application. Themodel predicted a polynya which responded dramatically to changing air temperature.The polynya radius varied from 0 to 4 metres over a 24 day simulation.Heat loss from a natural polynya in polar pack ice has been given a considerable amountof attention, but few studies have attempted to quantify rates of lateral ice growth in theopen water. Den Hartog et al. (1983) have predicted heat loss from an Arctic polynya inthe Canadian Archipelago. Good agreement was achieved when comparing net heat losscalculated using bulk aerodynamic formulae to heat flux tray measurements. No attempt,however, was made at predicting polynya size. Bauer & Martin (1983), however, haveexamined the growth of grease ice across polynyas in some detail. Grease ice, whichconsists of about 30% ice and 70% seawater by mass, grows in the upwind direction30under high wind conditions, from the downwind end of the lead. Ice forms as frazil discswhich are advected downwind to pile up against the lead edge to depths of 0.05 to 0.3 m.The ice growth model developed by Bauer & Martin (1983) predict high rates of growthbecause the exposed water is at or very close to its freezing point. The ice pile up depthis given as functions of wind speed parameters and fetch. For a constant wind speed of10 m/s and an initial fetch of 50 m, the average predicted grease ice thickness is 0.07 m.The predicted thickness increases with both wind speed and fetch.Lake AerationArtificial circulation of lakes and ponds in winter can be an effective method ofmaintaining dissolved oxygen concentrations which are required to support fish life inwinterkill lakes. This fact has been well established (see Ashley, 1987; Bandow, 1986;Lackey & Holmes, 1972), but little information regarding thermal impact, ice thickness,and polynya behaviour is available. It has frequently been observed, however, thataverage lake temperatures are lower than would be expected under natural conditions(Lackey & Holmes, 1972).2.4.2 Heat Transfer Between Water and AirThe heat flux at the air-water interface of the polynya is determined in the same manneras that at the snow-air interface, with To replaced by Tw , the water temperature, and thesurface albedo changed to a value appropriate for the open water. Once this flux is knownthe net effect of the polynya on the lake heat budget can be calculated. In addition toshort wave radiation, sensible and latent heat, and the balance of long wave radiation, thecontribution of free convection, and cooling due to snowfall are required in the balance.A further consideration is the turbulence induced by the circulator at the air-waterinterface. Sensible heat transfer measurements made in the ocean, however, have shownthat sensible heat transfer is a function of wind speed, but is independent of sea state31Tv — T+273 - 273 (2.26).378svnPatm )(Kraus, 1972). Given this finding, it is assumed that the turbulence generated bycirculator is only important insofar as it increases the heat transfer to the ice cover (see§2.5.3).Free convection due to unstable humidity profiles, is considered important in coolingponds, Arctic leads and polynyas, where the air temperature is on average much less thanthe surface water temperature, and wind speeds are low (Ryan et al., 1974; Den Hartog,1983). Its contribution to the evaporative heat loss is given as follows (Ryan et al.,1974):Qfc = X•(svp0 - svpd)•Twv - Tav^ (2.25)where, svp0 = saturation vapour pressure at the water temperature, Tw (mbars)svpd = vapour pressure of the air (mbars)Twv^= virtual surface water temperature (0C)Tav = virtual air temperature (0C)x = 2 . 7 Wm-2mbar 1(pc)-1/3Since water vapour is lighter than air, evaporation increases the driving buoyancy forcedriven by the temperature difference T w - Tom. This effect is accounted for by employingthe virtual temperature difference, Twv - Tay. The virtual temperatures are calculated asfollows:32where, Patm = atmospheric pressure (mbars)Reservoir operators have often observed a reduction in water temperature followingsnowfall over open water (Ashton, 1986). The following equation may be employed:Qsp = L Pw Is - Cpi Pw (Tair - Tm) Is^ (2.27)33where, Qsp^= heat removed from the water (W/m 2)Is^= snowfall (mm/hr, rain equivalent)The first term accounts for the latent heat of fusion, and the second for the heat used inwarming the snow to the melting temperature.The net heat flux, Hp , at the air-water interface is calculated as follows:Hp = Iwp + Rli - 12.10p - Qcp - Qep - Qef - Qsp (2.28)where the subscript p refers to the polynya.2.4.3 Heat Transfer Between Water and IceIn present sea ice models, the transfer of heat between water and ice is often ignored orset to a constant value, even in the vicinity of polynyas where horizontal variations canbe significant (Hamblin & Carmack, 1990). Gilpin et al. (1980) parameterized theturbulent vertical heat transfer using Newton's law of cooling:q = h (Tb - T.)^ (2.29)where, q = local heat transfer per unit areah = convection heat transfer coefficientTb = boundary temperatureT. = temperature of the free streamA more useful form of this law is the bulk aerodynamic formulation, in which heattransfer is given as a function of a dimensionless coefficient, C t, and the flow velocity, U,at some distance away from the boundary (Hamblin & Carmack, 1990):qt = Ct p cp U (Tw - Tm)^ (2.30)where all other terms are as previously defined.From Figure 13 in Gilpin et al. (1980) it can be found that, under experimentalconditions, C t varies from 0.6 x 10 -3 to 1.0 x 10-3 depending on the ice roughness. AsHamblin & Carmack (1990) point out, however, the scales of turbulent motion associatedwith laboratory flow fields are generally not representative of lakes or oceans.Furthermore, the use of cross-sectionally averaged velocities in Equation 2.30, is notapplicable to these water bodies.To determine the appropriate depth where velocity measurements should be taken,Hamblin & Carmack (1990) make reference to the standard height above ground used inmeteorological applications and scale the analogous distance below the ice-waterinterface by:(Pair )4). 5^ (2.31)Pwwhere, pair = air densitypw = water densityThis factor is based on the following requirements:a) stress continuity across the boundaryb) proportionality of boundary layer thickness to friction velocity.Given a standard above-ground height of 10 m, the appropriate below-ice depth is about1.0 m. Using under ice-current and temperature measurements in three large northernlakes with significant flow-through, Hamblin & Carmack (1990) found C t to be in therange of 0.5 x 10 -3 to 1.1 x 10-3 .34Boundary Layer DevelopmentThe details of boundary layer development may, in some applications, be of importanceto the net turbulent heat flux between the water and the ice. Hirata et al. (1979) examinedthese details in their laboratory investigation of the steady state profile of an ice layer in aforced convection flow of water. It was found that the Nusselt Numbers throughout theturbulent regime were in the order of 35% larger than those given by the von Karmancorrelation for turbulent flow over a flat plate. The von Karman relationship is given asfollows:Nux = 0.0296 Pr1 /3 Rex4/5^(2.32)where, Nux = Nusselt NumberPr^= Prandtl NumberRex = Reynold Numberhx x- KwCpw p.- KwU. x vx^= distance from leading edge of platehx^= convective heat transfer coefficient at location xU.^= free stream velocityTherefore, from the results of Hirata et al (1979), the equivalent expression for an icelayer would be:Nu x = 0.039 Pr 1/3 Rex4/5^ (2.33)In terms of the convective heat transfer coefficient, h x , the above formulation isexpressed as:KwxPr 1 /3 hx = 0.039^Rex4/5 (2.34)35Assuming that U.. is equivalent to U in Equation 2.30, then C t may be estimated usingthe above result:0.039 K Pr1/3 Rex4/5Ct —p Cpw U x(2.35)This formulation implies that C t is not constant unless U = f(x -1 ). For laminar flow overan ice sheet, Hirata et al (1979) found that the ice thickness varies as They did notprovide a formulation for the ice thickness in turbulent flow, but from their Figure 3, it isclear that the ice thickness increases much more slowly than it would in laminar flow.This is consistent with the field observations made by Stigebrandt (1978). This latterinvestigator found that, downstream of a river inflow into an ice-covered Norwegianlake, the ice thickness increases gradually towards the lake proper.Gilpin et al. (1980) explain that the heat transfer rates for forced convection over an icelayer are greater that those predicted by flat plate boundary layer theory because ofirregularities that will develop in the ice layer profile as a result of variability in the heattransfer coefficient. Since heat transfer coefficients increase with surface roughness(Holman, 1990), it is clear that greater heat transfer rates would be expected across anirregular ice surface, than across a smooth plate.2.4.4 Flow FieldThe flow field generated by the Air-o-lator ® is complicated by the following conditions:a) in lakes of sufficient depth, density stratification requires that the internalhydraulics be consideredb) fountain-like pumping of water into airc) irregular lake morphometryd) impingement of radial jet, in some directions, on the float used to support thecirculation device.36As described in Chapter 1, the first condition may be ignored owing to the small depth ofthe lakes considered. Given that the flow field model will be based on near-surfacevelocity measurements, the mechanics of the second problem may be ignored. It isassumed that the near-surface velocities are proportional to the pumping rate Q, andtherefore, an empirical formula based on measured velocities will apply for Air-o-lators®of various sizes. It is also assumed that the circulator operates in a semi-infinite body offluid, and that the radial jet velocity is not dependent on the direction of flow.The flow field generated by the circulator is modeled as two independent components.First, the radial flow to the intake at some distance below the water surface is assumed tobe governed by potential flow theory. Applying the semi-infinite body assumption andthe method of images, it can be shown that the intake velocities are small near the surfacewhere the outflowing jet velocities are important. The circulator is considered to beimportant to the heat budget only in terms of the size of polynya generated and thetransfer of sensible heat to the ice-sheet. This means that only the near surface velocitiesneed be considered in the model, and the intake flow field can therefore be ignored.This leaves the problem of modeling the outflowing radial jet produced by the circulator.The resulting flow field is assumed to be consistent with submerged, non-constrained,radial jet theory as described by Rajaratnam (1976) and Witze & Dwyer (1976). Thevelocity generated along the jet centerline varies as r 1 . This result is derived from thegoverning equations of motion. The standard semi-empirical relationships which givethe centerline velocity as a function of the initial flow rate, are also dependent on theoutlet geometry. Tolmeinn & Goertler have also derived solutions that require empiricalcoefficients for the non-centerline flow (Rajaratnam, 1976).372.5 PREVIOUS WORK AND THE MLI MODELThe work of Patterson & Hamblin (1988) constitutes a major step forward in the area ofice-covered lake modeling. Their model will be utilized in this study and extended toaddress the issues which have been introduced in this chapter:a) rainfall,b) heating of the snowpack due to flooding of the ice-cover,c) sediment heat transfer,d) variable albedo,e) variable snow density.In addition to these factors which may be important with regard to the heat budget of theice-covered lake in the natural state, the DYRESMI model will be extended to include theimpact of the Air-o-lator®. The heat transfer associated with snowfall on open water andfree convection (see §2.4.2) will be added to the list of bulk aerodynamic equationscommonly used in lake modeling. A more difficult problem is that of predicting thepolynya size and the turbulent heat transfer from the water to the ice. The predictivemodel developed by Ashton (1979) is particular to bubbler aeration systems and is basedon empirical results involving an air jet impinging on a flat plate and not on boundarylayer heat transfer considerations. Furthermore the ice melting formulation developed byAshton (1979) is much cruder than that of Patterson & Hamblin (1988).Given the lack of a satisfactory model for the prediction of polynya size, development ofa new model is required. Ice melt in the vicinity of the circulator is a function ofturbulent heat transfer from water to ice as described by Equation 2.30:38qt = Ct p Cp U (Tw - Tm)^ (2.30)The velocity field associated with the radial jet produced by the circulator will beassumed to provide a reasonable measure of the dependency of turbulent heat transfer onradial distance from the circulator. The heat transfer coefficient C t is expected to requirecalibration given the unique nature of the problem but should be of the same order as therange found by Hamblin & Carmack (1990) for river flow through ice-covered lakes.The details of the modifications to DYRESMI are dealt with in Chapter 4. Verificationof the new model and calibration of parameters requires a thorough field investigation ofice-covered lakes in a natural state and artificially circulated. This investigation isdescribed in detail in the following chapter.39Chapter 3FIELD INVESTIGATION3.1 CHARACTERISTICS OF STUDY LAKESThe two lakes which were selected for this study are both located in British Columbia'sSouthern Interior Plateau, about 300 km north-east of Vancouver, and 15 km south-eastof Merritt (Figure 3.1). Both are small, naturally eutrophic lakes which experiencesignificant rates of oxygen depletion below the thermocline in summer and under ice-cover in winter. Harmon Lake is the centrepiece of a small provincial park and is verypopular among recreational anglers. Menzies Lake, 5 km to the north of Harmon, hasbeen set aside by the B.C. Ministry of the Environment for Fisheries Research purposes.Harmon, which is left in its natural state over winter, was selected as a control lake, forcomparisons with Menzies Lake, which is artificially circulated.3.1.1 MorphometryAs reliable morphometric information was only available for Harmon Lake, a survey ofMenzies Lake was conducted in August 1991. This was accomplished by taking echo-soundings along 13 transects from a boat moving at constant velocity. Bathymetric mapsfor both lakes are shown in Figure 3.2. A summary of morphometric information derivedfrom these maps is given in Table 3.1. Hypsographs for both lakes are given in Figure3.2b. The wind sheltering coefficients were determined following USACE (1986). Atransit was used to determine the average tree height at the end of each lake from whichthe dominant winds blow. For abrupt changes in relief at the water's edge, it will take adistance of roughly 8 times the height of the relief before the wind field and resultant40BRITISHCOLUMBIA41Figure 3.1 Field Study LocationFigure 3.2 Bathymetry of Menzies and Harmon LakesMenzies Lake42Scale: 1:3620Depths given in metresHarmon Lake43stress reattach to the lake surface (USACE, 1986). The coefficient was calculated asfollows:1 8 illsw = 1. - Lf (3.1)where, Sw = wind sheltering coefficientH = average tree height at the water's edgeLf. = effective length of lake or fetch in dominant wind directionThe surrounding topography (along the dominant wind axes) is sufficiently flat as to havelittle influence on S w .Table 3.1 Morphometry of Menzies and Harmon LakesParameter Menzies Lake 1^Harmon LakeElevation (m) 1250 1150i (m) 8.2 8.0Volume (m3 ) 390 000 2 480 000Surface Area (m2) 47 500 311 000zmax(m) 16.5 22.0S w 0.68 0.82Maximum Fetch (m) 500 (NE) 1100 (N)3.1.2 Annual Temperature and Oxygen BudgetsIsopleths of oxygen and temperature for the period of June, 1991 to April, 1992 are givenin Figures 3.3 and 3.4. As described below, these isopleths exhibit typical characteristicsof dimictic eutrophic lakes. The maximum summer thermocline depths of about 9 m at44Aug 06 Oct 09 Dec 12 Feb 14 Apr 1822.0Jun 04HARMON LAKE0.04.48.8WrZ13.217.6MENZIES LAKEFeb 14 Apr 1^i Oct 09^Dec 12DATEFigure 3.3 Isotherms at Menzies and Harmon Lakes, June 1991 - April, 199245Oxygen (mg/1)22 0 I^^►Jun 04^Aug 06^Oct 09 Feb 14^Apr 18Dec 120.04.418.817.6Aug 06^Oct 09^Dec 12^Feb 14DATEHARMON LAKEMENZIES LAKE46Figure 3.4 Oxygen Isopleths at Menzies and Harmon Lakes, June, 1991 - April, 1992Harmon Lake and 6.5 m at Menzies Lake are evidence of greater fetch and windexposure at Harmon Lake (§ 3.1.1).The Oxygen BudgetBoth lakes exhibit the typical clinograde oxygen profiles of productive lakes: as shown inFigure 3.4, there is a rapid drop in oxygen concentrations across the thermocline, mainlyas a result of biological oxidation of organic material. This oxygen is not replacedbecause the hypolymnetic water is never exposed to the atmosphere over the course ofthe stratified period. By late June, an oxygen peak is observed above the thermocline.This is mainly due to the photosynthetic production of oxygen by phytoplankton.Phytoplankton commonly adapt to the low temperatures and dim light associated withthis depth, where they take advantage of relatively high nutrient concentrations (Wetzel,1983). Approaching the thermocline from above, there is an increase in oxygensolubility due to the drop in temperature which may contribute to this peak. The peak isfurther pronounced by the rapid decline in oxygen below the thermocline due to bacterialconsumption. The peak subsides by early September, at the end of the summer algaebloom.The oxygen budget at both lakes are similar up until a few weeks after the Menzies Lakecirculator was activated in early November, 1991. Once ice had formed later that month,oxygen depletion began near the bottom of both lakes. Diffusion of oxygen towards thelake bottom where it was being rapidly consumed resulted in significant depletionthroughout the entire water column at Harmon Lake. At Menzies Lake however, most ofthe water column was being re-circulated and replenished with oxygen by turbulententrainment and diffusion across the air-water interface within the ice-free area.47Interactions at the Mud-Water InterfaceFurther evidence of oxygen consumption is available from the measurement of dissolvedsolids. Dead plankton which has settled onto the lake bottom decompose, therebydecreasing the oxidation-reduction potential. Once the potential reaches zero at thesediment-water interface, dissolved metals are released as described by Mortimer(1971;1941;1942) and Stauffer (1987). Dissolved solids concentrations were consistentlyhigher in the hypolimnion than in the epilimnion over the summer period (suspendedsolids concentrations at both lakes were small compared with the dissolved fraction at alltimes). Over the winter, samples taken 1.0 m above the lake bottom also had a relativelyhigh concentration of dissolved solids, most notably at Menzies Lake (Figure 3.5). Thisis a noteworthy point since it provides further evidence that Menzies Lake experienced48440 —420 —400 —380 —o Top Menzies-AU BottomA-- Top HarmonBottom360 —340 —111111101111kb320 ii^IM ii300 —280 —•260Jun 11, 1991^Sep 19^Dec 28^Apr 6, 1992DateFigure 3.5 Total Dissolved Solids Near Lake Surface and Sedimentsincomplete circulation due to the Air-o-lator®: the small volume of water below 13 m(approximately 5% of the total lake volume [see Figure 3.2b]) remained stagnantthroughout the artificial circulation period. In addition to solids measurements (byevaporation of a known volume of lake water), a strong smell of hydrogen sulfide in thebottom samples indicated that reducing conditions prevailed at depth due to oxygenconsumption.3.2 DATA COLLECTION PROCEDURES AND OBSERVATIONS3.2.1 Data RequirementsIn order to utilize and verify MLI and MLI-C, a significant effort was required to acquirefield data.Input DataAs the model was designed for future compatibility with DYRESM, the inputrequirements are similar. The bulk aerodynamic equations suggested by the TVA (1972)and employed in MLI require the following meteorological data:a) air temperatureb) wind speedc) total daily solar radiationd) vapour pressuree) cloud coverf) precipitationAll of the above data are daily averages except for solar radiation and precipitation. Asdescribed in § 3.2.2, meteorological sensors were installed in order to collect all of thedata except precipitation and cloud cover. The latter variable enters into the calculation49of incoming long-wave radiation, Rli. No direct means was available of sensing Rli, butan empirical relationship exists which expresses cloud cover in terms of the observed andtheoretical clear-sky solar radiation. This formula is given along with the bulkaerodynamic equations in Chapter 4. Snowfall depth was measured by the Menzies Lakeresident using snow boards and a metre stick as described in § 3.2.2. Given that rainfallat the field location is rare over the winter months, no means for measuring rainfall wasprovided. (Less than 2 cm in total fall, on average, from December to March in theSouthern Interior Plateau. Even less would be expected at the field location due to itshigh elevation.) In accordance with Murphy's Law, however, rainfall was a significantfactor over the investigation period. Fortunately, the Menzies Lake resident noted alldays on which rainfall occurred. For these days, the rainfall data from the city of Merritt,500 m lower in elevation and about 20 km to the North, was used for snowmeltpredictions as described in Chapter 4. Although rainfall was observed at Merritt on everyday that it was noted at the lake, there was no significant correlation between theprecipitation measurements at Menzies Lake and those at Merritt (r = 0.37). Therefore,no adjustments to the Merritt rainfall data could be justified for use in the MLI model.Verification DataIn order to verify the model predictions, the following data were required:a) average temperature of the unfrozen lakeb) ice thicknessc) snow thicknessd) snow-ice thicknesse) polynya radiusAs described in § 3.2.2, sensors were installed to continuously measure the laketemperature at various depths. All of the remaining variables were only measured duringfield trips to the lakes.503.2.2 InstrumentationMeteorological data was collected using sensors connected to a multi-channel dataacquisition unit. This unit was housed in an aluminum buoy which was anchored in thelake at the location of maximum depth prior to the onset of ice cover. An anemometer,pyranometer, humidity sensor and thermistor were installed on a 2 m tower which wasfixed to the buoy (Figure 3.6). Light-gauge wire and surveyor pins frozen into the ice-cover were used to provide rigidity to the tower. In addition to this weather station, atemperature/dissolved oxygen meter was used to measure profiles of these twoparameters during field trips. All of this instrumentation is described below.DataloggerA Multidata GNOME datalogger, developed by Terrascience Systems Ltd. wasemployed. This unit scans up to fifteen channels simultaneously and logs data whileunattended. It runs under the control of a BASIC language interpreter that stores thecollected data in non-volatile random-access memory. It also has a time-out featurewhich conserves power by switching to low-power between acquisitions. An external12 V battery was used to power the unit, leaving the internal battery to act as a back-up. The scanning interval used was 0.5 hours. This interval was chosen to provideadequate data, while ensuring that enough memory would be available for field tripdelays of up to one month beyond the scheduled field trip interval. The unit was newat the time of installation.Portable ComputerA SHARP PC-6220 notebook computer was used to provide an interface for fieldcommunication with the datalogger via an RS232 serial port cable. The data acquiredby the datalogger was downloaded in hexadecimal form to the computer using51temperature sensorthermistors 12V batteryanemometerradiation shieldanchor lineIIpyranometer/•^ ^I levelingV■■■■ fixtureNilI)1humidity sensordataloggericesnowwaterFigure 3.6 Menzies Lake Weather Station52buoyinterfacing software provided by Terrascience. This procedure required 5 to 10minutes to complete. The software provided was also used to average and translatethe data to a file which could by imported to commercial spreadsheet packages.SensorsAll sensors provide reliable data under the field conditions observed at Menzies andHarmon Lakes over data collection periodTemperature Sensors: Thermistors were connected to eight of the analog channels.Two models were used (4 LM335 and 4 LM34), both producing an output voltage ofwhich temperature is a linear function. In order to ensure accurate data, the sensorswere field calibrated in the narrow temperature range characteristic of ice-coveredlakes using the Y.S.I. temperature sensor. In this manner, any non-linearity in thevoltage-temperature function would be minimized. One of the sensors was used tomeasure air temperature, and two were fixed at the same depth in the lake in order toprovide further assurance of sensor accuracy. A hole was augered into the ice and thesensors placed at 3, 7, 8, 11, 12.5 and 15 metres below the water surface. The cablessubsequently froze into the ice cover. Periodic checks following the installation of theweather station showed that the LM335/34 sensors were in agreement with the Y.S.I.sensor to within 0.1 °C. The LM34 sensors were new at the time of installation.Solar Radiation Sensor:  A LI-COR 200SA pyranometer was employed to measuredirect and diffuse solar radiation. Although this sensor does not cover the full rangeof the solar spectrum, the absolute error involved in the measurement of total solarradiation is only ±5%, and typically ±3% under daylight conditions. The sensitivityof the sensor is about 80 1.LA per 1000 W/m 2 , with a maximum linearity deviation of1% to 3000 W/m2 . A millivolt adapter (LI-COR model 2220S) was required toamplify the signal for use with the datalogger. The sensor was mounted on a LI-COR532003S leveling fixture which was bolted to a pipe that extended 0.3 m away from theanemometer post (see below) in order to avoid shading by other instruments.Stabilization of the anemometer post was required in order to ensure that thepyranometer remained level. Four lengths of light gauge wire were fixed to the postand tied back to surveyor pins frozen into the ice on the previous day. The unit wasnew at the time of installation.Wind Monitor: An RM Young model 05305 anemometer was used to measure windspeed and direction. It has a working range of 0-40 m/s, and a propeller threshold of0.4 m/s. The rotation of the propeller produces an AC sine wave signal withfrequency proportional to wind speed. According to factory calibrations, thepropeller is accurate to within 2% of the actual wind speed. The anemometer wasfixed to a post which was bolted to the buoy. The approximate height of the propellerwas approximately 2.0 m above the ice surface. The unit was new at the time ofinstallation.Relative Humidity Sensor: A VAISALA HMP 35A humidity and temperature probewas installed at a height of approximately 1.0 m above the ice surface and wasprotected from solar radiation using a shield made up of four aluminum pie plates.The HMP 35A probe contains a HUMICAP H-Sensor (model 0062). It is accurate to±1% relative humidity against factory references at 20°C. There is a temperaturedependence of ±0.04% / °C, and its long term stability is better than 1% per year.The unit was new at the time of installation.Temperature / Dissolved Oxygen ProbeTemperature and dissolved oxygen measurements were made using a Y.S.I. (YellowSprings Instruments) digital dissolved oxygen meter (model 58). The meter waschecked using a precision mercury thermometer (accurate to better than 0.1 °C) in the54laboratory. Calibration for dissolved oxygen was performed in the field by wrappingthe probe in a wet cloth and adjusting the reading to correspond to oxygen-saturatedwater at the local atmospheric pressure once the probe electrodes had polarized. Infreezing conditions, polarization required up to 20 minutes.Other equipmentPrecipitation Gauges: Two snow boards were placed in locations where windexposure was low and well away from tall obstructions such as trees and buildings. Alocal resident measured snow accumulation on the boards following each snow eventusing a metre stick. Following each measurement, the boards were cleared of snowand placed on top of the new snow. These measurements are considered accurate to ±1 cm. Neither the duration of the snowfall, nor the time of measurement relative tothe end of the snowfall event were recorded. Since much compaction can take placeover the first 24 hours following a snowfall (see Chapter 2), and since the snowboards may not have been ideally placed with respect to exposure, the errors inestimating the actual depth of snowfall is probably much greater.3.2.3 Field ProcedureField trips were scheduled every 2 to 3 weeks, depending on expected rates of change inice-cover, in addition to the availability of both equipment and assistants. The field tripschedule is given in Table 3.2. One full day (including travel time) was required tocomplete the observation procedure at Menzies and Harmon Lakes.Basic ObservationsAll observations were made near the location of maximum depth at Harmon Lake. AtMenzies Lake, additional measurements were made within the polynya and at a pointnear the opposite end of the lake. This was done to ensure that the circulation device was55never causing any significant horizontal density variations other than those caused byheat and mass transfer across the mud-water interface. Furthermore, examination of thetemperature data from within the polynya was required in order to determine if therewere any buoyancy effects in the behaviour of the radial jet.Table 32 Field Trip ScheduleDATEDec. 12/13, 1991Jan. 7, 1992Jan. 19Feb. 4Feb. 27Mar. 10Mar. 27Apr. 7meteorological sensors installedOne hole was augered into the ice at each station. Ice depths were measured in thefollowing manner: a metal bar was submerged into the hole, and pulled back up alongthe edge of the hole until a metal plate welded to its end made contact with the ice. Thedistance between the top of the plate and the point on the bar corresponding to the top ofthe ice sheet was then measured. Flooding due to excessive snowfall (see § 2.3.3) wasobserved on February 4th, 1992. Snow-ice thickness was only measured for the first timeon the following field trip, on February 27th. Total ice measurements are consideredaccurate to ±1 cm, whereas snow-ice measurements are not considered reliable for thereasons outlined in § 3.2.4.The Y.S.I probe was lowered through each hole, and measurements of both oxygen andtemperature were made at 1 m intervals below the water surface.Water samples were taken using a Van Doom sampling bottle at 2.0 m below the surfaceand 1.0 m off the lake bottom. These samples were analyzed in the laboratory for totalsolids using the standard evaporation technique. The mixed layer sample was taken at 256m in order to avoid the ice-water boundary layer stratification in which solidsconcentrations may be relatively high due to salt rejection. The deep sample was assumedto represent total solids in the stagnant layer, and was taken only 1 m above thesediments given the great depth of mixing induced by the circulator.Additional Observations Made at Menzies LakeAir temperature and relative humidity measurements were taken using a mercurythermometer and sling psychrometer. The psychrometer readings were only consideredreliable when the air temperature was above the freezing mark. These measurementswere compared with instantaneous readings from the humidity sensor on the weatherstation. The weather station thermistors were also checked irregularly to verify thatreadings were consistent with the Y.S.I. meter.On each trip, the data collected by the weather station was downloaded to the portablecomputer. The datalogger was subsequently cleared of memory and re-initialized toimmediately begin data collection.Oxygen and temperature profiles were also observed from a boat within the polynya, asclose to the circulator as possible. In addition, a surveyor's chain was used tosupplement visual estimates of the polynya radius. The end of the chain was tied to thecirculator. The chain was then permitted to unspool as the boat was moved to thepolynya perimeter, where measurements were make in 2 to 10 directions, depending onthe time available.Radial Jet Velocity MeasurementsIn order to quantify the sensible heat transfer to the polynya edge, it was necessary toformulate an equation describing the velocity field induced by the radial jet of waterproduced by the circulator. This required field measurements of velocity as a function of57radial distance away from the circulator. This was carried out using a propeller-drivenOTT® meter. Measurements were taken from the bow of a boat set inside the polynyasuch that the meter was always located upstream of the boat. This was done so as toavoid erroneous measurements due to disruption of the flow field. It was found that onlysurface measurements provided data which were reproducible and useful in terms ofproducing a reliable velocity field formula.A surveyor's chain was fixed to the circulator at one end, and a tree on shore at the other.This produced no disruption to either the operation of the circulator, nor the inducedvelocity field. The boat was pulled along the chain and three measurements were takenat each station. The sampling interval was as short as 10 cm close to the circulator,increasing with radial distance to a maximum of 2 metres. At about ten metres from thecirculator, the surface velocities were below the threshold of the current meter.In addition, current velocities were estimated by timing drogues placed in the water at theraft, and allowed to float toward the edge of the ice-free zone. Although much lessaccurate, the velocities calculated from the drogue data were of the same order as thedirect measurements. Wind speed was less than 0.5 m/s during all current measurements,and was therefore ignored.3.2.4 Field Trip ObservationsThe Menzies Lake PolynyaIt was found that the polynya radius varied only slightly over the entire winter period (seeChapter 5). Measurements were not of adequate quality or frequency to detect significantchanges in size until spring when the ice cover retreated rapidly to produce completelyice-free conditions within 2 weeks. The radius remained at about 20 m on average. Theshape of the polynya was elliptical rather than circular, with the upwind end tending to be58several metres shorter, and the downwind end several metres longer. When MenziesLake was visited in the early morning in cold weather, there tended to be a very thin icelayer, about 1 m wide, around the perimeter of the polynya. This was the only clearevidence of polynya contraction due to cold weather. Frazil ice was never observed inthe polynya because the average lake temperature never dropped below about 1.4 °C.Since isothermal conditions (to an accuracy of +/- 0.1 deg C) were measured within thepolynya over at least the top 3 m (the depth of the circulator intake) , buoyancy wasdismissed as an important factor in the circulation pattern and heat transfer.Lake TemperatureIsotherms based on the Y.S.I. probe data, for the period spanning December 13th andMarch 10th are given in Figure 3.7. The latter date is that of the last field trip duringwhich full ice-cover was observed at Hannon Lake. Clearly, the impact of the circulationdevice on the heat content of Menzies Lake is substantial. While the water below ice-cover at Harmon Lake gradually but consistently heats up over the winter due to solarradiation and sediment heat transfer, the circulator at Menzies Lake cause the water tocool down at a relatively high rate. As early as late February, however, average airtemperatures almost consistently above zero and increased solar heating rapidly drive thetemperature of both lakes up (Figure 3.8). By early April, both lakes were completely icefree.On March 10th, a reproducible unstable temperature profile was observed at HarmonLake. Throughout most of the water column on this date, the temperature was near thepoint of maximum density. The vertical variation in solids concentration was insufficientto explain the instability. It is hypothesized that the lake was undergoing vigorousconvective circulation. as solar radiation penetrated deep into the lake (the Secchi depthtransparency on this date was 11 m, the greatest ever observed at Harmon Lake). At 4 0C,59HARMON LAKEMENZIES LAKE0.0603.36.6ar49.913.216.5Dec 13 Jan 04 Jan 26 Feb 17 Mar 10DATEFigure 3.7 Isotherms at Menzies and Harmon Lakes, December 13th - March 10th61Feb 8^Feb 13^Feb 18 Feb 23^Feb 28^Mar 4 Mar 9DateFigure 3.8 Air Temperature and Solar RadiationObserved'Feb 8^Feb 13^Feb 18^Feb 23^Feb 28^Mar 4^Mar 9Date640E-Ua)cdg -2-4 —0.)-6 —-t>t.14x103 —12 —10 —286 —4 —c.c. = cloud cover (fraction of sky)1clear-sky radiationc.c=0.1c.c. = 0.1any local changes in temperature (except at the surface) results in instabilities whichcause convection. A re-circulation pattern may have developed as water heated by solarradiation rose up through the water column, only to be cooled near the thermal boundarylayer at the ice-water interface. This water would then begin to sink again as itapproached the temperature of maximum density.Another potentially important factor in the heat budget during the melting period is there-radiation of solar heat from the lake bottom in the littoral zone, which may haveproduced greater melting of the ice-cover in the near-shore area. In addition, the heattransfer assumed to be active at the ice-water interface may have increased due to anincreased thermal gradient at the boundary and/or the development of much higher ratesof turbulent heat transfer.Comparison of the temperature profiles at the three Menzies Lake stations has revealedthat there are no significant horizontal variations at any time over the winter period otherthan those caused by heat and mass transfer across the sediment-water interface (see§3.2.3).Ice and Snow ThicknessIce and snow thicknesses over the observation period are shown in Figure 3.9. Snow-icethicknesses may not be reliable for several reasons. On January 19th, a floodedsnowpack was observed, presumably as a result of a recent snowfall. (The most recentoccurred the day before, according to the records provided by the Menzies Lakeresident.) The thickness of flooded snow tended to increase towards the centre of thelake, with no observed flooding at the lake edge. The response of the ice sheet could becompared to that of a simply supported beam bending under a distributed load. Thedepth of flooding was observed, but snow-ice thickness was not, since the phenomenonwas not previously evident as an important factor in the heat transfer problem. The6230-Dec-91 19-Jan-92 19-Mar-92 8-Apr-920.45 —0.4 —0.35 —0.3 —0.25 —0.2 —0.15 —0.1 —0.05 —0 ^10-Dec-91—6— Menzies (stn 1)—o— Menzies (stn 2)---0-- Harmon8-Feb-92 28-Feb-92Date--A-- Menzies (stn 1)—G— Menzies (stn 2)—0— HarmonIce Thicknesses at Sampling Stations63Snow Thicknesses at Sampling Stations0.140.12E 0.1v,8 0.08...s4..c4g 0.063o= 0.04CD0.02010-Dec-91 30-Dec-91 19-Jan-92 8-Feb-92 28-Feb-92 19-Mar-92 8-Apr-92DateFigure 3.9 Ice and Snow Thicknessesactual total thickness of snow-ice was only observed for the first time on February 27th.Furthermore, there was a subjective transition from snow-ice to pure ice which may havebeen interpreted differently during the proceeding field trips. Since the modeling periodends on the third trip which included snow-ice observations (March 27th : see Chapter 4),no snow-ice observations will be used for model verification purposes.The Melting PeriodBetween the field trips of February 27th and March 10th, heating of both Menzies andHarmon Lakes increased significantly (Figure 3.7). Observations near the locations ofmaximum depth, however, indicate an insignificant decrease in total ice thickness atHarmon Lake, and even an increase of similar magnitude at Menzies Lake (Figure 3.9).Just over two weeks later, on March 27th, there was an ice-free zone 1 to 4 m widearound the main basin of Harmon Lake. This zone was substantially wider in the shallowbasin at the south end of the lake. The ice at the sampling station, however, was 13 cmthick. At Menzies Lake, the polynya had reached the east end of the lake, and care wasrequired to reach the weather station: it was clear that the maximum ice thickness in thelittoral zone was barely sufficient to support the weight of an adult (during field trips inthe late fall of 1991, it was found that an adult could be supported by as little as 5 cm ofice). However, the ice was 13 cm thick at Station 2, and 23 cm thick at Station 3. Thiswas by far the most substantial spatial variation in ice thickness observed at either lake.Judging by the changes in water column temperature, and the observed ice-free zone onMarch 27th, it is suspected that substantial thinning of the littoral ice would haveoccurred between February 27th and March 10th.64Boundary Effects at Menzies LakeA closer inspection of the Menzies Lake isotherms as generated by the weather stationsensors (Figure 3.10), indicates that the circulator's zone of influence extended down tothe lake bottom by about December 18th, 1991. Without regard to morphometry, itwould appear that a sudden increase in the rate of cooling occurred between December15th and 17th. There is, however, neither consistent information in the meteorologicaldata, nor is there conclusive evidence that a concurrent decline in the whole lake heatcontent occurred. The upper and lower heat content bounds using the weather stationthermistor data and applying density considerations are shown in Figure 3.11. Theinformation contained in this figure is insufficient to be conclusive, but when combinedwith the meteorological data, there is strong evidence that the sudden drop in theisotherms between 2.0 and 2.8 0C can be explained by examining the dynamics ofcirculation within the lake.Laboratory investigations have been performed in which a fluid was injected at a constantrate on top of another fluid of greater density. This lower fluid being withdrawn at thesame rate (see Wood, 1978; Jirka & Katavola, 1980). The interface which developsdescends through the water column and is drawn down at some critical point above thewithdrawal depth. The interface then continues to drop until it reaches a second criticalpoint below the intake, where the denser fluid will no longer be withdrawn. It ishypothesized that the equivalent interface at Menzies Lake dropped from about 8 to 13 mbetween December 15th and 17th. Hence the rate at which the isotherms shown in Figure3.10 drop is a function of pumping rate and lake morphometry, rather than the rate ofcooling.Given that the lake depth is approximately 13 m directly below the circulator and the lakeis almost isothermal above this depth by the 18th of December, it is concluded that the65YSI PROBE DATADATALOGGER THERM'S TORS13.216.5Dec 130.0Feb 17^Mar 10I^I^i^1^JDec 16^Dec 20^Dec 23^Dec 27^Dec 31DATE3.166Figure 3.10 Isotherms Generated from Datalogger Data at Menzies Lakemaximum withdrawal depth of the circulator is limited by the lake bottom. A thermalboundary layer caused by the diffusion of heat from the sediments would account for theobservation that purely isothermal conditions only existed above about 12 m after thisdate.Figure 3.11 Maximum and Minimum Likely Heat Content at Menzies Lake(based on density considerations [temperature effects only])Cloud CoverCloud cover was also observed from time to time during field trips. Estimates of dailyaverage cloud cover, expressed in tenths of the celestial hemisphere, are given in Table3.3. These observations will be compared with the estimates made in the model asdescribed in the following sub-section.67Weather stationIn the few weeks before and after the winter solstice, the tall trees on the south side ofMenzies Lake cast a shadow over the pyranometer all day except for 3 to 4 hours startingin the late morning. Although the solar radiation data during this period, may beappropriate in terms of the heat transfer between the lake and the atmosphere, twoproblems are created by the shading. First, as described below, comparisons cannot bemade with the theoretical clear-sky radiation. Second, the empirical formula expressingcloud cover in terms of the ratio of observed to clear-sky radiation cannot be consideredaccurate for use in the model (see Chapter 4). These problems apply from the start of theobservation period (December 13th) to some time between January 7th and January 19th.No shading of the pyranometer was observed on the latter date.Table 3.3 Observed Average Daily Cloud CoverDate I^Cloud coverDecember 13 5January 7 1January 19 4February 4 4February 27 1March 10 1March 27 9April 7 3During a cold spell in the data collection period, the glass which encases the pyranometerleveling bubble cracked, and the alcohol inside evaporated. As a result, accurate levelingof the pyranometer could not be ensured for a significant portion of the data collectionperiod.As described in §3.2.3, the meteorological tower was anchored to the ice cover usinglight gauge wire and four surveyor's pins frozen into the ice. As melting of the top of theice cover proceeded, however, the pins loosened, but at different rates depending on their68location and the depth and angle of penetration. In addition, melting of the ice around thecircumference of the buoy occurred due to solar heating. As result of these two meltingprocesses, the entire weather station became tilted to the west. This occurred at sometime between February 27th and March 10th. The angle of tilt was not measured, but isestimated at 100. Therefore, the deviation of the pyranometer from the horizontal mayhave been 100 or more for up to two weeks.The impact of this tilt may have been considerable. From the cloud cover observations,the solar radiation should have been nearly equal to clear-sky levels on February 27th andMarch 10th (the theoretical clear-sky radiation is shown in Figure 3.8, and was computedas given in Appendix A. It should be noted that the clear-sky radiation does not increasemonotonically after the winter solstice as a result of the influence of albedo on theindirect portion of the total incoming radiation.). On March 10th, however, the ratio ofobserved to clear-sky radiation is comparable to that on February 4th, when the observedaverage cloud cover was 0.4. No other observations of clear sky days were made duringthe period when there was insignificant topographic obstruction (tree shading) of theincoming solar radiation. Confidence in the clear-sky radiation calculations is limited bythe lack of reliable comparisons between the observed and calculated cloud cover overthe ice-covered period.Anthropogenic EffectsHarmon Lake is a popular ice-fishing lake over the winter period. In the early spring,photosynthetic activity resumes under the ice-cover as a result of increased levels of solarradiation and reduced snow-cover. In response, the resident fish move up the watercolumn to feed in the photic zone. Their presence near the lake surface makes the fishmore susceptible to angling (Ashley, pers. comm., 1992). The anglers use snowmobilesextensively at this time of year, resulting in a significant impact to the ice surface. Snow69and snow-ice were melted in the snowmobile tracks, which were estimated to coverbetween one third and one half of the lake surface area. In some areas, ponding wasobserved in the depressions caused by these vehicles. The important consideration to bemade here is the impact of this melting and ponding on the surface albedo.The Radial JetThe circulator is supported by a square raft and the jet strikes the raft at the corners (referto Figure 1.1). The jet, however, extends beyond the raft, entering the water directly,along the sides. This results in the current being strongest in the NNW,WSW,ENE andSSE directions. Measurements were made in these directions, and any dependencies ofvelocity on the direction of flow have been ignored.The surface velocities measured using the current meter (Figure 3.12) agree very wellwith turbulent radial jet theory (see §2.4.4). Measurements at 10, 15 and 20 cm belowthe water surface are, however, only qualitatively consistent with theory. The maximumdepth at which velocities were large enough to measure increased with radial distancefrom the circulator. This is consistent with the concept that the jet thickness, orentrainment zone should thicken with radial distance. Furthermore, the velocities at agiven depth reach a peak then decay, as would be expected as the jet continues to thickenwith increasing radius. The vertical velocity gradient reduces rapidly and at a radialdistance of 7 m, the velocity is nearly uniform over the top 20 cm.The motion of the drogues was observed to be qualitatively similar to the velocityrelationship shown in Figure 3.12. The average speed in the first 3 to 5 metres wasestimated to be in the order of 0.25 m/s. These estimates are consistent with the directmeasurements.7071Figure 3.12 Velocity Measurements in Vicinity of CirculatorChapter 4NUMERICAL MODELIn this chapter, the development of MLI and MLI-C is given, along with a discussion onthe choice of the many parameters required to run the model.The ice cover component of the model developed by Patterson & Hamblin (1988) hasbeen extended and adapted to a bulk representation of the lake heat budget (in otherwords, it is assumed that the small lakes considered are of uniform temperature under ice-cover). The formulation, as described below, is based on the same basic assumptions andequations as those presented in Chapter 2. Improvements to the model include theaddition of the following features:a) a snow-ice component in the heat conduction equationb) a heat source distributed over the top component of the cover to account for theheat associated with snow-ice formationc) an extra term in the surface heat budget to account for rainfalld) sediment heat transfere) a variable albedo.0 a variable snow density and conductivityIn addition, routines were added to account for the impact of artificial circulation in theMLI-C model. These routines feature the following algorithms:a) an estimate of polynya size based on the balance struck by the meteorologicalconditions and turbulent heat transfer as caused by the radial jetb) an estimate of the net turbulent heat transfer between the water and the ice-coverc) free convection between the atmosphere and the water in the ice-free aread) cooling due to snowfall on the open water.72Also considered in this chapter are the following issues and parameters:a) thermal characteristics of the lake sediments,b) the turbulent heat transfer coefficient, CLc) the area over which turbulent heat transfer to the ice cover is influenced by theradial jet,d) the thermal gradient at the ice-water interface,e) the attenuation of solar radiation through the ice and snow cover4.1 MLI FORMULATIONPatterson & Hamblin (1988) applied a two-band radiation absorption law to the steadystate conduction equations for an ice and snow cover. These equations can be extendedto include snow-ice as a third component in the cover:a2TsKs a z2 + A 1XslicreXP{ -XSi(hi+he+hs -z)] + A2Xs2IoexP[ -Xs2(hi+he+hs-z)] + Qsi = 0,hi+he+hs z hi+hea2TeKe azz + A1XelloexP[-Xs1hs-Xel(hi+he-z)] + A2Xe21oexP[As2hs -Xe2(hi+he-z)]=0,hi+he z ?_ hia2TiKi az2 + AIXilloexp[Xsihs -Xelhe-Xil(hi-z)] + A2A42I0exp[-Xs2hs-Xe2he-Xi2(hi -z)]=0,hi z 0,^ (4.1)where the subscript e refers to the snow-ice layer and all other terms are as previously73defined.The heat source, Qsi, accounts for the heat supplied by lake water flooding the ice coverfollowing a heavy snowfall (see §4.3.2). This source of heat, which is assumed to bedistributed over the entire snow thickness, includes the sensible heat given up to the snowlayer as it cools to the freezing mark, and the latent heat which is produced when thiswater freezes to form snow-ice. Although it would be more appropriate to add a fourthmedium to the lake cover to reflect that this heat is not distributed over the entire snowthickness, the data available are not of a sufficient quality to justify further complicatingthe governing equations. This process is described in greater detail in Chapter 2.The appropriate boundary conditions for the above equations are:74Ti = TfTi = TeK. aTi . Ke aTe^K , az^ azTe = TsaTe _Ke— = ^az^azTs:--- Toz = 0z = hiz = hi + hez= hi + he + h s^(4.2)Again, qo is defined as the conductive heat flux between the ice-cover and theatmosphere:aTsC10 = - .1T, 5 az z= hi + he + hs (4.3) The solution can be written as:Chl(1.2411e_F hsKe K s ) (c10 -10) = Tf - To Ago { 1-expeXsihs) exp(-Xsihsill-exP(-Xeihe)]KsXsi^KeXeiexp(-Xsihs-Xeihe)[1-exp(-ailhi)l }KiXii_A2101 ^ics1-exp(21:22 s2hs) exp(-Xs2hs)[1 -exP(-Xe2he)]KeXe2exp(-Xs2hs-Xe2he)[1-exp(-Xi2hi)]lKiXi2+Qsi hs^+^hs Qsi hs2Ki Ke KO - 2 Ks(4.4)The extension to any number of components is obvious. The surface condition may beexpressed as given by Equation (2.9).Ablation and accretion of ice at the ice-water interface is modeled as given in Chapter 2.The heat flux, qf, is again obtained from the solution of the heat conduction equation:qf = -Kiaz= go - Al lo( 1 - exp-(?sihs+Xeihe+Xi1))- A2I0 ( 1 - exp-(Xsihs -I-Xeihe+Xii)) - Qsi hs^ (4.5)4.2 BULK AERODYNAMIC FORMULAEThe heat balance at the interface between the ice cover and the atmosphere is dependenton meteorological conditions. The components of the balance may be computed usingbulk aerodynamic formulae such as those given in TVA (1972), Imberger & Patterson75aTiz=076(1981) and the DYRESM User's Manual. These formulae are given here. Adjustment ofthe meteorological data for use in the bulk aerodynamic formulae is discussed in §4.5.The net incoming meteorological flux is given as:H = Rh - Rio(To) + Qc(ro) + Qe(To) + Qs(To)^ (4.6)The heat supplied by rainfall, Q r, includes sensible heat as described in §2.4.4 and latentheat as described in §4.3.1. All other terms are as defined in §2.1.3. It is important tonote that the heat generated both by solar radiation and ice cover flooding are notincluded in the balance since they appear in Equation 4.4.Long-wave radiation is governed by black body radiation principles, and is thereforedependent on the temperature and emissivity of the body in question. The incomingradiation, R11, depends on the characteristics of the sky and is thus more difficult todetermine than the outgoing radiation from the lake surface, R1 0 . Rh may be estimatedusing the Swinbank (Equation 3.7) or Anderson (Equation 3.8) relationships (TVA,1972):Rh = 9.37 x 10 -6 6 (Tair + 273)6 (1+0.17C2)^ (4.7)Rh = (0.74+0.0049 svpd) 6 (Tair + 273)4 (1+0.17C2) (4.8)where, RH^= incoming long-wave radiationC^= cloud cover (fraction of sky)svpd = vapour pressure at Taira^= Stefan - Boltzmann ConstantTair^= air temperatureHamblin (1990), in a study involving a meteorological data set for the Yukon Riverbasin, found that the former relationship was more accurate when the air temperature,Tair , is above freezing, while the latter is superior below freezing.The vapour pressure of the air is calculated from the relative humidity (RH) data and theair temperature (TVA, 1972):airsvpd = -RHi-F0 exp (2.303 ( TaaT.---r_ f-Fb c)) (4.9)where the constants a,b,c depend on whether open water or ice-covered conditionsprevail.The cloud cover, C, was not directly observed on a daily basis but may be estimatedusing the solar radiation data. The following relationship is one of three recommendedby the TVA (1972) and was found to provide reasonable estimates in Hamblin's (1990)study:.\1 (1 - SW C =^SWCS )0.65 (4.10)where, SW^= observed total daily solar radiationSWCS = theoretical total daily solar radiation under clear sky conditions.The theoretical clear sky total includes indirect radiation and depends on the Julian day,latitude, atmospheric attenuation and the local albedo. Details of its calculation are givenin Appendix A.The outgoing radiation, Rh, is given as:Rio = C a ( To + 273.15 )4^(4.11)where, E = emissivity of the lake surface = 0.97The sensible heat flux between the lake surface and the atmosphere is governed by forcedconvection heat transfer principles, and is thus a function of the temperature differenceand the wind speed:7778Qc = Pair Cp Ch Uw ( To - Tair)^ (4.12)where, pair= density of airCp = heat capacity of waterCh = dimensionless empirical coefficientUw = wind speedFor wind speed measured at 10 m above the lake surface, Imberger & Patterson (1981)recommend Ch = 1.5 x 10-3 . Equation 4.12 should be multiplied by the shelteringcoefficient described in Chapter 3.The evaporative heat flux, Qe , is dependent on the difference between the vapourpressure of the air, svpd, and the saturated pressure, svp0, at the temperature of the lakesurface, To . In a vacuum, vaporization from the surface would continue until the vacuumreaches saturation pressure. Under atmospheric conditions, however, a vapour blanket atthe surface approaches saturation and turbulently diffuses into the air above. The rate ofevaporation is therefore also dependent on wind speed:Qe = Ce Uw ( svp0 - svpd)^ (4.13)where, Ce = a dimensionless empirical coefficient. The saturation vapour pressure, svp0,is given by Equation 4.9, with RH = 100% and Tair replaced by To, the lake surfacetemperature. Equation (4.13) also applies in the case of condensation when svpd > svp0.The evaporative flux should also be reduced by the wind sheltering coefficient.The only variable appearing in the above formulae which was not measured directly orestimated using techniques found in the literature, is To, which is provided in the solutionto Equation MLI: LAKE IN THE NATURAL STATE4.3.1 MLI Main RoutineThe MLI flow chart is given in Figure 4.1. Once the initial conditions and parametershave been set, the daily loop commences. Within this loop, the daily averages of themeasured meteorological variables are read into memory. If snow is present on the icecover, control is passed to SNOWDENS, the subroutine in which the snow density andconductivity are calculated, and then to SNOWICE where it is determined if flooding ofthe ice cover occurs.Subroutine SOLAR is then called, where both SWCS and C are estimated. Since SWCSdepends on the local albedo, the ALBEDO subroutine is called from SOLAR.Control is returned to the main programme where the sub-daily time step begins. Sincethe internal hydraulics, which can cause rapid changes to thermal structure are ignored,the short sub-daily time steps used in DYRESM are not required. Without resorting tothe requirement of sub-daily averages of the meteorological data, a sub-daily time step isonly considered appropriate to reflect that solar heating is not distributed over the entireday. The time over which solar heating is active will vary over the ice-covered periodbut, for simplicity, two 12 hour time steps have been incorporated into MLI, with all ofthe observed solar radiation applied to the first sub-daily step.Within the sub-daily time step, the surface temperature, To , and the heat fluxes at thesurface-atmosphere interface are calculated. Newton's iterative method is used to find T ofrom Equation 4.4, with Tair used as the initial estimate of To . If it is found that T o isgreater than Tm the melting temperature, it is set equal to T m and subroutine MELT iscalled to quantify the amount of surface melt. In this case qo , the flux in the surface7980Figure 4.1 Mill Flow Chartcomponent at the atmospheric interface, is determined from Equation 4.4, with To = Tm .If To is less than Tm , then qo is set equal to the net meteorological flux, H.The heat source, Qs;, is required in the calculation of To described above. It is assumedthat the heat is added to the cover at a uniform rate over the entire day. Q si is calculatedin Subroutine SNOWICE (§4.3.2).It has been found that snowmelt due to rainfall is only important compared withcondensation melt when total daily rainfall is greater than about 13 cm/day (Harr, 1981).Rainfalls of this magnitude were not observed at Merritt at any time over the winter of1991-92 (rainfall was observed but not measured at Menzies Lake). However, the rainon snow events of late January may have been significant in increasing the snow-covertemperature from sub-freezing temperatures due to the latent heat of fusion, Qri,produced upon freezing of the rainwater:Qii = L pw p , (4.14)where all terms are as previously defined. This heat will increase the temperature of thesnow pack and may result in melting if the surface temperature reaches 0 °C. Both thesensible (Equation 2.22) and latent heat (Equation 4.14) associated with rainfall areincluded in the surface meteorological balance as given by Equation 4.6. It would likelybe more realistic to treat the latent heat as an additional internal source of heat in thesnow cover. The use of daily meteorological averages, however, results in under-estimation of both condensation melt and sensible heat transfer from the rainfall to thesnow. Ideally, the average humidity and air temperature associated with the rainfallperiod would be used in calculating Qrs and Qe , not the average daily values. Withoutthe latent heat term in the surface balance, the model would incorrectly predict very littlesnowmelt. This issue is discussed further in Chapter 5.81From the main programme, subroutine ICEWATER is invoked, where the fluxes of andqw are calculated using Equations 4.5 and 2.12 respectively. These fluxes are used inICEWATER to determine the change in ice thickness at the ice-water interface as givenby Equation 2.11. Finally the net heat transfer to the unfrozen lake is computed insubroutine TEMP. The change in lake heat content is determined by the balance ofresidual solar energy, sediment heat transfer, qsed, and the conductive heat flux fromthe water to the ice, qw . The formulation of these variables, as applied to MLI, isdescribed in §4.3.2 and §4.3.3. Once the new lake temperature is computed, the timestep is incremented. If the new time step is the second of the day, then the procedure isrepeated starting from the calculation of To and the net meteorological flux, this timewith the solar radiation set to zero in Equation MLI SubroutinesSNOWDENS The SNOWDENS flowchart is given in Figure 4.2 First, the number of days since thelast snowfall, is incremented. If no snow has fallen that day, then densification of thesnowpack will proceed as follows:If there are two layers of snow (that is, if there have been two or more snowfalls over thesimulation period), then the lower layer is assumed to be at a pre-set maximum. Thismaximum density is increased to a second pre-set value if the surface temperature T o isgreater than 00C. A maximum density which is higher still is assigned if rainfall occurs.Once any accumulated snow has been assigned to a higher maximum density, only newlyfallen snow may be characterized by a lesser maximum. The upper layer of snow isassumed to settle according to an exponential decay function:8283Figure 42 SNOWDENS FLOW CHARTPs = Pnwhere,(Pm - Pn)k^- In(1 - e - kd)iPm - P1 1(4.15)Pm - PnPn = density of fresh snowP1 = density of snow 24 hours after fallingpm = maximum density of snowd = number of days since the last snowfallIf there is one layer of snow, it is assumed to settle in accordance with the aboveequation. Although somewhat arbitrary in form, this function will predict densitieswithin the appropriate ranges given in the literature (see §2.3.1), provided that reasonablevalues for pn , p1 and pm are chosen.In the case where snowfall occurs, then all of the underlying snow is assumed to be atPm, the fresh snow is set to pn , and d is reset to 1.The empirical equation given by Grant & Rhea (1973) was specifically considered inchoosing an initial density, pn , of 80 kg/m3 because it is based on observations at analtitude which is comparable to that of Menzies and Harmon Lakes. The density after 24hours, pi, was set to 100 kg/m 3 , and the maximum density, p m, was set to 200 kg/m3 forTo < 0°C, or 300 kg/m3 for To > 0°C. These values may be altered somewhat in thetesting of MLI in order to produce results which conform with the observations.It should be noted that the ice cover observed on December 13th, 1991 at Menzies Lakewould have been insufficient to support the observed fresh snow cover if the snowdensity had been greater than about 90 kg/m 3 . Hence the assumed fresh snow density of80 kg/m3 is unlikely to be too low. On January 19th, for the observed 7.5 cm of snowwhich was saturated with flooded water, the snow density should have been in the orderof 250 kg/m3 . Although fresh snow had fallen on the previous day to cause flooding, at84least two-thirds of the snow cover was old snow. According to the model describedabove, the average density of the entire snowpack would have been in the order of 230kg/m3 which is in good agreement with the first estimate (the maximum density of theunderlying snow layer would have been at the highest maximum density on this date dueto the rain which fell about a week earlier).The final steps in SNOWDENS include the calculation of the average density of theentire snow thickness, which is used in Equation 2.16 to determine the conductivity.SNOWICE The SNOWICE flowchart is given in Figure 4.3. Equation 2.21, which provides anestimate of snow-ice production, must be modified for the case where a snow-ice layeralready exits in the cover:hsm — hi (p w - Pi) + he Ow - Pe) (4.16)Pswhere the subscript e refers to the snow-ice.The density of snow-ice is generally in the range of 880 to 900 kg/m 3 . Lower densitiesresult from finer snow grains (Ashton, 1986). A constant value of 890 kg/m 3 is assumedhere.If hsm is less than the total snow thickness, h s , including the most recent snowfall, thenflooding will occur. The calculations are complicated by the possible existence of twosnow layers of different densities. First, Equation 4.16 may be applied using the averagesnow density to determine if flooding will occur. If so, then the depth of flooding may becalculated in one of several ways depending on the thickness of each snow layer. If twosnow layers exist then Equation 4.16 is modified so that the depth of flooding iscalculated based on the existence of 4 mediums. Otherwise Equation 4.16 is applied as85Calculate hen , he,^Pset 'Is = hsmshown to determine the maximum snow thickness which can be supported by the ice.The new snow-ice produced, hen , is equal to the difference between hs and hsm , and isadded to the total snow-ice thickness, h e . Qsi, the sensible and latent heat added to thesnow layer by virtue of the flooding, is calculated as follows:Qsi = (Tw cpw + L) pw S^pw^ (4.17)where, Qsi = heat generated per unit volume of snowTw = water temperatureL = latent heat of fusionFigure 4.3 SNOWICE FLOW CHARTThe factor hen/hs appears in Equation 4.16 to account for the fact that, although it isassumed that Q si is uniformly distributed over the new value of h s, in reality it is onlydistributed over hen . The (1-ps/pw) term is added in recognition of the fact that the water86which seeps into the snow layer must be distributed throughout the snow pore structure(i.e. the thickness h en is composed of both snow and water, but only the water constitutesthe heat source). Q si is assumed to be distributed over the entire day.SOLARThe SOLAR flowchart is given in Figure 4.4. This subroutine begins with the calculationof, %, the total solar radiation at the top of the atmosphere for the given day andatmospheric attenuation. The ALBEDO subroutine is called at this stage, as it isimportant to the total amount of reflected or indirect solar radiation which could bemeasured near the lake surface. Hence it is required for the calculation of the theoreticalshort wave radiation under clear sky conditions, SWCS. The cloud cover, C, is thenestimated using Equation 4.10. Since albedo is a function of C however, reiteration ofthe above steps starting from the call to ALBEDO is required until convergence of C isachieved. Once the first iteration is complete, then a flagging variable is set equal to 1 sothat only the cloud dependent calculations are performed in Subroutine ALBEDO. Thisflag is reset to zero after convergence is achieved. The details involved in the calculationof SWCS is given in Appendix A.ALBEDOThe ALBEDO flowchart is given in Figure 4.5. In the case where there is no snow layer,the albedo, a, is determined as given in Table 4.1Table 4.1 Snow-Ice and Pure Ice Albedo (Henderson-Sellars, 1984)Surface Condition apure ice 0.02snow-ice 0.4587Figure 4.4 SOLAR FLOW CHARTGiven the lack of grain size and snow quality data, it would appear that the USACEfunctions and Petzold's (1977) equations are the only quantitative means of altering thesnow albedo for the Menzies and Harmon Lake data set. Tabulated lake albedoinformation (see §2.3.2) were referenced, in order to determine the applicability of themodel.If new snow has fallen, then d, the number of days since the last snowfall, is set to zero,and the base albedo, ab, is set to 0.84 as is customary when using the USACE albedodecay functions. This base value is modified in accordance with the sky conditions toproduce a. If no snow has fallen, ab is modified using the aging functions (Equations8889Figure 4.5 ALBEDO FLOW CHART2.16) prior to the sky condition adjustments. In the case of rain, 0.05 is subtracted fromab provided that it is not already less than 0.8. These latter adjustments, to account forsolar angle and cloud cover are given by Equations 2.19 and 2.20.Given that Menzies and Harmon Lakes are both remote and well over 1000 m above sealevel, snow purity should be high. The solar angle is also low over much of thesimulation period. It would therefore be reasonable to suspect that the albedo was oftenin excess of 0.9 following snow accumulation events. Although the reduction in albedoin the case of rain is arbitrary, albedos which are consistent with this expectation areproduced. In addition, results for the full spectrum of observed conditions are within therange given by Henderson-Sellars (1984). Furthermore, a significant reduction in albedofollowing rain events is justifiable given that an increase in average grain size will resultdue to fine flake destruction and metamorphosis to granular snow (see Chapter 2). BothPetzold (1977) and Strickland (1982) attributed lower than predicted albedo observationsto rainfall.MELTThe flowchart for MELT is given in Figure 4.6. Since surface melting is taking place,the heat fluxes at the surface, H and q 0 , are calculated using T o = Tm . All meltingproceeds according to Equation 2.9, with the density set to the value which ischaracteristic of the surface component whether it be snow, snow-ice or ice. Heatderived from rainfall is included in the surface balance as described in §4.3.1. If theentire top component is melted then the excess heat is used to melt the layer(s) below.ICEWATERThe flowchart for ICEWATER is given in Figure 4.7. The fluxes, of and qw , at the ice-water interface are calculated using Equations 4.5 and 2.12 respectively. The melt or90hS = 0 'Calculatenew heNNhe = 0Calculatenew hi 1^ MELTo = T4.Calculate H, goYCalculatenew heYCalculatenew hsCalculatenew hi,hs = he = 0RETURNFigure 4.6 MELT FLOW CHART91Calculatehi>I Calculatenew heyCalculatenew heV >^(, RETURN Figure 4.7 ICEWATER FLOW CHART92accretion of ice proceeds in a manner similar to that of subroutine MELT using Equation2.11. Although melting of either pure or snow-ice may occur, only pure ice may beformed at the ice-water interface. Turbulent heat transfer is assumed to be of importanceonly in the case of artificial circulation since there is no significant flow-through at thelakes considered. This process is treated in §4.4.In order to calculate qw, information regarding the thermal boundary layer at the interfaceis required. Detailed measurements are not available, but a maximum layer thicknessmay be estimated from the field data. The shape and thickness of this layer is assumed tobe constant over the entire ice-covered period for both Menzies and Harmon Lakes.Temperature measurements made 1 m below the water surface at both lakes (usuallyabout 0.8 to 0.9 m below the ice-water interface, although this distance was notmeasured) indicated that the boundary layer did not extend to this depth. On March 10th,measurements were made at 0.5 m. At Menzies Lake there was an indication that theboundary layer was just intersected, but the data are inconclusive since the temperature atthis depth was only 0.1 0C less than at 1.0 m. At Harmon Lake, however, it is clear thatthe boundary layer had been intersected. Nonetheless, the temperature at 0.5 m (3.3 0C)was still much closer to the temperature at 1.0 m (4.10C) than Tm . In light of this verylimited data, a linear decrease from Tw to Tm over a thickness of 0.5 m is considered areasonable model. The equation for qw is therefore:93, aTqw = -Kw az z=0TW " TM (4.18) TEMPThe water temperature, Tw , is calculated in TEMP and control is returned to the mainprogramme where the time step is incremented. Intermediate steps include thecalculation of Iw , the energy available to heat the water following reflection at the surfaceand absorption in the cover, and Qsed, the sediment heat transfer. The net heat is givenby:clnet = Iw Qsed qw,and the change in temperature in one time step, At, is:AT — tlnet AL AtV Cpw Pwwhere, AL = lake areaV = lake volumeFrom §4.1, it can be seen that Iw is calculated as:Iw = Al 10 exp (-Xs ihs - ?eihe - Xil ) + A2 10 exp (-Xs2hs - Ae2he -The term Qsed is considered in detail in the following section.4.3.3 Sediment Heat Transfer Estimates(4.19)(4.20)(4.21)Using Falkenmark's empirical relationship (Ashton, 1986; see Chapter 2) as a guide, thesediment heat budget at both Menzies and Harmon Lakes should constitute roughly 25%of the heat budget of the water compared with only about 8% at Lake Mendota. Sincethe effect of latitude and altitude may be considered to be reflected in the annual heatbudget, average sediment heat transfer rates for both Menzies and Harmon Lakes may beestimated using data available for Lake Mendota and the following formula:(25 1 (clsed)m - (00H (clsed)H^8 '^(OL)Hwhere, qsed^= sediment heat transfer rate= annual heat budgetH^= Harmon and Menzies Lake(4.22)9495M^= Lake MendotaThe inherent assumption made here is that the sediment heat transfer is proportional tothe sediment heat budget.First it is necessary to determine the appropriate period over which qsed applies at LakeMendota. All three lakes reach maximum epilimnetic temperatures in mid-August. Datafrom 1908 (Hutchinson, 1957), however, indicate that fall turnover at Mendota maysometimes occur several weeks before turnover at Menzies and Harmon. With the onsetof ice-cover in January, Mendota may often turn over for up to 3 months compared withabout 1 month for Harmon and Menzies as observed in 1991. This information suggeststhat a greater proportion of stored heat is released during turnover at Mendota, and thatthe flux of 2.9 W/m2 calculated by Birge et al. (1927) may be appropriate in spite of thetwo weeks of ice-free conditions that the associated period includes. Furthermore, datafor use in the MLI model is available starting only two days before the beginning of theperiod considered by Birge et al. (1927).Estimates of the annual heat budget, 9L, are also required. For small lakes it is importantto include the latent heat required to melt the maximum observed ice thickness, inaddition to the heat needed to bring the ice to the melting temperature (Hutchinson, 1957)Excluding heat absorbed by the sediments, (OL)M is 924 MJ/m 2 (Birge et al., 1927).Using the data collected in 1991-1992, (OL)H is estimated at 562 MJ/m 2 for MenziesLake, and 549 MJ/m2 for Harmon Lake. It has been assumed that the maximum heatcontent associated with the hottest summer is 110% of the maximum heat contentobserved in the summer of 1991. Similarly, a factor of 90% was applied to the minimumheat content observed over the winter period. It is interesting to note that in spite of thesevere winter cooling due to artificial circulation at Menzies Lake, the annual heatbudget, and therefore the expected sediment heat transfer rate is almost equal to that atHarmon Lake. This is a reasonable result, however, given that the two lakes are about 5km apart, are at approximately the same elevation, and are similar with regard tochemistry and sediment type. Using Equation 4.22, the expected average sediment heattransfer rate is 5 W/m2.For comparison, the sediment heat transfer rate can also be estimated using Equation 2.13and an appropriate sediment conductivity. Ty was estimated using the temperature datacollected from June 4, 1991 to April 18, 1992. The temperatures measured over thedepth of the water column were weighted in accordance with the lake morphometry inorder to compute the whole lake average temperature, Tavg, on the date of each field visit(field visits varied in frequency from once a week to once every three weeks). Theannual average Ty was approximated by weighting Tavg in accordance with the samplingfrequency. From the literature, an average conductivity of 1.2 W/m°C produces a heattransfer rate of about 2.8 W/m 2 at Harmon Lake, if an average Tv, of 3 0C is assumed.This is only about half the rate estimated using Falkenmark's (Ashton, 1986) empiricalresult.The average rate could very well be even lower than the second estimate, however, giventhe qualitative similarity between the lake sediments of this study and the WisconsinLake sediments for which Likens & Johnson (1969) reported conductivities equal to thatof water (about 0.6 W/m/°C). The deep sediments at Menzies and Harmon Lakes are softand highly organic. A spherical 25 kg anchor lowered from the water surface wasobserved to penetrate over 0.5 m into the sediments. No qualitative observations,however, were made of the shallow (epilimnetic) sediments from which the bulk of theheat transfer would be expected to originate. Since finer grained material would not beexpected to settle out in the littoral zone as readily as in the relatively quiescenthypolimnion, it would be reasonable to expect coarser, more conductive sediments in theshallow regions of the lakes. Furthermore, as suggested in Chapter 2, low level96turbulence under the ice cover may lead to greater heat transfer rates than predicted bythe conduction equation.Clearly, further work is required to better quantify sediment heat transfer. As a firstestimate, however, the second method (Equation 2.14), using an average conductivity of1.2 W/m°C, was selected for use in MLI since it is based on theoretical considerations.4.4 MLI-C: ARTIFICIALLY CIRCULATED LAKEThe MLI model is a one-dimensional balance of heat flux across an ice and snow cover.All fluxes are assumed to be uniformly distributed over the entire lake area. It has beenmodified, however, in order to account for the non-distributed effect of the Air-o-lator®device. The procedure for calculating the net heat transfer across both the polynya andthe ice-covered area is similar to that used by Patterson & Hamblin (1988) to account forthe effect of a partial ice cover. The modified version, MLI-C, is described here.The new key variables which must be calculated at each time step is the radius of thepolynya generated by the circulator, and the turbulent heat flux from the water to the icein the vicinity of the polynya. It is assumed that the polynya is circular at all times, withthe circulator at the centre. Beyond its importance with respect to sensible andevaporative heat loss, the wind is neglected as far as the shape of the polynya and thesurface currents are concerned. This is reasonable given that, upwind of the circulator,surface currents are impeded while downwind they are assisted by the wind. The lakemorphometry may also play a role in the polynya shape and size, but this possiblity hasalso been ignored.Little information is available regarding the shape and other characteristics of thetransition region between the open water and full ice and snow cover. This region may97HPOLYNYA ICE - COVERED REGIONHP Hpe airsnow(qo)ppolynya ice(qf)pesnow-ice1qt(r) + qwsubscript 'p' refers to the polynyasubscript 'pe' refers to the polynya edgef qt(r) + qwiceI qfbe of significant importance to the flow regime, and hence the heat transfer at the ice-water interface. The lack of information described above, however, is such thatprediction of the shape of the transition region is not justified. Instead, a simple model,based on the assumption that there is an abrupt change from open water conditions tosome minimum ice thickness, will be employed to calculate the polynya radius, and anaverage turbulent heat transfer from the water to the ice.All heat fluxes across the lake surface, including the polynya and the ice at the polynyaedge, are shown in Figure 4.8. Details are given in the sections that follow.98Figure 4.8 Heat Fluxes Across Ice-covered and Ice-Free Zones of LakeJust as MLI requires the input of initial non-zero ice cover conditions, MLI-C requiresinput of an initial polynya radius (although a radius equal to the average lake radiusshould work as well, under ice-free conditions). The second input parameter required isthe flow rate, Q, of the Air-o-lator®. Here it is assumed that the surface velocitiesproduced are proportional to Q and inversely proportional to r, the radial distance awayfrom the centre of the device. Detailed measurements were made of the surface velocityat Menzies Lake, which clearly varies as r-1 as shown in Figure 3.12. Thesemeasurements are consistent with the theoretical centerline velocity for an unconstrainedradial jet (see §2.4.4). From these measurements the following relationship was derived:6.13 Q U —^ (4.23)It should be noted that the above relationship is only appropriate for the Air-o-lator®device complete with flow deflector, as shown in Figure 1.1 (item #1).4.4.1 MLI-C SubroutinesThree new subroutines are added to MLI to account for the effect of the circulator. ThePOLYNYA subroutine is called just before the average lake temperature is calculated insubroutine TEMP. From POLYNYA, where the heat fluxes between the atmosphere andthe open water are calculated, subroutine NEWAREA is called to compute any change inthe polynya radius. If, in this latter routine, it is determined that surface melting of theice should occur at the polynya edge, the quantity of melt is determined in subroutineMELTPE.99POLYNYASubroutine POLYNYA is responsible for the calculation of the heat fluxes across the air-water interface and the turbulent heat flux, Q t, from the water to the ice. In general,Equations 4.6 to 4.13 are employed with T o replaced by Tw , the mean water temperature.Further modifications include recalculation of the non-reflected solar radiation to accountfor the much lower albedo of the open water and re-evaluation of the saturated vapourpressure (Equation 4.9) using the appropriate coefficients for water. In addition to thesemodifications, two additional fluxes are included in the surface heat balance to accountfor cooling due to snowfall on the open water (Equation 2.27) and the increase in theevaporative heat flux due to free convection (Equations 2.25 and 2.26). This latter flux isonly applied when the water temperature is greater than the air temperature. The net heattransfer across the polynya, H p , is given by Equation 2.28.The heat loss from the lake across the polynya area is insufficient in explaining the lowwater temperatures observed at Menzies Lake. In addition to conduction of heat throughthe ice cover, the remaining heat loss is due to the turbulent heat transfer, Q t, from thewater to the ice in the vicinity of the polynya. The extent of the zone in which turbulentheat transfer is important is unknown. This significant heat transfer mechanism may onlybe of importance over a small area in the vicinity of the ice-free patch. An importantconsideration which has not been dealt with is the profile of the ice edge. The steadystate profile may be such that the flow is deflected downwards to such an extent that,beyond a short transition zone, the velocities at the interface may be negligible. At theother extreme, the transition to full ice thickness may be so smooth that, outside a thinboundary layer the interface velocities continue to vary as r -1 with no discontinuity at theice edge. The average heat transfer,Qt, will be estimated by integrating qt as defined by100Equation 2.30 (see Subroutine NEWAREA) from the ice edge to some radial distance, r.. :,away from the edge:^27c^I.'^12.26n ct Pw cp Q (Tw - Tm) (r00 - rp)^Qt — AL _ Ap  r dr — AL - ApP rp(4.24)where Ap = polynya arearp = polynya radiusThe radial distance, r., will be determined by calibration of the model. It is assumed thatthis parameter varies with the size of the polynya. If friction dissipates the velocity to anegligible value over the same distance under the ice, regardless of rp, then thecalibration parameter becomes Or. Therefore r o. = rp + Or. The value of ct, the turbulentheat transfer coefficient, is discussed in the section on Subroutine NEWAREA whichfollows.NEWAREAThe flowchart for Subroutine NEWAREA is given in Figure 4.9 (note that most variablesshown in the figure are modified by the subscript 'pe' in the model to indicate valuesassociated with the polynya edge). The polynya size is determined at each time step inNEWAREA in a manner similar to the calculation of the ablation and accretion of ice atthe ice-water interface for a lake with full ice cover. In this case, the change in icethickness at the polynya edge is determined as follows:dh (qf)pe - qw - qt. dt —^p e.L(4.25)where, (qf)pe= the flux of heat through the ice at the ice-water interface of the polynyaedge,and all other terms are as previously defined.101Figure 4.9 NEWAREA FLOW CHART* hi, To, Qc, Qe, Rlo, Qr, Qsi, go, Hp and of should be modified bythe subscript 'pe' to indicate values associated with polynya edge102The molecular heat flux from water to ice, q ty , is calculated using Equation 4.18. Moredifficult is the evaluation of q t which is expected to decrease rapidly with distance fromthe circulator, possibly in accordance with the empirical surface velocity relationship(Equation 4.23). The turbulent heat flux, q t, will be assumed to be dependent on thesurface velocity, U, and the difference between the water temperature and the meltingtemperature. The appropriate equation will have the same form as that given in Hamblin& Carmack (1990):qt = Ct pw Cp U(Tw - Tm) (2.30)As described in Chapter 2, the velocity used in Equation 2.30 should be measured 1.0 mbelow the ice-water interface. In all probability, a velocity at 1.0 m is inappropriategiven that the sensible heat transfer equation was developed for a flow which can beconsidered uniform with depth as in the case of an inflowing or outflowing river (seeHamblin & Carmack, 1990). Similarly, the coefficient C t is possibly only appropriate forthese types of situations. Here C t must account for the turbulence generated by thecirculator and the shape of the boundary layer at the ice-water interface. This value of Ctmay require adjustment in order to achieve good agreement with observations. However,a unique value for C t is conveniently available using the information given in §2.4.3 andthe flow field observations of §3.2.4. As can be seen in Figure 3.12, the velocity over thetop 20 cm tends to uniformity as the jet moves towards the polynya edge. Right at theedge it is likely that uniformity exists over a much greater thickness. It may therefore bereasonable to treat the extrapolated radial velocity given by Equation 4.23 as the free-stream velocity, U, in Equations 2.30 and 2.35. Given that the radial jet velocity variesas r', Equation 2.35 may be used to calculate a constant value of C t, if x is set equal to r.By inserting Equation 4.23 into Equation 2.35, this value is 5.9 x 10 -4 . This number is atthe low end, but within the range given by Hamblin & Carmack (1990).103This value of Ct, however, cannot be considered reliable since it is based on experimentalwork using a laminar, uniform and uni-directional free velocity field. It is more likelythat C t is at the high end or above the range given by Hamblin & Carmack (1990) giventhat the radial surface jet produces maximum velocities which are less than 1 m from theice-water interface. In the river flow examples by Hamblin & Carmack (1990), thevelocities decreased monotonically from 1 m to meet the no-slip condition at theinterface. In any event, it is reasonable to suspect that the velocity profile within 1 m ofthe interface is unique for this problem and Ct must be determined by calibration.The heat flux through the ice at the ice-water interface, qf, as calculated for the ice-covered area of the lake, is inappropriate for the polynya edge. First, to be consistentwith field observations, the snow layer is removed in the vicinity of the edge.Furthermore, the minimum ice thickness, h min , concept employed by Patterson &Hamblin (1988) to account for partial ice cover is assumed to apply to the ice at thepolynya edge. This minimum thickness should, based on observations, be less than the10 cm employed by Patterson & Hamblin (1988). It will, however, be used as anadditional calibration constant for the model (see Chapter 5). In addition, this thicknesswill be assumed to consist entirely of white ice, as field observations have indicated.This latter assumption is supported by the fact that lake surface agitation (such as thatproduced by the circulator) will lead to granular white ice during freeze-up (Gray &Male, 1981). (It is assumed that the properties of white ice are the same as snow-ice). Anew value of To is calculated for this white ice layer in order to determine if melting willoccur at the surface. For the polynya size to increase, the entire ice layer must be meltedat the radius defined by the location of the polynya edge. Here, an average value of ofcould be employed to reflect the ice thickness over the course of the time step as itreduces to zero. However, given that the model requires calibration using three constants(Or, C t, and hm i n), calculation of an average value using an integration technique is104unjustified. It is therefore only important that the value of of reflect the relative day today changes in meteorological fluxes. This will be achieved for of defined at anyreasonable value of knit' .In the above procedure it is assumed that there is an immediate transition from openwater to full ice cover. This being said, it should be noted that it is also assumed that theice cover has no effect on the surface velocity as given by Equation 4.23.The computational procedure begins with calculation of the heat fluxes at the surface inthe same manner as for the fully covered ice region. T o , here the temperature of thewhite ice surface at the polynya edge, is again determined using Newton's iterativeprocedure. The heat associated with flooding, Qsi, is included in the polynya ice model,and, in this case, is assumed to be distributed over the entire ice thickness. Melting mayoccur at the surface just as it does for the rest of the lake. This process is quantified inSubroutine MELTPE.At this stage, the change in ice thickness at the polynya edge is calculated using Equation4.25. If the new ice thickness is less than zero, then the polynya is assumed to beundergoing expansion, and Equations 4.23, 4.25 and 2.30 are solved for r, the radius atwhich the polynya edge ice thickness is completely melted over the time step. IfEquation 4.25 produces an ice thickness which is greater than zero, contraction of thepolynya occurs and these equations are again solved simultaneously. In this casehowever, of and qw are replaced by Hp, the net meteorological flux between the openwater and the atmosphere, in Equation 4.23:dh -Hp - qtdt — Pe L(4.26)105This is done to reflect the fact that accretion is assumed to be occurring on the face of theice edge, where it is exposed to Hp . Melting, on the other hand, is assumed to occur atthe bottom of the ice cover. This point is discussed further in Chapter 5.MELTPESubroutine MELTPE is equivalent to Subroutine MELT, except that it is only applicableto the polynya edge. It is much simpler than MELT because of the assumption that onlywhite ice exists in this region. The surface temperature, To , is set to Tm , and the netmeteorological flux, H, and the flux through the ice at the air-ice interface, q0 , arecalculated. The imbalance between these fluxes determines the amount of ice which ismelted (Equation 2.9). The modified ice thickness at the polynya edge is calculated andcontrol is returned to Subroutine NEWAREA.Modifications to Subroutine TEMPSubroutine TEMP must be modified to include the heat transfer terms related to thecirculator. To Equation 4.17 must be added the turbulent heat flux from the water to theice, Q t and the net heat transfer across the polynya, Hp:clnet = Hp + 6 + qsed - qw - Ot^ (4.27)The components of net are scaled appropriately to equivalent values which all applyover the same area.4.5 DATA PREPARATIONSome of the data retrieved from the datalogger must be adjusted and transformed in unitsof measurement required by MLI-C. Wind data is required for a point 10 metres above106the lake surface. Following the National Research Council of Canada (1985), a powerrelationship describing wind speed profiles was used to transform wind speed at 2.0 m towind speed at 10.0 m. The appropriate correction factor is 1.25. Although the airtemperature was not measured at the correct height for strictly proper use in the bulkaerodynamic equations, a simple correction factor is not available. The development ofthermal inversions is extremely important in determining the nature of the boundary layerabove the surface. Strong inversions would have occurred and persisted over more thanhalf the average day in mid-winter. These inversions would become weaker and shorterin duration as melting conditions began to predominate (Geiger, 1965). Therefore, acharacteristic daily average thermal profile which is applicable to the entire ice-coveredperiod does not exist. Furthermore, the coefficient Ch in Equation 4.12 was determinedusing data for ice-free conditions and may not be appropriate for a lake with ice cover.As these difficulties are beyond the scope of this investigation, the raw temperature datawill remain unmodified for use in the model.The solar radiation data required integration to give daily totals. The radiation isassumed to remain constant over the half hour period surrounding the time ofmeasurement. The total, therefore, is equal to the sum of the products of eachmeasurement and the sampling period. For all other meteorological data, daily averageswere calculated using the 48 half hour readings.4.6 SUMMARY OF PARAMETERSA summary of the values chosen for the parameters required for MLI and/or MLI-C aregiven in Table 4.2. Some of these parameters will be varied to determine the sensitivityof the model to their values or adjusted for calibration purposes (Chapter 5). With regardto solar attenuation, it is clear from the literature (and intuition) that 2e should lie107between 24 and Xs. Therefore an average of Xs and Xi as given by Patterson & Hamblin(1988) was selected as an appropriate first estimate of A. Similarly, the conductivity ofsnow-ice should lie in between that of snow and that of ice. Since the voids in the snoware filled with water and subsequently frozen, it is reasonable to expect K e to be closer toKi than Ks. A value of 2.0 W/m0C was selected.Table 4.2 Summary of ParametersParameter Value CommentXit 1.50 m-1 see Patterson & Hamblin (1988)2si 6.00 m-1 see Patterson & Hamblin (1988)Xei 3.75 m-1 see text above2i2 20 m-1 see Patterson & Hamblin (1988)Xs2 20 m-1 see Patterson & Hamblin (1988)Xe2 20 m-1 see text aboveKs Equation 2.16 see §2.4.1Ki 2.3 W/m/°C see Patterson & Hamblin (1988)Ks 2.0 W/m/°C see text aboveKs 1.2 W/m/°C see §4.3.3Pi 917 kg/m3 see Patterson & Hamblin (1988)Pe 890 kg/m3 see Ashton (1986)Pi Equation 4.15 see §4.3.2Ct *** *** calibration constant (see §5)hmin <0.10 m see §4.4.1Or *** *** calibration constant (see §5)as variable see §4.3.2ai 0.02 see §4.3.2asi 0.45 see §4.3.2aw 0.06 see §2.4.2Pn 80 kg/m3 see §4.3.2P1 100 kg/m3 see §4.3.2Pm 200 - 300 kg/m 3 see §4.3.2108Chapter 5RESULTS AND DISCUSSIONThe models described in Chapter 4 were written in BASIC and run on an 80386 DX (25MHz) PC equipped with an 80387 math co-processor. For 105 days of data, less than 10seconds were required to produce results for either MU or MLI-C.5.1 LAKE IN NATURAL STATEThe results produced by MLI for Harmon Lake using the default parameters given in §4.6are shown in Figure 5.1. Given the approximations made in the development of themodel, the results are very good. The estimated errors in the observations account formost of the difference between the observed and predicted variable. It will be shown thatreasonable variations in some of the default parameters provide an improved fit to theobservations. The model continues to provide good predictions even in late March whenan ice-zone 1 to 4 metres in width around the perimeter of Harmon Lake was observed.This observation emphasizes the fact that, although assumed to be uniform, the icethickness can be highly variable over the lake surface. Results using the DYRESMIformulation are also given in Figure 5.1. A detailed analysis of the results follow theproceeding discussion on the observation errors, which are needed in order to assess themodel's success in deciphering nature's work.1090.0Dec 13, 1991^ snow (observed)• ice (observed)— default parameters^ P & H, 1988 , . . ... ,. . . ..... .. . . .. ........... . ..'.. ••....Jan 19, 1992^Feb 8I5.5 —A observed temperature— default parameters---- Patterson & Hamblin, 1988........................................................ 0. ...... •••••3.0 —12.5^ II^I 15.0 —4.5 —4.0 —3.5 —DateDec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.1 Ice and Snow Thickness and Whole Lake Temperatureat Harmon Lake: A Comparison of Results Using MLI andPatterson & Hamblin's (1988) Model1105.1.1 Observation ErrorsThe error in the ice and snow observations were estimated using all of the thicknessesobserved at the two Menzies Lake stations over the simulation period. The thicknessesobserved at Station 1 were plotted as a function of the thicknesses at Station 2, and alinear regression was performed. The error bars shown in Figure 5.1 represents ± 2standard deviations from the regression line. With regard to the ice thickness, this errorwas assumed to be applicable to both lakes. At Harmon Lake, the expected error wouldbe larger given that measurements were only taken at one station. However, there is alsoless expected variability over the entire lake surface because there is no circulation devicein operation which may be responsible for horizontal variability at a scale that iscomparable to the lake length. It is assumed that the impact of the first factor on thestandard error negates the second. Turning to snow thickness, the standard error for themeasurements at Harmon Lake must be greater than that at Menzies Lake given thatthese measurements were only made at one station. If it is assumed that the variance ofthe regression is equal to the variance of the sample mean from the population mean (thisis reasonable for the almost uniform distribution of snow thickness), and that the truepopulation variance is the same at both lakes, then the sample variance for the HarmonLake snow thickness data may be estimated by the Law of Large Numbers:a 2 _ 02-x Nwhere, ax2 = variance of the sample mean from the population meana = population meanN = number of samplesSince N=1 at Harmon Lake, and N=2 at Menzies Lake, a reasonable estimate of thesample variance at Harmon Lake would be twice that at Menzies Lake.111With regard to whole lake temperature observations, the temperature probe is consideredaccurate to ±0.3°C by the manufacturers. This error may be considered too large giventhe following reasons:a) periodic checks with a precision mercury thermometer indicated ±0.1°C accuracyb) whole lake temperature data is based on 16 measurementsHowever, given the following additional sources of error, an accuracy of ±0.3 0C isconsidered reasonable:a) there is error in the morphometric data used to calculate the whole laketemperaturesb) interpolation between measurement depths was requiredc) horizontal variations in temperature due to sediment heat transfer were notconsidered.5.1.2 Analysis of MLI PredictionsA Comparison with DYRESMIThe MLI model provides an improved prediction of a mid-latitude lake winter heatbudget over the DYRESMI model. Considering the latter model first, the initialpredicted total ice and snow thicknesses do not match with observations because theinitial ice thickness could not support the weight of the observed snow (Figure 5.1). Thisis explained by the high snow density assumed to be appropriate for the lakes consideredby Patterson & Hamblin (1988). Hence DYRESMI predicts instant flooding of the icecover from the initial conditions. Furthermore, DYRESMI predicts ice thicknessesalmost 50% greater than the observations throughout most of the modeling period. Thisis primarily due to the relatively high snow conductivity associated with the assumeddensity. Using Equation 2.15 and the range of densities considered in this investigation,112MLI predicts snow conductivities of about 20% to 60% of the constant value of 0.31W/m°C assumed in DYRESMI.In addition, the DYRESMI model generally over-predicts the snow thickness for twoprincipal reasons. First, it is assumed that no settling or wind transport of snow occurs,and secondly, rain melt is ignored. Some surface melt is predicted by DYRESMI at theend of January, when most of the rain fell, because the meteorological flux balanceproduces a surface temperature of 00C on several consecutive days. Neglect of heatreleased to the snow cover due to flooding may also contribute to the inadequatesnowmelt predicted by DYRESMI.It should be noted that, in spite of the deficiencies of DYRESMI, the trend in icethickness is very well predicted from mid-January until the end of the modeling period.Although the rate of ice melt is somewhat under-predicted in March, the predicted ratesare only marginally less than MLI predictions, and furthermore, explanations involvingthe quality of the meteorological input data are available (see Chapter 3). This is in spiteof significant discrepancies in snow thickness prediction. This is, for the most part, dueto the fact that excess snow is converted to pure ice in DYRESMI, not snow-ice. Thisresults in much less attenuation of solar energy through the total ice thickness, therebyoffsetting the effect of the thick snow layer. Also, as stated earlier, the conductivity ofthe DYRESMI snow is much greater than that of the relatively thin MLI snow layer,which tends to increase the heat flux into the lake during the melting phase. The samemay be said of the pure ice which DYRESMI produces as a result of ice cover flooding.Furthermore, DYRESMI predicts a much lower albedo (a = 0.6) than MLI for a briefperiod due to above zero air temperatures immediately following the mid-Marchsnowfall.113Turning to lake temperature, DYRESMI predicts cooling over the first two month of thesimulation because of the neglect of sediment heat transfer. Less solar heating is alsopredicted during some periods when over-prediction of snow and ice thicknesses occur.DYRESMI also under-predicts the rapid increase in lake temperature in March. This canonly be explained by an under-prediction of solar penetration. There are two majorfactors which lead to this result. First, MLI, on average, predicts a lower albedo thanDYRESMI over the latter period of the simulation. Secondly, there is high solarattenuation through the thick snow and ice cover which is not completely offset by thelack of snow-ice as described above.Where Improvement is NeededThe MLI predictions are almost always within the observations of total ice and snowthickness, as well as temperature. Consideration will now be given to where possibleproblems lie, and how results could be improved. The most notable deviation fromobservation is the snow thickness prediction in early January. As described in Chapter 4,the initial snow density could have been no more than about 90 kg/m 3 to be supported bythe observed ice thickness. The Merritt precipitation data indicates that most of thissnow was probably deposited over the previous two days. Unless snowfall events werenot recorded between the first and second observation dates, or the two early Januarysnowfalls were grossly under-estimated by the Menzies Lake resident, then very littledensification and wind transport of the existing snow must have occurred. Clearly, morework is required in order to derive a physically based model of snow densification andwind transport. In spite of this problem, all of the remaining predictions are within theobservation errors. This may be partially a result of the snow load restriction, but it isalso possible that a meteorological event such as the freezing rain of January 9th mayhave increased the average snow density to within the range of the snow density formulaused in the model. Furthermore, following the rainy period at the end of January (see114Figure 5.2), the snow is so thin that the model is not well tested as far as snow coverpredictions are concerned.It would also appear that the rate of snowmelt due to rain is over-predicted. MLIproduces completely snow-free conditions by the end of January and a net reduction inice thickness over four consecutive days of rain. On February 4th, however, up to 2 cmof snow still existed on the lake. The predicted reduction in ice thickness due to rain isdue to melting both at the surface and at the ice-water interface. The latter melting isprobably well predicted since its only dependence on the surface balance is through thevalue of To , which, for melting conditions, must be 00C. The amount of surface melt,however, depends directly on H, of which Qr is an important part. As described inChapter 4, the latent heat was included in the surface balance because both sensible andcondensation heat is under-predicted due to the use of average daily meteorological data.The trend in ice thickness between the 18th of January and the 4th of February does notseem to be consistent with the observations. Ice growth appears to be under-estimatedbetween these dates. There are two likely explanations for this. First, the snowconductivity may be under-estimated using Equation 2.15. Secondly, there isconsiderable rainfall and snow-ice production over the interval. The incorrect predictionof complete snowmelt, and a further melting of ice from the top has been discussed, but itis also possible that ice melt from the bottom due to snow-ice production has been over-estimated. This would have occurred if the prediction of heat generated by the snow-iceis excessive. Such a problem, however, would be offset by an over-estimation of snow-ice production. Finally, with the accumulation of a significant amount of snow-ice in thelatter part of January, the snow-ice conductivity increases in importance (although it wasassumed that the initial ice cover was 50% snow-ice). If snow-ice conductivity is moreaccurately represented by the higher value of pure ice, then greater ice thicknesses wouldbe expected during subsequent cold weather.115116Surface Heat Budget Components•A.r^4", .•,-•, , r.4.0.^a?. „A^a.^ .A, 6At/ "' NA- liby v u sz 2kotorh, vv.,/300-200-10o-ET-100 —0 ^,%1 ^ Ri (in) ^ Ri (out)Qc (out) ^ Qe (out)- - - Io (in)^o Qsi (in)Qrain (n)r%^I •^• ;;•-;71.1`^****************************^• ... -4/^1 1 ^.. II ^►^1► 1^►^1► 1 ► 1^i^ri^11^III^IN ► 11^/^1^1r-- •I „ii -- tl^1 1 I II^1I ►^I1 I11 Iit^I^1I IDec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateHeat Fluxes to and from Lake Water^100(left hand axis: Qsed, qw ; right hand axis: I,„)Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 88(-4E60 4—• .=20^ qsedqw—806040 ,N)200Feb 28^Mar 19• •••^• •• •^•^• .• • „^•DateFigure 5.2 MU: Surface Heat Budget Componentsand Heat Fluxes to and from Lake WaterThe rate of ice melt in March is also under-predicted. Several explanations are available.The difficulty of accurately predicting the presence (or absence) of very thin layers ofsnow, which may have a huge effect on the albedo, is recognized. The weather stationtilt likely resulted in some error in the radiation data. The open ice area around theperimeter of the lake which developed in March probably served to encourage convectivecirculation in the lake, resulting in an increased transfer of heat from the water to the ice.Convective circulation is discussed further in Chapter 3. It is probable that this lastmechanism was most important in increasing ice melt since it does not necessarilyinvolve an increase in lake temperature, the trend of which is quite well predicted overthe last month of the simulation.The lake temperature is actually very well predicted on the whole, aside from theimprobable surge in mid-February. The Menzies Lake weather station data (see §5.2)indicate that a minor surge in temperature was likely but, with the incorrect prediction ofcomplete snowmelt, solar penetration would have been over-estimated, thereby causingtoo large an increase in the rate of lake warming.The effect of varying the most influential and least well established parameters on snow,ice and temperature will be studied in §5.1.4. In the section which follows, the resultswill be interpreted in terms of the individual heat flux terms, and the surface heat budget.5.1.3 Components of the Heat BudgetThe Surface Heat BalanceThe components of the surface heat budget and the direct heat fluxes to and from the lakewater are shown in Figure 5.2. The meteorological balance results in a net loss of heatfrom the lake surface over the first half of the simulation, facilitating a period ofsustained ice growth. Most important in the early winter balance is a net emittance of117long-wave radiation. The evaporative flux is generally negligible over the entiremodeling period. The sensible heat transfer almost always constitutes a source of heat tothe lake surface, not a sink. This latter result is due to the fact that the solution to theconduction equation yields surface temperatures, T o , which are consistently below airtemperatures, Tair (more on this later). Both Ashton (1986) and the U.S. Army Corps ofEngineers (1956) provide thermal profile measurements across the air-snow interfacewhich are consistent with these results. No examples were found in the literature whereTo > Tair for a snow covered surface.Although there is, on average, an increase in the net emittance of long-wave radiationtowards the end of the simulation period, both solar radiation and sensible heat transferincrease substantially. This results in a reverse in the flux imbalance at the ice-waterinterface thereby causing the ice to melt.Although most of the transmitted solar radiation is absorbed in the cover (about 90% inDecember, decreasing to about 75% in March), the heat is distributed throughout thecover, albeit not uniformly, thereby diminishing the tendency to counter the heat flux outof the lake. Therefore, when solar radiation appears to tip the balance in favour of a netgain of heat in early February, there is actually a continued period of heat loss and iceformation until well into the latter half of the month.It can be seen from Figure 5.2 that the heat due to rainfall is in most instances, sufficientwhen combined with high rates of sensible heat transfer into the lake to reverse thedirection of the surface heat flux, and cause surface melting. The rainfall of January 9th,was described as freezing rain by the Menzies Lake resident, which explains therelatively low sensible heat transfer on this date compared with the other rainy days.(Although the latent portion of the heat produced by the rain was not removed on thisdate, the total rainfall was light, and surface melt was not predicted. The reduction in118snow thickness on this date is due to the assumed increase in maximum snow densitywhich occurs when rain falls).Since almost all of the heat generated by the rainfall is assumed to be latent heat, therewould not be sufficient melt predicted if this heat were not included in the surfacemeteorological balance. According to Harr (1981), there should have been significantamounts of heat derived from condensation (as quantified by Qe) to produce surface melt.The warm temperatures and high humidity which likely persisted over the rainfall periodsare not reflected in the daily averages. Although daily average meteorological data areuseful for computing daily average heat fluxes which are active over the entire day, theyare clearly not appropriate for determining heat transfer due to precipitation which isusually only active over much shorter periods, when the meteorological variables are atextreme values. It is therefore concluded that a purely theoretical representation ofsnowmelt due to rain will not produce accurate results unless sub-daily averages ofmeteorological data are used.The heat produced by ice cover flooding is substantial (again, the bulk of the heat islatent). The impact of this term on the heat budget is considered in §5.1.4.The heat fluxes to and from the unfrozen lake are also shown in Figure 5.2. Althoughsolar radiation is responsible for almost all of the heat absorbed by the lake water overthe modeling period, there are a few significant periods over which sediment heat transferrates exceeded the rate of solar heating. In fact, from December 13th until the end ofJanuary, the average rate of solar heating was less than the average 2.8 W/m2 deliveredfrom the sediments to lake water. As described in §5.1.4, however, the model resultsagree more with the field observations if the sediment heat conductivity is reduced by50%. In that case, the solar heating over this period would be slightly greater than119sediment heating. By February, however, solar heating dominates strongly as thesediment heat transfer begins to fall in response to a reduced thermal gradient.Over most of the period, sediment heating is only slightly more than offset by heatconduction through the ice. The gap widens, however, starting in February, as the lakebegins to heat up more rapidly. By the end of the simulation period, conduction throughthe ice cover is more than four times the diminishing sediment heat transfer rate. At thistime, however, conduction of heat to the ice is only 5 to 10% of the daytime averages ofsolar heating.Comparisons with DYRESMIWith a few notable exceptions, there is little difference between the DYRESMI and MLIsurface heat balance. Of course, Q and Q si do not appear in the DYRESMI heat balance,but the Rli term is identical, and there is no significant difference in Qe which is againnegligible. The remaining terms for both models are shown in Figure 5.3. The differentalbedo representations account for the mismatch in I o. The difference is most significantin mid-February, and during some periods in March, when MLI predicts snow-freeconditions, allowing a great increase in the quantity of solar heat penetration into thelake. This comparison emphasizes the importance of a good albedo estimate in latewinter when solar radiation dominates the heat budget.Both Rio and Q correspond much more closely. The most significant deviations occur atthe start of the simulation, when immediate snow-ice formation is predicted usingDYRESMI, and in early February when MLI predicts complete snowmelt. Thedifferences in these two terms are due to To. At the start of the simulation, the immediatesnow thickness discrepancy, coupled with the relatively high conductivity of DYRESMIsnow produces a greater conduction of heat to the surface, thereby increasing To . Theprediction of a snow layer in mid-February, according to DYRESMI, means that there is120100 —^ R1(out)- Qc (out) DYRESMI- - - Io (out)— MLI Results200-100Heat Fluxes Through Ice Cover — 60-20-150 ^I^ I— goof- qwAir and Surface TemperaturesToComparison of Surface Heat Budget Components: MLI and DYRESMI300 —Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateMLI ResultsDec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.3 a) MLI & DRESMI: Comparison of Surface Heat Budget Componentsb) MLI: Heat Fluxes Through the Cover, Air and Surface Temperatures121reduced heat conduction to the surface, resulting in an under-estimation of T o. The lackof solar heat in the surface layer due to the high snow albedo also contributes to this latterresult. Hence the sensible heat transfer into the lake is over-estimated and the long-waveemission is under-estimated in this case, while the reverse is true at the start of thesimulation.As described earlier in explaining the rate of ice melt predicted by DYRESMI, relativelyhigh conductivity leads to slightly more sensible heat transfer into the lake and somewhatless outgoing long-wave radiation during the melting phase to offset the relatively lowpredictions of solar heating.Heat Flux Through the CoverThe heat fluxes qo and of are shown in Figure 5.3. These two terms are equal over thenight time step, when there is no solar heat distribution in the cover. The extremedaytime deviation from this equilibrium illustrates the importance of solar radiation incausing melt at the ice-water interface. The exception to this rule occurs during snow-iceproduction in January, since the associated heat is distributed over both daily time steps.A consistent net upward flux of heat through the ice, which nearly always exceeds therate of heat conduction to the ice, results in the almost steady gain in ice thickness upuntil mid-February. The model predicts that this flux is reversed for a significant part ofthe latter half of January, resulting in a small net loss of about 1 cm of ice. This is due toice cover flooding, rainfall, and warm air temperatures. The actual ice melt predicted atthe ice-water interface is much higher, since there are significant accumulations of snow-ice during this period. Evidence from the field observations, however, suggest that therewas little ice melt to offset the snow-ice growth. Further work is required to betterquantify both Qr and In the latter half of the simulation, daytime solar heatingresults in significant melting at the ice-water interface. Even in late March, however,122there is some ice growth at night, although it is much less significant than the daytimemelt.Both Tair and To are plotted for the simulation period in Figure 5.3. As described earlier,it is predicted that To is almost always less than Tair, which explains the sensible heattransfer into the lake. The temperature data indicate a sustained surface melt over severaldays in late January. This was made possible by including Qn in the surface heat budget.During this brief period, surface melt is predicted for both time steps. In contrast, surfacemelt only occurs in the first time step in late March due to high levels of solar heat in thecover.5.1.4 Impact of MLI Parameters on PredictionsThe impact of change to some of the default parameters and functions are examined here.The most likely combination of changes is considered in Chapter 6.Effect of Ice Cover FloodingThe impact of flooding on the Harmon Lake heat budget is illustrated in Figure 5.4.Here, results are given based on the expected snow-ice formation, but with no associatedheat addition to the snow (Qsi = 0). Several snowfalls occur in mid-January whichexceed the bearing capacity of the ice cover. Figure 5.3 indicates that the flooding causesa significant reversal in the balance between of and qw , which translates into substantialmelting at the ice-water interface. The primary source of this heat, however, is the latentheat of fusion resulting from the formation of snow-ice. The net effect of this process isan increase in the total ice thickness. In fact, the net rate of growth is substantially higherover the period of these snowfalls than over the preceding month.If the formation of snow-ice produced no heat at all, the total ice thickness would havebeen expected to be about 2 cm greater over most of the remainder of the modeling123^ snow (observed)• ice (observed)— default parameters^ Qsi = 0DateWhole Lake TemperatureA observed— predicted^ Qsi = 0Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar^ 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.4 Effect of Heat Associated with the Formation of Snow-iceon Ice and Snow Thickness and Whole Lake Temperature124period. This heat, however, appears to have virtually no impact on either the expectedsnow thickness, or the whole lake temperature.More work is needed to quantify both the snow-ice thickness and the associated heat.One likely possibility is that high densification of the flooded snow occurs before itactually freezes. This would mean both less snow-ice and heat production. The expectedreduction in flooding depth may be offset by capillary rise of flood water in the snow.Effect of Snow-Ice Solar Attenuation CoefficientThe impact of doubling the attenuation coefficient associated with the visible band ofsolar radiation for snow-ice is shown in Figure 5.5. With regard to ice cover, there islittle dependency on these coefficients in the first two months of the simulation when thenon-reflected solar energy is well below daily averages of 100 W/m 2. Despite these lowlevels of radiation, the impact on the lake temperature is more notable because radiationis absorbed twice as quickly through the snow-ice (which is assumed to constitute 50% ofthe ice cover at the start of the simulation) before reaching the water. There is asignificant reduction in lake warming due to the higher attenuation coefficient by earlyFebruary. By the end of March, however, there is effectively no difference in the rate ofwarming because the higher attenuation results in a greater rate of melting, and thereforeless material through which absorption may take place.Regarding ice cover, the maximum impact over the entire modeling period is only lessthan half the expected error in the measurements. By the end of the simulation period,The lake temperature is less than the prediction based on the default coefficient by morethan the expected error.It seems reasonable to assume that solar attenuation may increase with snow density. Inthe early winter, when the average density is low, there would presumably be less125^ snow (observed)• ice (observed)— default parameters6.0 m 1I8 0.2 -U5.5 —A observed temperature— default parameterse = 6.0 m-12.55.0 —4.5 —4.0 —3.5 —3.0 —Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateDec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.5 Effect of Increased Solar Attenuation Through Snow-iceon Ice and Snow Thickness and Whole Lake Temperature126attenuation than in late winter when the average density can be up to four times higher.This possibility is beyond the scope of this project but should be considered in futureimprovements to MLI.Effect of Snow Density and ConductivityThe important effect of snow density on the lake condition is shown in Figure 5.6. Themaximum snow density in Equation 4.15 was increased to 300 kg/m 3 for To < 0 0C, and400 kg/m3 for To > 0 °C. This results in a substantial increase in predicted ice thicknessover almost all of the simulation period. The increased density translates into a greaterheat flux out of the cover because of the corresponding increase in snow conductivity asquantified by Equation 2.15, and a reduction in snow thickness. This increases theimbalance between of and qw, thereby creating more ice. Almost all of the increase overthe default condition occurs in the first 3 weeks of the simulation. The ice thicknessthereafter changes at essentially the same rate as the default case for the followingreasons:a) When rain falls, the maximum density is assumed in both cases to rise to 400kg/m3 . Therefore, following the first rainfall in early January, there is much lessdisagreement in the average snow density.b) The increase in of due to the higher snow conductivity is reduced due to theincreased insulation of the additional ice previously produced.c) The extra ice provides greater snow bearing capacity, thereby increasing the massof snow which could be supported. In spite of the higher rate of snowdensification, the net effect is slightly increased snow thicknesses over the defaultcase in the latter part of January. The extra snow provides a little more insulationto further temper qf. (The greater snow thickness in late January can also beattributed to a reduced rate of surface snowmelt because this rate is inverselyproportional the density.)In summary, the increased snow density facilitates a substantial increase in ice thickness,but has little impact on the whole lake temperature. A slight decrease in the latter1275.55.0 128^ snow (observed)• ice (observed)— default parameters^ (ps ). = 300 - 400 kg/m3i2 0.2 —Ug0.0Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateA observed— default parameters^ (P s )max= 300 - 400 kg/m33.02.5Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.6 Effect of Increased Maximum Snow Density on Iceand Snow Thickness and Whole Lake Temperaturevariable occurs due to the increased solar attenuation caused by the extra ice in the latterpart of winter.The impact of increased maximum snow density on snow thickness is most substantialwhen there is a completely fresh, low density cover and a subsequent period free ofsnowfall and surface melting.Patterson & Hamblin (1988) assume a constant snow density of 330 kg/m3 , and aconstant snow conductivity of 0.31 W/m°C. If this density is substituted in Equation2.15, the result is a conductivity of only 0.24 W/m°C. To examine the effect of snowconductivity while maintaining the original density function, Equation 2.15 wasmultiplied by the ratio of these two conductivities, 0.31/0.24. The result is given inFigure 5.7. The result is almost identical to the density adjustment described above. Theincrease in ice thickness is slightly less, but the predicted lake temperature is almostidentical to the previous case. The only substantial difference is the lack of a reducedsnow thickness (as caused by increased densification) in the early stages of thesimulation.Effect of Sediment ConductivityIn the majority of the investigations regarding sediment heat transfer which werereviewed for this project, the sediment conductivity, Ksed, was found to equal that ofwater (about 0.6 W/m°C). A substantially higher conductivity was assumed for thedefault case because values of up to 3.5 times that of water have been reported in theliterature (see Chapter 2). If Kj were reduced from the assumed value of 1.2 W/m°Cback down to 0.6 W/m°C, the predicted agrees with the observed lake temperature towithin the expected observation error on all dates as shown in Figure 5.8. Thecumulative deviation from the default case, however, is only about 0.2 OC, which is alsowithin the expected error.129^ snow (observed)• ice (observed)— default parameters^ Increased snowconductivityDec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19Date5.5A observed temperature— default parameters^ increased snowconductivity3.02.5Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.7 Effect of Increased Snow Conductivity on Iceand Snow Thickness and Whole Lake Temperature1305.04.5o snow (observed)• ice (observed)— default parameters- • - — K sal = 0.6 W/m°C0.1 -0.0Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19Date5.5A observed temperature— default parameters•- - - . K sed = 0.6 W/m°C................................... ..2.5Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.8 Effect of Decreased Sediment Conductivity on Iceand Snow Thickness and Whole Lake Temperature3.01310.4 -0.3 -5.04.5There is effectively no impact on the ice and snow cover aside from a negligible increasein ice thickness late in the period owing to the reduction in heat conducted from the waterto the ice.Effect of Albedo Decay RateThe result of multiplying the USACE decay functions (Equations 2.16) by a factor of 0.5is presented in Figure 5.9. The higher albedos which are produced result in a substantialincrease in ice thickness over the default case during the period spanning from lateFebruary to early March. At that time, there is a light snow cover and solar radiation ismore than an order of magnitude more intense than over the first month of the simulation.Although solar heating is as important in early February, there is no snow cover duringthat period and therefore no change in the results. The only change in the predicted snowcover occurs in mid-March. The reduced absorption of solar energy at that time causes adelay in the complete loss of snow due to melting. Clearly, the albedo is of extremeimportance to the whole lake temperature in late winter, especially if the meteorologicalconditions produce a situation where an extended snow-free period may or may notoccur, depending on how much solar heat is reflected at the surface. Although this is nota significant factor for the two albedo functions compared here, careful examination ofFigure 5.9 reveals that the impact on lake temperature is still quite substantial over thefew days in March when this condition occurs using the two decay rates. The averagerate of change in whole lake temperature over the last month is only about 70% of thedefault case, producing results which are more consistent with observations. Assuggested below, however, it is more likely that the deviations from the observed aremore likely due to other parameters.These results suggest that the predicted albedo should have been lower over most ofMarch in order to match the rapid ice melt. Complete snow melt actually occurred earlier132^ snow^I observed• ice• half a decay ,— full a decay .1 predictedDec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19Date5.5observedhalf a decay— full a decay3.02.5Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.9 Effect of Reduced Albedo Decay Rate on Ice andSnow Thickness and Whole Lake Temperature1335.04.5than either decay rate predicted in the first half of the month. This may be partly due tothe heavy snowmobile traffic on the lake at this time. Whatever the cause, there wouldhave resulted a much reduced albedo over several days due to snow-ice exposure.More research is required in order to establish a more reliable relationship betweenalbedo and meteorological conditions. A more likely solution which would producemore accurate results is direct albedo measurement.Effect of the Thermal Gradient at the Ice-Water InterfaceThe effect of reducing the thermal gradient at the ice-water interface by 62.5% isillustrated in Figure 5.10. (The assumed distance, dz, over which the lake temperaturelinearly decreases to 0 0C at the interface is increased from 0.5 to 0.8 m.) As far as iceand snow thicknesses are concerned, the reduced conduction from the lake waterproduces results which are similar to those obtained from increasing the snowconductivity (Figure 5.7). This is due to the fact that the reduction in q w in the first caseis of similar magnitude to the increase in of in the latter case. The result is similar fluximbalances which leads to similar rates of ice growth.While the lake temperature effectively does not change with increased snowconductivity, it undergoes a substantial increase as a result of the reduced thermalgradient. By the end of the simulation period, the lake temperature is more than twice theexpected error greater than the observed value. Although this result suggests that thelower gradient is unlikely, the probability of decreased lake temperature predictions dueto reasonable variations in sediment conductivity, solar attenuation and albedo must beassessed. Furthermore, it is also quite possible that the gradient is significantly variableover the course of the winter. Convective mixing, as described in Chapter 3, may haveresulted in turbulent heat transfer to the ice in late winter, thereby greatly reducing therate of increase in lake temperature caused by solar heating.134o snow• ice^) observed— default I predicted^ dz = 0.80.2 —,.. —0.3 —5.5A observed— default----. dz = 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateDec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.10 Effect of Reduced Thermal Gradient at Ice-Water Interfaceon Ice and Snow Thickness and Whole Lake Temperature1354a"5.2 ARTIFICIALLY CIRCULATED LAKE5.2.1 Analysis of MLI-C PredictionsResults from MLI-C using the Menzies Lake data set are given in Figure 5.11. Allparameters used in MLI were set to their default values. Although reasonable variationsin some of these parameters have shown to produce better predictions, the correctcombination of variations remains unknown. The most likely changes to theseparameters are considered in Chapter 6. One of the key recommendations which stemsfrom this work is that more effort be made to narrow the range of possible parametervalues for a given location and set of (easily measurable) conditions.Model AdjustmentsSome modifications to the subroutines which are specific to MLI-C were required inorder to avoid unlikely polynya radius predictions. It was found that the incorporation ofsurface melt into the polynya edge heat balance produced improbable radii on those dayswhen surface melt was predicted. This is attributed to the small value assigned to h mjn(the value is discussed later in this section). Significant surface melt results in aproportionally large reduction in this parameter. This significantly increased thesensitivity of the polynya radius to the warm weather conditions associated with surfacemelt, resulting in very large radius predictions. The predictions were much improvedwhen the MELTPE subroutine was removed. This is justified by the fact that the snow-free zone around the polynya is limited in extent (see Chapter 4). Therefore, if therewere to be a large increase in polynya radius, then some snow melt must occur. Sincethis snow is assumed not to exist, the incorporation of surface melt results in under-estimation of the thermal energy required to expand the polynya.136137• ice (observed)^— mli-c predictionso snow (observed) ^ mli predictions^ Harmon Lake Results (default Parameters) 1 0.3 -0.2 -0.0Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateA Temperature (Y.S.I data)— MLI-C predictionsMLI predictions^ Maximum & Minimum ProbableTemperature (weatherstation data)Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.11 Snow and Ice Thickness and Whole Lake Temperature Predctions:A Comparison of MLI and MLI-C Ouput for Menzies LakeTesting MLI-C at extreme values of Ar resulted in water temperatures which dippedbelow 1°C (see §5.2.3). Even on warm days with Hp > 0, (and therefore a small value of(qf)pe), Equation 2.25 did not yield a melting rate which would satisfy the polynyaexpansion criterion. This is due to the low value of qt associated with the very cold laketemperature. Hence Equations 4.26, 4.23 and 2.30 would be solved to calculate a new,presumably contracted, radius. With H p > 0 however, the solution gives an extremelylarge polynya radius. In response to this result, the model was changed so that when Hpis greater than zero, expansion is assumed and Equations 4.25 (which is based on (Ope,not Hp), 4.23 and 2.30 are solved for the new radius.Ice and Snow PredictionsFigure 5.11 includes ice and snow predictions for Harmon Lake, as well as results forMenzies Lake in the case of no artificial circulation. All three results are nearly identicalfor snow cover, and only minor differences in ice thickness are predicted over the lattertwo-thirds of the simulation. All differences are accounted for by q, which is lower for alower T. The artificially circulated lake, being the coolest, therefore results in thegreatest ice thicknesses.Since both the snow cover results and observations are almost identical to those atHarmon Lake, it is unlikely that artificial circulation has any direct impact on the snowcover. There are sound explanations, however, for the small differences observed. Theobserved ice thickness was significantly greater at Menzies than at Harmon Lakebetween late December and mid-January. This translates into a higher snow bearingcapacity at Menzies Lake which explains the greater snow thickness observed at the lakeon January 19th. Following the late January rain the snow thicknesses are similarbecause the snow bearing capacity at both lakes greatly exceed snow loading. The smallobserved differences may be, to a lesser degree, also due to dissimilar sheltering138characteristics and slight variances in the rates of surface melt due to the indirectinfluence of heat conduction through the ice.Unfortunately, the MLI-C does not predict the Menzies Lake ice cover satisfactorily. Itseems unlikely that errors were made in the observations given the simplicity of themeasurement, as well as the fact that significant differences from the Harmon Lakeobservations were noted on almost half of the field trip dates. As described in §5.2.4,however, improper sampling may have led to poor estimates of average lake icethickness. The greatest discrepancies between the predicted and observed ice thicknessesoccur in the first month of the simulation when the lake water is experiencing its highestrate of cooling. If the observations are reliable, MLI-C must significantly under-predictthe conduction of heat through the ice and snow cover. There are two possibleexplanations for this. First, the conductivity of ice and snow at Menzies Lake may bemuch greater than what is assumed. Second, the meteorological balance is improperlyrepresented. These possibilities are considered in §5.2.4. It should be noted that, asexplained in §5.2.3, although the calibration parameters Ct , hm i n, and Ar have asignificant impact on lake temperature and polynya radius, they have virtually no effecton the ice cover.With regard to the latter half of the winter, the explanations for any differences inobserved ice and snow thickness at Hannon Lake using MLI apply equally well here.Water Temperature PredictionsThe predicted water temperature at Menzies Lake is also given in Figure 5.11. The bestresults were achieved for the following calibration parameter values:Ct 1.1 x 10-3hm in 0.02 mAr 13 m139These values produced reasonable estimates of whole lake temperature, although the rateof warming in March is under-estimated. As described in §5.1, this may be due to lowradiation measurements resulting from weather station tilt or poor albedo representation.Furthermore, factors which are particular to Menzies Lake include under-estimation ofthe polynya radius and boundary effects . This latter factor may be of great importancegiven that the polynya edge was only two to three metres from the North shore of thelake, in very shallow water.The above value of Ct is slightly above the range given by Hamblin & Carmack (1990).In their study, C t is based on velocities measured 1 m below the ice-water interface. Inthis study, a relatively high value of c, is expected since the velocity distribution used isassumed to describe the water motion immediately below the interface. For a moreaccurate model, ct should be determined for the Air-o-lator® in a manner like that whichis described by these authors.Initially ro. was set to the average lake radius, but this produced much too high a rate oflake water cooling. The average lake radius would be expected to be an appropriatechoice provided that:a) Equation 4.23 is based on the true velocity field below the ice cover. It is likelythat the ice roughness accelerates the velocity decay, thereby reducing theturbulent heat transfer to below the rates predicted using extrapolated polynyasurface measurements.b) The polynya is located closer to the centre of the lake, well away from shore.Minimum ice thicknesses of approximately the same magnitude as the selected h min havebeen observed at both Harmon and Menzies Lakes. The value of h min should drop withreduced fetch length because the maximum wind speed and wave height development islimited, resulting in less turbulent energy available to melt ice. Although the circulatorcreates turbulence it is quickly dissipated as the radial jet spreads, and the wave140amplitudes produced are small. Therefore a relatively thin ice sheet is able to persistaround the polynya edge. Although no direct measurements were made, ice thicknessesof order 0.02 m were observed around the perimeter of the polynya. Directmeasurements made when both Menzies and Harmon Lakes were beginning to freezealso indicated that this minimum ice thickness is appropriate.Further fine-tuning of both h m in and Ar are unlikely to produce more reliable resultswithout making the model more sophisticated, or by providing better velocity field data.In order to improve the model, much attention must be given to the shape andcharacteristics of the ice edge at the polynya boundary and to the hydrodynamic andthermal boundary layers.If the reason for the under-prediction of ice thickness is adequately addressed, it isnonetheless possible that the parameter values given may produce reasonable results forother artificially circulated lakes (any changes to the MLI parameters examined in §5.1.4may result in minor alterations to the MLI-C calibration parameters). Indeed the fact thatthe weather station thermistor data (Figure 5.11) confirms that all of the trends intemperature are predicted by MLI-C (including the slight rise and subsequent fall in mid-February) is encouraging. It should be pointed out, however, that the value of Ar usedhere may be dependent on the lake morphometry given the proximity of the polynya tothe lake edge.The effect of varying the three calibration parameters is described in §5.2.3.Polynya RadiusAs suggested in Chapter 3, the errors involved in estimating the average polynya radiusgenerally exceed its variability over the observation period. The most detailedobservations were made on March 10th, 1992 when ten measurements of polynya radius141were made each with an estimated accuracy of ±0.3 m. Given the great local variabilityhowever, the average radius is considered only to be accurate to ±2 m. Ninemeasurements were made on February 27th, and only 2 on all the remaining field tripdates. The errors in all cases were estimated by scaling the estimated error of March 10thin the manner described in §5.1.2. An error of ±10 m is assigned to the final observation,on March 28th since no actual measurements were made and the polynya was very largeand irregular in shape. The location of the polynya edge relative to distinct features weremerely noted in this case for an off-site size estimate (the edge had almost reached thenorth shore of the lake by this time).MLI-C correctly predicts a narrow range in the polynya size as shown in Figure 5.12.There is a sharp diurnal trend due to the effect of solar heating in the first time step. Thetrend cannot be confirmed given the available data since all observations were made nearmid-day. Aside from this diurnal variation, there are, at times, some significant day today changes notable, especially during the January rainy period and near the end of thesimulation. The sampling was insufficiently frequent to confirm any rapid day to daychanges.MLI-C produces results which agree quite well with the observations but the lack of longterm variations which exceed the expected error does not provide an adequate test of themodel. Given the asymmetries in the jet and the potential impact of the lake boundarieson the hydrodynamics and the heat transfer processes, better results are not likely withoutincreasing both the sophistication of the model and the input data requirements.5.2.2 Components of the Heat BudgetThe surface heat budget components for the ice-covered portion of Menzies Lake is givenin Figure 5.13. The indirect effect of much cooler lake temperatures has virtually noeffect on any of the components in relation to the Harmon Lake simulation. If the ice14270 —60—E 50 —,t,.—cid 40 —a4>, 30—al. 20 —10 —— polynya radius (predicted)^ polynya radius (observed)Dec 13, 1991 Dec 30 Jan 19, 1992^Feb 8^Feb 28^Mar 19Datethickness were more accurately predicted, some minor differences may appear, the mostimportant of which may be a reduction in Q si due to the greater snow bearing capacity.Figure 5.12 Predicted and Observed Polynya RadiusThe heat fluxes to and from the lake water are also given in Figure 5.13. All of the fluxesshown have been adjusted to apply over the entire lake surface to facilitate directcomparisons. The only flux which is virtually identical to the Harmon Lake simulation isthe solar radiation which penetrates the ice cover, I. This is explained by the similarityof the snow and ice predictions. Initially the sediment heat transfer rate is almost thesame because the reduced thermal reserve in the sediments at Menzies Lake is offset by alower lake temperature. The sediment heat transfer rate, q sed, continues to increase,however, as the lake cools before a gradual reduction which begins in earnest at the endof February. The opposite trend is observed in q w . As the lake cools, there is a reducedthermal gradient between the lake water and the ice. As a result, the heat transfer rate143-20-151015Ched ^ Hp(IW(Qt )avgHeat Fluxes to and from Lake Water(left hand axis applies to all components except k„all fluxes are adjusted to entire lake area)02080100Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.13 MLI-C: Surface Heat Budget Components and Heat Fluxesto and from Lake Water at Menzies Lake144goes down until a steady recovery begins as the lake warms up at the end of February. Ingeneral both qsed and qw vary little compared with the remaining sources and sinks ofheat.Menzies Lake cooled steadily down to a minimum temperature of less than 1.4 0C. Fromthe start of the simulation until the end of January, when this minimum was reached, Qtaveraged about 7.6 W/m2 over the entire lake area. This rate was over three times therate of cooling across the polynya (averaged over the lake area) over the same period.From the start of February until the end of March, there was, on average, a slight nettransfer of heat into the lake across the polynya, while Qt increased to an average rate of8.6 W/m2 due to elevated water temperature. It is noteworthy that there is only about a 4W/m2 range in Qt over most of the simulation, mainly due to gradual changes in laketemperature. This is explained by the opposing effects of lake temperature and polynyaradius. Over the course of the winter no dramatic variations are seen in these keyvariables. A sharp diurnal trend in polynya size is predicted due to the lack of solarheating in the second time step (see §5.2.1). This trend is only mildly reflected in Q tbecause Ot is more dependent on Ar (which is assumed constant) than the actual polynyaarea (see Equation 4.24). There is no diurnal trend in Tw, so as r is reduced in the secondtime step, there is only a small corresponding increase in Q t. By mid-March, the polynyabegins its rapid expansion and there is an equally rapid increase in temperature. Qt isshown to steadily increase on average over this time. The area over which O  is assumedactive must increase for the given Or when the polynya area increases. This factor isoffset, however, by the inverse relationship between q t and r (Equation 2.30). So as rpincreases, the effective area of Qt also increases but q t decreases, resulting in little changein Qt . Since Qt, however, is proportional to Tw, it must follow the same trend as the laketemperature. It is likely that there is a dependency of Or on rp. That is, the distance overwhich the energy associated with the jet is dissipated decreases as the polynya radius145increases. If this is so, then MLI-C would tend to over-predict the heat loss due to Qtnear the end of the simulation when the polynya is large. This provides another possibleexplanation for the under-prediction in the rate of temperature increase in March (seeFigure 5.11). In addition it provides further impetus for future study of thehydrodynamics at the polynya edge.Hp varies strongly due to diurnal changes in the meteorological balance, but average heatloss over the first two months is only in the order of about 3 W/m2. In terms of heat lossfrom the entire lake, this means that H p is only slightly more important than simpleconduction of heat from quiescent water to ice. Furthermore, when warmer conditionsprevail, as early as the late January rainy period, the average daily heat exchange acrossthe polynya is close to zero. Following a brief cooling period in mid-February, whenaverage daily heat losses were only about 1 W/m2 , the direction of net heat transfer isreversed. Greatly increased solar energy penetrating through the ice-covered area morethan offsets the consistent turbulent heat transfer, Q t, to produce warming of the lakewater. The range in Hp is shown to increase dramatically at the end of March becauseboth solar heating and polynya area also increase in the same manner. Hence for eachtime step, the importance of Hp increases greatly with respect to the other sources andsinks of heat in the lake. In terms of average daily heat exchange, however, Hp remains amodest term compared with both I w and Qt .The components of the polynya heat budget are plotted in Figure 5.14. (They have notbeen adjusted to apply over the entire lake area.) As the lake surface is relatively warmcompared with the snow surface, net emission of long-wave radiation is much greaterfrom the polynya than the rest of the lake. Net  long-wave heat loss ranges from about 25to 100 W/m2 over the entire period. The greater temperature difference between thewater and the atmosphere also translates into greater evaporative and sensible heattransfers. In this case, however, the sensible component removes heat from the polynya146150 —----- Q (out) ^ Qe (out)Q1 (0110 • Q (01100—50100 —147Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19Date350 —,Ne 300 —at4.250 —200 —[L-1.3.) 150 —=100 —50 —Dec 13, 1991 Dec 30 Jan 19, 1992^Feb 8^Feb 28^Mar 19 13DateFigure 5.14 MLI-C: Surface Heat Budget ComponentsAcross Polynya at Menzies Lakesurface over most of the simulation (up to an occasional maximum of over 100 W/m 2).This is in spite of wind speeds which rarely exceed about a daily average of about 3 m/s,but are more commonly in the range of 1 to 2 m/s. Although the lake cools more inresponse to turbulent heat transfer to the ice cover, increased exposure to high windsduring cold periods may lead to more significant heat loss across the polynya surface.From Figure 5.14, it is also clear that heat loss due to free convection, Qfc, is much moreimportant than the evaporative loss as quantified by the standard aerodynamic formula.Maximum heat transfer rates due to this mechanism are in the range of 25 to 35 W/m2.Heat loss due to the cooling and melting of snow, Q sp, falling into the polynya is alsoimportant when the total daily accumulation is greater than about 5 cm. For this amountof precipitation, about 10 W/m 2 of heat is lost.By early February, solar heating begins to overcome all of the heat sinks describedabove. Following another brief cooling period later that month, only the long-wave heatloss term is significant, but is no match for the solar heating which ranges from about 100W/m2 at the beginning of March to over 400 W/m 2 at the end of the simulation.5.2.3 Impact of Varying MLI-C Calibration ParametersAny reasonable change in the calibration parameters has virtually no impact on the lakeice and snow cover. The impact of each on lake temperature and polynya radius,however, is substantial. The effect of changing the parameters on these two key variablesare studied in this section.Radial Extent of Turbulent Heat TransferThe impact of increasing Ar from 13 to 26 m is shown in Figure 5.15. The temperaturepredictions must first be examined in order to explain the effect on the polynya radius.By doubling Or, the area over which turbulent heat transfer more than doubles.14860 —— default parameters^ r... = rp + 26 m■ radius (observed)S.. • ..S..^.^0.• / •.••• "S i ' ss's .•: # `... ... :^*-.• .^...,050 —40 —30 —20 —10 —Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateA Temperature (observed)— default parameters^ r,,,, = rp + 26 m 1--..,...- . ...... ••0Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.15 Effect of Increased Radius of Turbulent Heat Transferon Polynya Radius and Whole Lake Temperature149Therefore, there is greater heat loss from the lake on each and every day of the simulationby virtue of increased Q t. The divergence in the temperature predictions is gradualbecause there is limited day to day variation in Q t for a given Or. The divergenceproceeds at a somewhat higher rate in December because of the relatively high thermalreserve from which heat can be drawn. This rate of divergence is not regained when thelake heats up again because el by then is relatively unimportant term compared withsolar heating. Hp is relatively unimportant throughout much of the simulation. Thereforethe slight reduction in polynya size (see below) results in less heat loss across thepolynya but this is insignificant compared with the increased Q t .In order to produce a clear picture of the effect on the polynya radius, only the first timestep predictions of this variable are shown (this gives an improved look at how well themodel agrees with observations since all of the observations were made during daylighthours). About a 15 m difference from the default case gradually builds up over thesimulation period. This is due primarily to the difference in lake temperature. For coolerwater, there is less heat available to keep the ice edge from encroaching on the polynyaarea. It is not Qt that is important to the polynya radius, but the value of qt right at the iceedge. As qt is proportional to Tw , the decrease in temperature results in a shrinkingpolynya.Minimum Ice ThicknessThe effect of increasing the minimum possible ice thickness from 2 to the 10 cm used inDYRESMI for much larger lakes is shown in Figure 5.16. Since this entire thicknessmust be melted in order for the polynya to increase, then q t must be larger for a largerhm in in order for expansion to occur. Since q t is proportional to r 1 , but not dependent onhm in , a balance will only be struck at a much smaller rp for the larger hm in . Hence thepredicted radii are less than half those for the default case. In addition, since a great deal150— default parameters^ hind, = 10 cm■ radius (observed)''''' - ----------------------- .. ............................................. .--- ...... /...0Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19Date4A Temperature (observed)— default parameters^ h.= 10 cm15150 -40 -30 -20 -10 -3 1... . .... . ... ...•••■•••• ........ .... ..... ...1Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19DateFigure 5.16 Effect of Increased Minimum Possible Ice Thicknesson Polynya Radius and Whole Lake Temperaturemore heat is required to melt 10 cm of ice than 2 cm, there is almost no day to dayvariation in the predicted radii. Over the long term, there is a very slight trend whichfollows that of the lake temperature. Since the thermal reserve is important indetermining the heat balance at the polynya edge, the size of the polynya decreases as thelake temperature goes down, and increases as it goes up.It is only with this strong variation in hmin that a significant change in water temperaturedue to the change in Hp occurs. With less than half the radius, the area over which Hp isactive is less than one quarter of the default area. Hence there is a significant reductionin the heat loss from the lake. Essentially none of this change can be attributed to Q tbecause it is only very weakly dependent on r p. The temperature curves stop diverging inlate January when the daily average flux approaches zero. With further cooling inFebruary, there is little further widening in the temperature gap because H p isunimportant relative to both Q t and Iw . The gap closes almost completely as the end ofthe simulation is approached, because the smaller polynya now restricts the heat transferback into the lake.Turbulent Heat Transfer CoefficientThe effect of reducing the turbulent heat transfer coefficient from 1.1 x 10 -3 to0.8 x 10-3 , the average value found by Hamblin & Carmack (1990), is shown in Figure5.17. The impact on the polynya radius is immediate since, from the start, there is onlyabout 70% of the turbulent heat transfer to the ice edge. This contraction of about 7 mreduces Hp due to the smaller area across which the heat transfer can take place. Therelatively warmer temperatures, however, are mainly due to the 70% reduction in Q twhich has been shown to be more important to the lake temperature than Hp. As warmertemperatures are produced, however, there is more heat available to melt the ice at thepolynya edge, and the gap between the polynya radius for C t = 1.1 x 10 -3 and Ct = 0.8 x152— default parameters.----. ct = 0.0008■ observed radius40 -30 -20 -10-070 -60 -50 -1.0^Dec 13, 19911^1Dec 30^Jan 19, 1992 Feb 8^Feb 28^Mar 19Date1Dec 13, 1991^Dec 30^Jan 19, 1992^Feb 8^Feb 28^Mar 19Date4.5 — A Temperature (observed)— default parameters---- ct = 0.00084.0 —3.5 —3.0 —153Figure 5.17 Effect of Decreased Turbulent Heat Transfer Coefficienton Polynya Radius and Whole Lake Temperature10 -3 is reduced. As Q t reduces in relative importance with increased solar heating, thetemperature no longer increases over the default case. With the temperatures risingtogether, the difference in ct values increases in importance again with regards to theheat balance at the polynya edge, and the radius must therefore decrease relative to thedefault case.5.2.4 Improving Ice Thickness PredictionsIn §5.2.1, it was shown that the MLI-C ice thickness predictions are somewhat less thansatisfactory due to the under-prediction of heat conduction through the cover during thecooling period of the simulation. Three possible explanations have been identified:a) under-estimation of ice and snow conductivityc) improper selection of ice thickness measurement locationb) improper representation of the meteorological balanceGiven the good results for the Harmon Lake simulation using the same snow and iceproperties, the first explanation is unlikely. One possibility, however, is the tendency toform white instead of black ice at Menzies lake due to the turbulence induced by thecirculator. Although, given the lower density of white ice, it is more likely that thiswould tend to decrease conduction, not increase it. The second explanation deservesmore consideration. In retrospect, ice measurements were often made too close to theweather station, which probably served as a heat sink. An increased flux of heat from thewater up through the weather station would have resulted in greater local ice thicknesses.This possibility, which brings the ice cover data into question, emphasizes the need for amore thorough sampling procedure, so that statistically reliable estimates can becalculated for each observation date.154The third explanation deserves detailed examination, as it provides a starting point forfurther research. To illustrate the potential impact of the polynya on the conduction ofheat through the ice cover in the rest of the lake, the analogy of parallel electricalresistors may be considered. Typically, the heat flux through a material which ischaracterized by variable conductivity over its surface is calculated using this concept. Inthis manner, reasonable results may be obtained for many heat transfer problems withoutresorting to a two-dimensional model. The best results are obtained if the variability iswell distributed over the surface of the material, such as in the case of equally spacedsteel bolts passing through an insulated wall. Although this is clearly not the case atMenzies Lake, the resistance concept could nonetheless be applied in future work to seeif the results improve. It is likely that increased melt from the underside of the ice due toturbulent heat transfer in the vicinity of the polynya would be offset by the tendency ofthe heat flux through the ice cover to increase in this region due to the effect of the lowrelatively thermal resistance associated with the polynya. Again this suggests the needfor more thorough ice thickness sampling in the future.The thermal resistance approach provides explanations which are specific to the poorrepresentation of the rate of ice growth in early winter, and the rate of melt late in theseason. Both of these rates are under-estimated using MLI-C, especially the former.Referring back to Figure 5.13, it can be seen that the importance of heat transfer acrossthe polynya, Hp is most significant during these two periods. Hence the impact of thepolynya on the ice thickness would be most important at these times. Although the rangein Hp is greatest in late winter, the average daily value per unit area of the polynya ismuch less significant. (Recall that the polynya at this time is large, and so Hp, per unitarea of the entire lake surface, as given in Figure 5.13, would also tend to increase inimportance.) It follows that, as the results indicate, the impact of the thermal resistanceof the polynya on ice growth should be most significant in early winter.155Chapter 6CONCLUSIONS AND RECOMMENDATIONS6.1 GENERAL REMARKSThe field observations at Menzies and Harmon Lakes coupled with the modeling resultsprovide conclusive evidence of the severe impact that artificial circulation has on thetemperature of a small mid-latitude lake. It has been confirmed that a gradual increase intemperature would likely take place at Menzies Lake if it were left in a natural state.Instead, significant cooling takes place, bringing the lake dangerously close to atemperature which is intolerable to the aquatic life. (As described in Chapter 1,deleterious effects on aquatic life may occur below about 1 0C). Over theuncharacteristically warm winter of 1991-92, the average temperature in Menzies Lakedipped to 1.40C. This fact alone suggests that the peril of harming the biota of lakessubjected to artificial circulation should not be under-estimated.MLI constitutes a first attempt to fine-tune an ice and snow cover model for small mid-latitude lakes left in a natural state. It has been shown in Chapter 5 that the extensions toDYRESMI provide significant improvements to ice, snow and temperature predictions.Although highly simplified, MLI-C also provides good predictions of water temperature.Given the questions which are arisen regarding the quality of the ice cover observations,it is possible that the model may also provide better ice cover predictions than theobservation suggests. Certainly it would appear that the trend in ice thickness is wellpredicted once the ice thickness approaches its maximum. Furthermore it should beemphasized that the ultimate goal of this project was to determine if there would be a156deleterious effect on the aquatic life in an artificially circulate lake due to excessively lowlake temperatures. To this end, the project has been successful. As described later in thischapter, however, more work is required before reliable predictions can be made basedon a hypothetical meteorological data set.An important feature of MLI and MLI-C is the reduced computing power requirementcompared with DYRESMI. It has been shown that good results are achieved for smalllakes even when the internal hydrodynamics are ignored. Hence the short time stepsemployed by DYRESMI can be avoided, resulting in significant savings in computingtime. For this project MLI was run on a 80386 (25 MHz) PC. Less than 10 seconds wererequired for both MLI and MLI-C to produce results for the 105 days of simulation.In this chapter, the results of both the Menzies and Harmon Lakes simulations will besummarized and conclusions drawn regarding the DYRESMI extensions, the actual heattransfer processes, and the parameter values.6.2 MLI6.2.1 Summary of ResultsMLI provides good predictions of ice thickness and water temperature, and reasonablepredictions of snow thickness. It has been found that solar heating is the main cause oflake warming over the ice-covered period but that sediment heat transfer may at leastcounter conductive heat loss to the ice cover to prevent a reduction in temperature inearly winter. In addition, rain events may lead to a significant increase in lake warmingdue to reduced albedo. Ice thicknesses may be reduced as a result of the heat associatedwith rain, but it may eventually cause an increase in ice thickness if the insulationprovided by the snow cover is lost and cold conditions follow. The impact of snow-ice157formation is similar in that mass is added to the total ice thickness, but the associatedlatent heat may result in melting at the ice-water interface.The following extensions to DYRESMI are most important in achieving the results asdescribed in Chapter 5:a) rain meltb) solar attenuation through snow-icec) a snow density function appropriate for mid-latitudesd) sediment heat transfere) variable albedoIt is clear that the water which floods the snowpack when the bearing capacity of ice isexceeded constitutes a source of heat within the cover. As described in §5.1.4, this heat,as represented by Qsi, may be over-estimated since further snow densification likelyoccurs prior to the snow-ice formation. In addition, capillary rise may contribute to thesnow-ice thickness.The latent heat associated with rainfall was needed to produce the order of observedsurface melting. Although the results are reasonable, the mechanism is unlikely. Asdescribed in Chapter 2, condensation melt should be the main driving force behindsnowmelt for the given amounts of precipitation. It is suspected that the use of dailyaverage meteorological data leads to an under-representation of condensation melt.Results indicate that the surface melt due to rainfall is over-estimated. Accuratepredictions of snowmelt is essential, especially when the snow cover is thin. Even slightinaccuracies may lead to much more significant inaccuracies in ice and temperaturepredictions because of the significant drop in albedo (and to a lesser extent, the increasein conduction) associated with ice exposure when the snow is completely melted. Thisexplains why the unlikely surge in lake temperature due to solar heating following theincorrect prediction of complete snowmelt in mid-February.158The correct combination of parameter values which would provide more accuratepredictions is not completely evident from this study. Much work is required to quantifynarrow ranges of possible values for mid-latitude locations. The relationships betweensnow density, conductivity, solar attenuation and albedo are of the greatest importance.A physically based model of snow densification is also needed.From §5.1.4 however, it is hypothesized that the following combination of parameterchanges are most likely to produce a more accurate representation of the heat budget:a) Increased K s : Figure 5.7 illustrates the significant impact of snow conductivity onice thickness. A similar impact was only achieved by increasing maximum snowdensities, but this leads to less likely snow cover predictions. The initial densitycannot be increased by much to improve this latter result because of the restrictedload bearing capacity of ice.b) Increased L.: Figure 5.5 shows the impact of increasing the visible radiationattenuation coefficient of snow-ice from 3.75 m -1 to 6.0 m-1 . The result is a morelikely rate of ice melt in March due to the absorption of more heat. This extra heatdoes not penetrate to the water, thereby producing improved temperaturepredictions.c) Reduced a in March: The impact of reduced albedo in March is suggested byFigure 5.9. More than a simple parameter change is required here, since theresults are very much dependent on whether or not a thin snow layer is predicted.From the few observations, made, it would appear that snow free conditionspersist for a longer period of time than what is suggest by the predictions. Theeffect of heavy snowmobile traffic is suspected a being important in reducing thealbedo to below those values predicted by MU. The result would be a faster rateof ice melt and temperature increase.d) reduced Kj: It is suspected that the conductivity of the sediments is close to thatof water as suggested by most of the studies available in the literature. Included inthis majority is a study of gelatinous, organic sediments which are qualitativelysimilar to the sediments of Menzies and Harmon Lakes (see Chapter 2). It isexpected that this change, in combination with the increase in Ae would159appropriately offset the lake temperature increase due to a reduced albedo inMarch.6.2.2 Recommendations for Future ResearchFurther research should concentrate on confirming that the parameter changes describedin §6.2.1 are appropriate. This would involve including direct measurements of snowdensity, conductivity, solar attenuation, albedo, and sediment conductivity, in addition torepeating the data collection and modeling process over a full winter at Harmon Lake. Inorder to show that there are no other parameters of importance to the results, additionalsteps will be required to prove the accuracy of the model. These include a betterunderstanding of Qr, Qsi, qw and qsed as described below. Furthermore, the observationerrors must be reduced. For the ice and snow cover, many more observations should bemade over the entire lake area. With regard to lake temperature, more measurementsover the thermal boundary layer at the ice-water interface may improve estimates (thiswould also help to better quantify q w). In addition, instrumentation error may be reducedslightly by measuring the lake temperature over several profiles. If field measurementsare restricted by time, however, the emphasis should be on making several ice and snowmeasurements. Temperature estimates would be improved to a greater extent by a highlyaccurate bathymetric survey and by better quantifying boundary effects.More research is required to assess the importance of convective stirring near thetemperature of maximum density. This condition probably occurred in March to result inincreased rates of ice melt due to turbulent (convective) heat transfer from the water tothe ice. It is possible that this mechanism may entirely explain the difference betweenthe observed and predicted rates of ice melt. If this were the case, then the changes to Xeand a described in §6.2.1 would be less likely.160Of great importance is the need to develop a more accurate rain melt model withoutresorting to sub-daily meteorological data. The apparent importance of condensationmelt, which is not reflected in the MILI model, is iterated in Chapters 2, 4 and 5. Relatedissues include the impact on both the snow albedo and density. Furthermore, the role oflatent heat due to the freezing of rain water as it approached the ice layer should be betterunderstood.The nature of the thermal boundary layers both at the ice-water interface and the lakebottom should be explored further. Of particular importance are the thermal gradients atboth interfaces. In a global sense, the heat flux from the water to the ice is treatedindependently of the heat flux through the ice and snow cover in the MLI formulation. Ifthis is so, further research may produce a model which predicts the thermal gradientbased on the water temperature outside the boundary layer alone if the water is still. Anyturbulent motion in the water below the ice, such as that created by convective mixingshould also be considered. With regard to sediment heat transfer, a relationship that thethermal gradient may have with the annual average lake temperature and time sinceturnover would be useful for modeling purposes.Qsi should be better quantified. This could be accomplished by making several fieldmeasurements of snow density, flooding depth, and snow, ice, and snow-ice thicknessesover the period in which snow falls, flooding occurs and snow-ice formation takes place.The spatial variation of these quantities should also be assessed by sampling over theentire lake area. It is likely that bending of the ice cover under the weight of snow resultsin significantly more flooding and snow-ice production near the centre of the lake, asobservations have indicated. Establishing the importance of capillary rise of waterthrough snow would also be useful.161Finally, the cloud cover estimate algorithm should be tested thoroughly by comparingactual long-wave radiation measurements to model predictions. Given the importance ofcloud cover to both the long-wave radiation balance and albedo predictions, the need tobetter quantify this variable is substantial.6.3 MLI-C6.3.1 Summary of MLI-C ResultsThe MLI-C model provides good predictions of temperature at Menzies Lake throughoutthe entire simulation period. The deviations from the observations can be attributed toinappropriate choices of parameters which are independent of the artificial circulationprocess. This issue is dealt with in §6.2. However, the reliability of the results arebrought into question by the apparently poor ice predictions in early winter. Throughoutmost of the period, the predicted trend in ice thickness is good, and the deviations maysimply be explained in terms of the MLI parameters. Further field work is required inorder to determine if the lack of agreement in early winter was due to an impropersampling procedure as described in §5.2.4.The following values for the calibration parameters have been found to give the bestpredictions:hmm 0.02 mAr 13 mCt 1.1 x 10-3Although determined through calibration, the values of hmin and C t make good physicalsense as described in §5.2.1. Only Ar has no physical meaning except that it is the radialdistance over which the turbulent heat transfer for frictionless flow is equal to that162produced by a velocity field which is dissipated by friction at the ice-water interface overthe entire ice-covered area. According to visual observations through ice holes, thevelocities became negligible below the ice at much less than the average lake radius awayfrom the polynya edge.In spite of an apparently low value of Or, the turbulent heat transfer from the water to theice, Qt, at an average of about 7.6 W/m2 , is over three times as important as the heattransfer across the polynya, H p, over the lake cooling period. After the lake had reachedits minimum temperature, there is, on average, a slight net transfer of heat across thepolynya until the end of the simulation. Meanwhile, Q t increased in importance in thelatter half of the simulation, with a lake average heat transfer rate of about 8.6 W/m 2 .Furthermore, in spite of such a small value of h min , there is little day to day variation inrp until late March when ice-free conditions are rapidly produced. In this regard, themodel simulates field conditions well.The results described in Chapter 5 indicate that there is a delicate heat flux balance whichresults in lake temperatures which are just adequate in supporting aquatic life. Forexample, the heat which is provided to the lake by the sediments could, over coldwinters, often constitute the difference between fish survival and fish mortality.Similarly, if overcast skies prevailed, the lake may be deprived of the solar energy itneeds (fortunately, extremely cold weather is often associated with clear skies).Extremely windy conditions could also increase the importance of Hp through substantialgains in sensible heat loss. Furthermore, although snow cover does insulate the lake, itsmain impacts are to reduce the growth in ice thickness and to prevent solar radiation frompenetrating into the water. Hence a winter characterized by consistently thick snowcover may deprive an artificially circulated lake of the solar heat it needs over the winter.163Finally, if a large lake, characterized by a larger value of hm in , were equipped with aproportionately large circulation system, a proportionately large polynya would not beexpected. Similarly, a polynya which is larger than expected may develop at a very smalland highly sheltered lake.6.3.2 Recommendations for Artificial Circulator DesignFrom the results of this study, it is not the polynya size which must be controlled to avoidover-cooling a lake, but the turbulence beneath the ice. If aeration is improved merely byincreasing the polynya size (without increasing the pumping rate) then a lake whichpersistently experiences winterkill due to inadequate aeration could be rehabilitatedwithout increasing O. Unfortunately, for the given system, a larger polynya could onlybe created by increasing the radial jet velocity. As both q t and Qt are proportional to thisvelocity, then the polynya cannot be increased in size without increasing Q t .The cost associated with the option of virtually eliminating Q t by installing large bafflesin a suitable geometry at some distance away from the circulator would, in possibly allcases, be prohibitive. A more viable option may involve the installation of two or moresmaller units in such a way that an adequate polynya surface area (and adequate totalpumping rate if required) is maintained. In this manner, Qt is reduced since, although thesame surface velocities must exist along the perimeter of the polynya for a given polynyaradius and any system configuration, the use of several units should mean that the rate atwhich this velocity decreases beyond the polynya edge should be higher. This isexplained by the inverse relationship between velocity and the radial distance from eachindividual circulator. The net effect is a reduction in Or. Figure 5.15 suggests that boththe average lake temperature and the average polynya radius would increase. Such astrategy may therefore actually improve re-aeration over the use of one large unit, whilediminishing the impact on lake temperature.164A third option involves the use of a distributed air bubble diffuser system. In mostsystems described in the literature (see, for example, Ashton, 1979), a point sourcebubble diffuser creates a rising plume which entrains water, strikes the water surfacewhere it is redirected in the form of a radial jet. Such systems would create a similarvelocity decay function as the Air-o-lator®. An alternative would be to install adistributed source diffuser which would create a grid of smaller plumes impinging at thesurface. The principal regarding reduced Qt is the same as that described above. Thisalternative may be less costly, since only one air compressor would be required.A final suggestion involves the use of a pump which could produce a relatively fine sprayof water. In this case the jet would be projected high into the air, such that trajectory ofthe falling drops of water are nearly vertical. The larger the required polynya, the greaterthe maximum height of the fountain would have to be. This system would almostcompletely avoid any significant turbulent heat transfer to the ice. The design wouldinvolve calculation of the minimum rate of constant rainfall required to maintain ice-freeconditions, and the amount of cooling that the drops of water would experience. If thedrops were to cool substantially, or even freeze, the result could be very high rates of lakecooling.6.3.3 Recommendations for Future ResearchThe first step in improving MLI-C must be to re-evaluate the MLI parameters asdescribed in §6.2.2. Once MLI is fine-tuned, attention should be focused on thecharacteristics of the polynya edge and the hydrodynamic boundary layer which becomesestablished at the ice-water interface. Once the hydro/thermodynamics are betterunderstood, it may be possible to abandon the hmin approach in favour of a moretheoretical boundary layer heat transfer model in which the growth of ice from thepolynya edge out the full ice thickness is also predicted. If this could be done, then there165would also be no further need for the Ar calibration constant. The remaining parameter,Ct, would best be assessed in the field specifically for the Air-o-lator®. Alternatively, C tcould be determined with better confidence by means of calibration once thehydrodynamic model is established.The model described above would ideally be verified in the field. The use of under-icedrogues such as those described by Hamblin (1990) would allow for measurement of thevelocity field below the ice. In addition, large sections of the ice could be cut away for adirect measurement of the ice profile in the vicinity of the polynya. Furthermore,detailed measurements of the thermal boundary would provide critical information on theturbulent heat transfer rates.More thorough ice thickness observations is also recommended for any future studiesinvolving artificially circulated lakes. Several measurements would be useful both formeaningful values of average thickness, and for providing more insight on spatialvariations which may be dependent on the circulator. This may allow for an indirectevaluation of the velocity decay below the ice, and O.Finally the impact the polynya has on the heat flux through the ice covering theremainder of the lake should be well established and quantified. A first attempt atsolving this problem could involve the thermal resistance approach described in Chapter5.166REFERENCESAshley, K.I., 1987, Artificial Circulation in B.C.. Review and Evaluation, Province ofBritish Columbia, Ministry of the Environment and Parks, Fisheries TechnicalCircular no. 78, 34 pp.Ashley, K.I., 1991, personal communication, Fish and Wildlife Branch, B.C. Ministry ofthe Environment, Vancouver, B.C.Ashton, G.D., 1974, Air Bubbler Systems to Suppress Ice, Cold Regions ResearchEngineering Laboratory, Special Report 210, AD-A008867.Ashton, G.D., 1979, Point Source Bubbler Systems to Suppress Ice, Cold RegionsScience and Technology, no. 1, pp. 93-100.Ashton, G.D., 1986, ed., River and Lake Ice Engineering, Water Resources Publications,Colorado.Bandow, F., 1986, Evaluation of Winter Lake Aeration Techniques in Minnesota, Pub.No. 386, Minnesota Dept. of Natural Resources, Division of Fish and Wildlife, St.Paul, Minnesota.Bauer J., Martin, S., 1983, A Model of Grease Ice Growth in Small Leads, J. Geophys.Res., 88(C5), pp. 2917-2925.Bergen, J.D., Hutchinson, B.A., McMillan, R.T., Ozment, A.D., Gottfried, G.J., 1983,Short-wave Reflectivity of Recent Snow, Journal of Applied Meteorology, 22(2),pp.193-200.Birge, E.A., Juday, C., March, H.W., 1927, The Temperature of the Bottom Deposits ofLake Mendota: a Chapter in the Heat Exchanges of the Lake, Trans. Wisc. Acad. Sci.Arts Lett. 23, pp. 187-231.Biier, K.W., 1980, The Terrestrial Solar Spectrum, Solar Energy Technology Handbook,Part A (W.C. Dickenson & P.N. Cheremisinoff [eds.]), Dekker, pp. 65-87.Bohren, C.F., and Beschta, 1979, Snowpack Albedo and Snow Density, Cold RegionsScience and Technology, 1(1), pp.47-50.Chen, C.W., Orlob, G.T., 1975, Ecological Simulation for Aquatic Environments, SystemAnalysis and Simulations in Ecology, B. Patton, ed., Vol. 3, Academic Press, NewYork.Choudhury, B.J., 1982, Spectral Albedos of Mid-latitude Snowpacks, Cold RegionsScience and Technology, 6(2), pp.123-139.Den Hartog, G., Smith, S.D., Anderson, R.J., Topham, D.R., Perkin, R.G., 1983, AnInvestigation of a Polynya in the Canadian Archepelago, 3. Surface Heat Flux,Journal of Geophysical Research, 88(C5), pp. 2911-2916.167Eranti, E., Leppanen, E., Penttinen, M., 1983, Ice Control in Finnish Harbours,Proceedings, 7th International Conference on Port and Ocean Engineering UnderArctic Conditions, Vol. 1, Espoo, Finland, pp. 370-380.Gardon, R., Akrifirat, J.C., 1966, Heat Transfer Characteristics of Two-Dimensional AirJets, Trans. Am. Soc. Mech. Eng., series C, J. Heat Transfer, 88, pp. 101-108.Geiger, R., 1965, The Climate Near the Ground, Cambridge, Harvard, University Press,611 pp.Gilpin, R.R., Hirata, T., Cheng, K.C., 1980, Wave Formation and Heat Transfer at an Ice-Water Interface in the Presence of a Turbulent Flow, Journal of Fluid Mechanics, pp.619-640.Gosink, J.P., 1987, Northern Lake and Reservoir Modeling, Cold Regions Science andTechnology, no. 13, pp. 281-300.Grant, L.O., Rhea, J.0., 1973, Elevation and Meteorological Controls on the Density ofNew Snow, Advanced Concepts and Techniques in the Study of Snow & IceResources, U.S. Committee of the International Hydrological Decade, pp. 169-181.Gray, D.M., 1970, Handbook on the Principles of Hydrology, Canadian NationalCommittee for Hydrological Decade, Ottawa.Gray, D.M., and Male, D.H., 1981, Handbook of Snow: Principles. Processes,Management and Use, Pergamon Press.Grenfell, T.C., Maykut G.A., 1977, Optical Properties of Ice and Snow in the ArcticBasin, Journal of Glaciology, 18, pp. 445-463.Grenfell, T.C., Perovich, D.K., Ogren, J.A., 1981, Spectral Albedos of an AlpineSnowpack, Cold Regions Science and Technology 4(2), pp. 121-127.Hamblin, P.F., 1990, Yukon River Headwater Lakes Study. 1983 and 1985: Observationsand Analysi s, Environment Canada, Scientific Series No. 175, Lakes ResearchBranch, Canadian Centre for Inland Waters.Hamblin, P.F., 1992, personal communication, Lakes Research Branch, Canadian Centrefor Inland Waters (National Water Research Institute), Burlington, Ont.Hamblin, P.F., Carmack, E.C., 1990, On the Rate of Heat Transfer Between a Lake andan Ice Sheet, Cold Regions Science and Technology, 18, pp. 173-182.Harr, R.D., 1981, Some Characteristics and Consequences of Snowmelt During Rainfallin Western Oregon, Journal of Hydrology (53), pp. 277-304.Henderson-Sellers, B., 1984, Engineering Limnology, Pitman Publishing Ltd., 356 pp.Hill, J.M., Kucera, A., 1983, Freezing a Saturated Liquid Inside a Sphere, Int. J. HeatMass Transfer, 26, pp. 1631-1638.168Hiratia, T., Gilpin, R.R., Cheng, K.C., 1979, The Steady State Ice Layer Profile on aConstant Temperature Plate in a Forced Convection Flow. II. Turbulent Regime,International Journal of Heat and Mass Transfer, 22, pp. 1435-1443.Holman, J.P., ed., 1990, Heat Transfer, 7th ed., McGraw-Hill Publishing.Hutchinson, G.E., 1957, A Treatise on Limnology. Vol. 1: Geography. Physics andChemistry, Wiley, New York.Imberger, J., & Patterson, J.C., 1981, A Dynamic Reservoir Simulation Model-DYRESM:5, in Transport Models for Inland & Coastal Waters: Proceedings of asymposium on Predictive Ability ( H.B. Fischer, ed.), Academic Press, New York.Jirka, G.H., Katavola, D.S., 1980, Supercritical Withdrawal from Two-Layered FluidSystems. Part 2: Three-Dimensional Flow into Round Intake, Journal of HydraulicResearch 17, No.1, pp. 53-62.Kirk, J.T.O., 1983, Light and Photosynthesis in Aquatic Ecosystems, Cambridge.Kraus, E.B., 1972, Atmosphere-Ocean Interaction, Clarendon Press, Oxford, 275 pp.Lackey, R.T., Holmes, D.W., 1972, Evaluation of Two Methods of Aeration to PreventWinterkill, Prog. Fish Cult., 34(3), pp. 175-178.Likens, G.E., Johnson, N.M., 1969, Measurement and Analysis of the Annual HeatBudget for the Sediments in Two Wisconsin Lakes, Limnology & Oceanography,114(1), pp. 115-135.Maykut, G.N., Untersteiner, N., 1971, Some Results from a Time Dependent.Thermodynamic Model of Sea Ice, J. Geophys. Res. (83), pp. 1550-1575.McKay, G.A., 1968, Problems of Measuring and Evaluating Snowcover, SnowHydrology: Proceedings of the Canadian National Committee of the InternationalHydrological DecadeMortimer, C.H., 1941-42, The Exchange of Dissolved Substances Between Mud andWater in Lakes, Journal of Ecology, 29, pp. 280-329; 30, pp. 147-201.Mortimer, C.H., 1971, Chemical Exchanges in the Great Lakes - Speculations onProbable Regulatory Mechanisms, Limnology & Oceanography, 16(2), pp. 387-404.Mortimer, C.H., and Mackareth, 1958, Convection and its Consequences in Ice-coveredLakes, Verh. internat. Ver. Limnol., pp. 923-932.National Research Council of Canada, 1985, National Building Code of Canada, 1985,Supplement, Report No. 23178.Neumann, J., 1953, Energy Balance and Evaporation from Sweet-Water Lakes of theJordan Rift, Bull. Res. Coun. Israel, 2, pp. 337-357.O'Neill, K., Ashton, G.D., 1981, Bottom Heat Transfer to Water Bodies in Winter, U.S.Army Corps of Engineers, Special Report 81-18.169Patterson J.C., & Hamblin, P.F., 1988, Thermal Simulation of Lakes with Winter IceCover, Limnol. Oceanogr. 33(3): 328-338.Petzold, D.E., 1977. An Estimation Technique for Snow Surface Albedo, ClimatologicalBulletin, McGill University, Dept. of Geography, 21, pp.1-11.Rajaratnam, N., 1976, Developments in Water Science: Turbulent Jets, ElsevierScientific Publishing Co., Amsterdam, 304 pp.Raphael, J.M., 1962, Predictions of Temperature in Rivers and Reservoirs, Journal of thePower Division, A.S.C.E., 88(P02), Proc. paper 3200, pp. 157-182.Robinson, D.A., and Kukla, G., 1984, Albedo of a Dissipating Snow Cover, Journal ofClimate and Applied Meteorology, 23(12), pp.1626-1634.Ryan, P.J., Harleman, D.R.F., Stolzenbach, K.D., 1974, Surface Heat Loss from CoolingPonds, Water Resources Research, 10(5), pp. 930-938.Semter, A.J., 1976, A model for the Thermodynamic Growth of Sea Ice in Numerical Investigations of Climate, J. Phys. Oceanogr., 6, pp. 379-389.Stauffer, R.E., 1987, A Comparative Analysis of Iron, Manganese, Silica, Phosphorusand Sulpher in the Hypolimnia of Calcareous Lakes. Water Research 9(21), pp.1009-1022.Stigebrandt, A., 1978, Dynamics of an Ice-Covered Lake with Through Flow, NordicHydrology, 9, pp. 219-244.Strickland, N., 1982, Factors Affecting Temporal Variations in the Albedo of SnowCover, Western Snow Conference, Proceedings, pp. 186-187.Tennessee Valley Authority (TVA), 1972, Heat and Mass Transfer Between a WaterSurface and the Atmosphere, Division of Water Control Planning, Report No. 0-6803, Water Resources Research Laboratory Report No. 14, Norris, Tennessee.United States Army Corps of Engineers (U.S.A.C.E.), 1956 , Snow Hydrology, SummaryReport of the Snow Investigations, North Pacific Division, Portland, Oregon, 437 pp.United States Army Corps of Engineers (U.S.A.C.E.), 1986, CE-QUAL-R1: ANumerical One-Dimensional Model of Reservoir Water Quality (User's Manual),Instructional Report E-82-1, Environmental Laboratory, Vicsburg, Mississippi.Urban, T.C., Diment, W.H., 1988, Thermal Regime of the Shallow Sediments of ClearLake, California, Late Quaternary Climate, Tectonism, Sedimentation in Clear Lake,Northern California Coasts, Geol. Soc. of America, Boulder, CO., pp. 207-221.Welch, H. E., and Bergmann, M.., 1985, Water Circulation in Small Arctic Lakes inWinter, Can. J. Fish. Aquat. Sci. 42, pp. 506-520.Wetzel, R.G., 1983, Limnology, 2nd Edition, Sanders College Publishing.Wiscomb, W.J., and Warren, S.G., 1980, Model for the Spectral Albedo of Snow, Journalof the Atmospheric Sciences, 37(12), pp. 2712-2733.170Witze, P.O., Dwyer, H.A., 1976, The Turbulent Radial Jet,  J. Fluid Mech., 75(3), pp.401-417.Wood, I.R., 1978, Selective Withdrawal from Two-Layered Fluid, Journal of theHydraulics Division, ASCE, pp. 1647-1659.171cos(hss f) - cos(0) cos(8)sinNss) - sin(0) sin(8)APPENDIX ATheoretical Solar Radiation Under Clear SkiesThe Tennessee Valley Authority (1972) has published a comprehensive report coveringempirical formulae which permit calculations of heat and mass transfer between at watersurface and the atmosphere. Included in the report are the formulae required to calculatethe total amount of solar radiation, SWCS, which could reach the earth's surface underclear skies at a given location and at any time of the year. SWCS requires estimates ofatmospheric attenuation which is controlled by the optical air mass, precipitable watercontent and atmospheric dust. In addition, since the total daily radiation includes bothdirect and indirect energy, an estimate of the local albedo is also required. First,however, the solar radiation at the top of the atmosphere must be calculated as follows:Qs - 24-Ic [ hss sin(0) sin(6) + cos(0) cos(6) sin(hss) ]nrd2where, Ic = solar constant= 4871 kJ / m2 / hrrd = ratio of actual to mean distance of the earth from the sun2n= 1+ 0.017 cos(36-5  (186-jd))jd = Julian day0 = latitude of geographical location6 = solar declinationhssThe hour angle at sunset, h ss, is given by:172= hour angle at sunsetwhere, lliss = solar altitude at sunset (-.=, 15° at Menzies and Harmon Lakes due totopographic obstructions; assumed to be equal to the solar altitude atsunrise)If cos(hss ') = 0, then hss = 71/2If cos(hss ') < 0, then h ss = TC - cos-1 1cos(hss ')IIf cos(hss ') < 0, then h ss = cos -11cos(hss1 )1An accurate formulation for the sun's declination, 8, is given as:sin(8) = sin(Smax) sin(a)2it^where, 6^= s66 [279.9348 + 360 d + 1.914827 sin(d) -0.079525 cos(d) + 0.01993827csin(d) - 0.001620 cos(2d)]d^= angular fraction of a year2ir^,,,, ,— 365.242 u" - I)7C8max = 23.445 180Two transmission coefficients are required for the calculation of SWCS:a' = exp[-(0.465 + 0.134 w) (0.129 + 0.171 exp(-0.88 mp)) mp], anda" = exp[-(0.465 + 0.134 w) (0.179 + 0.421 exp(-0.721 mp)) mp]where, a' = mean atmospheric transmission coefficient after scatteringa" = mean atmospheric transmission coefficient after scattering and absorptionw = precipitable water content in the atmosphere= 0.85 exp(0.110 + 0.0614 Td)Td = mean dew point temperature (°C)265.5 (log(svpd) - 0.7858) — 9.5 + 0.7858 - log(svpd)mp = optical air mass[(288 - 0.0065 z)/288]5256 _sin(W) + 0.15(tv 180/n + 3.885) -1 .253z = local elevation173tV = solar altitude (the altitude at solar noon is used in MLI)= sin(0) sin(8) + cos(0) cos(8)Although the equation for the precipitable water content, w, has been calibrated for amonthly average Td (TVA, 1972), daily averages are used here since they are available.SWCS may be computed as:SWCS – a" + 0.5(1 - a' - dust) — 1- 0.5 Rg (1 - a' + dust)Where, dust^= depletion coefficient of the direct solar beam by dust scatteringRg^= local reflectivity of the ground in the vicinity of the siteA value of 0.13 for dust was selected. Few data are available but the range given in TVA(1972) for the winter season is 0.06 to 0.13, depending on the air mass, m p . It was foundthat a value chosen from the high end of the range produced the most reasonable cloudcover predictions. The local reflectivity is dependent on the surrounding forest, the smallarea of snow-covered fields adjacent to the lake, and the lake surface itself. The lakesurface must only have an indirect effect on the reflected solar radiation measured at thepyranometer since it is entirely below the elevation of the instrument. The tall treessurrounding the lake, however, have a more direct impact. Rg was estimated by adding70% of the forest reflectivity to 30% of the lake surface albedo. An average value of0.07 for forest cover was selected with reference to TVA (1972). This is much less thanthe lake surface albedo with either snow or snow-ice cover. Since the lake surface albedomay vary substantially from day to day, then SWCS will not necessarily be expected tosteadily increase each day as time progresses beyond the winter solstice.174


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items