UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Heat budget for an oil sands pit lake Chang, Sarah 2020

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

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata


24-ubc_2020_november_chang_sarah.pdf [ 4.96MB ]
JSON: 24-1.0392976.json
JSON-LD: 24-1.0392976-ld.json
RDF/XML (Pretty): 24-1.0392976-rdf.xml
RDF/JSON: 24-1.0392976-rdf.json
Turtle: 24-1.0392976-turtle.txt
N-Triples: 24-1.0392976-rdf-ntriples.txt
Original Record: 24-1.0392976-source.json
Full Text

Full Text

Heat Budget for an Oil Sands Pit Lake  by  Sarah Chang  B.A.Sc, University of Waterloo, 2013  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Civil Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  August 2020  © Sarah Chang, 2020    ii  The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis entitled:  Heat Budget for an Oil Sands Pit Lake  submitted by Sarah Chang in partial fulfillment of the requirements for the degree of Master of Applied Science in Civil Engineering  Examining Committee: Gregory Lawrence, Professor, Department of Civil Engineering, UBC Co-supervisor Edmund Tedford, Research Associate, Department of Civil Engineering, UBC Co-supervisor  Roger Pieters, Research Associate, Department of Earth, Ocean and Atmospheric Sciences, UBC Additional Examiner iii  Abstract At Syncrude’s Base Mine Lake, a hydrocarbon sheen arising from a tailings deposit below the water column modifies the heat exchange between the lake and its surroundings. A heat budget was implemented to determine how heat fluxes are affected, and to assess resultant impacts to lake dynamics. The period of record spans May 2015 to April 2016, separated into open-water and under-ice periods.  In the open-water period, individual heat fluxes were estimated using an existing parameterized model of the physical processes, the Air-Sea Toolbox, and validated using eddy covariance measurements and other methods to increase confidence in the results. The hydrocarbon sheen introduces resistance to evaporation at the water surface. This was modelled as a reduction in the relative humidity at the base of the air-water boundary layer. Latent heat flux due to evaporation dropped as a result by 10 % overall, or 20 % in low winds; in windy conditions, surface waves form, temporarily dispersing the sheen.  The hydrocarbon sheen fundamentally impacts many of the physical processes of the lake. Reduced latent heat flux means a directly proportional drop in evaporation, which must be reflected in the water budget. As evaporation is a key cooling mechanism, elevated surface temperatures also result. Increasing the surface temperature also impacts other fluxes, partially offsetting this effect. Nevertheless, the lake remains warmer than would otherwise be expected, with potential consequences to biological activity. Reduced cooling also means a decrease in the buoyancy flux and in buoyancy-induced mixing.  In winter, individual fluxes are again estimated using physical models, and a rough heat budget is developed. The hydrocarbon sheen does not play a major role in winter. However, the brackish water and heat exchange with the FFT deposit below the water cap introduce interesting under-ice dynamics. Rapid warming and occasional buoyancy-driven mixing are observed at the lake bottom, while at the surface, the majority of the fluxes relate to the growth and loss of the ice cap. These observations contribute to the understanding of heat fluxes and their impacts on water dynamics in tailings-affected systems. iv  Lay Summary Environmentally conscious mining practices are essential for balancing a healthy environment with society’s economic needs. Syncrude’s Base Mine Lake is a unique opportunity for researching the effectiveness of current techniques and policies. At Base Mine Lake, mine waste has been deposited and isolated from the surface by a water cap. This has resulted in a thin oil film on the water surface. The film forms a barrier between the atmosphere and water, changing the behaviour and characteristics of the lake. This study modelled the heat exchange between the lake and its surroundings, and found a reduction in evaporation at the water surface, leading to changes in the processes governing lake behaviour. This research estimates the strength and timing of key heat sources and sinks in a tailings-affected lake over a full year, offering insight into whether water-capped mine waste deposits are sustainable solutions for disposing of mine waste.   v  Preface All of the data analysis for this work is my own, with guidance and supervision from Dr. Gregory Lawrence and Dr. Edmund Tedford. The data collected for this work was done by Dr. Tedford, David Hurley, and others who worked at the site.  Partial meteorology observations presented in Chapter 2 have previously been published in Wind Waves and Internal Waves in Base Mine Lake by D.L. Hurley (MASc thesis). Further partial observations of meteorology and temperatures presented in Chapter 2 and Chapter 5 have been published in Temporal variations in turbidity in an oil sands pit lake by Tedford, E., Halferdahl, G., Pieters, R., & Lawrence, G. A., which was published in Environmental Fluid Mechanics (volume 19, Issue 2). Partial observations and results for the winter under-ice period were presented at the 2018 Ocean Sciences Meeting for the American Geophysical Union. This work was completed by me.  vi  Table of Contents Abstract ......................................................................................................................................... iii Lay Summary ............................................................................................................................... iv Preface .............................................................................................................................................v Table of Contents ......................................................................................................................... vi List of Tables ..................................................................................................................................x List of Figures ............................................................................................................................... xi List of Symbols ........................................................................................................................... xiii List of Abbreviations ................................................................................................................. xvi Acknowledgements ................................................................................................................... xvii Chapter 1: Introduction ................................................................................................................1 1.1 Surface Fluxes ................................................................................................................. 1 1.2 Evaporation Models and Inputs ...................................................................................... 2 1.2.1 Model Formulation ..................................................................................................... 3 1.2.2 Humidity ..................................................................................................................... 5 1.2.3 Near-Surface Temperatures: Stratification and Skin Effects ...................................... 7 1.2.4 Effect of Hydrocarbon Sheen on Evaporation ............................................................ 8 1.3 Summary ....................................................................................................................... 11 Chapter 2: Field Methods and Observations ............................................................................12 2.1 Site Description ............................................................................................................. 12 2.1.1 Water Budget and Hypsometry ................................................................................. 14 2.2 Field Methods ............................................................................................................... 16 2.2.1 Temperature Measurements ...................................................................................... 17 vii  2.2.2 Meteorological Data.................................................................................................. 17 2.2.3 Eddy Covariance System .......................................................................................... 18 2.3 Field Observations ........................................................................................................ 18 2.3.1 Water Temperatures .................................................................................................. 18 2.3.2 Meteorological Conditions........................................................................................ 19 2.3.3 Eddy Covariance Data .............................................................................................. 23 2.3.4 Focus Period.............................................................................................................. 23 Chapter 3: Modelling ...................................................................................................................27 3.1 Model Setup .................................................................................................................. 27 3.1.1 Estimation of Heat Budget Components................................................................... 27 3.1.2 Cumulative Heat Budget ........................................................................................... 29 3.2 Air-Sea Toolbox Approach ........................................................................................... 30 3.2.1 Determining TS .......................................................................................................... 31 3.2.2 Standard AST Results ............................................................................................... 32 3.2.3 Comparison with Eddy Covariance Methods ........................................................... 32 3.3 Extended Air-Sea Toolbox ........................................................................................... 35 3.3.1 Adjustments to the Interfacial Relative Humidity .................................................... 35 3.3.2 Extended Air-Sea Toolbox Results ........................................................................... 38 3.4 Examining the Overall Heat Budget ............................................................................. 38 3.5 Heat Flux Comparisons................................................................................................. 40 3.5.1 Radiative and Sensible Heat Fluxes .......................................................................... 40 3.5.2 Latent Heat: Overnight Cooling................................................................................ 42 3.6 Summary ....................................................................................................................... 42 viii  Chapter 4: Open-Water Period Discussion ...............................................................................44 4.1 Physical Processes and Modelling ................................................................................ 44 4.1.1 Skin Temperature of the Water ................................................................................. 44 4.1.2 Effects of Hydrocarbon Sheen .................................................................................. 46 4.1.3 Reduced Evaporation ................................................................................................ 48 Chapter 5: Under-Ice Period ......................................................................................................51 5.1 Field Observations ........................................................................................................ 51 5.1.1 Thermistor Data ........................................................................................................ 51 5.1.2 Meteorological Data.................................................................................................. 55 5.2 Heat Flux Calculation Methods .................................................................................... 56 5.2.1 Ice Cap Formation and Loss ..................................................................................... 56 5.2.2 Heat Flux from Snow ................................................................................................ 63 5.2.3 Basal Heat Fluxes ..................................................................................................... 63 5.3 Winter Heat Budget Results.......................................................................................... 65 5.3.1 Upper Water Column ................................................................................................ 65 5.3.2 Lower Water Column ............................................................................................... 70 5.4 Winter Summary ........................................................................................................... 73 Chapter 6: Conclusion .................................................................................................................75 6.1 Open-Water Period ....................................................................................................... 75 6.2 Under-Ice Period ........................................................................................................... 76    ix  Bibliography .................................................................................................................................78 Appendix A : Parameterization of the Air-Sea Toolbox .......................................................... 86 A.1 Theory and Previous Models .................................................................................... 86 A.2 Theory of Similitude ................................................................................................. 89 A.3 Parameterization ....................................................................................................... 92 A.4 Surface Temperature Estimation............................................................................... 99 A.5 List of Assumptions .................................................................................................. 99 Appendix B : Determining the Near-Surface Temperature .................................................... 101  x  List of Tables Table 1.1 Equations for heat transfer from momentum, sensible heat and latent heat. ...................4 Table 2.1 Water balance data in mm for April 2015 – April 2016 ................................................16 Table 3.1 Methods used to estimate each component of the heat budget ......................................28 Table 5.1 Parameter values used to estimate ice growth ...............................................................61 Table A.1 Variables and parameters used in COARE. ..................................................................87 Table A.2 Formulae for calculating ψ in different atmospheric conditions ..................................95 Table A.3 Parameter values as assigned in the Air-Sea toolbox .................................................100 Table B.1 Stratification parameter values in various wind and net radiation conditions ............103  xi  List of Figures Figure 1.1 Specific humidity values as a function of relative humidity and temperature ...............6 Figure 1.2 Typical diurnal patterns in water and skin temperatures ................................................9 Figure 2.1 Location of Base Mine Lake in Alberta, Canada .........................................................13 Figure 2.2 A photo of the hydrocarbon sheen during the open-water season ................................14 Figure 2.3 Hypsometric curve of horizontal lake area at each depth.............................................15 Figure 2.4 Atmospheric and limnological conditions for the 2015 open-water period .................20 Figure 2.5 Meteorological conditions binned by time of day ........................................................22 Figure 2.6 Relative and specific humidity at BML in the 2015 open water season ......................24 Figure 2.7 Sensible and latent heat fluxes from eddy covariance instruments at P1 .....................25 Figure 2.8 Atmospheric and limnological conditions at BML during the focus period ................26 Figure 3.1 Comparison between measured and modelled BML heat content ...............................33 Figure 3.2 Comparison of eddy covariance methods and model output for HS nd HL ..................34 Figure 3.3 Visual representation of sheen effects on RH at the air-water interface ......................36 Figure 3.4 RMSE between the Extended AST and the TBE estimates of HL................................37 Figure 3.5 Comparison between measured and modelled lake heat content .................................39 Figure 3.6 Averaged hourly heat fluxes through autumn unstratified overturn periods ................41 Figure 4.1 Percentage difference in HL from Standard and Extended AST results .......................47 Figure 4.2 Transfer coefficients for momentum, heat, and moisture.............................................49 Figure 5.1 Field conditions during the 2015-2016 under-ice period .............................................52 Figure 5.2 Lower water column and FFT temperature measurements ..........................................54 Figure 5.3 Winter meteorological conditions binned by time of day ............................................57 Figure 5.5 Measured hourly average snow depth ..........................................................................58 xii  Figure 5.5 Modelled ice thickness over winter ..............................................................................62 Figure 5.6 Winter heat fluxes at the water surface ........................................................................67 Figure 5.7 Heat fluxes in the lower half of the water column .......................................................72 Figure B.1 Influence of wind and net radiation on stratification during fall overturn .................102 Figure B.2 Scatter plot of corrected and uncorrected surface temperature estimates ..................104  xiii  List of Symbols A Lake surface area 𝑏𝑜 Radiative losses to the atmosphere during freezing 𝑏1 Turbulent losses to the atmosphere during freezing 𝐶𝑑 Transfer coefficient  for shear stress 𝐶𝑒 Transfer coefficient for latent heat (moisture) 𝐶ℎ Transfer coefficient for sensible heat 𝑐𝑝 Specific heat of air or water at constant pressure 𝑑 Lake depth 𝑒𝑤 Water vapor pressure 𝑘 Thermal conductivity of ice 𝑘𝑠 Thermal conductivity of the sediments 𝑘1 Coefficient for turbulent exchange of heat with the atmosphere in melting ℎ𝑝 Water equivalent of precipitation ℎ𝑠 Height of snow cover 𝐻𝑓 Latent heat flux due to ice cap freezing 𝐻𝐹𝐹𝑇 Heat flux from advection of water from underlying sediments 𝐻𝐿 Latent heat flux due to evaporation 𝐻𝑚 Latent heat flux due to ice cap melting 𝐻𝑝 Latent heat flux due to precipitation 𝐻𝑆 Sensible heat flux 𝐻𝑠𝑒𝑑 Heat flux due to conduction with sediments 𝐻𝑠𝑢𝑏 Latent heat flux due to sublimation xiv  𝐿𝐸 Latent heat of vaporization 𝐿𝑓 Latent heat of freezing and melting 𝐿𝑠 Latent heat of sublimation 𝐿𝑊𝑜𝑢𝑡 Outgoing longwave radiation 𝑃𝐴 Air pressure 𝑃 Pressure at the freezing surface 𝑄 Net heat flux 𝑄𝑆 Incident solar radiation  𝑄𝑇 Solar radiation transmitted through the surface layer 𝑞 Specific humidity, or water vapor mixing ratio 𝑞∗ Monin-Obukhov similarity scaling parameter for humidity 𝑞𝑠 Saturation specific humidity RH Relative humidity 𝑅𝐻𝑊 Relative humidity at the air-water interface 𝑅𝑛𝑒𝑡 Net radiation 𝑆 Water salinity 𝑇∗ Monin-Obukhov similarity scaling parameter for temperature 𝑇𝐴 Air temperature 𝑇𝑓 Temperature of freezing 𝑇𝑚(𝑧) Water temperature measured at depth 𝑧 𝑇𝑝 Temperature of precipitation 𝑇𝑆 Water surface temperature 𝑇0 Reference temperature  xv  𝑢∗ Shear velocity 𝑢𝑧 Measured wind speed at height z  𝑢10 Wind speed normalized to standard height of 10 m 𝑉 Volume 𝑧 Measurement depth/height (𝑑𝑇𝑠𝑑𝑧)𝑤𝑠 Sediment temperature gradient at the sediment-water interface 𝛥ℎ𝑓 Change in height of ice due to freezing 𝛥ℎ𝑚 Change in height of ice due to melting 𝛥ℎ𝑠 Change in height of ice due to sublimation 𝛥𝑡 Time step duration 𝛥?̅? Change in mean water column temperature 𝛥𝑇𝑐 Temperature correction for cool-skin effects 𝛥𝑇𝑤 Temperature correction for the warm layer 𝛼 Surface albedo 𝛼𝑡 Parameter describing near-surface stratification 𝛽 Parameter describing near-surface heating 𝛽𝑠 Insulation efficiency of snow 𝜃 Potential temperature of air brought adiabatically to a standard pressure 𝜌𝑎 Density of air 𝜌𝑖𝑐𝑒 Density of ice 𝜌𝑤 Density of water 𝜏 Shear stress 𝛾 Fraction of incident solar radiation transmitted through the surface xvi  List of Abbreviations  AST Air-Sea Toolbox  BML Base Mine Lake  COARE Coupled Ocean Air Response Experiment  COSIA Canada’s Oil Sands Innovation Alliance  DOY Day of year  EC Eddy covariance  FFT Fluid Fine Tailings  masl Metres above mean sea level   NTU Nephelometric Turbidity Units  ppt Parts per thousand  RMSE Root mean squared error  TBE Temperature-Based Estimate  TVA Tennessee Valley Authority xvii  Acknowledgements I owe tremendous thanks for the help and support of many individuals over the course of my research.  Many, many thanks must go to Gregory Lawrence and Ted Tedford for their tremendous patience and guidance over the evolution of this research. Their encouragement and belief in my results and research abilities helped me to persist in my investigations, and their insights were invaluable contributions to this thesis.  Additional thanks must go to Ted Tedford for the time spent at Base Mine Lake – the humour and insights, scientific and otherwise, that you provided in Fort McMurray regardless of conditions were most enjoyable. To the friendly, knowledgeable, and capable team at Building 578: thank you for facilitating this work and working to advance reclamation and research day in and day out, rain or shine. This would not have been possible without you, and without your friendly cheer, the work would not have been nearly as enjoyable.  The National Sciences and Engineering Research Council of Canada (NSERC), Syncrude Canada Ltd., and the University of British Columbia must be thanked for providing the funding for this research.  To Kelly Graves, Lydia Smith, and Cristi Pinela, many thanks for the moral support, for looking over my document to help keep me on track, and for otherwise helping me to get the job done.  Thanks to the UBC EFM team – the camaraderie, encouragement, and overall team spirit shown by everybody on the research team made the entire experience more enjoyable. Thank you to my family and friends for always being there for me. 1  Chapter 1: Introduction Safe and environmentally sustainable tailings management is a key part of a holistic mining plan, and requires an approach that is tailored to the type of waste material being produced. As a prototype for tailings sequestration in the closure landscape of the Canadian oil sands, the characteristics and processes of Base Mine Lake (BML) are the subject of significant research (Canada’s Oil Sands Innovation Alliance (COSIA), n.d.; Prakash et al., 2011). Investigations are underway on a variety of topics including biochemistry and transport processes of the tailings, meteorology and water budget of the site, and lake physics. This research investigates the impact of hydrocarbons on the heat budget of Base Mine Lake during the open water season. This section will provide a background on heat budgets and their component fluxes, followed by an overview of different methods used to estimate these fluxes, with a focus on latent heat flux due to evaporation.  The bulk of this thesis focuses on the open water period, during which no ice is present on the lake surface, with an aim of developing an understanding of the extent to which heat fluxes are affected by the hydrocarbon sheen on the water surface at the study site. In Chapter 5, a simple model of fluxes across the ice-water and the water-sediment interfaces is presented in order to assess the relative magnitudes of each heat flux. The impacts of these heat fluxes on the under-ice dynamics in the lake are then discussed. 1.1 Surface Fluxes Heat budgets, which estimate the sources and quantities of heat entering and leaving a system, provide insight into the extent and duration of stratification, the influence of buoyancy-driven mixing, the timing for the onset of ice cover, and the amount of evaporation from the water surface (Momii & Ito, 2008). These are key parameters influencing the biological productivity, water quality, and water budget of the lake.   2  In a basin without significant inflows or outflows, the heat budget is dominated by the surface fluxes across the air-water interface:   𝑄 = 𝑅𝑛𝑒𝑡 + 𝐻𝐿 + 𝐻𝑆 (1.1) where 𝑄 is the net heat flux, 𝑅𝑛𝑒𝑡 is the net radiation including shortwave and longwave components,  𝐻𝐿 is latent heat flux and 𝐻𝑆 is sensible heat flux (Imboden & Wuest, 1995). Radiation at the site was assessed using direct measurements of downward shortwave, upward shortwave, downward longwave and upward longwave radiation collected at P1. The presence of a hydrocarbon sheen on a water surface impacts radiation and its interactions with the water body sufficiently that it can be detected by remote sensing (Jha et al., 2008). However, by measuring radiation directly at the site rather than modelling, as is done here, the uncertainty in this variable is reduced, allowing for more accurate assessment of the effects of the hydrocarbon sheen on other fluxes.  Sensible heat fluxes account for the conduction of heat between the air and water, and are a function of the temperature difference between the water and atmosphere at the lake surface (refer to Section 1.2.1). Although a hydrocarbon sheen may change the temperature of the water due to changes in radiation and evaporation, here it is assumed that the physical mechanism causing sensible heat fluxes is not changed with the addition of a hydrocarbon sheen.   By comparison, surface sheens have long been observed to suppress evaporation of water (Rideal, 1925; Roberts, 1962; refer to Section 1.2.4), and in so doing alter the heat budget by reducing the latent flux from evaporation. The contribution of the latent heat flux to the overall heat budget of a lake is determined by the lake characteristics and setting, sometimes comprising 40 % or more of the annual heat budget (Momii & Ito, 2008; Xing et al., 2012). This potentially large contribution makes accurate determination of latent heat flux critical to closing the energy budget.  1.2 Evaporation Models and Inputs A report by the Tennessee Valley Authority (TVA) (Wunderlich, 1972) identifies 4 methods for estimating evaporation rates: the water budget method, the energy budget method, the empirical method, and model-based estimates. The first two are simple mass or temperature budget 3  methods, which can offer insights given sufficiently detailed measurements of water quantity or temperature changes in the lake, respectively. However, these are challenging to obtain, and reflect only the net heat flux, with individual component fluxes difficult to discern. Empirical estimates are site-specific, using field results from evaporation pans, and this approach cannot offer insights more detailed than the long-term means at the site (Wunderlich, 1972). Model-based methods, described in greater detail in Section 1.2.1 and in Appendix A, take readily measured site conditions, such as wind speed, air temperature, and water temperature, and input them into models developed from the principles of heat transfer physics. Finally, heat fluxes can be determined directly, via eddy covariance (EC) measurements. Eddy covariance instruments measure exchanges of heat, mass, and momentum between a flat, homogeneous surface and the overlying atmosphere (Aubinet et al., 2012). By estimating the covariance between the turbulent fluctuations in vertical velocity and the parameter of interest, measurements of many variables can be obtained over extensive spatial and temporal scales (Smith & Jones, 1986; Businger, 1986; Gu et al., 2012).  Eddy covariance methods are limited, however, by a number of factors. When local flows around the instrument are not representative of site conditions, inaccurate readings can result. Data losses can also occur due to instrument errors caused by poor adjustment, electronic noise, contamination from external factors, and other sources, or due to limitations in instrument design or signal processing (Aubinet et al., 2012; Massman & Lee, 2002). Although EC methods cannot be used to close the energy budget (Wilson et al., 2001), they can provide some insights into lake behaviour. 1.2.1 Model Formulation Model-based estimates use the atmospheric and water surface conditions to compute heat fluxes. The formulation and parametrization of these equations vary from model to model (refer to Liu et al., 1979; Stull, 1988; Godfrey & Beljaars, 1991; Garratt, 1994). However, the general form and inputs are consistent between models: transfer of momentum, heat, and moisture are governed by the difference between the atmospheric and boundary layer conditions, and subject to a transfer function that describes the extent to which this forcing is transmitted from the boundary layer to the atmosphere.  4  As presented in Fairall et al. (1996b), analogous equations estimate transfer of momentum, heat, and moisture between the lake and the atmosphere. Examining the momentum transfer equations to understand the principles involved, shear stress at the water surface (𝜏), can be expressed as:  𝝉 = −𝝆𝒂𝒖∗ 𝟐 =  𝝆𝒂𝑪𝒅𝒖𝟐 (1.2) where 𝜌𝑎 is the density of air; 𝑢∗ is the shear velocity, alternately known as the Monin-Obukhov similarity scaling parameter for momentum; 𝐶𝑑  is a momentum transfer coefficient; and 𝑢 is the average value of wind speed relative to the water surface.  Shear velocity, 𝑢∗, is a term developed from the assumption that the theory of similitude applies to momentum within the surface layer (Von Kármán, 1931; Monin & Obukhov, 1954). Similitude implies that a universal relationship exists between wind speed and shear stress, which can be used to describe the behaviour of a system at different scales than the one being measured (Hutter et al., 2014). Analagous equations for 𝐻𝐿 and 𝐻𝑆 from three different sources are introduced in Table 1.1. Equations have been modified for consistency in variable usage. In the TVA report (Wunderlich, 1972), latent heat was presented as a function of the drag coefficient on the water surface. This approach was possible because the transfer coefficients for momentum and moisture were assumed to be equal in value. Though expedient, this assumption is not accurate. Furthermore, inaccurate estimation of transfer coefficients can lead to significant Table 1.1: Equations for heat transfer from momentum, sensible heat and latent heat. Variables are defined in the List of Symbols.   Monin-Obukhov (1954) TVA (1972) Fairall et al. (1996b) Momentum 𝜏 =  −𝜌𝑎𝑢∗2 - 𝜏 =  𝜌𝑎𝐶𝑑𝑢2 Sensible Heat 𝐻𝑆 =  −𝜌𝑎𝑐𝑝𝑢∗𝑇∗ 𝐻𝑆 = −𝜌𝑎𝑐𝑝𝐶𝑑𝑢(𝑇𝑠 − 𝑇𝑎) 𝐻𝑆 = −𝜌𝑎𝑐𝑝𝐶ℎ𝑢(𝑇𝑆 − 𝜃) Latent  Heat 𝐻𝐿 =  −𝜌𝑎𝐿𝐸𝑢∗𝑞∗ 𝐻𝐿 = −𝜌𝑎𝐿𝐸𝐶𝑑 𝑢(𝑞𝑠 − 𝑞) 𝐻𝐿 = −𝜌𝑎𝐿𝐸𝐶𝑒𝑢 (𝑞𝑠 − 𝑞)  5  errors in flux calculations (Blanc, 1985), and subsequent research has developed methods for estimating transfer coefficients more accurately. Fairall et al. (1996b) use an iterative approach to estimate the transfer coefficients. This is necessary because the transfer coefficients are a function of atmospheric stability and therefore of the fluxes themselves (refer to Appendix A for further details). Despite the relative complexity of this approach, however, the key inputs that affect results are straightforward to measure: air temperature, water surface temperature, the humidity of the atmosphere, and the wind speed. Humidity and water surface temperature considerations are discussed further below. 1.2.2 Humidity Relative humidity is the amount of moisture in the air relative to the moisture content at saturation for the same atmospheric conditions:  𝑅𝐻 = 100 ∗ 𝑞𝑞𝑠 (1.3) where 𝑅𝐻 is relative humidity (%), 𝑞 is specific humidity (kg water per kg moist air), and 𝑞𝑠 is the saturation specific humidity, i.e. the specific humidity for air at 100 % relative humidity. Relative humidity does not by itself offer any information about the quantity of water in the air.  Specific humidity, 𝑞, is the mass of water that is contained in a given mass of moist air, and is calculated using air temperature (𝑇𝐴), and air pressure (𝑃𝐴), and relative humidity. Many formulae exist for this calculation; here, we follow the approach of the Coupled Ocean Atmosphere Response Experiment (COARE) as described in Fairall et al. (1996b), which uses the formulae first presented in Buck (1981):  𝑞 =𝑅𝐻100∗  0.622 𝑒𝑤𝑃𝐴−0.378 𝑒𝑤 (1.4)  𝑒𝑤  = 6.112 (1.001 + 3.46 x 10−6 𝑃𝐴 17.502 𝑇𝐴240.97+𝑇𝐴) (1.5) where 𝑒𝑤 is the vapor pressure of water in the air. 𝑞𝑠 is calculated by setting 𝑅𝐻 to a value of 100 % and setting 𝑇𝐴 to the same value as 𝑇𝑠, the temperature of the water surface. 6  The nonlinear influence of 𝑇𝐴 means that while cool air at a low 𝑅𝐻 holds slightly less water than cool air at 100 % 𝑅𝐻, the amount of moisture that is held in a parcel of air increases dramatically when the air is both warm and humid (Figure 1.1). Thus, 𝑞 provides a measure of the actual quantity of water in a given parcel of air and so can be used to track the exchange of moisture between the lake and the atmosphere. In warm conditions, the air can hold more moisture, increasing the maximum specific humidity possible, but because the relative humidity is a ratio of the measured to saturation specific humidity, short-term increases in air temperature without corresponding increases in the actual moisture content of the air will lower the relative humidity.  Figure 1.1 Specific humidity values as a function of relative humidity and temperature.     7  1.2.3 Near-Surface Temperatures: Stratification and Skin Effects The strong influence of the water surface temperature (𝑇𝑆) on humidity means that an accurate estimation of this variable is required for estimating both 𝐻𝐿 and 𝐻𝑆 (Table 1.1; Equation 1.4). Infrared measurements of the skin temperature can be used as a direct estimate of 𝑇𝑆; however, when these are unavailable, near-surface bulk temperature measurements must be used. As outlined in Fairall et al. (1996a), the temperature used for input into dynamic models is:  𝑻𝑺 = 𝑻𝑴(𝒛) + 𝜟𝑻𝑾(𝒛) − 𝜟𝑻𝑪 (1.6)  where 𝑇𝑆 is the water temperature at the air-water interface, 𝑇𝑀(𝑧) is the water temperature measured at depth 𝑧, 𝛥𝑇𝑊 is the temperature correction for the warm layer between depth 𝑧 and the surface, and 𝛥𝑇𝐶 is the temperature correction for the cool skin of the lake. In order to correctly estimate 𝑇𝑆, both the warm layer and the cool skin effects must be evaluated (Alappattu et al., 2017).  Skin effects arise due to heat and moisture fluxes at the water surface. The boundary layer of the water varies in thickness as a function of wind speed (Wu, 1985); however, it is generally <1 mm in thickness (Handler et al., 2001), with the skin itself on the order of 10 μm (Kawai & Wada, 2007).  Further from the boundary layer, the thermal gradient is destroyed by turbulent mixing (Fairall et al., 1996b). Within the interfacial boundary layer, heating and cooling of the water via solar insolation and surface fluxes create a thermal gradient (Saunders, 1967). These processes typically cool the skin of the water surface; Donlon et al. (1999) found the mean 𝛥𝑇𝐶 was -0.17 K cooler than the bulk water temperature measured directly (~5 cm) beneath the water surface. Skin effects have strong diurnal trends, with the cool-skin effect being elevated in morning and in periods with high radiation, and becoming less significant towards night (Wilson et al., 2013).  Convective cooling at night causes warm water near the surface to mix with water below it, creating a diurnal mixed layer in the upper few metres (Imberger, 1985). This can be seen in Figure 1.2a, where the water column above 𝑇𝑚(𝑧) is near-uniform up to the water surface. Cool-skin effects are still observed at night (𝛥𝑇𝐶), leading to a slightly lower skin temperature than what is measured by the thermistors. In contrast, during the day, heating from insolation and 8  other factors causes the near-surface water to stratify (Figure 1.2b). The depth at which the bulk temperature is measured will therefore significantly impact resulting observations, and consequent warm-layer (𝛥𝑇𝑊) adjustments, especially on days with high insolation and limited mixing. The presence of hydrocarbons on a water surface introduces further complications in determining 𝑇𝑆. Monolayers on water surfaces have been observed to reduce evaporation, thus reducing cooling from the latent heat of evaporation and reducing cool-skin effects (Barnes, 1993). However, it has also been hypothesized to contribute to immobilization of the near-surface water (Saunders, 1967), which could have impacts to both warm-layer and cool-skin temperature adjustments. 1.2.4 Effect of Hydrocarbon Sheen on Evaporation Researchers, starting with Rideal (1925), have previously noted the influence of surface sheens on evaporation, and their consequent potential to decrease the latent heat flux. Much of the research on the effects of hydrocarbons on lakes has focused on developing water conservation techniques in water-scarce regions. Field tests performed in the 1950s involving spreading emulsifiers and alkanols onto pilot-scale tanks, found reductions in evaporation as high as 27 % over the course of a two-week period (Koberg et al., 1963). Similar effects have been observed for hydrocarbon sheens; Anikiev et al. (1988) observed that water evaporation from a system with a hydrocarbon sheen was reduced by 15-33 %. The reduction in the evaporation rate, expressed as resistance to evaporation, is caused by properties of the sheen that resist diffusion of water from the liquid phase to the vapour phase (Langmuir & Schaefer, 1943). Resistance occurs because of the additional energy required for the water molecules to permeate the sheen, and is a function of temperature and hydrocarbon molecule chain length (Henry et al., 2010). Film thickness, and therefore the spreading properties of the compound, is also important, with thicker films observed to cause greater reductions in evaporation rates (Barnes, 1993).  Typically, moisture transfer is assumed to be air-limited, meaning that the moisture in the air at the air-water interface is replenished instantaneously and moisture transfer into the atmosphere is 9   Figure 1.2 Typical diurnal patterns in water and skin temperature. Note the scale change at 0.5 m. The thickness of 𝜟𝑻𝒄 is lower than what is shown here (typically order of 1 mm thick). Data below 0.05 m depth are from on-site thermistors on August 28, 2015; above 0.05 m, modelled cool-skin effects are presented.  The topmost sensor measuring water temperature over the full period of record was at 0.5 m, which is used for 𝑻𝒎(𝒛) in Equation 1.6. At night, radiative cooling combined with elevated winds create a mixed layer near the water surface, and the bulk water temperature near the surface is close to 𝑻𝒎(𝒛). During the day, typically lower winds and high insolation stratify the epilimnion, inducing a warm layer (𝜟𝑻𝑾). Both periods are subject to surface skin effects (𝜟𝑻𝑪).  10  limited by the rate at which moisture leaves the atmospheric boundary layer (Haghighi et al., 2013). For this reason, relative humidity at the boundary layer on the water side of the air-water interface is typically maintained at 100 %, though above sea water, the presence of 20 ppt salt reduces saturation to 98 %, as dissolved salts block the remaining 2 % of the water surface (Katsaros, 2001). However, because the sheen introduces resistance to evaporation at the interface, moisture in the air at the water surface can no longer be replenished instantaneously. Here, the presence of the sheen is expressed as a reduction in the relative humidity at the boundary layer (𝑅𝐻𝑊). Determining the parameters under which this effect is applicable requires an examination of the behaviour of the hydrocarbon sheen, which appears to have a strong wind dependency. Hurley (2017) observed that at low wind speeds (<5 m∙s-1), there is significant dampening of the <0.1 m waves that would be expected to form. This was postulated to be due to the presence of the hydrocarbons on the water surface, which dampen the development and propagation of waves due to the changes in surface tension at the air-water interface. Fitzgerald (1963) found that waves were suppressed at mean wind velocities below 5 m∙s-1. Fitzgerald’s experimental setup yielded a critical shear velocity of 25 cm∙s-1, a threshold that was also observed in Gottifredi & Jameson (1968) to be associated with wave suppression.  At low wind speeds, the dampening of the capillary waves observed at Base Mine Lake by Hurley (2017) and further documented in Hurley et al. (2019) suggests that a persistent hydrocarbon layer is present at the lake surface and is limiting air-water interaction. The reduction in wave action and air-water interaction may also lead to a reduction in the transfer of momentum from the atmosphere to the water, resulting in reduced mixing in the epilimnion and changes to the stratification in this portion of the water column. In contrast, during storm events with persistent wind speeds above 5 m∙s-1, measured wave magnitudes corresponded well with model predictions (Hurley, 2017). This change in behaviour is likely from the breakup of the hydrocarbon layer that occurs once sufficient stress is applied to the water surface. The effects of the hydrocarbon sheen were therefore assumed to approximate a step function. At wind speeds below 5 m∙s-1, where the sheen is expected to be relatively intact across the water 11  surface, 𝑅𝐻𝑊 should be reduced (refer to Figure 3.3). At wind speeds greater than 5 m∙s-1, where the sheen is expected to be absent, 𝑅𝐻𝑊 is set to 100 %. 1.3 Summary Closing the heat budget for a lake requires accurate determination of input parameters for modelling individual fluxes, most notably water temperature, air temperature, and humidity levels. The presence of a hydrocarbon sheen on the water surface introduces additional challenges in determining the water surface temperature, and potentially changes the amount of moisture available at the air-water interface, limiting evaporation. The impacts of the sheen on water availability are expressed herein as a reduction in the relative humidity at the air-water interface. 12  Chapter 2: Field Methods and Observations This chapter outlines the site history and conditions. The field methods used to record site conditions are then described and the resulting field measurements are presented. This research focuses on the open-water period (May 6 – Nov 18, 2015), with particular focus on the autumn turnover period (Sept 4 – Nov 18, 2015). Discussion of the winter inverse stratification period (Nov 19, 2015 – April 19, 2016) is presented in Chapter 5. 2.1 Site Description Base Mine Lake is situated north of Fort McMurray, Alberta, Canada at Syncrude’s Mildred Lake Mine (57°1’N, 111°37’W). The lake has a surface area of roughly 8 km2 (Figure 2.1). BML is a dimictic lake, with partial or full overturn observed in spring and autumn each year, and inverse stratification observed in winter. The lake is located on the footprint of a former oil sands pit, previously known as the West In-Pit. Upon completion of mining in 1994, fluid fine tailings (FFT) were dredged from the nearby Mildred Lake. FFT is a fluid suspension of bitumen (1-6% w/w) and mineral solids (initially 25-35% w/w, dominated by silt and clay (Dompierre et al., 2017)). The FFT is experiencing either hindered or self-weight settlement (COSIA, 2012). These phenomena occur when the solid matrix lacks the strength to support the mass above it (Guo et al., 2015). This causes upward flux of pore water through the matrix and out into the lake, leading to a decrease in the FFT deposit thickness. Because settlement is governed by matrix strength and surrounding stresses, the uneven FFT depth and placement history has resulted in spatially non-uniform settlement. On average, the FFT deposit initially settled by 1 m∙year-1, a rate that will decrease exponentially over 30 years to a rate of 0.1 m∙year-1 (Dompierre et al., 2017).  Residual hydrocarbons are released from the FFT deposit and result in sheens that are readily observed on the surface of BML in calm conditions (Figure 2.2). The strength, duration, and direction of winds on the lake have a strong influence on the behaviour of the hydrocarbon sheen. The hydrocarbon sheen is therefore spatially and temporally non-uniform in character.  13   Figure 2.1 Map of Base Mine Lake in Alberta, Canada (Geoff Halferdahl, Syncrude Canada Ltd. personal communication 2016), showing lake bathymetry and platform locations. Meteorological data and near-surface water temperatures are taken at Platform 1 (P1); moorings with temperature chains, pressure sensors, and turbidity sensors are located at Platform 2 and 3 (P2 and P3).   14   Figure 2.2 A photograph of the hydrocarbon sheen during the open-water season (Tedford, 2016). Note the smooth lake surface, presenting few capillary waves. 2.1.1 Water Budget and Hypsometry The water cap volume at the start of the 2015 open water period was 60 Mm3, corresponding with an average depth of 7.9 m (Carey, 2017). Pumping of water into and out of the lake occurs seasonally due to process water requirements, and to ensure that the lake elevation is maintained close to 308 masl (metres above mean sea level). This, in conjunction with FFT settlement, has resulted in a deepening water cap, with an average depth of 8.1 m and a volume of 62 Mm3 as of November 2015. The hypsometry of the basin was surveyed in October 2013 (Figure 2.3). The water level has not changed significantly in that time, remaining at approximately 308 masl. However, due to the spatially variable settlement of the FFT surface, the bathymetry of the FFT-water interface for  15   Figure 2.3 Hypsometric curve of horizontal lake area at each depth (Geoff Halferdahl, Syncrude Canada Ltd. personal communication 2017). Survey taken in October 2013 when the lake level was 308.54 m. The surface of the FFT, and therefore the surface area of the lake at depth, has continued to drop since the survey was taken. 2015 is difficult to estimate with accuracy. The water balance estimates for porewater input into the lake and a rough estimate of the settlement of the basin since the hypsometric survey occurred were used to create approximate volume estimates every 0.2 m over the full depth of the lake.  The water budget, prepared by researchers at McMaster University (Carey, 2017), was also used for estimates of rainwater input, runoff into the lake, and volumes of water pumped into and out of the lake (Table 2.1).  16  The water budget also includes evaporation estimates; however, the reliability of these estimates has not been established, and so they were used primarily as an initial guess and for comparison purposes.  2.2 Field Methods Data from the 2015 open water period and the winter of 2015-2016 were reviewed. Sensors are located at three large floating platforms anchored within BML (Figure 2.1), referred to as P1, P2 and P3.  At P2 and P3, permanent moorings record water temperature, pressure, and turbidity across the water column. Meteorological instruments are located at P1, near the centre of the  Table 2.1 Water balance data in mm for April 2015 – April 2016.  (Carey, 2017).  Month Volume (Mm3) Mean Depth (m) Monthly Water Inputs (mm) ET Pump OUT Pump IN Rain Snow Runoff FFT Apr-15 60.9 7.91 10.9 129.9 39.9 5.0 50.0 11.5 42.9 May-15 62.5 8.12 30.0 79.6 243.4 29.1 0.0 0.1 42.9 Jun-15 62.5 8.12 58.5 247.9 233.7 34.4 0.0 0.0 42.9 Jul-15 62.6 8.13 78.4 226.8 190.0 85.1 0.0 0.0 42.9 Aug-15 61.8 8.03 94.8 98.1 17.6 26.7 0.0 0.0 42.9 Sep-15 62.1 8.07 45.3 51.4 40.0 56.0 0.0 0.0 42.9 Oct-15 61.6 8.00 31.7 118.9 15.4 6.6 0.0 0.0 55.8 Nov-15 62.0 8.05 0.8 0.0 0.0 0.8 0.0 0.0 55.8 Dec-15 62.4 8.11 0.0 0.0 0.0 0.0 0.0 0.0 55.8 Jan-16 62.9 8.16 0.0 0.0 0.0 0.0 0.0 0.0 55.8 Feb-16 63.3 8.22 0.0 0.0 0.0 0.0 0.0 0.0 55.8 Mar-16 63.7 8.28 0.0 0.0 0.0 0.0 0.0 0.0 55.8 Apr-16 64.1 8.33 6.5 69.7 9.7 15.3 50.0 0.2 55.8 17  lake. Wind speed measured at a meteorological station located 1.5 km to the northeast of the lake, Sandhill Fen, was used to fill gaps in the data from P1. McMaster University Hydrology Group (McMaster University) operated eddy covariance instruments at P1 measuring sensible and latent heat flux.  2.2.1 Temperature Measurements RBR SoloT sensors were installed at P2 and P3 at 1 m intervals through the water column. At P2, this spanned 1 to 10 m; P3 is slightly shallower than P2, recording water temperatures at 1-9 m depth. Sensors also recorded temperatures at 0.5 m at P2 and P3, and at 1.5 m depth at P3. The thermistors recorded water temperatures every 10 seconds, and a 30-minute averaging period was used for the analysis. The data from the P2 and P3 stations were averaged for this work to reduce spatial effects, such as internal waves. The instruments have a resolution of <0.00005 °C and an accuracy of ± 0.002 °C (RBR Ltd., 2019). Due to fluctuations in the platform position and the depth of the FFT-water interface, the deepest thermistors (at 11 m at P2 and 10 m at P3) are below the FFT-water interface during most of the period of record.  The sensors are serviced twice annually: during fall turnover and after ice-off in spring. In 2015, the moorings were serviced in early May and late September. This process involves removing the chains from the lake, downloading the data that has been stored since the previous servicing, and resetting the mooring with reprogrammed instruments. The servicing process creates a gap in the data sets, typically on the order of a day.  These instruments were supplemented by Campbell Scientific Canada Model 107 temperature probe thermistors moored at P1 (tolerance of 0.02 °C). These sensors are less accurate than those at P2 and P3; however, they include near-surface placements at 0.05 m, 0.15 m, and 0.3 m depths, offering insights into the behaviour of the top half metre of the water column. The P1 thermistors recorded data between 29 July 2015 (day of year (DOY) 210) and 14 Oct 2015 (DOY 287), encompassing late summer into fall overturn.  2.2.2 Meteorological Data Hourly time-averaged measurements of air temperature, relative humidity, barometric pressure, shortwave and longwave radiation, and wind speed and direction are available from the P1 18  station. The mechanical anemometer at P1 is inaccurate at very low wind speeds. For this reason, the wind data was supplemented by data from the Sandhill Fen meteorology station.  The Sandhill Fen station has the wind sensor installed at a height 2.7 m above the ground. Equation 2.1 converts wind speed measurements taken at height 𝑧 to the standard height used for most meteorological calculations, 10 m, and includes a factor of 1.2 to convert from over-land to over-water wind speeds:  𝑢10 = 1.2 ∗ 𝑢 (10𝑧)17 (2.1) where 𝑢10 is the wind speed at a height of 10 m and 𝑧 is the measurement height of the instrument (Manwell et al., 2010). 2.2.3 Eddy Covariance System Eddy covariance instruments, measuring sensible and latent heat, were operated by researchers at McMaster University. Measurements are available at half-hour intervals between 26 Aug 2015 (DOY 238) and 13 Oct 2015 (DOY 286), spanning the end of summer stratification and a portion of fall turnover.  There are gaps in the data following quality control by the McMaster University Hydrology Group; these are observed regardless of wind speed. Flux measurements were used for comparison with and validation of the model results. 2.3 Field Observations Meteorological and thermistor data from the 2015 open water period are presented in Figure 2.4. This period extends from 6 May 2015 (DOY 126) until 18 Nov 2015 (DOY 322), the date of ice-on.  A sub-period from September 21-28, 2015 (DOY 264-271), highlighted in Figure 2.4 and presented in Figure 2.8, is included to provide detail on the daily patterns in the data. 2.3.1 Water Temperatures As of 6 May 2015 (DOY 126), the start of the period of record, the water column was well-mixed, with a mean temperature of 5.5 °C (Figure 2.4a). The water column then rapidly stratified, with the 0.5 m water reaching 14 °C by 12 May 2015 (DOY 132). The temperature of 19  the lake continued to warm through spring, with periodic storms cooling the epilimnion and mixing heat downwards through the water column; notable storm events occurred on 16 May, 28 May, and 13 June 2015 (DOY 136, 148, and 164). These storms mixed heat downwards through the water column and increased the thickness of the surface mixed layer. When wind speeds were persistently above 5 m∙s-1, the temperature difference between the 0.5 m and 1 m sensors decreased and at times approached zero. Stratification of the water column became more consistent after 14 June (DOY 165). The surface sensors warmed to the peak seasonal temperatures on 30 June (DOY 181), with the 0.5 m sensor climbing to 24 °C on this date. The epilimnion did not warm significantly beyond this date; instead, heat was transferred downwards through the water column as the season progressed (Figure 2.4a). The overall heat content of the lake reached a peak in mid-August (DOY 225) before declining. Fall overturn began 4 Sept (DOY 247) and continued until 18 Nov 2015 (DOY 322), the date of ice onset. Brief stratification events occurred through this period, particularly in the beginning of overturn. The instantaneous temperature difference rose as high as +1.8 °C between 0.5 m and 1 m depth during overturn. This phenomenon was highly dependent on atmospheric and lake conditions. Towards the later portion of fall overturn, the upper water column was consistently well-mixed, but occasional warming of the bottom thermistors was observed. The 11 m thermistor, located at P2 in the deepest part of the lake, was initially within the water cap; the thermistor entered the FFT layer by 19 June (DOY 170) and remained there through the rest of the open water period. This change in placement caused recorded temperatures to be less variable, changing slowly with changes to the FFT heat content. The thermal peak of the FFT occurs later than that of the water cap, lagging by 41 days, with the surface of the FFT continuing to warm until it crosses above the cooling lake temperature on 23 Sept (DOY 266), at which point it begins to cool (Figure 2.4a).  2.3.2 Meteorological Conditions Air temperatures were at their highest in summer and during the day, with daily maximums of 27-29 °C recorded for much of July before declining in autumn (Figure 2.4b). The air temperature at BML peaks in the evening and drops overnight, reaching a minimum around 20     Figure 2.4 Atmospheric and limnological conditions at BML during the open-water period (May 6 to November 18, 2015). Data presented are 1-h averages. Panel (a) presents the average of the water temperatures recorded at P2 and P3. The 11 m sensor is located below the water-FFT interface at P2 and is recording FFT temperatures at a point no more than 1 m below this interface. Panels (b)-(e) display the air temperature, relative humidity, wind speed, and net radiation, respectively. The shaded region highlights the period presented in Figure 2.8. The vertical dashed lines at DOY 165 indicates the start of the stable summer stratification period. A horizontal dashed line in panel (e) shows where net radiation is zero.  21  dawn before beginning to increase again when the sun rises.  The relative humidity levels also had a seasonal pattern, with peak values sometimes reaching 100 % in autumn, significantly higher than the daily peak of closer to 60-80 % in spring and summer (Figure 2.4c). The daily lows were also typically lower in spring and summer, with typical values of 10-50 % near the start of the open water period, much lower than the 40-60 % daily minimums observed in late summer and autumn.  Wind speed, 𝑢, had median and maximum hourly values of 4.5 m∙s-1 and 16.5 m∙s-1 respectively (Figure 2.4d). Seasonal trends in 𝑢 were not observed. Net radiation, 𝑅𝑛𝑒𝑡, had hourly averages reaching peaks of 700 W∙m-2 daily for the majority of the summer before declining steadily through autumn (Figure 2.4e).  Meteorological conditions, binned by time of day and averaged to show the typical values that occur each hour across a 24-hour day, are shown in Figure 2.5. The left column includes all data, while the right column shows the average when wind was below 5 m∙s-1. There is more scatter in the latter case due to the lower number of data points. The period starts at 06:00, and noon and midnight are marked with black lines. Net radiation rose during the day from a near-constant minimum at night of just over -100 W∙m-2 (Figure 2.5a). The peak daily net radiation was slightly higher in low-wind periods (Figure 2.5b). Wind speeds increased slightly over the day, dropping again in the evening (Figure 2.5c). 𝑅𝐻 reaches a minimum around midday, rising overnight to a peak near dawn (Figure 2.5e-f). This phenomenon is due to its inverse relationship with air temperature, which displays the opposite trend (Figure 2.5 g-h). The mean water temperature of the lake is near-constant throughout the day, due in large part to the high heat capacity of water and the influence of the deep parts of the water column, which are less influenced by meteorological conditions than the epilimnion. Due to the inverse relationship of relative humidity and temperature, the seasonal patterns in humidity are not always clear. Evaluating the specific humidity over the season shows more clearly that the moisture in the air increases through summer and peaks when air temperatures are at their highest, dropping again through autumn (Figure 2.6).   22   Figure 2.5 Meteorological conditions binned by time of day and averaged, for 6 May to 18 Nov, 2015.  Dots on each figure are sized to indicate the relative number of data points averaged in that period. Vertical solid lines mark noon and midnight. The left column includes all data; the right column includes only periods with winds < 5 ms-1. Panels (a) and (b) display net radiation; (c) and (d) show wind speeds; 𝑹𝑯 is presented in (e) and (f); and water and air temperatures are in (g) and (h). The increased variability in the low-wind periods is due to the reduced number of data points available. There is little variability in 𝑻𝑾 due to the high heat capacity of water. Dashed lines in panels (a)-(b) and (c)-(d) mark 0 W∙m-2 and 5 m∙s-1, respectively. 23  2.3.3 Eddy Covariance Data Eddy covariance instruments measured 𝐻𝑆 and 𝐻𝐿 between Aug 26, 2015 (DOY 238) and Oct 13, 2015 (DOY 286), spanning the end of summer stratification and a portion of fall overturn.  Sensible and latent heat measurements for a 25-day period are presented in Figure 2.7.  The sensible heat flux out of the lake decreases during the day as the air temperature warms, even becoming positive at times when warm air heats the water. In general, sensible heat loss increases in magnitude at night when the air cools relative to the water temperature. Minimum and mean sensible heat flux values are -192 W∙m-2 and -23 W∙m-2, respectively (Figure 2.7a). Latent heat flux measurements, shown in Figure 2.7b, indicate that this parameter is consistently negative or close to zero; the minimum recorded value was -257 W∙m-2. The mean heat loss was -47 W∙m-2. Note that due to quality control measures that removed many invalid or suspect data points, these means do not represent the true mean fluxes over the season. However, based on these estimates, latent heat flux from evaporation contributes more to lake cooling than does the sensible heat flux. 2.3.4 Focus Period To highlight the daily patterns in the data, a focus period between 21-28 Sept 2015 (DOY 265-271) was selected for closer examination. The near-surface temperatures and eddy covariance measurements taken during the focus period are presented in Figure 2.8. The data presented in panels (b)-(d) are identical to that in Figure 2.4. In panel (a), temperature measurements from select near-surface temperature sensors at P1 are included (the 0.05 m and 0.15 m data series). The lake is undergoing fall overturn; however, partial or full stratification of the water column occurs in certain conditions. Eddy covariance measurements are presented in panels (e)-(f). The inverse relationship between air temperature and RH is readily apparent in Figure 2.8(b), with RH dropping when temperature rises due to the corresponding increase in the saturation specific humidity, 𝑞𝑠 (see Equations 1.3 and 1.4). On DOY 265 and 267, low winds and high incoming radiation warm and stratify the water surface, temporarily pausing overturn, with nearly the full water column restratifying on DOY 267. Elevated winds and high net radiation on DOY 269 and 270 correspond with increased latent heat flux out of the lake. 24   Figure 2.6 Relative and specific humidity at BML in the 2015 open water season. Data presented is a 24-h moving average to reduce noise. Relative humidity is shown by the solid blue line and specific humidity corresponds with the dashed red line.  25   Figure 2.7  (a) Sensible and (b) latent heat flux measurements from eddy covariance instruments at P1, Base Mine Lake, Sept 7 to Oct 2, 2015 (DOY 250-275). Grey bars indicate nighttime. The dashed lines indicates zero heat flux.   26    Figure 2.8 Atmospheric and limnological conditions at BML during the focus period (Sept 21-28, 2015). The temperature data is the mean of P2 and P3, panel (a). Air temperature and relative humidity, wind speed, and net radiation are in panels (b)-(d), with sensible and latent heat flux measurements in panels (e) and (f). During this period, data from P1 are used. Grey bars correspond with night. The dashed line in panel (c) highlights where wind drops below 5 m∙s-1, while horizontal dashed lines in panels (d)-(f) on 0 W∙m-2 indicate where the heat flux changes sign.  27  Chapter 3: Modelling The methods used to close the heat budget are addressed in this chapter, with particular focus on the methods used to estimate the latent heat fluxes into the lake.  The results from four methods for estimating latent heat flux are compared, and cumulative heat budgets over the open-water season are used to evaluate how the hydrocarbon sheen impacts heat fluxes at the lake. This study will not attempt to model stratification within BML, and will instead focus on how heat flux processes and calculations vary from standard models. This is to decouple heat flux calculations from errors that may arise when estimating mixing.  3.1 Model Setup 3.1.1 Estimation of Heat Budget Components As discussed in Section 1, the heat budget in the open-water period in a basin with limited inflows and outflows consists primarily of fluxes across the water surface: 𝑄 = 𝑅𝑛𝑒𝑡 + 𝐻𝐿 + 𝐻𝑆 (see Equation 1.1). The methods used to estimate each component of the heat budget are listed in Table 3.1 and described further below. Because in-situ instruments measuring incident shortwave radiation, reflected shortwave radiation, and longwave radiation are located at P1,  providing a more reliable estimate of net radiation than could be achieved through modelling, particularly in the absence of detailed cloud cover data (Luo et al., 2010), the net radiation data from this source was assumed to be reliable. The remaining terms of the heat budget (𝑄, 𝐻𝑆, and 𝐻𝐿) were estimated using four different methods: the Air-Sea toolbox with standard parameterization (Standard AST); the extended Air-Sea toolbox with modifications for addressing the effects of the hydrocarbon sheen (Extended AST); in-situ instrumentation, namely EC measurements; and by using the measured lake water temperatures to estimate the change in heat content of the lake (temperature-based estimate, TBE).  Both AST approaches use the MATLAB® Air-Sea Toolbox (Beardsley et al., 1999), code that was adapted from COARE (Fairall et al., 1996b), to calculate 𝐻𝑆 and 𝐻𝐿. Eddy covariance measurements for 𝐻𝑆 and 𝐻𝐿, shown in Figure 2.7, provide a comparison against which model 28  Table 3.1 Methods used to estimate each component of the heat budget, where the Standard and Extended Air-Sea Toolbox (AST) use models to simulate physical processes; eddy covariance (EC) techniques use instruments to measure turbulent fluxes; and a temperature-based estimate (TBE) uses data from the thermistor records. Heat Flux Estimation approaches Net Radiation In-situ instrumentation Sensible Heat Standard AST,  EC Latent Heat Standard AST, Extended AST, EC, TBE Net Heat Flux Standard AST, Extended AST, TBE results can be evaluated; however, they cannot by themselves close the heat budget due to the limited time spanned and gaps in the data.  The temperature-based estimate for the heat fluxes (Equation 3.2) used spatial averages of the measured lake temperature. The spatially averaged mean lake temperatures were obtained by averaging the P2 and P3 sensors at each depth. The average temperature for each depth is then multiplied by the corresponding lake volumes indicated by the hypsometric curve (Figure 2.3); these are summed and divided by the total lake volume for a vertically and horizontally averaged lake temperature.  This approach can provide insight into the net heat flux at the lake surface, but not the magnitudes of each component of the heat budget. However, it can be used in combination with the other approaches to obtain estimates of individual components. For example, by rearranging Equation 1.1, the net heat flux not accounted for by net radiation or sensible heat flux was assumed to be from the latent heat of evaporation (Equation 3.1).     29  Thus, a temperature-based estimate (TBE) of latent heat flux can be obtained if the values of 𝑅𝑛𝑒𝑡 and 𝐻𝑆 can be ascertained using other methods:  𝐻𝐿 𝑇𝐵𝐸 = 𝑄𝑇𝐵𝐸 − 𝐻𝑆𝐴𝑆𝑇 − 𝑅𝑛𝑒𝑡 (3.1) where  𝑄𝑇𝐵𝐸 = 𝜌𝑤𝑐𝑝𝑑𝛥?̅?𝛥𝑡 (3.2) where 𝜌𝑤 is the water density, 𝑐𝑝 is the heat capacity of water, 𝑑 is the water depth, 𝛥?̅? is the change in the mean water column temperature (spatially averaged) over the time step, and 𝛥𝑡 is the duration of each time step. Sensible heat, 𝐻𝑆𝐴𝑆𝑇 is obtained from the Standard AST, and 𝐻𝐿 𝑇𝐵𝐸 is the expected contribution of latent heat to the change in lake heat content as estimated by the thermistor records.  Thus, by comparing the heat flux estimates from each of these sources, and using them in combination in the case of the TBE estimate for 𝐻𝐿, a clearer picture of the heat fluxes in BML is obtained.  At P2 and P3, the temperature sensors are well-distributed in the water column, capturing the vertical temperature gradient of the lake. However, particularly during periods with high stratification, horizontal movement of water in the lake results in uncertainty in the estimated instantaneous lake heat content due to wind mixing and other factors. The analysis was therefore restricted to fall turn over due to the relative homogeneity of the lake in this time.  3.1.2 Cumulative Heat Budget For the Standard AST and Extended AST approaches, each component of the surface heat budget (Equation 1.1) is estimated and summed to yield a net heat flux, 𝑄, for each time step in the 2015 open water period, where each time step is 30 minutes in duration. Springtime up to 15 June 2015 (DOY 166) were omitted from the analysis due to the challenges involved in accurately estimating the heat content of the lake when the system is not yet fully stratified.   30  This net heat flux is used to calculate a mean lake heat content in each time step:  𝐻𝐴𝑆𝑇(𝑡) = 𝐻𝐴𝑆𝑇(𝑡 − 1) + 𝑄𝐴𝑆𝑇𝐴𝑑𝑡 (3.3) which can be compared against the actual mean lake heat content estimated from the thermistors  𝐻𝑇𝐵𝐸 = 𝑐𝑝 ∑ 𝜌(𝑧)(𝑇(𝑧) − 𝑇0)𝐴(𝑧)𝑑𝑧 (3.4) where 𝐻𝐴𝑆𝑇 is the modelled heat content based on the estimated surface heat fluxes, 𝑄𝐴𝑆𝑇; 𝐻𝑇𝐵𝐸 is the heat content based on the thermistor measurements; 𝐴 is the surface area of the lake; and 𝑇0 is the reference temperature (the temperature of maximum density is used here).  This is a 1D model, meaning that the vertical and horizontal variation in the lake will not be modelled, and only a single average lake temperature is estimated for each time step. As the AST models require surface temperatures as an input in the calculations, a surface temperature for each time step is estimated using the thermistor data (refer to Section 3.2.1 for details). As mentioned above, the surface temperature and stratification will not be modelled, in order to decouple errors in this parameter from the main research questions; instead, the temperature is estimated using lake and atmospheric conditions at each time step.  The model volume is bounded by the lake sidewalls, the FFT-water interface at the base, and the air-sea interface at the top. The lake surface was assumed to remain at a constant elevation of 308 masl. The lower boundary of the model, the FFT-surface interface, dropped gradually over the duration of the simulation due to FFT settlement.  3.2 Air-Sea Toolbox Approach The Air-Sea Toolbox (AST) uses meteorological and limnological conditions at the lake surface to iteratively determine atmospheric stability and transfer coefficients at each time step, and thereby estimate the surface fluxes at the lake surface. Refer to Fairall et al. (1996b) and Appendix A for further details. Inputs to the AST include water temperature at the air-water interface, 𝑇𝑆, and meteorological data. As reliable measurements are available for most key input parameters and fluxes, including 𝑇𝐴, 𝑅𝐻, 𝑢, and longwave and shortwave radiation, the main source of uncertainty in the heat budget calculations was 𝑇𝑆. Warm-layer effects causing an increase in temperature beyond what is measured by the shallowest thermistor (𝛥𝑇𝑊, Equation 1.6) must be taken into 31  consideration when there is stratification in the system; cool-skin effects must be accounted for as well (refer to Figure 1.2). The Air-Sea Toolbox includes an algorithm to calculate the temperature difference across the cool skin of the water surface; however, the warm-layer effects must still be accounted for, as the Air-Sea Toolbox does not include an algorithm for these calculations. Here, the shallowest thermistor available throughout the study period is at a depth of 0.5 m.  3.2.1 Determining TS Near-surface temperature models are readily available within the literature (Alappattu, 2017; Noh et al., 2011), and a warm layer model is detailed in Fairall et al. (1996a) and updated in Wick et al. (2005). However, the complexities of these algorithms make them challenging to implement for a full season. The applicability of existing algorithms may also be reduced due to elevated turbidity levels in the lake, as well as the presence of the hydrocarbon sheen on the lake surface. Near-surface temperatures are therefore estimated using a data-driven approach, as outlined in Appendix B.  A skin effect model is included in the Air-Sea toolbox and is used to estimate the difference between the bulk water column temperature and the skin temperature, 𝛥𝑇𝐶 (Equation 1.6). The cool-skin temperature difference is a function of the heat flux through the thin skin layer, which is in turn dependent on the temperature of the water surface. Cool-skin effects therefore require iterative calculation. The Air-Sea toolbox does not consider warm-skin effects; however, this occurs only when there is a reversal in the usual sign of the surface fluxes, usually when the atmosphere is warm and moist (Murray et al., 2000). Over the entire open water period, BML had only 30 hours with high 𝑅𝐻 and air temperature greater than the temperature of the upper water column; for this reason, warm-skin effects were disregarded.  The cool-skin and warm-layer phenomena are of opposite direction and are similar in average magnitude, as described below. However, the timing and magnitude of each of these effects make them important for accurate determination of the heat fluxes across the water surface. The surface of BML is not consistently stratified, meaning the warm-layer adjustment is not always high, with an average value of 0.15 °C (being zero during unstratified periods and averaging 0.5 °C when there is an adjustment). However, the correction increases at elevated water 32  temperatures (Figure B.2). These warmer periods are also when more evaporation would be expected to occur, as the air can hold more moisture in these periods. These adjustments lead to a 0.6 Wm-2 and 3.9 Wm-2 increase in the magnitude of the sensible and latent heat fluxes predicted by the Standard AST.  By comparison, the cool-skin effect overall averages -0.3 °C, peaking at -0.6 °C. Using the Standard AST, accounting for the cool-skin of the water decreases the average magnitude of the sensible heat flux over the period of record by 2.4 Wm-2 (or a 13 % reduction in sensible heat lost to the atmosphere), and decreases the average magnitude of the latent heat flux by 4.4 Wm-2 (or a 5 % reduction in latent heat lost to the atmosphere).  3.2.2 Standard AST Results After correcting the Standard AST inputs to account for warm-layer and cool-skin temperature effects, the Air-Sea toolbox calculated the sensible and latent heat fluxes for each time step. These estimates were used to calculate the heat content of the lake over the study period based on the surface heat fluxes (Equations 1.1 and 3.3) and the thermistor record (Equation 3.4). The results are presented in Figure 3.1a, which indicates that the Standard AST predicts a cooler system than is observed. The difference between the two estimates increased at an average rate of 5.4x109 kJ∙d-1 or 9 W∙m-2 over the study period (dotted line, Figure 3.1b). Note that periods where thermistor records are unavailable were omitted from the data set; the difference between the measured heat content of the lake at the start and end of these periods was also applied to the simulated cumulative heat content in order to account for heat fluxes that occurred during these periods. 3.2.3 Comparison with Eddy Covariance Methods Data from the EC instruments were compared with the Standard AST results with the intent of identifying which fluxes or periods predicted by the Standard AST correspond poorly with measurements.  The agreement between the 𝐻𝑆 estimates from the Standard AST and EC methods is high, with an average difference of -4.0 W∙m-2 and an R2 of 0.6 (Figure 3.2a). Strongest agreement was  33   Figure 3.1 Comparison of BML heat content estimates from thermistor measurements and cumulative fluxes estimated using the Standard AST approach over the open-water period, 6 May to 18 Nov, 2015 (DOY 166-322). Panel (a) presents the cumulative heat budget for the lake, while panel (b) presents a time series of the difference between the Standard AST and the measured heat content. The dashed line on panel (b) indicates where there is no difference between these estimates, while the dotted line shows the mean daily increase over the reporting period.   34   Figure 3.2 Comparison of EC measurements and Standard AST estimates for 𝑯𝒔 (panel a) and 𝑯𝑳(panel b), BML, 26 Aug to 14 Oct 2015 (DOY 238-287). The R2 values indicating goodness-of-fit are 0.6 and 0.4 for 𝑯𝑺 and 𝑯𝑳 respectively, with R2 dropping to 0.1 when considering 𝑯𝑳 values in low wind periods only. Blue data points correspond to low winds (𝒖<5 ms-1), and red points correspond to high-wind periods. The straight line represents a 1:1 fit.  35  observed in periods with moderate wind speeds, with the difference between the two estimates dropping to -0.8 W∙m-2 when wind speed exceeded 5 m∙s-1. In low winds, accord decreased (e.g. DOY 266; DOY 267.5), likely due both to the reduced certainty in the surface temperatures in these conditions, as well as due to the challenges in measuring turbulent fluxes when winds are low (Barr et al., 2013). Despite the reduced correlation in low wind conditions, the sensible heat flux output from the model generally represented the measured fluxes well and was therefore assumed to accurately represent the actual sensible heat fluxes from the lake surface.  Unlike 𝐻𝑆, estimates of 𝐻𝐿 using the Standard AST do not correlate well with measured values (Figure 3.2b). The difference between the EC measurements and the Standard AST estimates is significant for 𝐻𝐿, with an average difference of 46 W∙m-2 and an R2 of 0.4. The RMSE for 𝐻𝐿 was more than double that of the sensible heat fluxes. This difference is likely due in part to the effects of the hydrocarbon sheen on the latent heat flux; however, it is also likely due in part to the additional complexity that arises when estimating latent heat fluxes, as this involves the measurement of both atmospheric turbulence and water vapor fluctuations. For these reasons, and because the eddy covariance data has been updated since this analysis was performed (Sean Carey, McMaster University, personal communication 2020), the latent heat estimates from the EC and other methods will not be compared further herein. 3.3 Extended Air-Sea Toolbox 3.3.1 Adjustments to the Interfacial Relative Humidity Based on the results described in Section 3.2.2, it was determined that the Standard AST does not accurately characterize the lake heat budget. In addition to potentially affecting near-surface lake temperatures by reducing mixing in the upper portion of the water column, the hydrocarbon sheen is expected to cause resistance to evaporation at low wind speeds. In low wind conditions, the hydrocarbon layer on the water surface introduces a resistance to evaporation; above a critical wind speed, gravity waves can form and break up the hydrocarbon layer, reducing the impact of the hydrocarbon sheen, as outlined in Section 1.2.4.   Here, the effects of the hydrocarbon sheen are approximated by a step function with a critical wind speed of 5 m∙s-1. At wind speeds below 5 m∙s-1, where the sheen is expected to be relatively intact across the water surface, 𝑅𝐻𝑊 should be reduced (Figure 3.3).  The Standard AST with 36   Figure 3.3  (Line) Proposed effects of sheen on the 𝑹𝑯𝑾 value at the air-water interface, applied in the Extended AST. (Bar) The hourly wind speed distribution, Sandhill Fen, 15 June to 18 November, 2015 (right axis).  Winds skew towards lower values, only rarely falling below 1 m∙s-1 or exceeding 10 m∙s-1. this reduced 𝑅𝐻𝑊 is here referred to as the Extended AST. The Extended AST was run and cumulative heat fluxes calculated, varying the value of 𝑅𝐻𝑊 each time (full calculations provided in Appendix A). The daily net heat flux out of the lake was estimated using both the thermistor record and the modelled heat flux estimates, and the root mean squared error (RMSE) between these two data sets was calculated over the autumn turnover period (Sept 4 – Nov 18, 2015, i.e. DOY 247-322). The autumn turnover period was selected due to the relative confidence in 𝑇𝑆 and in the average lake temperature during this period. The RMSE values were plotted (Figure 3.4) and the 𝑅𝐻𝑊 value that yielded the minimum RMSE was selected as representative of the reduction in relative humidity at the water interface due to the influence of  37   Figure 3.4 RMSE calculated between the daily Extended AST and TBE estimates of 𝑯𝑳 with varying 𝑹𝑯𝑾 values. Only the fall turnover data set was used for this analysis due to the relative confidence in the surface temperatures during this time.   hydrocarbons. A final value of 92 % RH at the air-water interface during low-wind periods was thereby established for the Extended AST.  Low wind conditions (< 5 m/s) occur 60 % of the time during the period of record (Figure 3.3). At higher wind speeds, the hydrocarbon sheen was assumed to have dispersed such that it no longer affects large sections of the lake, and 𝑅𝐻𝑊 was set at the standard value of 100 % in the Extended AST.  38  3.3.2 Extended Air-Sea Toolbox Results The key outputs from the Air-Sea Toolbox pertinent to this study are the cool-skin effect, sensible heat flux, and latent heat flux. For both the Standard and Extended AST, the cool-skin effects and sensible heat fluxes averaged -0.3 °C and -18 W∙m-2, respectively. The average latent heat fluxes estimated by the Standard AST and Extended AST parameterization were -94 W∙m-2 and -86 W∙m-2, respectively. 3.4 Examining the Overall Heat Budget When the sensible and latent heat estimates from the Standard AST are inputted into the 1D model and heat fluxes calculated cumulatively, a cooler system is predicted than what was observed by the temperature sensors (Figure 3.5a).  Figure 3.5b presents a time series of the difference between the measured (TBE) and modelled (Standard and Extended AST) heat contents of the lake, and indicates that the error between the Standard AST and measured heat contents increases steadily over the course of the study. Once the adjustment to 𝑅𝐻𝑊 is applied in the Extended AST, this seasonal error is reduced and the cumulative heat budget follows the measured lake heat content well. After adjusting for the effect of the hydrocarbon sheen on 𝑅𝐻𝑊, the mean magnitude of 𝐻𝐿 calculated by the Extended AST was reduced by 8 W∙m-2 on average, going from a mean of -94 W∙m-2 to -86 W∙m-2. In wind speeds below 5 m∙s-1, the Extended AST reduced the mean 𝐻𝐿 by -12 W∙m-2 compared with the Standard AST results, reducing it from a mean of -71 W∙m-2 to -59 W∙m-2.  Evaporation is calculated by summing the latent heat fluxes over the period of record, as latent heat flux is directly proportional to evaporation. Over the course of the 2015 open water period, this reduction in 𝐻𝐿 calculated by the Extended AST translated to a reduction in evaporation of 22 % when considering only low wind periods. Because winds are low during 60% of the period of record, the effect during these periods results in an overall decrease in evaporative water loss of 9 % over the entire period of record.  Examining the magnitudes of each heat flux over the open water period highlights the importance of the hydrocarbon sheen on the heat budget. The magnitude of the difference  39   Figure 3.5 Comparison between measured (TBE) and modelled (Standard and Extended AST) BML heat content, 15 June to 18 Nov, 2015. Panel (a) presents the cumulative heat budget for the lake, while panel (b) presents the difference between the modelled and measured lake heat content. The dashed line on panel (b) indicates where there is no difference between the modelled and measured heat contents. 40  between the Standard and Extended 𝐻𝐿 estimates in low wind is nearly as high as the total magnitude of the sensible heat flux out of the lake, which varies from -160 W∙m-2 to +214 W∙m-2 over the period of record but has an average value of -18 W∙m-2, which is on the same order of magnitude as the modelled effect of the hydrocarbon sheen on the latent heat flux. Furthermore, although net radiation varies drastically both seasonally and on a daily basis, the mean net radiation over the period under analysis is +71 W∙m-2; this magnitude is close to the total magnitude of the latent heat flux estimated after accounting for the hydrocarbon sheen. The relative contribution of the latent heat flux is therefore quite high at BML, and the importance of the hydrocarbon sheen is correspondingly high on a seasonal basis, being nearly the same magnitude as the total contribution of the sensible heat flux. 3.5 Heat Flux Comparisons To compare the magnitudes and trends of the component heat fluxes from various estimates, data from the fall turnover period were binned by time of day, then averaged to obtain 24-h flux curves (Figure 3.6). The figure includes the measured air temperature and average water column temperatures, 𝑅𝑛𝑒𝑡 estimates from in-situ longwave and shortwave radiation sensors, 𝐻𝑆 fluxes calculated using the Standard AST, and 𝐻𝐿 flux estimates from the Standard and Extended AST, as well as the TBE method.  The data sets were differentiated into periods of low and high winds (defined as an hourly average wind speed below and above 5 m∙s-1), and restricted to periods where the mean temperatures of the 0.5 m and 1 m thermistors at P2 and P3 were near-uniform (absolute temperature difference no greater than 0.02 °C). The latter measure was taken due to the tendency of the lake surface to stratify, even during turnover (refer to Figure 2.8).  3.5.1 Radiative and Sensible Heat Fluxes Incoming solar radiation peaks shortly after noon and the average net radiation is near-constant at night (Figure 3.6a). Throughout fall overturn, the nights get longer and incoming radiation progressively decreases. The range of hours over which sunset and sunrise occur over the course of the measurement period are shaded on Figure 3.6 to help visualize this process. Sensible heat flux is below zero on average throughout the day, reaching a minimum near dawn when the difference between the water and air temperature is greatest (Figure 3.5 e-f and g-h). In the  41   Figure 3.6 Averaged hourly heat fluxes in autumn overturn in unstratified periods, BML, 4 Sept to 18 Nov, 2015. Dots on each figure are sized to indicate the relative number of data points averaged in that period (fewer in daytime as there are more stratified periods when net radiation is high). Vertical solid lines mark noon and midnight; dashed horizontal lines on panels (a)-(f) indicate where the heat fluxes change sign. Light grey shading indicates the range of times for sunrise and sunset, as this changes rapidly in autumn. Panels a-b present mean 𝑹𝒏𝒆𝒕values. Panels c-d show latent heat flux estimates from the TBE method (Equation 3.1) and the Extended and Standard AST. Panels e-f show the sensible heat flux estimate from the Standard AST. Panels g-h show the mean water column and air temperatures over the autumn overturn period. The figures on the left contains only data where the average wind speed is below 5 m∙s-1.  42   afternoon, the temperature difference between the air and the water approaches zero, causing the sensible heat flux to decrease. 3.5.2 Latent Heat: Overnight Cooling The TBE 𝐻𝐿 estimate suggests that in the first few hours after sunrise, the sign of the latent heat flux is reversed and condensation is occurring. An examination of the lake temperatures suggests that the average measured water temperature remains nearly constant throughout the night, followed by rapid cooling in the morning after sunrise. This trend is contrary to what is expected in this period, as night is typically associated with cooling due to longwave radiation. It is not clear why this happens. One possibility is that overnight mixing in the lake is weak, causing the cooling to not register in the thermistors until morning. The moorings are located in the deepest parts of the lake where the heat content to surface area ratio is high, and they may experience lag in registering the changes in temperature that occur in shallower portions of the lake.   The phenomenon is amplified when considering only low-wind periods and effectively vanishes in high winds, suggesting that the delay in response is indeed caused by insufficient overnight convective mixing to fully mix the lake (Figure 3.6c-d). This implies that the cooling of the lake should occur gradually overnight, rather than appearing around sunrise.  3.6 Summary By measuring the atmospheric conditions and the lake temperature over the open water period, and using these measurements to estimate the heat fluxes using several methods, including in-situ measurements, parameterized models of physical processes, and a combination of the two methods, it was possible to estimate the contributions of each component of the heat budget and the impacts of the hydrocarbon sheen on evaporation. Net radiation and sensible heat flux were determined through in-situ measurements and the Standard AST, respectively, and the latent heat flux was determined to be impacted by the hydrocarbon sheen. The impacts of the hydrocarbon sheen were modelled in the Extended AST by reducing the relative humidity on the water side of the air-water interface.  43  The result of this phenomenon is that at BML in the 2015 open-water season, the latent heat flux is reduced by an average of -12 W∙m-2 in low-wind periods.  This corresponds to a reduction in evaporation of 22 % in these periods, or 9 % over the total analysis period, as low winds occur 60 % of the time across the period of record. The predominant heat fluxes in the open-water season were found to be longwave and shortwave radiation, quantified here as net radiation (mean +71 W∙m-2), balanced by latent heat flux (mean -86 W∙m-2), and sensible heat flux (mean -18 W∙m-2). 44  Chapter 4: Open-Water Period Discussion This section will discuss the results and highlight some of the challenges in closing the heat budget, as well as interesting phenomena, among them the consequences of elevated lake temperatures and reductions in evaporation arising because of the hydrocarbon sheen on BML. The sheen acts as a barrier on the water surface through which water must permeate before it can be evaporated. This was represented here as a reduction in the relative humidity at the air-water interface. This reduction occurs during periods of low wind, as at higher wind speeds, gravity waves form, dispersing the hydrocarbon sheen. 4.1 Physical Processes and Modelling The hydrocarbon sheen introduces uncertainties in several of the physical processes at BML, which in turn impact the heat budget. 4.1.1 Skin Temperature of the Water The presence of a hydrocarbon sheen on the water surface introduces challenges in determining the temperature of the water surface, a critical variable in determining the quantity of moisture that saturated air can hold. The value of 𝑇𝑆 is used for calculations of both 𝐻𝑆 and 𝐻𝐿 and is a key parameter in the overall heat budget, particularly in periods where the air can retain high quantities of moisture.  Surface stratification is typically pronounced in summer, and indeed that is the case at BML. However, it was also observed in periods where surface mixing would be expected, such as during fall turn over. Measurements taken at 0.15 m varied from the 0.6 m sensor by as much as +4.1 °C during turn over, a period where mixing is typically assumed to result in uniform temperatures throughout the water column. The frequent near-surface stratification events observed in this lake in autumn may also be due in part to the effects of the hydrocarbon layer. Suppression of surface waves reported by Gottifredi and Jameson (1967), among others, may also correspond with reductions in momentum transfer across the air-water interface during low-wind periods.  Turbidity in the lake may play a role; Tedford et al. (2019) reported high turbidity levels in BML in 2015, ranging from a minimum of 45 NTU (Nephelometric Turbidity Units) at the beginning 45  of turnover up to a peak of 300 NTU. In experiments conducted on small pools, elevated turbidity levels within this range caused an increase in the temperature of the surface (top 10 mm) layer on both cloudy and sunny days, with the upper layer of the most turbid water found to be 2.8 °C warmer on average than the clear water at the peak of solar radiation (Paaijmans et al., 2008). Similarly, Lee and Rast (1997) found light extinction coefficients to decrease with increased turbidity, indicating that for higher-turbidity water, more light and therefore more heat is absorbed in the upper water column.  The consequences of elevated surface temperatures to the heat budget would impact the sensible heat, latent heat, and longwave radiation terms. Elevated surface temperatures increase the outgoing longwave radiation exponentially: 𝐿𝑊𝑜𝑢𝑡 α  𝑇𝑆4 where 𝐿𝑊𝑜𝑢𝑡 is the outgoing longwave radiation (Imboden & Wuest, 1995). Impacts to sensible and latent heat flux are less straightforward; the sensible heat is dependent on the difference between the air and water surface temperatures. Increasing the lake surface temperature will reduce the temperature difference between the water and air when the air is hot, reducing sensible heat input to the water, but will also increase the temperature difference when the air is cooler than the water, increasing sensible heat loss during these times. Over the period of record, the air is slightly cooler on average than the mean lake temperature, at 13 °C and 15 °C respectively. However, changes to the surface temperature would be most significant during the day when there is higher incoming solar radiation to drive stratification. These dynamics mean the impacts of the hydrocarbon sheen on the water surface cannot be captured with daily estimates of the heat budget; more frequent calculations are needed to capture the daily cycle of temperature fluctuations and resultant impacts to sensible heat flux.   Meanwhile, a reduction in the latent heat flux of evaporation also leads to a reduction in cooling at the lake surface, as latent heat is in most cases negative out of the lake, which may reduce the cool-skin effect or lead the surface layer to warm more than would be expected in a lake of equivalent size and meteorology without a hydrocarbon sheen. Having said that, an overall increase in the surface temperature (warm-skin effect) that the hydrocarbon layer may induce would increase the amount of water that can be held in a given parcel of air, meaning that more moisture is available to be evaporated and more latent heat flux can occur (Equations 1.4 and 46  1.5). This implies by definition an increase in cooling surface flux at the water surface and so an increase in the cool-skin effect, counteracting some of the effects of the increased surface temperature in increasing the latent heat flux. These self-insulating phenomena mean that the impacts to latent heat flux require iteration in any given time step in order to accurately estimate the surface temperature, as described in Fairall et al. (1996a).  4.1.2 Effects of Hydrocarbon Sheen  In the model formulations, it is the difference between the specific humidity in the bulk air and the specific humidity at the water interface (𝑅𝐻𝑊) that controls the amount and direction of latent heat flux (Fairall et al., 1996b; Wunderlich, 1972), and it is postulated here that the hydrocarbon sheen causes a reduction in 𝑅𝐻𝑊, resulting in a decrease in evaporation. A reduction in 𝑅𝐻𝑊  in periods with low winds can be interpreted as a reduction in the rate at which moisture in the air-water interface is replenished, with instantaneous recovery to 100 % relative humidity prevented due to the need to diffuse water through the hydrocarbon layer that occurs at a rate slower than the transfer of moisture through the boundary layer. Note that the reduction in 𝑅𝐻𝑊 could be interpreted as a bulk value across the entire lake. The hydrocarbon sheen is not always present across the entire lake surface, but rather appears in patches that drift around BML, and so regions with more hydrocarbon sheen might experience lower 𝑅𝐻𝑊, with less affected areas closer to typical values.  The magnitude of the latent heat flux increases with increasing wind speed, as would be expected by the governing equations presented in Table 1.1; increasing the wind speed linearly increases the transfer of moisture away from the boundary layer. Accordingly, the magnitude of the hydrocarbon-induced reduction of latent heat flux also increases with increasing wind speed (Figure 4.1a). Thus, the reduction in latent heat flux of -12 W∙m-2 in wind speeds below 5 m∙s-1 translates to a difference of -22 % in low winds when normalized as a percent change from standard estimates (Figure 4.1b). Above 5 m∙s-1 there is no difference due to the hydrocarbon sheen dispersing at elevated wind speeds.   There are outliers in the data set where the percent difference is far higher in magnitude. These generally occur when the estimated latent heat flux is very close to zero.  47   Figure 4.1 Absolute and percent difference (panels and and b, respectively) in 𝑯𝑳 outputs from 15 June to 18 Nov, 2015 from the Standard and Extended AST estimates.  Above 5 m∙s-1, no difference is observed due to wind dispersal of the hydrocarbon layer. Below 5 m∙s-1, the median percent difference between the two estimates is -22 %. The dashed line at 0 % in panel (b) indicates no difference between the standard and extended 𝑯𝑳 estimates.  48  The transfer coefficients for momentum (𝐶𝑑), sensible heat (𝐶ℎ), and moisture (𝐶𝑒) for the BML open water period are presented as a function of wind speed in Figure 4.2(a-c), respectively. These transfer coefficients are used in models to describe the degree to which forcings from the atmosphere are imparted to the water (refer to Appendix A). Early models assumed constant values, which evolved into assumptions that values of 𝐶ℎ and 𝐶𝑒 were identical to 𝐶𝑑 (Wunderlich, 1972; Panofsky & Dutton, 1984).  In the Air-Sea Toolbox, each transfer coefficient is calculated separately, although the profile functions, a related concept, are assumed to be the same for temperature and humidity (Appendix A). As seen in Figure 4.2, the transfer coefficients calculated by the Extended AST are a function of wind speed, among other things. Based on these plots, it appears that 6 m∙s-1 is an inflection point for heat and moisture; above this point, the transfer coefficient decreases. By comparison, the momentum transfer coefficient is at a minimum in moderate wind speeds, increasing in magnitude at low and high winds. This phenomenon is not induced by modelling of the hydrocarbon sheen, being present in the Standard AST as well, and may be instead due to impacts of increased wind speed on the skin layer or other factors, which may reduce the degree of transfer between the atmosphere and water at high wind speeds.  4.1.3 Reduced Evaporation The overall reduction in evaporation estimated using this approach (~10 % over the full period) is lower than that estimated by the water budget, which indicated a reduction in evaporation by approximately 30 % (Carey & Drewitt, 2017). This may be due in part to the uncertainty in the TS estimates, which strongly influence the results (refer to Section 3.2.1). The data set used to determine the effects of near-surface stratification also started in late summer on 29 July (DOY 210), meaning the data used to establish the warm-layer adjustment was biased towards late summer and autumn overturn when the water column was less stratified. Thus, over the full period, the warm-layer adjustment is likely more significant than what is suggested by the above data. Furthermore, the adjustment was based on data obtained from the 0.15 m sensor, as this is the shallowest sensor that provided reliable data. However, in the absence of mixing, warming of the surface would continue up to the depth where the cool-skin effect begins to be observed (Figure 1.2), and so the warm-skin effect in high radiation, low-wind periods would almost  49   Figure 4.2 Transfer coefficients for momentum, heat, and moisture as a function of wind speed. Coefficients were obtained from the Extended AST.    50  certainly be greater than what was predicted here. An increase in the warm-layer effect would alter the resultant calculations of 𝐻𝑆 and 𝐻𝐿 and may therefore lead to a greater reduction in evaporation predicted by the hydrocarbon sheen. Reducing the extent of hydrocarbon impact to the lake surface will have a number of effects on the lake physics. Notably, the evaporative loss of water into the atmosphere would be expected to increase, particularly from mid-summer to early autumn. This process has already begun, with the water budget indicating gradual recovery to expected evaporation rates (Carey & Drewitt, 2017). This may be due to reductions in coverage of the hydrocarbon sheen. Skimming of the sheen at the water surface and dredging of the FFT regions producing high quantities of bitumen in recent years may have led to reductions in the hydrocarbon presence at the water surface. Natural slowing of the hydrocarbon release from the FFT or breakdown of the surface sheen may have also played a role.  The increase in evaporation resulting from reduced hydrocarbon coverage would also result in increased cooling of the system. This change in the heat budget would also be reflected in changes in the lake dynamics. Increased cooling would lead to reduced temperatures, possibly changing the biological activity within the lake, as well as an increase in convective mixing near the surface. Depending on the depth of mixing and the geometry of the basin, this may also mean increases in entrainment from the bottom, as well as strengthened turnover.  It can thereby be seen that the presence of the hydrocarbon sheen acts not only to change the heat flux across the lake surface but affects the overall dynamics of the lake in the open-water period.  51  Chapter 5: Under-Ice Period This chapter presents the field data from the under-ice period and presents a preliminary heat budget for this period. The goal of this work is to estimate the magnitudes of the predominant heat transfer mechanisms at BML in winter, and to discuss the implications of these heat sources and sinks on the under-ice dynamics. The heat budget for the lake is split into two parts, looking first at the top half of the water column, along with ice formation and loss, then examining the exchange of heat with the FFT and the temperature trends in the bottom half of the water column. Measurements and field methods from the 2015-2016 ice-on period, as well as early spring of 2016, will be presented. An overview of the calculation methods and sources of uncertainty for the main heat budget components in winter is then provided. Finally, the results of these calculations and a description of the impact of the heat budget components on the actual behaviour of the lake over winter are addressed.  5.1 Field Observations The thermistors moored at P1 and P2 and the meteorological instruments installed at P1, described in detail in Sections 2.2.1 and 2.2.2, recorded data throughout the ice-on period.  Meteorological and lake temperature data are presented in Figure 5.1. Ice-on occurs on 19 Nov 2015 (denoted as DOY -42, 2016), while ice-off occurs around 19 April 2016 (DOY 110), marked with a dashed line in  Figure 5.1a. Spring conditions will not be considered in the analysis, as this section is focused on winter processes; however, observations of general trends in the spring field data are included below. 5.1.1 Thermistor Data At the end of fall turnover on DOY, the overturn of the water column ceases and inverse stratification sets up rapidly. The water column temperature at ice-on is 1.4 °C; however, once inverse stratification begins in late November 2015 (DOY -42 2016), the bottom sensor rapidly warms to over 3 °C, while the near-surface sensors cool (Figure 5.1a).  The near-surface thermistors (up to 2 m depth) cool from the near-uniform lake temperature of 1.4 °C at the end of turnover to below 0.7 °C over the course of 24 h, and the top 4 m cools to  52   Figure 5.1 BML field data from the 2015-2016 under-ice and spring period. Data presented are 1-h averages. Panel (a) presents the averaged water temperatures recorded at P2 and P3. The 11 m sensor is recording FFT temperatures at P2 at a depth up to 1 m below the water-FFT interface. Panels (b)-(e) display the air temperature, relative humidity, wind speed, and net radiation, respectively. The vertical dashed line at DOY 110 indicates the end of the winter ice-covered period, while the horizontal dashed lines on panels (a) and (b) indicate 0 °C and the horizontal dashed line on panel (e) highlights where net radiation changes sign. 53  below 0.3 °C by 8 Dec 2015 (DOY -23). Prior to the 0.5 m sensor freezing into the ice, the 0.5 m and 1 m thermistors track closely, with a mean temperature difference of 0.07 °C through this period, rising only infrequently above 0.2 °C. Cooling of the upper water column then slows, and the thermistor records at 1 m, 2 m, and 3 m each fluctuate by less than 0.2 °C through the remainder of the ice-on period, with all three remaining below 0.4 °C through this period. The water at the bottom of the lake exhibits very different behaviour from that at the surface; below 6 m depth, the thermistors warm or remain near the water temperature at the beginning of winter (highlighted in Figure 5.2). At the start of the winter record on 19 Nov 2015 (i.e. DOY 323 2015), the thermistors at 6-10 m depth are close in temperature, ranging from 1.3 to 1.4 °C. For the 9 m and 10 m sensors, this is their thermal minimum for the season; the 10 m sensor at P2 rapidly warms, reaching 3.3 °C 9 days later (28 Nov 2015) and reaching a thermal peak of 3.8 °C by 27 Dec 2016. It then cools gradually to a temperature of 2.9 °C by the end of the winter. This thermal peak is above the temperature of maximum density (𝑇𝑚𝑑) at BML; in pure water at standard atmospheric pressure, 𝑇𝑚𝑑 is 3.98 °C. However, this decreases by about 0.2 °C for every part per thousand (ppt) salt content (Chen & Millero, 1986). Given the salinity at BML of approximately 2 ppt, the expected 𝑇𝑚𝑑 of the bottom waters is closer to 3.5 °C. Consequently, the maximum water temperature observed by the 10 m sensor is higher than the temperature of maximum density. Temperatures above 𝑇𝑚𝑑 are recorded by the 10 m sensor sporadically between 16 Dec 2015 and 25 Jan 2016.  The 9 m thermistor warms through the winter as well, starting shortly after the 10 m sensor starts to warm. The warming rate is more gradual than that of the 10 m sensor, with the thermistor rising from a temperature of 1.5 °C to 3.1 °C between 26 Nov 2015 and 8 Jan 2016. Warming of the 8 m and 7 m thermistors is also observed; each warms more slowly and reaches a cooler temperature than the thermistor below it. The 6 m thermistor ranges between 0.5 °C and 1.3 °C, cooling slightly in the early half of the winter and warming again to close to its original temperature in the later half. The 10 m thermistor cools through the latter half of the winter, but the thermistors between 7 m and 9 m do not display a trend between when they reach their thermal maximum and the end of the winter: once they finish warming (on 10 Jan, 20 Jan, and   54   Figure 5.2 Lower lake and FFT temperatures at BML for winter 2015-2016. Thermistor data partially presented in Figure 5.1. The dashed line at 3.5 °C represents 𝑻𝒎𝒅. The dashed line at DOY 110 indicates ice-off. 55  19 Feb 2016 at 9 m, 8 m, and 7 m, respectively), they remain at roughly that temperature through the remainder of the winter (Figure 5.2). The temperature trends of this underlying layer are captured by the 11 m sensor, located up to a metre into the surface of the FFT. This sensor cools rapidly in the early stages of winter from an initial temperature of 8.5 °C to 6.6 °C by the time that the 10 m sensor first crosses 𝑇𝑚𝑑 on 16 Dec 2015. Cooling then slows to a near-linear rate, reaching 4 °C by ice-off. After 13 Jan 2016, the sensor installed at 0.5 m depth undergoes significant variations in temperature. This occurs because by this point in the winter, the ice cover thickness exceeds 0.5 m and so the sensor is frozen into the ice and no longer captures the water column temperature. This sensor offers the clearest insights into the ice temperature and how it changes through the season with shifts in air temperature. Most notably, over the course of 48 hours, the air temperature falls from above freezing to below -30 °C by 29 Feb 2016. Over this time span, the temperature of the 0.5 m sensor also drops sharply, falling from -0.2 °C to -4.4 °C, reaching a minima 6 hours after the air temperature reached its coolest point.  This sensor also captures the daily melting and freezing cycle towards the end of winter; starting on 28 March 2016, the recorded temperature oscillates each day, spiking above 0 °C during periods of high net radiation and dropping again at night. 5.1.2 Meteorological Data Air temperatures in winter are cold, rarely exceeding 0 °C until late March (Figure 5.1b); this warming coincides with increases in net radiation warming the site (Figure 5.1e). Prior to this, net radiation is very low in magnitude, rarely exceeding ±100 W∙m-2. From the start of the ice-on period until 28 March 2016 (DOY 88), when spikes of the 0.5 m sensor begin, indicating daily freezing and cooling of the surface ice, net radiation had an average value of -8 W∙m-2, meaning that solar heating was too low to balance outgoing longwave radiation from the lake surface. After this date, the mean net radiation rose drastically to +100 W∙m-2, rising as high as +700 W∙m-2 by late spring. Relative humidity values were high throughout winter, rarely dropping below 70 % between mid-December and late March (Figure 5.1c). As air temperatures warmed in spring, causing the 56  saturation moisture content of the air to increase (Equation 1.5), the average value of 𝑅𝐻 began to drop, falling as low as 20 % 𝑅𝐻 on May 1 2016 (beyond scope of investigation). Daily patterns are highlighted in Figure 5.3, with under-ice conditions shown in the left panels and spring measurements shown in the right-hand panels. Panels (a)-(b) highlight the very low incoming solar radiation in winter, when net radiation is almost always negative. In winter, wind speeds and relative humidity are near-uniform regardless of time of day (Figure 5.3c and e). In winter, the water column is warmer than the atmosphere, while in spring this is reversed (Figure 5.3g and h). The average water temperature does not change significantly between these two periods due to the relatively short period of warming included in the spring data set (Figure 5.1). The warm air temperatures during the day in spring correspond with decreases in relative humidity (panel f).  Hourly average snow depth data are available from in-situ instrumentation (Figure 5.4). This instrument is located at P1 and records distance from the sensor to the nearest surface; the recorded data are therefore susceptible to error from snow drift and differential melting.  5.2 Heat Flux Calculation Methods Estimation methods for the primary heat fluxes in winter are presented below, with the aim of obtaining a preliminary estimate of the general order of magnitude and timing from each of the various heat fluxes at the lake, starting with the surface fluxes and then addressing processes at the base of the water column. These regions are treated separately due to the distinct behaviours observed in each, as discussed in Section 5.1.1. Surface conduction across the snow-ice cap is an important process in winter. In this model, it is accounted for in Section 5.2.1, as the ice formation model provided by Lepparanta (2015) includes conduction of heat away from the water to the cold atmosphere, as well as conduction of the latent heat of freezing released during ice formation in the model. 5.2.1 Ice Cap Formation and Loss At the latitude of BML, there is very low incoming radiation in winter, and an ice and snow cover on the lake surface reduces the amount of heat lost due to convection, as well as the amount of incoming radiation interacting with the lake. The dominant heat flux is therefore 57    Figure 5.3 Meteorological conditions binned by time of day and averaged. Vertical solid lines mark noon and midnight. The column on the left shows data from the under-ice period (23 Nov 2015-19 April 2016), with the column on the right showing early spring conditions (20 April – May 1 2016). Panels (a) and (b) display net radiation; (c) and (d) show wind speeds; 𝑹𝑯 is presented in (e) and (f); and water and air temperatures are in (g) and (h). Dashed lines in panels (a)-(b) and (g)-(h) mark 0 W∙m-2 and 0 °C, respectively.58   Figure 5.4  Measured hourly average snow depth (m) in 2015-2016. Drops in depth are associated with temporary rises in air temperature above freezing.  related to the growth and decay of the ice cover. An extensive discussion of ice formation and loss mechanisms is presented by Lepparanta (2015). This source presents empirical equations and parameter values for ice growth, sublimation, and melting using meteorological data. The change of state of the ice sheen was then used to estimate the associated heat fluxes.  Ice Growth The formation of ice at the water surface is associated with a heat flux; for every kilogram of ice that forms, energy must be inputted to form the molecular bonds in the ice. The latent heat of fusion of ice (𝐿𝑓) in pure water is 334 kJ∙kg-1, although the latent heat of fusion decreases slightly in saline solutions (Kumano et al., 2007).   59  Lepparanta (2015) provides a model for ice growth under snow:  𝛥ℎ𝑓𝛥𝑡=𝑘𝐿𝑓𝜌𝑖𝑐𝑒𝑇𝑓−𝑇𝐴−𝑏𝑜𝛽𝑠ℎ𝑠+𝑏1 (5.1) where 𝛥ℎ𝑓 is the growth of the thickness of the ice due to freezing for every time step, 𝛥𝑡; 𝑘 is the thermal conductivity of ice; 𝜌𝑖𝑐𝑒 is the density of ice; 𝛽𝑠 describes the insulation efficiency of snow; ℎ𝑠 is the height of the snow; 𝑏𝑜 and 𝑏1 account for radiative and turbulent losses to the atmosphere; and 𝑇𝑓 is the temperature of freezing. The temperature of freezing is a function of the salinity of the water and the water pressure; here, we use the equation of state given by Fofonoff and Millard (1983):   𝑇𝑓 = −0.0575 ∗ 𝑆 + 1.710523𝑥10−3 ∗ 𝑆32 − 2.154996𝑥10−4 ∗ 𝑆2 − 7.53𝑥10−3𝑃 (5.2) where 𝑆 is salinity in parts per thousand and 𝑃 is the pressure in bars. The mean temperature of freezing at BML was -0.11 °C in winter 2015-2016; 𝑇𝑓 increased gradually over the under-ice period as a result of increases in salinity caused by salt exclusion that occurs during ice formation (Bluteau et al., 2017). The effects of salt exclusion are not quantified explicitly here.  Sublimation Sublimation, the process whereby ice vaporizes, or changes to a gaseous state without first passing through a liquid state, requires a very high energy input. This is represented by the latent heat of sublimation, 𝐿𝑠, which has a value of 2,838 kJ∙kg-1 (Datt, 2011). This energy requirement, an order of magnitude greater than the latent heat of freezing (Table 5.1), means that although sublimation may be small in magnitude, the associated energy can be high.  Ice loss due to sublimation is calculated as:  𝛥ℎ𝑠𝛥𝑡=𝜌𝑎𝜌𝑤𝐶𝑒𝑢(𝑞 − 𝑞𝑠) (5.3) where 𝛥ℎ𝑠 is the change in the thickness of the ice due to sublimation and 𝐶𝑒 is the transfer coefficient for moisture (refer to Appendix A). In neutral atmospheric conditions,  𝐶𝑒 is approximately 1.2x10-3 (Lepparanta, 2015).  60  Melting When there is solar radiation or air temperatures sufficiently elevated to cause warming, melting of the ice cap can occur. Melting can be calculated as a function of both of these factors:  𝛥ℎ𝑚𝛥𝑡= −𝑄𝑇+𝑘1𝑇𝑎𝜌𝑖𝑐𝑒𝐿𝑓 (5.4) where 𝛥ℎ𝑚 is the loss of ice thickness due to melting, 𝑄𝑇 is the solar radiation transmitted through the snow and ice cap, and 𝑘1 is a factor describing turbulent exchange. Full details of the derivation are presented in Lepparanta (2015).   The solar radiation transmitted to the surface, 𝑄𝑇, is calculated as:  𝑄𝑇 = 𝛾𝑄𝑠(1 − 𝛼) (5.5) Where 𝑄𝑠 is the incident solar radiation, 𝛾 is an empirical fraction of non-reflected incident solar radiation that is transmitted through the surface, and 𝛼 represents albedo, which is a measure of how much light is reflected from the surface instead of absorbed. Albedo drops during the melting period as the snow cover is lost, and is estimated here using a ratio of the incident and reflected solar radiation.  Ice Cover Evolution Using the equations given above, the snow depth data recorded on site (Figure 5.4), and parameter values on the order of those given by Lepparanta (Table 5.1), the growth and loss of the ice cover over winter from each individual mechanism was estimated, along with the net ice thickness (Figure 5.5). The figure also includes the data point for when the 0.5 m sensor fell below 0 °C, indicating that the ice was approximately 0.5 m thick, for comparison of the model against site observations.  As seen in Figure 5.5, ice growth is most rapid in the start of the season, slowing as the winter progresses. This is due to the insulation that the growing ice provides itself, as the increased ice thickness causes a corresponding decrease in the thermal gradient driving ice growth; this effect is amplified by the snow layer capping the ice, which also typically grows throughout the winter period and is a very effective insulator.   61  Table 5.1 Parameter values used to estimate ice cover growth and loss. Values are taken from Lepparanta (2015), with some modifications. Symbol Definition Value Units 𝐶𝑒 Transfer coefficient for moisture 1.2x10-3 - 𝑏𝑜 Parameter representing radiative losses during freezing -3 - 𝑏1 Parameter representing turbulent losses during freezing 0.1 - 𝑘 Thermal conductivity of ice 2.18 W∙m-1∙°C-1 𝑘1 Coefficient for turbulent exchange of heat with the atmosphere during melting 15 W∙m-2∙°C-1 𝐿𝑓 Latent heat of freezing 333.5 kJ∙kg-1 𝐿𝑆 Latent heat of sublimation 2,838 kJ∙kg-1 𝛽𝑠 Relative insulation efficiency of snow 4 - 𝜌𝑖𝑐𝑒 Density of ice 917 kg∙m-3 𝛾 Fraction of non-reflected incident solar radiation transmitted through the surface 0.5 -  62   Figure 5.5  Modelled cumulative ice growth and loss at BML over winter 2016. The black line indicates the total ice thickness, while the blue, magenta, and red lines represent ice growth, sublimation, and melting, respectively. The star marks the date at which the 0.5 m sensor fell below 0 °C, indicating that the ice was approximately 0.5 m thick.  Here, loss of ice thickness due to sublimation and melting was assumed to be consistent with the models described above. In practice, sublimation would act first on snow covering the ice cap, and melting is a very complex process that does not necessarily involve a linear relationship between heat input and loss of ice thickness. However, these calculations were deemed to be sufficient for the purposes of estimating the heat fluxes associated with these phenomena.  63  Heat Fluxes due to Change of State Having estimated the rate of formation and loss of ice, it is possible to evaluate the heat flux due to freezing, sublimation and melting (𝐻𝑓, 𝐻𝑠𝑢𝑏, and 𝐻𝑚):   𝐻𝑓 = 𝜌𝑖𝑐𝑒𝐿𝑓𝛥ℎ𝑓𝛥𝑡 (5.6)  𝐻𝑠𝑢𝑏 = 𝜌𝑖𝑐𝑒𝐿𝑠Δℎ𝑠Δ𝑡 (5.7)  𝐻𝑚 = 𝜌𝑖𝑐𝑒𝐿𝑓Δℎ𝑚Δ𝑡 (5.8) These equations use the parameters as given in Table 5.1. 5.2.2 Heat Flux from Snow  Precipitation in winter, typically consisting of snowfall, adds mass to the snow and ice cap, and therefore is a factor in the heat budget. The heat associated with precipitation, 𝐻𝑝, is analogous to that presented by Equation 5.11:  𝐻𝑝 = 𝑐𝑝 ∑ 𝜌𝑤(𝑇𝑝 − 𝑇0) ℎ𝑝𝐴𝑑𝑡 (5.9) where 𝑇𝑝 is the temperature of the precipitation and ℎ𝑝 is the snow-water equivalent of precipitation. The precipitation temperature can be assumed to be the same as 𝑇𝐴; however, the precipitation volume is more difficult to ascertain, due to the lack of snowfall density data. Here, the daily snowfall was converted into an equivalent water volume by dividing by 12, the mean snow-water equivalent ratio in the climate zone where BML is located (Sturm et al., 2010). The heat flux and associated impacts on the ice-snow cap due to rain are not incorporated into this model. 5.2.3 Basal Heat Fluxes Conduction and advection of heat into and out of the lake occurs at the base of the water cap, as described below.   64  FFT Heat Source At the bottom of the water cap is the FFT deposit that was placed in the pit prior to commissioning. This layer, described in detail in Section 2, is releasing water into the water cap. Although having an influx of fluid to a water body from the base of the basin is not dissimilar to having a net groundwater input, the elevated temperature of the FFT when it was deposited (at 11-19 °C) relative to the average annual temperature (~4 °C in 2015-2016), combined with the high heat capacity of the FFT, means that this may be a lingering source of heat to the water body. The surface of the FFT also warms through summer and early autumn until the temperature of the bottom of the water column drops below the temperature of the surface of the FFT. Advection from FFT The heat flux from advection of warm water into the basin can be expressed as:  𝐻𝐹𝐹𝑇 = 𝑐𝑝 ∑ 𝜌(𝑡)(𝑇𝐹𝐹𝑇 − 𝑇0) 𝑉𝐹𝐹𝑇𝑑𝑡 (5.10) where 𝑉𝐹𝐹𝑇 is the incoming fluid volume from the FFT; and 𝑇𝐹𝐹𝑇 is the temperature of the water released from the surface of the FFT. Here, the 11 m thermistor is used as it is located within the surface of the FFT. The volume of water entering the lake from the FFT is based on the water budget (Table 2.1).   Conduction from FFT Due to the relatively high temperature difference between the FFT and the lake water, particularly near the beginning of winter, conduction of heat from the warm FFT layer is also an important factor. The heat flux associated with conduction from sediments below the lake basin, 𝐻𝑠𝑒𝑑, can be expressed as:  𝐻𝑠𝑒𝑑 = −𝑘𝑠 (𝑑𝑇𝑠𝑑𝑧)𝑤𝑠 (5.11) where 𝑘𝑠 is the thermal conductivity of the sediments, and (𝑑𝑇𝑠𝑑𝑧)𝑤𝑠 is the sediment temperature gradient at the basin-water interface (Fang & Stefan, 1994). Here, the sediment temperature gradient is not precisely known, as there is only one thermistor within the FFT body. However, 65  an estimate of the gradient can be obtained by taking the difference between the deepest thermistor in the water column and the thermistor located in the FFT.  In a study described in Dompierre and Barbour (2016), the thermal conductivity of the FFT material in a laboratory setting ranged between 0.7 W∙m-1∙K-1 and 1.2 W∙m-1∙K-1. In the field, it was predicted that the thermal conductivity near the FFT surface would decrease somewhat to be closer to the thermal conductivity of water (Dompierre and Barbour, 2016). Here, a value at the bottom end of the laboratory range, 0.7 W∙m-1∙K-1, was selected. Due to uncertainty in both the temperature gradient within the sediments and the thermal conductivity of the FFT, this calculation may not be highly accurate, however, it can provide order-of-magnitude estimates. 5.3 Winter Heat Budget Results In winter, heat fluxes acting at the top and bottom of the water column lead to mixing, stratification, and other physical processes in each of these regions, which will be described further below.   5.3.1 Upper Water Column Governing Processes Near the water surface, the water body cools during the main ice formation period. Ice formation is initially rapid but slows after January 1, 2016 due to insulation of the ice formation surface from the atmosphere as the ice-and-snow cap thickens (Figure 5.5). During the rapid ice-formation period, the shallowest sensors at 0.5 and 1 m are close in temperature, as discussed in Section 5.1. Once ice growth slows, the thermistors remain cool and do not vary significantly in temperature. This steady temperature in the upper water column results from the heat fluxes and associated physical processes and properties of water. Because water density is nonlinear, reaching a maximum near 4 °C and decreasing above and below this point, lakes in regions with cold winters become inversely stratified (Whipple, 1898). Thus, cold water near 0 °C is buoyant, stabilizing the water column and floating near the cold ice surface, thereby reducing both heat loss to the atmosphere and temperature fluctuations in the upper water column.  Other density effects related to freezing also play a role in under-ice dynamics and heat transfer: salt exclusion during ice formation may help to mix the water surface, as the resulting saline 66  waters are denser, and therefore heavier, than the surrounding waters, causing them to fall through the water column. This may contribute to the low temperature difference between the 0.5 m and 1 m sensors during the main ice formation period.  The presence of the ice cap, in addition to insulating the water below it, is also by definition a solid mass and so unable to mix with the fluid below it. Fluctuations in the temperature of the ice cap, observed by the 0.5 m sensor after 13 Jan 2016, therefore cannot immediately be transferred to the water below via mixing, like they would be in the open water period. Rather, they must be transferred into the water column by convection, conduction, or a change of state. Similarly, snow falling on the frozen lake surface, despite being associated with a heat flux, is isolated from the water body by the ice cap and so the main body of water does not register the heat flux associated with snow accumulation until melting begins in spring. Thus, in the upper portion of the lake, the majority of the heat fluxes act on the ice cap; it is here that most of the temperature fluctuations through winter are observed, and the size of the ice cap results directly from changes of state due to sublimation, melting, and freezing described by the heat budget.  Heat Flux Estimates The cumulative heat fluxes from each source are presented in Figure 5.6, below. As discussed above, the most significant heat fluxes relate to the formation and loss of the ice cap. The freezing of water into ice, caused by the cold atmosphere above, releases heat, yielding the positive sign in 𝐻𝑓. This is partially offset by latent heat of sublimation through the winter period. Although the mass of matter lost by sublimation is a fraction of the ice cap thickness, the heat flux associated with it is disproportionately higher due to the high energy input required to sublimate frozen water into water vapour. Here, the model estimates a loss of 8.4 cm of ice over the winter period, corresponding with a cumulative heat loss of -9.3x104 kJ∙m-2. In contrast, although the ice sheet was modelled to grow to nearly ten times that thickness, with a total growth of 0.72 m, the magnitude of the heat flux was not correspondingly large, with a cumulative flux of 20x104 kJ∙m-2 by the end of winter, only slightly more than double the heat loss from sublimation. 67   Figure 5.6  Winter heat fluxes at the water surface, including the latent heat of freezing (𝑯𝒇), latent heat due to sublimation (𝑯𝒔𝒖𝒃), latent heat of melting (𝑯𝒎), heat flux from precipitation (𝑯𝒑), and heat flux from net radiation (𝑹𝒏𝒆𝒕). Panel (a) shows the cumulative heat gained and lost from each source through the winter, while panel (b) shows the daily average heat fluxes. The black lines are the sum of the individual heat fluxes at the water surface, and the grey lines are the heat fluxes as estimated by the thermistors. Fluxes are smoothed over 5 days to reduce noise.   68  Although the direct impacts of snowfall on the heat budget are quite low, with the cumulative heat flux due to snowfall below 1x104 kJ∙m-2 over the winter (Figure 5.6a), the most substantive effect of snow on the winter heat budget lies in its insulation properties. This is impactful for several of the fluxes. Sublimation, discussed above, is calculated in terms of the ice thickness equivalent. However, in practice, sublimation of the snow cover would occur before sublimation of the ice cap itself, as this is the surface in contact with the atmosphere. The magnitude of 𝐶𝑒 (Equation 5.3) may be affected due to differences in the properties of ice and snow, such as surface roughness. Furthermore, in windy conditions, sublimation may occur on blowing snow above the ground surface, potentially impacting local humidity and temperatures and therefore reducing sublimation at the surface of the snow-and-ice cap (Taylor, 1998). In addition to affecting sublimation, snow traps air, reducing convective heat transfer with the atmosphere. This slows freezing, reduces the amount of radiation penetrating the ice and snow cap, and otherwise forms a barrier between the atmosphere and the lake. For these reasons, the snow cap introduces significant uncertainty into the winter heat fluxes. These impacts are difficult to quantify in the absence of detailed snow data, as the density and thickness of the snow layer play a large role in the amount of heat transfer (Armstrong & Brun, 2008; Equation 5.1), and these properties can vary spatially and temporally, being affected by wind, air temperature, incoming solar radiation, and other atmospheric factors. Solar and longwave radiation result in a net loss of energy through the majority of the ice-on period, becoming a net contributor of heat only at the end of the season. The heat flux due to melting begins to appear at around the same time as net radiation is increasing, in part due to the importance of incoming solar radiation on melting (Equation 5.4). Net radiation is relatively straightforward to capture with instrumentation; in contrast, melting is a very complicated phenomenon. Here, the heat flux due to melting is estimated using a parameterized model to calculate loss of thickness in the ice sheet due to melting, and calculating the associated heat flux (Equations 5.4, 5.8). However, in practice, melting is not such a straightforward process. Heat input causing melting does not generally yield a proportionate loss in ice thickness; rather, it may go first towards warming the ice sheet to a temperature of zero degrees. Subsequent ice melting releases liquid water which can move and percolate through the ice cap, and this in concert with the crystalline structure of the remaining frozen ice can produce ice candling, a common ice 69  decay structure (Grenfell & Perovich, 2004). However, although estimating melting accurately is challenging, the overall magnitude and general timing can be assessed using the equations provided above.  Overall, the net surface heat flux calculated from the individual surface fluxes, indicated by the thick black line in Figure 5.6a, indicates very little heat flux in the winter. Some overestimation of the total cumulative heat flux is observed, particularly in early winter. During this time, the P2 and P3 thermistors at the top of the water column (0.5 m – 5 m) are continuing to cool, and these measurements can be used to obtain a temperature-based estimate of the net heat flux at the water surface using Equations 3.1 and 3.3.  The results are shown in Figure 5.6a (grey line). This estimate is imprecise due to a lack of temperature measurements in the ice cap; however, general observations can still be drawn.  The positive heat flux predicted by the modelled surface fluxes at the beginning of the ice-on period is not observed by the thermistors. This discrepancy could be due to a number of factors. Exchange of heat with the bottom waters due to falling salt plumes from salt exclusion may play a role, particularly in early winter when the water column is less stratified and the ice sheet is freezing rapidly. Alternatively, errors in estimating the heat flux from the freezing of the ice sheet may help explain this discrepancy. The freezing of the ice sheet, as presented by Lepparanta (2015), considers conduction of heat through the ice, which both drives ice formation and releases the resultant latent heat of freezing back into the atmosphere. The model assumes the ice temperature profile is linear, ignoring the thermal inertia of the ice. The ice growth model also requires the use of a number of parameters (Table 5.1), which in this model did not vary in time; however, in reality, many of these may evolve over the course of winter. Inaccurate estimation of the ice formation and loss may also lead to further inaccuracies, as thinner ice would be less effective at insulating the water from the atmosphere, meaning higher heat flux across the ice layer than what is predicted (or vice-versa for thicker ice estimated by the model than what is observed).    70  5.3.2 Lower Water Column Governing Processes The warming of the bottom of the water column is driven by the underlying warm FFT which conducts and advects heat into the lake. Due to inverse stratification in the lake, which stabilizes the water column, surface heat fluxes may be somewhat restricted to the upper half of the lake with more limited heat flux between the waters above and below 6 m depth after stratification sets up. Inflows and outflows of water in the winter are also minimal, as pumping ceases during the ice-on period (Table 2.1). Conduction and advection of heat from the FFT are therefore likely the only sources of heat in the lower water column.  The heat from the bottom acts first on the water closest to the FFT, as evidenced by the rapid warming in the bottom thermistor. As this water warms, it approaches 𝑇𝑚𝑑, becoming denser, and so introduces vertical structure into the water body and limits convective mixing with the water above it. However, once the base of the water column reaches 𝑇𝑚𝑑 and continuing heat influx from the bottom warms it above 𝑇𝑚𝑑, making it more buoyant than the water above, increased mixing is expected to occur as the warm, buoyant waters rise to a location of neutral density. This process is observed in the 10 m sensor, which rises above 𝑇𝑚𝑑 for brief periods after 16 Dec 2015 (DOY -15) but does not remain above 𝑇𝑚𝑑 overall, with warming of the thermistors above it observed instead. Warming of these thermistors likely occurs both by convective mixing and molecular diffusion; the precise mechanism is beyond the scope of this thesis. Heat Flux Estimates On average, conduction of heat imparts 1.9 W∙m-2 of heat through the winter, and advection of water from the FFT introduces a further 0.2 W∙m-2. These fluxes, although small in magnitude, have substantial effects on the water column temperature in winter as described above due to the limited volume on which they act and the lack of other fluxes in this time.  In the 10-day period immediately following ice-on, the incoming flux of heat acts mostly on the very bottom waters; as discussed above, as they warm to 𝑇𝑚𝑑, they become denser and so mixing of this heat upwards in the water column is limited until the bottom reaches 𝑇𝑚𝑑. During this 71  period, the mean flux from the bottom is higher due to the higher temperature gradient between the FFT and the water column, with an initial flux of 5 W∙m-2 and average flux of 3.1 W∙m-2. From back-of-the-envelope calculations, this translates to a warming of 1.4 K in the bottom 0.5 m of the water column, the same warming as is seen in the 10 m thermistor.  After the bottom water reaches 𝑇𝑚𝑑, further increases in temperature lead to more buoyant water, leading to mixing upwards through the water column and an effective increase in the volume of water being affected by heat fluxes from the bottom. This is reflected in Figure 5.7a, which shows the cumulative heat flux in the water bottom over the winter period both using the 10 m thermistor, and using the average of the bottom 4 m of the water column (P2 thermistors at 6-10 m). It also shows the incoming conductive and advective heat fluxes from the FFT. Although the cumulative incoming convective heat flux over the season is greater than the measured change in heat in even the warmest sensor (i.e. 10 m), the cumulative fluxes do not reflect the rapid warming observed at the bottom of the water column in the early period, suggesting instead a slow but steady heat input. This is due to the volume across which the fluxes are acting.  Applying the incoming fluxes to only the bottom 0.5 m of the water column in the first 10 days of winter when the 10 m thermistor is warming most rapidly, the resultant calculated change in temperature is comparable to that measured by the 10 m thermistor (estimated warming of 1.6 °C and 1.4 °C, respectively). It can thus be inferred that the effective volume being warmed by the bottom fluxes is limited to the deepest 0.5 m of the water column. As these waters approach 𝑇𝑚𝑑 and increased warming leads to increased buoyancy and therefore upwards mixing of heat in the water column, the volume of water being affected by incoming heat from the FFT increases.  This principle is further confirmed by the black line on Figure 5.7a; this shows the mean change in heat content in the bottom 4 m of the water column. This estimate shows far slower warming than that estimated for the 10 m water, reflecting the time needed for the heat from the bottom to reach the sensors further from the FFT.  Across the course of the winter, the two estimates converge as heat is transferred upwards, though the thermistor closest to the FFT remains warmer throughout, reflective of winter inverse stratification. The cumulative estimated change in heat content in the bottom 4 m is notably less  72   Figure 5.7 Heat fluxes into the lower water column, and heat fluxes in the bottom based on the measured change in heat content. The point of zero flux is indicated by a solid black line. The daily measured fluxes in panel (b) are smoothed over 10 days to reduce noise.  73  than the incoming heat fluxes. Measurements and Equation 3.1 were used to calculate an effective incoming heat flux of 0.3 W∙m-2. As mentioned above, the incoming heat flux from the FFT was 2.1 W∙m-2 (average daily fluxes shown in Figure 5.7b). As there are no other notable fluxes across the bottom boundary, this discrepancy is reflective of heat transfer with the upper half of the water column. The heat from the bottom could partially account for the difference between the calculated surface heat fluxes, which did not account for heat transfer with the bottom lake waters, and the measurements. The calculated fluxes at the ice-atmosphere interface indicated an average heat flux of -1.2 W∙m-2, cooler than the measured change in heat content in the upper 6 m of BML, which indicated a flux of -0.4 Wm-2 across the winter ice-on period. Incorporating the heat flux from the bottom half of the water column would therefore reduce the amount of cooling predicted by the model; however, the higher magnitudes of the individual lake-atmosphere heat fluxes, as well as uncertainty regarding the impact of the surface ice-and-snow cap, means that the surface heat budget cannot be closed with this heat flux alone. 5.4 Winter Summary Overall, the calculated heat fluxes across the lake surface and bottom are very low over the winter ice-on period, averaging at -1.2 W∙m-2 and +2.1 W∙m-2, respectively. This is reflective of the insulating ice cap on the water surface which limits the amount of heat exchange with the atmosphere. However, individual fluxes, particularly related to the ice cap formation and loss, are much higher in magnitude. This, in conjunction with the unique properties of water and its  resulting changes in density with salinity and temperature, means that under-ice mixing and heat transfer are more dynamic than would be predicted by the average heat fluxes.  The role of salt exclusion in mixing the water column and transferring heat was not quantified here. However, as detailed in Olsthoorn et al. (2020), ice formation from low-salinity water can lead to localized downward transport of salt, and these plumes would also transport heat as they fall. The consequences of this mixing mechanism on the under-ice heat budget warrants investigation, in particular the effect of falling plumes on homogenizing the near-surface waters and the depth to which heat may be transferred downward. 74  Although the FFT is a persistent source of heat to the overlying water column, the magnitudes of the fluxes are quite low, peaking at 5 W∙m-2 when the temperature difference between the water column and the FFT is highest, and averaging +2.1 W∙m-2. Despite the relatively low magnitude of the bottom fluxes, the warming of the bottom of the water column yields rapid onset of inverse stratification, which persists throughout the ice-on period, helping to stabilize the water column.  75  Chapter 6: Conclusion This document presents the methods used to develop a heat budget for a tailings-affected lake in the open-water and under-ice periods, as well as results and analysis.  6.1 Open-Water Period This study examined the heat budget of BML, which was of particular interest due to the presence of a hydrocarbon sheen at the water surface, the source of which is an FFT deposit at the base of the water column. This sheen introduces unique challenges in closing the heat budget and is of interest because of its possible implications on lake dynamics and management. Field observations of atmospheric conditions and water made at BML in the open-water period are presented in Chapter 2. This is followed by a description of the heat budget component estimation methods and the impacts of the hydrocarbon sheen in Chapter 3 and a discussion of key open-water period influences and considerations in Chapter 4.  The sheen introduces a barrier to evaporation on the water surface. In this thesis, this barrier was incorporated into an existing parameterized physical model of the system, the Air-Sea Toolbox, by reducing the relative humidity on the water side of the air-water interface (𝑅𝐻𝑊). This parameter is typically set at 100 % in atmospheric exchange models, indicating a ready source of water that can be transported across the boundary layer and into the atmosphere. In the presence of a hydrocarbon barrier, available moisture is reduced. Here, closing the heat budget demanded a reduction in 𝑅𝐻𝑊 to 92 % in low-wind periods. At higher wind speeds, 𝑅𝐻𝑊 is not expected to be affected, as waves form on the water surface and the sheen is dispersed, leading to 100 % saturation at the air-water interface in these conditions. The reduction in 𝑅𝐻𝑊 led to a reduction in the magnitude of the latent heat flux by an average of 12 W∙m-2 in low-wind periods.  Low winds occur 60 % of the time across the period of record. This is a substantial portion of the surface heat fluxes, with seasonal contributions from each heat source averaging at +71 W∙m-2 for net radiation, balanced by latent heat flux (mean -86 W∙m-2), and sensible heat flux (mean -18 W∙m-2).  The reduction in the magnitude of the latent heat flux has a number of implications on the heat budget and the dynamics of the lake. The reduction in available moisture at the air-water 76  interface translated to a reduction in evaporation of ~20 % in low wind speeds, or ~10 % overall. This change would be reflected in the water budget at the lake, with more of the lake volume retained each season than would otherwise be expected. Reduced latent heat flux also implies reduced cooling across the water surface, leading to a warmer lake on average, particularly in the epilimnion. This further affects the heat budget, as several of the major heat fluxes are dependent on the surface temperature of the lake. The mixing dynamics and biological activity of the lake are also influenced. The impacts of hydrocarbon sheens on the heat budgets and mixing dynamics of bodies of water are rarely studied and not fully understood. Here, the impacts to surface temperatures and resultant impacts on mixing are decoupled from the latent heat flux estimates by estimating surface temperature separately using empirical observations made at the lake. Further research on the impacts of a hydrocarbon sheen on the water surface would benefit from measurements of the lake surface temperature to confirm the impacts of the hydrocarbon sheen on near-surface stratification, as this was a major source of uncertainty in this work. Research to model stratification and mixing of a water column below a hydrocarbon sheen should take the reduction in latent heat flux into consideration, as this was found to be substantially reduced at BML in the 2015 open-water season.  6.2 Under-Ice Period Field data collected during the under-ice period are presented in Chapter 5, along with the methods and results for a winter heat budget. The calculations for the winter heat budget offer a rough estimate of the importance of different heat fluxes in winter, and provide a first glimpse at the under-ice dynamics within a tailings-affected lake. During this time, the hydrocarbon sheen is frozen into the ice surface and covered by snow and so does not affect the heat budget as it does in the open-water period; however, the tailings deposit at the base of the water column is a major source of heat to the lake. At the base of the water column, conduction of heat from the warmer underlying FFT and advection of warm water from settling of the FFT layer inputs 2.1 W∙m-2 on average to the water column. This flux is greatest at the beginning of winter, when the temperature difference between the two layers is at its peak. Although the magnitude of the flux is low relative to what is seen in summer, it causes a 77  rapid warming of the overlying waters, with the warm layer gradually expanding as the bottom reaches the temperature of maximum density and heat is transferred upwards.  In the upper water column, the changes of state creating and degrading the ice cap, whether by freezing, sublimation, or melting, provide the bulk of the heat fluxes in winter, as net radiation is low. Heat transfer with the cold overlying atmosphere causes the top surface of the water column to freeze, a process that both releases heat and insulates the water column, limiting further cooling of the underlying waters.  Surface heat fluxes were estimated using existing parameterized models of ice growth and formation. Although these are useful for providing initial estimates, the presence of snow cover greatly increases uncertainty in the model, as the insulation imparted by this layer varies by orders of magnitude depending on its condition, which can change on a daily or even hourly basis in winter. Thus, even occasional measurements of snow cover thickness and density during the ice-on period would enable increased accuracy in the parameters used in the winter flux models.  This work highlights interesting under-ice dynamics brought about due to heat fluxes, which should be examined in greater detail. The moderate salinity of the lake waters and its potential effects on mixing from salt exclusion should be examined in further detail to verify its significance. Meanwhile, warming from the FFT layer that brings the bottom of the water column to, and occasionally above, the temperature of maximum density would be of interest in examining instabilities in the water column and contribute to the understanding of under-ice lake dynamics.     78  Bibliography Anikiev, V. V., Mishukov, V. F., Moiseevsky, G. N., & Tkalin, A. V. (1988). The effect of oil films on water evaporation and oxygen content in sea water. Geojournal, 16(1), 19-24. doi:10.1007/BF02626367 Alappattu, D. P., Wang, Q., Yamaguchi, R., Lind, R. J., Reynolds, M., & Christman, A. J. (2017). Warm layer and cool skin corrections for bulk water temperature measurements for air‐sea interaction studies. Journal of Geophysical Research: Oceans, 122(8), 6470-6481. doi:10.1002/2017JC012688 L. Armstrong, R., & Brun, E. (2008). Snow and climate : Physical processes, surface energy exchange and modeling. Cambridge University Press. Aubinet, M., Vesala, T., & Papale, D. (Eds.). (2012). Eddy Covariance: A Practical Guide to Measurement and Data Analysis. Dordrecht, the Netherlands: Springer Atmospheric Sciences.  Beardsley, B., Dever, E., Lentz, S., Edson, J., Payne, D (1999). AIR_SEA Toolbox (Version 2.0) [MATLAB®  toolbox]. https://www.mathworks.com/ Barnes, G. T. (1993). Optimum conditions for evaporation control by monolayers. Journal of Hydrology, 145(1), 165-173. doi:10.1016/0022-1694(93)90225-X  Barr, A. G., Richardson, A. D., Hollinger, D. Y., Papale, D., Arain, M. A., Black, T. A., ... & Law, B. E. (2013). Use of change-point detection for friction–velocity threshold evaluation in eddy-covariance studies. Agricultural and Forest Meteorology, 171, 31-45. Blanc, T. V. (1985). Variation of bulk-derived surface flux, stability, and roughness results due to the use of different transfer coefficient schemes. Journal of Physical Oceanography, 15(6), 650-669. Bluteau, C. E. (2006). Mixing in Brackish Lakes due to Surface Ice (University of British Columbia).  Bluteau, C. E., Pieters, R., & Lawrence, G. A. (2017). The effects of salt exclusion during ice formation on circulation in lakes. Environmental Fluid Mechanics, 17(3), 579-590. Buck, A. L. (1981). New equations for computing vapor pressure and enhancement factor. Journal of Applied Meteorology, 20(12), 1527-1532. doi:10.1175/1520-0450(1981)020<1527:NEFCVP>2.0.CO;2  Businger, J. (1986). Evaluation of the accuracy with which dry deposition can be measured with current micrometeorological techniques. Journal of Climate and Applied Meteorology, 25(8), 1100-1124. Businger, J. A. (1988). A note on the Businger-Dyer profiles. Boundary-Layer Meteorology, 42(1-2), 145-151. doi:10.1007/BF00119880  Campbell Scientific. (2019, August 27). 107 Temperature Probe Product Brochure. Carey, S. (2017). A 2016 Water Balance for Base Mine Lake and Revised 2014 and 2015 Water Balances. N.p.: McMaster University. 79  Carey S. K., & Drewitt, G.B. (2017). Multi-year water balance and gas exchanges from an end pit lake, Fort McMurray, Alberta. Presented at the Canadian Geophysical Union Annual General Meeting, Vancouver, Canada.  Chang, S., Tedford, E. W., Lawrence, G., & Olsthoorn, J. (2018, February). Winter stratification and turbulent mixing in Base Mine Lake. In 2018 Ocean Sciences Meeting. AGU. Charnock, H. (1955). Wind stress on a water surface. Quarterly Journal of the Royal Meteorological Society, 81(350), 639-640. Chen, C. A., & Millero, F. J. (1986, May). Precise Thermodynamic Properties for Natural Waters Covering Only the Limnological Range. Limnology and Oceanography, 31(3), 657-662. COSIA. (2012). Technical Guide for Fluid Fine Tailings Management. Retrieved April 10, 2020, from http://www.cosia.ca/initiatives/project-research COSIA. (n.d.). Pit Lake Research. Retrieved April 10, 2020, from https://www.cosia.ca/initiatives/water/projects/pit-lake-research Datt P. (2011) Latent Heat of Sublimation. In: Singh V.P., Singh P., Haritashya U.K. (eds) Encyclopedia of Snow, Ice and Glaciers. Encyclopedia of Earth Sciences Series. Springer, Dordrecht Dompierre, K. A.  (2016). Controls on mass and thermal loading to an oil sands end pit lake from underlying fluid fine tailings (University of Saskatchewan Saskatoon Canada). Dompierre, K. A. (2016). Controls on mass and thermal loading to an oil sands end pit lake from underlying fluid fine tailings (Master's thesis).  Dompierre, K. A., & Barbour, S. L. (2016). Characterization of physical mass transport through oil sands fluid fine tailings in an end pit lake: a multi-tracer study. Journal of Contaminant Hydrology, 189, 12-26.  Dompierre, K. A., Barbour, S. L., North, R. L., Carey, S. K., & Lindsay, M. (2017, May 17). Chemical mass transport between fluid fine tailings and the overlying water cover of an oil sands end pit lake. Water Resources Research, 53(6), 4725-4740. Donlon, C. J., Minnett, P. J., Gentemann, C., Nightingale, T. J., Barton, I. J., Ward, B., & Murray, M. J. (2002). Toward improved validation of satellite sea surface skin temperature measurements for climate research. Journal of Climate, 15(4), 353–369. https://doi.org/10.1175/1520-0442(2002)015<0353:TIVOSS>2.0.CO;2 Donlon, C. J., Nightingale, T. J., Sheasby, T., Turner, J., Robinson, I. S., & Emergy, W. J. (1999). Implications of the oceanic thermal skin temperature deviation at high wind speed. Geophysical Research Letters, 26(16), 2505-2508. doi:10.1029/1999GL900547  Donlon, C. J., & Nightingale, T. J. (2000). Effect of atmospheric radiance errors in radiometric sea-surface skin temperature measurements. Applied Optics, 39(15), 2387. https://doi.org/10.1364/ao.39.002387 80  Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., Edson, J. B., & Young, G. S. (1996a). Cool‐skin and warm‐layer effects on sea surface temperature. Journal of Geophysical Research: Oceans, 101(C1), 1295-1308. doi:10.1029/95JC03190 Fairall, C. W., Bradley, E. F., Rogers, D. P., Edson, J. B., & Young, G. S. (1996b). Bulk parameterization of air‐sea fluxes for tropical Ocean‐Global atmosphere Coupled‐Ocean atmosphere response experiment. Journal of Geophysical Research: Oceans, 101(C2), 3747-3764. doi:10.1029/95JC03205 Fang, X., & Stefan, H. G. (1996). Dynamics of heat exchange between sediment and water in a lake. Water Resources Research, 32(6), 1719-1727. doi:10.1029/96WR00274 Fischer, H. B.,. (1980). Mixing in inland and coastal waters. London, England; San Diego, California;: Academic Press, Inc.  Fitzgerald, L. (1963). Wind-induced stresses on water surfaces: A wind-tunnel study. Australian Journal of Physics, 16(4), 475. doi:10.1071/PH630475 Fofonoff, N. P., & Millard Jr, R. C. (1983). Algorithms for the computation of fundamental properties of seawater. Foken, T. (2008). The Energy Balance Closure Problem : An Overview. Ecological Applications, 18(6), 1351–1367. Francey, R. J., & Garratt, J. R. (1978). Eddy flux measurements over the ocean and related transfer coefficients. Boundary-Layer Meteorology, 14(2), 153-166. doi:10.1007/BF00122616  Garratt, J. R. (1994). The atmospheric boundary layer. Earth-Science Reviews, 37(1-2), 89-134. Geernaert, G. L., & Plant, W. J. (1990). Surface waves and fluxes. Norwell, MA, U.S.A;Dordrecht;Boston: Kluwer Academic Publishers.  Godfrey, J. S., & Beljaars, A. C. M. (1991). On the turbulent fluxes of buoyancy, heat and moisture at the air‐sea interface at low wind speeds. Journal of Geophysical Research: Oceans, 96(C12), 22043-22048. doi:10.1029/91JC02015  Goldacre, R. J. (1949, May). Surface Films on Natural Bodies of Water. Journal of Animal Ecology, 18(1), 36-39.  Gottifredi, J. C., & Jameson, G. J. (1968). The suppression of wind-generated waves by a surface film. Journal of Fluid Mechanics, 32(3), 609-618. doi:10.1017/S0022112068000911 Grenfell, T. C., and Perovich, D. K. ( 2004), Seasonal and spatial evolution of albedo in a snow‐ice‐land‐ocean environment, Journal of Geophysical Research, 109, C01001. Gu, L., Massman, W. J., Leuning, R., Pallardy, S. G., Meyers, T., Hanson, P. J., ... & Yang, B. (2012). The fundamental equation of eddy covariance and its application in flux measurements. Agricultural and Forest Meteorology, 152, 135-148. Guo, S. J., Zhang, F. H., Song, X. G., & Wang, B. T. (2015). Deposited sediment settlement and consolidation mechanisms. Water Science and Engineering, 8(4), 335-344. 81  Haghighi, E., Shahraeeni, E., Lehmann, P., & Or, D. (2013). Evaporation rates across a convective air boundary layer are dominated by diffusion. Water Resources Research, 49(3), 1602-1610.  Handler, R. A., Smith, G. B., & Leighton, R. I. (2001). The thermal structure of an air-water interface at low wind speeds Taylor & Francis. doi:10.3402/tellusa.v53i2.12187 Henry, D. J., Dewan, V. I., Prime, E. L., Qiao, G. G., Solomon, D. H., & Yarovsky, I. (2010). Monolayer structure and evaporation resistance: A molecular dynamics study of octadecanol on water. The Journal of Physical Chemistry B, 114(11), 3869-3878. Högström, U. (1988). Non-dimensional wind and temperature profiles in the atmospheric surface layer: A re-evaluation. Boundary-Layer Meteorology, 42(1-2), 55-78. doi:10.1007/BF00119875 Hutter, K., Wang, Y., & Chubarenko, I. P. (2014). Dimensional Analysis, Similitude and Model Experiments. In Physics of Lakes (pp. 307-396). Springer, Cham.  Hurley, D. L. (2017). Wind Waves and Internal Waves in Base Mine Lake. University of British Columbia. doi:10.14288/1.0351993  Hurley, D., Lawrence, G., & Tedford, E. (2020). Effects of Hydrocarbons on Wind Waves in a Mine Pit Lake. Mine Water and the Environment, 1-9. Imberger, J. (1985). The diurnal mixed layer. Limnology and Oceanography, 30(4), 737-770. Imboden, D. M., & Wüest, A. (1995). Mixing mechanisms in lakes. In Physics and Chemistry of Lakes (pp. 83-138). Springer, Berlin, Heidelberg.  Jha, M. N., Levy, J., & Gao, Y. (2008). Advances in Remote Sensing for Oil Spill Disaster Management: State-of-the-Art Sensors Technology for Oil Spill Surveillance. Sensors (Basel, Switzerland), 8(1), 236–255. https://doi.org/10.3390/s8010236 Katsaros, K. (2001). Evaporation and Humidity. In Encyclopedia of Ocean Sciences (Vol. 2, pp. 870-877). Kawai, Y., & Wada, A. (2007). Diurnal sea surface temperature variation and its impact on the atmosphere and ocean: A review. Journal of Oceanography, 63(5), 721-744.  Koberg, G. E., Cruse, R. R., & Shrewsbury, C. L. (1963). Evaporation Control Research, 1959-1960. Washington, D.C.: US Department of the Interior. Kondo, J. (1975). Air-sea bulk transfer coefficients in diabatic conditions. Boundary-Layer Meteorology, 9(1), 91-112. doi:10.1007/BF00232256  Kumano, H., Asaoka, T., Saito, A., & Okawa, S. (2007). Study on latent heat of fusion of ice in aqueous solutions. International Journal Rof efrigeration, 30(2), 267-273. Lawrence, G. A., Tedford, E. W., & Pieters, R. (2016). Suspended solids in an end pit lake: potential mixing mechanisms. Canadian Journal of Civil Engineering, 43(3), 211-217. Langmuir, I., & Langmuir, D. B. (1927). The effect of monomolecular films on the evaporation of ether solutions. The Journal of Physical Chemistry, 31(11), 1719-1731. doi:10.1021/j150281a011  82  Langmuir, I., & Schaefer, V. J. (1943). Rates of evaporation of water through compressed monolayers on water. Journal of the Franklin Institute, 235(2), 119-162. doi:10.1016/S0016-0032(43)90904-4  Lee, R. W., & Rast, W. (1997). Light Attenuation in a Shallow, Turbid Reservoir, Lake Houston, Texas. TX: U.S. Dept. of the Interior, U.S. Geological Survey ; Branch of Information Services [distributor], 1997. Leppäranta, M. (2014). Freezing of lakes and the evolution of their ice cover. Springer Science & Business Media.  Liu, W. T., Katsaros, K. B., & Businger, J. A. (1979). Bulk parameterization of air-sea exchanges of heat and water vapor including the molecular constraints at the interface. Journal of the Atmospheric Sciences, 36(9), 1722-1735. doi:10.1175/1520-0469(1979)036<1722:BPOASE>2.0.CO;2 Luo, L., Hamilton, D., & Han, B. (2010). Estimation of total cloud cover from solar radiation observations at Lake Rotorua, New Zealand. Solar Energy, 84(3), 501-506. Manwell, J. F., McGowan, J. G., Rogers, A. L., & Wiley InterScience (Online service). (2009;2010;). Wind energy explained: Theory, design and application (2nd;2. Aufl.; ed.). Chichester, U.K: Wiley.  Massman, W. J., & Lee, X. (2002). Eddy covariance flux corrections and uncertainties in long-term studies of carbon and energy exchanges. Agricultural and Forest Meteorology, 113(1), 121-144. doi:10.1016/S0168-1923(02)00105-3  Mobasheri, M. R. (2006). Reformation time for the thermal skin layer of the ocean. International Journal of Remote Sensing, 27(23), 5285-5299. doi:10.1080/01431160600836026  Momii, K., & Ito, Y. (2008). Heat budget estimates for lake Ikeda, Japan. Journal of Hydrology, 361(3-4), 362-370.  Monin, A. S., & Obukhov, A. M. (1954). Basic laws of turbulent mixing in the surface layer of the atmosphere. Contrib. Geophys. Inst. Acad. Sci. USSR, 151(163), e187. Murray, M. J., Allen, M. R., Merchant, C. J., Harris, A. R., & Donlon, C. J. (2000). Direct observations of skin‐bulk SST variability. Geophysical Research Letters, 27(8), 1171-1174. Noh, Y., Lee, E., Kim, D., Hong, S., Kim, M., & Ou, M. (2011). Prediction of the diurnal warming of sea surface temperature using an atmosphere‐ocean mixed layer coupled model. Journal of Geophysical Research: Oceans, 116(C11), n/a. doi:10.1029/2011JC006970  Olsthoorn, J., Bluteau, C. E., & Lawrence, G. A. (2019;2020;). Under‐ice salinity transport in low‐salinity waterbodies. Limnology and Oceanography, 65(2), 247-259. doi:10.1002/lno.11295 Paaijmans, K. P., Takken, W., Githeko, A. K., & Jacobs, A. F. G. (2008). The effect of water turbidity on the near-surface water temperature of larval habitats of the malaria mosquito 83  anopheles gambiae. International Journal of Biometeorology, 52(8), 747-753. doi:10.1007/s00484-008-0167-2  Panofsky, H. A., & Dutton, J. A. (1984). Atmospheric turbulence: Models and methods for engineering applications. New York: Wiley.  Prakash, S., J. A. Vendenberg, and E. Buchak (2011), The oil sands pit lake model—Sediment diagenesis module, in 19th International Congress on Modelling and Simulation, Modelling and Simulation Society of Australia and New Zealand, Perth, Australia, edited by F. Chan, D. Marinova, and R. S. Anderssen, pp. 3761–3767, Modell. and Simul. Soc. of Aust. and N. Z., Canberra. Roberts, W. J. (1962). Reducing water vapor transport with monolayers. In Retardation of Evaporation by Monolayers (pp. 193-201). Academic Press. RBR. (2019, October). Small Temperature Recorder Datasheet. In RBRsolo³ T | Temperature Logger. Rideal, E. K. (1925). On the influence of thin surface films on the evaporation of water. The Journal of Physical Chemistry, 29(12), 1585-1588. doi:10.1021/j150258a011  Saunders, P. M. (1967). The temperature at the ocean-air interface. Journal of the Atmospheric Sciences, 24(3), 269-273. doi:10.1175/1520-0469(1967)024<0269:TTATOA>2.0.CO;2  Schwerdtfeger, P. (1963). Theoretical derivation of the thermal conductivity and diffusivity of snow. Int. Assoc. Sci. Hydrol. Publ, 61, 75-81. Smith, S. D. (1988), Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature, Journal of Geophysical Research, 93( C12), 15467– 15472, doi:10.1029/JC093iC12p15467.  Smith, S. D., & Jones, E. P. (1985). Evidence for wind‐pumping of air‐sea gas exchange based on direct measurements of CO2 fluxes. Journal of Geophysical Research: Oceans, 90(C1), 869-875. Soloviev, A. V., & Schlüssel, P. (1994). Parameterization of the cool skin of the ocean and of the air-ocean gas transfer on the basis of modeling surface renewal. Journal of Physical Oceanography, 24(6), 1339-1346. doi:10.1175/1520-0485(1994)024<1339:POTCSO>2.0.CO;2 Stull, R. B. (1988). An Introduction to Boundary Layer Meteorology. Dordrecht, The Netherlands: Kluwer Academic Publishers. Sturm, M., Holmgren, J., König, M., & Morris, K. (1997). The thermal conductivity of seasonal snow. Journal of Glaciology, 43(143), 26-41. doi:10.3189/S0022143000002781 Sturm, M., Perovich, D. K., & Holmgren, J., Thermal conductivity and heat transfer through the snow on the ice of the Beaufort Sea, Journal of Geophysical Research, 107( C10), 8043, doi:10.1029/2000JC000409, 2002. Sturm, M., Taras, B., Liston, G. E., Derksen, C., Jonas, T., & Lea, J. (2010). Estimating snow water equivalent using snow depth data and climate classes. Journal of Hydrometeorology, 11(6), 1380-1394. 84  Taylor, P. (1998). The thermodynamic effects of sublimating, blowing snow in the atmospheric boundary layer. Boundary-Layer Meteorology, 89(2), 251-283. Tedford, E. (Photographer). (2016). Base Mine Lake Hydrocarbons [digital image].  Tedford, E., Halferdahl, G., Pieters, R., & Lawrence, G. A. (2019). Temporal variations in turbidity in an oil sands pit lake. Environmental Fluid Mechanics, 19(2), 457-473. doi:http://dx.doi.org.ezproxy.library.ubc.ca/10.1007/s10652-018-9632-6 UNESCO, I. (1981). The practical salinity scale 1978 and the international equation of state of seawater 1980. Tenth Report of the Joint Panel on Oceanographic Tables and Standards (JPOTS), 25. Von Karman, T. (1931). Mechanical similitude and turbulence. National Advisory Committee for Aeronautics; Washington, DC, United States Webb, E. K., Pearman, G. I., & Leuning, R. (1980). Correction of flux measurements for density effects due to heat and water vapour transfer. Quarterly Journal of the Royal Meteorological Society, 106(447), 85-100. Wei, Z., Miyano, A., & Sugita, M. (2016). Drag and bulk transfer coefficients over water surfaces in light winds. Boundary-layer meteorology, 160(2), 319-346. Wells, A. J., Cenedese, C., Farrar, J. T., & Zappa, C. J. (2009). Variations in ocean surface temperature due to near-surface flow: Straining the cool skin layer. Journal of Physical Oceanography, 39(11), 2685-2710. doi:10.1175/2009JPO3980.1  Wheeler, J. R. (1975). Formation and collapse of surface films. Limnology and Oceanography, 20(3), 338-342. doi:10.4319/lo.1975.20.3.0338  Whipple, G. C. (1898). Classification of lakes according to temperature. The American Naturalist, 32(373), 25-33. Wick, G. A., Ohlmann, J. C., Fairall, C. W., & Jessup, A. T. (2005). Improved oceanic cool-skin corrections using a refined solar penetration model. Journal of Physical Oceanography, 35(11), 1986–1996. https://doi.org/10.1175/JPO2803.1 Wilson, K. B., Hanson, P. J., Mulholland, P. J., Baldocchi, D. D., & Wullschleger, S. D. (2001). A comparison of methods for determining forest evapotranspiration and its components: sap-flow, soil water budget, eddy covariance and catchment water balance. Agricultural and Forest Meteorology, 106(2), 153-168. Wilson, R. C., Hook, S. J., Schneider, P., & Schladow, S. G. (2013). Skin and bulk temperature difference at Lake Tahoe: A case study on lake skin effect. Journal of Geophysical Research Atmospheres, 118(18), 10332–10346. https://doi.org/10.1002/jgrd.50786  Wu, J. (1985). On the cool skin of the ocean. Boundary-layer meteorology, 31(2), 203-207. Wunderlich, W. O. (1972). Heat and mass transfer between a water surface and the atmosphere. Norris, Tenn: Tennessee Valley Authority, Office of Natural Resources and Economic Development, Division of Air and Water Resources, Water Systems Development Branch. 85  Xing, Z., Fong, D. A., Tan, K. M., Lo, E. Y. M., & Monismith, S. G. (2012). Water and heat budgets of a shallow tropical reservoir. Water Resources Research, 48(6). Zhi, H., & Harshvardhan. (1993). A hybrid technique for computing the monthly mean net longwave surface radiation over oceanic areas. Theoretical and Applied Climatology, 47(2), 65-79. doi:10.1007/BF00866182     86  Appendix A: Parameterization of the Air-Sea Toolbox The Air-Sea Toolbox is a commonly used set of MATLAB® functions used for estimating surface fluxes into the water surface, using readily measured inputs. The code was developed using the algorithms presented in Fairall et al. (1996b), also referred to herein as COARE (Coupled Ocean Atmosphere Response Experiment).  COARE presents a method for estimating surface fluxes using readily measured bulk variables. The general modelling approach taken by COARE is based on an approach outlined in Liu et al. (1979); however, the parameterization and inputs were updated in Fairall et al. (1996b) to account for processes that had not previously been considered. Modifications made by Fairall et al. (1996b) to the parameter estimation methods, as well as background information, are described below. The symbols used herein are presented in Table A.1. A.1 Theory and Previous Models Bulk Equations for Stress, Latent Heat, and Sensible Heat The bulk expressions summarized by Fairall et al. (1996b) for the momentum, sensible heat, and latent heat fluxes from the atmosphere into the water surface are presented in Equations A.1-A.3. The momentum flux due to wind shear at the water surface, also referred to as turbulent stress (𝜏), can be expressed as:  𝜏 =  𝜌𝑎𝐶𝑑𝑆2 (A.1) where 𝜌𝑎 is the density of air; 𝐶𝑑   is the momentum transfer coefficient; and 𝑆 is the average value of wind speed relative to the water surface at 𝑧𝑟. 𝑧𝑟 is a reference height at which parameters are measured or to which they are scaled, and typically has a value of 10 m. The equation for sensible heat flux (𝐻𝑠) has a similar form:  𝐻𝑆 =  𝜌𝑎𝑐𝑝𝑎𝐶ℎ𝑆(𝑇𝑠 − 𝜃) (A.2) where 𝑐𝑝𝑎 is the specific heat of air; 𝑇𝑆 is the surface interface temperature; 𝜃 is the potential temperature of the air at the air-water interface; and 𝐶ℎ is the heat transfer coefficient.   87  Table A.1 Primary variables and parameters used in COARE. Some minor terms omitted for brevity. Symbol Definition Units 𝐶𝑑 Transfer coefficient for stress Unitless 𝐶𝑒 Transfer coefficient for latent heat Unitless 𝐶ℎ Transfer coefficient for sensible heat Unitless 𝑐𝑑, 𝑐𝑑𝑛 Profile component for stress; 𝑐𝑑𝑛 applies in neutral conditions Unitless 𝑐𝑇, 𝑐𝑇𝑛 Profile component for temperature; 𝑐𝑇𝑛 applies in neutral conditions  Unitless 𝑐𝑞, 𝑐𝑞𝑛 Profile component for moisture; 𝑐𝑞𝑛 applies in neutral conditions Unitless 𝑐𝑝𝑎 Specific heat of air at constant pressure kJ∙kg-1∙K-1 𝑔 Standard gravity m∙s-2 𝐻𝐿 Latent heat flux W∙m-2 𝐻𝑆 Sensible heat flux W∙m-2 𝐿𝐸 Latent heat of vaporization J∙kg-1 𝑞 Specific humidity, or water vapor mixing ratio g∙kg-1 𝑞′ Turbulent fluctuation in water vapor mixing ratio kg water per kg moist air 𝑞∗ Monin-Obukhov similarity scaling parameter for humidity kg water per kg moist air 𝑞𝑠 Saturation specific humidity at 𝑇 = 𝑇𝑆 g∙kg-1 𝑅𝑖 Gradient Richardson number Unitless 𝑅𝑟 Roughness Reynolds number Unitless 𝑆 Average wind speed relative to surface at height 𝑧𝑟 m∙s-1 𝑇∗ Monin-Obukhov similarity scaling parameter for temperature K 88  Symbol Definition Units 𝑇′ Turbulent fluctuation in temperature T 𝑇𝑆 Sea surface temperature K 𝑢′ Turbulent fluctuation in streamwise component of horizontal wind m∙s-1 𝑢∗ Shear velocity m∙s-1 𝑢 Measured wind speed  m∙s-1 𝑤′,?̅? Turbulent fluctuations and mean vertical wind components m∙s-1 𝑊∗ Convective scaling velocity m∙s-1 𝑧𝑜 Velocity roughness length m 𝑧𝑜𝑇 Scalar roughness length for temperature transfer m 𝑧𝑜𝑞 Scalar roughness length for humidity transfer m 𝑧𝑟 Reference height m 𝛼 Charnock coefficient Unitless 𝛽 Empirical constant for convective scaling velocity Unitless 𝜅 Von Karmán constant Unitless 𝜃 Potential temperature of air brought adiabatically to a standard pressure K 𝜈 Kinematic viscosity of air m2∙s-1 𝜌𝑎 Density of air kg∙m-3 𝜓𝑐 Profile function near the free convective limit Unitless 𝜓𝑢 Profile function in stable conditions Unitless 𝜓𝑥𝐾 Profile function in near-neutral conditions Unitless 𝜏 Turbulent stress N∙m-2 89  Finally, latent heat flux (𝐻𝐿) from phase changes is calculated using:  𝐻𝐿 =  𝜌𝑎𝐿𝐸𝐶𝑒𝑆(𝑞𝑠 − 𝑞)  (A.3) where 𝐿𝐸 is the latent heat of vaporization; 𝑞𝑠 is the saturation specific humidity at the air-water interface; 𝑞 is the actual specific humidity; and 𝐶𝑒 is the moisture transfer coefficient. The terms in the bulk equations are readily determined, with the exception of 𝐶𝑑, 𝐶ℎ  and 𝐶𝑒.  Governing Equations for Momentum, Latent Heat, and Sensible Heat  The bulk expressions were developed based on Equations A.1-A.3, which define fluxes in terms of Reynolds averages. Formulating the equations in terms of turbulent fluctuations offers a good representation of the physical processes involved; however, the measurement of turbulent fluxes is not a straightforward process.   𝜏 =  𝜌𝑎𝑤′𝑢′̅̅ ̅̅ ̅̅ =  −𝜌𝑎𝑢∗2 (A.4)  𝐻𝑆 =  𝜌𝑎𝑐𝑝𝑎𝑤′𝑇′̅̅ ̅̅ ̅̅ =  −𝜌𝑎𝑐𝑝𝑎𝑢∗𝑇∗ (A.5)  𝐻𝐿 =  𝜌𝑎𝐿𝐸𝑤′𝑞′̅̅ ̅̅ ̅̅ =  −𝜌𝑎𝐿𝐸𝑢∗𝑞∗ (A.6) where 𝑤′, 𝑢′, 𝑇′, and 𝑞′ represent the turbulent fluctuations of vertical wind, streamwise component of horizontal wind, temperature, and water vapor mixing ratio, respectively, and 𝑢∗, 𝑇∗, and  𝑞∗ are the related Monin-Obukhov similarity scaling parameters (Monin & Obukhov, 1954), described in a following section. The overbar over the turbulent fluctuations indicates an ensemble average; in practice, a space or time average is used (Fairall et al., 1996b). A.2 Theory of Similitude The terms 𝑢∗, 𝑇∗ and 𝑞∗ are the Monin-Obukhov similarity scaling parameters, which describe the transfer of momentum, heat, and moisture in the atmosphere and were developed based on the assumption that these processes obey the theory of similitude. The theory of similitude states that there are universal relationships between variables, which can be established using dimensional analysis (Garratt, 1994; Stull, 1988). These universal relationships can be used to describe the behaviour of a system at different scales than the one being measured, or in different settings than the original experiment. Consequently, if flow at a given height can be measured and the transformation function to scale the atmosphere to different heights can be determined, the flow 90  throughout the region of interest can be predicted. This assumes a homogeneous, semi-infinite plane. Monin-Obukhov Similarity Scaling Parameters Applying the theory of similitude, we first consider u∗, the similarity scaling parameter for velocity (also known as friction velocity or shear velocity). This parameter is defined as:  𝑢∗ =  √𝜏𝜌𝑎  (A.7) where 𝑢∗ is defined at the interface between the atmosphere and the surface it is interacting with, in this case the water.  However, assuming a steady-state horizontally homogeneous boundary layer, within the surface layer (approximately the bottom 50 m of the atmosphere), 𝑢∗ and 𝜏 are independent of height z  (Monin & Obukhov, 1954; Garratt, 1994). Because 𝑢∗ is the parameter used to scale the flow, with the assumption that it is self-similar with height, a function can be defined relating the friction velocity to how the wind speed changes with height in neutral conditions:  ?̅?𝑧2−?̅?𝑧1𝑢∗= 𝑓 (𝑧2𝑧1) (A.8) where the only solution to this function is:  ?̅?𝑧2−?̅?𝑧1𝑢∗=1𝜅𝑙𝑛 (𝑧2𝑧1) (A.9) where κ is the Von Kárman constant, and ?̅?𝑧2and ?̅?𝑧1are the mean wind speeds recorded at heights 𝑧2 and 𝑧1, respectively. The wind speed at any height above the characteristic roughness scale of the surface in question can thereby be estimated, given a known wind speed and measurement height (Monin & Obukhov, 1954; Garratt, 1994). 𝑢∗ is thus the velocity scaling parameter. However, in non-neutral conditions, additional complications appear. Due to roughness elements near the surface and their resultant turbulence, as well as buoyancy effects that arise in an unstable atmosphere, the near-surface 𝑢∗ needs to be further defined (Monin & Obukhov, 1954).   91  Profile Components In COARE, the scaling parameters are defined using profile components:  𝑢∗ = −𝑐𝑑1/2𝑆 (A.10)  𝑇∗ = −𝑐𝑇12 (𝑇𝑠 − 𝜃)   (A.11)  𝑞∗ = −𝑐𝑞12(𝑞𝑠 − 𝑞)  (A.12) where 𝑐𝑑 is a profile component used to describe the effects of the boundary layer on the flow, and 𝑐𝑇 and 𝑐𝑞 are the temperature and moisture profile components . The velocity profile component is calculated using Equations A.13 and A.17, below; the temperature and moisture profile component calculations follow analogous forms.  Transfer Coefficients for Stress, Latent Heat, and Sensible Heat Transfer coefficients, which describe the extent to which stresses on the system are imparted into the water, have been the subject of considerable study (O’Connor, 1983; Panofsky & Dutton, 1984; Geernaert & Plant, 1990; Kondo, 1975; Francey & Garratt, 1977). Transfer coefficients can be defined in terms of the profile components: 𝐶𝑑 = 𝑐𝑑1/2𝑐𝑑1/2; 𝐶ℎ = 𝑐𝑇1/2𝑐𝑑1/2; and 𝐶𝑒 =𝑐𝑞1/2𝑐𝑑1/2 (note the use of capitalization to distinguish the different parameters).  Calculating the profile components requires determining the neutral profile components (𝑐𝑑1/2, 𝑐𝑇1/2 and 𝑐𝑞1/2) and the Monin-Obukhov profile function (𝜓) for the atmospheric conditions under evaluation. The former are functions of the roughness length scales, which are dependent on 𝑢∗ (Equations A.13-A.17); and the latter is a function of temperature, 𝑢∗, and 𝑞∗ (Equations A.19-A.22; Equation A.25). In short, the variables used to calculate heat fluxes are also a function of those fluxes, due to their effects on the boundary layer profile.    92  A.3 Parameterization Roughness-Stress Relationship Neutral Transfer Coefficient Calculations In neutrally stratified atmospheric conditions, meaning that the lapse rate is adiabatic and that there is no convection, the profile components are a function of the roughness lengths for each parameter:  𝑐𝑑𝑛1/2=𝜅𝑙𝑜𝑔(𝑧𝑟𝑧𝑜) (A.13) Where 𝑐𝑑𝑛 is the neutral profile component for momentum, 𝜅 is the von Kármán constant (value of 0.4, and 𝑧𝑜 is the roughness length for momentum. Velocity Roughness Length: Smooth Surface and Rough Flow Despite the definition given above, the velocity roughness length is generally related to the physical roughness of the surface or the flow regime due to the impracticality of attempting to determine the height at which wind speed becomes zero for every measurement period.   The roughness length is a function of the properties of the roughness elements on the surface, predominantly their height (Garratt, 1994).  The roughness length is not the same as the height of the roughness elements; however, as the prominence (in height, density, and roughness) of the drag elements increase, roughness length increases as well. The roughness length for a flow regime over a smooth surface is presented in Equation A.14:  𝑧𝑜 =𝑅𝑟𝜈𝑢∗ (A.14) where ν is the kinematic viscosity of air and 𝑅𝑟 is the roughness Reynolds number, which approaches a value of 0.11 at low wind speeds where the flow is assumed to be laminar.   93  Charnock (1955) offers an alternate formulation for the velocity roughness length that does not consider the flow regime directly, but takes into account the physical roughness of the water surface:  𝑧𝑜 =𝛼𝑢∗2𝑔 (A.15) where 𝑔 is acceleration due to gravity and α is the Charnock coefficient with values ranging from 0.011 to 0.035 (Garratt, 1994). The value for 𝛼 is determined in part by the sea state, which describes the physical roughness of the system due to surface gravity waves, which are the dominant source of roughness on the ocean surface. In the open ocean, 𝛼 is typically close to 0.011, while close to the coast, 𝛼 has a value closer to 0.018. In lakes, which have limited fetch and space for the generation of swell, the value of α is assumed to be 0.035 (Xing et al., 2012).  Ultimately, COARE takes the approach given by Smith (1988), which considers each of the calculations given above to be a component of the total roughness length:  𝑧𝑜 =0.11𝜈𝑢∗+𝛼𝑢∗2𝑔 (A.16) Despite assuming a value of 𝑅𝑟 of 0.11, this formulation does not assume smooth flow; rather, it assumes that the first component of the equation describes the baseline roughness under smooth conditions, and the component of the equation with the Charnock coefficient will describe the additional roughness associated with a rough surface and turbulent flow.  Scalar Roughness Lengths and Iteration There are three roughness lengths: 𝑧𝑜, 𝑧𝑜𝑇, and 𝑧𝑜𝑞. The first is the velocity roughness length, while 𝑧𝑜𝑇 and 𝑧𝑜𝑞 are scalar roughness lengths, used to calculate the neutral transfer coefficients for heat and moisture fluxes, respectively. 𝑧𝑜𝑇 and 𝑧𝑜𝑞 are calculated based on the estimate for 𝑧𝑜 (refer to Liu et al. (1979) for details). In the COARE algorithm, an initial estimate for the roughness length is given; this is used to calculate an initial estimate for 𝑢∗, which is then used to recalculate 𝑧𝑜 until the values converge to a final estimate.  94  Convective Limit Profile Components in Non-Neutral Conditions Having determined the neutral profile component, the actual profile component for the atmospheric conditions must be determined. This is a function of the neutral profile component and a parameter describing the atmospheric conditions:  𝒄𝒅𝟏/𝟐=𝒄𝒅𝒏𝟏/𝟐𝟏−𝑪𝒅𝒏𝟏/𝟐𝝍𝜿 (A.17) Where ψ is the Monin-Obukhov scaling profile function describing the atmospheric conditions. Monin-Obukhov Profile Functions The equation determining ψ varies depending on the stability of the atmosphere. It is a parameter that arises due to the fact that in non-neutral conditions, the wind speed does not follow the log-linear relationship shown in Figure A.1, but rather varies due to buoyancy effects. In the constant-flux layer, the wind is represented by Equation A.9, also written as:   𝒅𝒖𝒅𝒛=𝒖∗𝒌𝒛 (A.18) However, in non-neutral conditions, this equation is modified to:  𝑑𝑢𝑑𝑧=𝑢∗𝑘𝑧𝜓 (𝑧𝐿) (A.19) And corresponding equations exist for temperature and moisture fluxes:  𝒅𝑻𝒅𝒛=𝑻∗𝒛𝝍 (𝒛𝑳) (A.20)   𝑑𝑞𝑑𝑧=𝑞∗𝑧𝜓 (𝑧𝐿) (A.21) Where 𝐿 is the Monin-Obukhov length scale, which relates the characteristic parameters of the atmosphere. Refer to Monin and Obukhov (1954) for the full derivation of the above equations. Monin-Obukhov Length Scale L is a parameter that arises in Monin-Obukhov similarity scaling theory and links the 3 main variable groups that define the surface-layer turbulence: 𝑔/𝑇𝑎, 𝑢∗, and 𝐻𝑠/(𝑐𝑝𝜌), where 𝑇𝑎 is the 95  air temperature. The term 𝑔/𝑇𝑎 is related to convection and buoyancy; the ratio 𝐻𝑠/(𝑐𝑝𝜌) is constant under similarity theory, indicating that the vertical fluctuations in temperature and wind maintain the same ratio with height in the boundary layer; and 𝑢∗ describes the same for shear stress.  L is therefore a scale relating these characteristic parameters to each other:  𝐿 =  −𝑢∗3𝜅𝑔𝑇𝑎𝐻𝑠𝑐𝑝𝜌  = (𝜅𝑔𝑇𝑎(𝑇∗+0.61𝑇𝑎𝑞∗)𝑢∗2 )−1 (A.22) It is related to 𝑧 to create a non-dimensionalized parameter:  𝜉 = 𝑧/𝐿  (A.23) 𝜉 is used for the calculation of 𝜓 in various atmospheric conditions (Refer to Table A.2). Table A.2 Formulae for calculating 𝝍 in different atmospheric conditions. 𝜸 is an empirical constant with a value of 12.87 (Fairall et al., 1996b). Atmospheric Conditions Corresponding equation for atmospheric profile functions Stable 𝜓𝑢 =  −4.7 𝜉  Near-neutral 𝜓𝑥𝐾 = 2 𝑙𝑛 (1 + √1 − 16𝜉2)  Unstable, near the free convective limit 𝜓𝐶 = 1.5 𝑙𝑛 (𝑦2 + 𝑦 + 13) − √3 𝑎𝑡𝑎𝑛 (2𝑦 + 1√3) +𝜋√3 𝑦 = √1 − 𝛾𝜉3  The profile functions given in Table A.2 are modifications to the general form of Businger-Dyer profiles (Businger, 1988). Although there are a variety of equations available in the literature, Högström (1988) found that much of the difference appeared to be due to inadequacy in the instrumentation.  96  Monin-Obukhov Length Scale: Interpretation and Implications In practical terms, 𝐿 describes the balance between shear and buoyancy fluxes. 𝐿 is related to the Richardson number, which describes the temperature and wind gradient of a fluid. The traditional formulation for the gradient Richardson number is:  𝑹𝒊 =  −𝒈 (𝒅𝑻𝒅𝒛)𝑻𝒂(𝒅𝒖𝒅𝒛)𝟐 (A.24) Where the numerator describes the production/consumption of turbulence due to buoyancy, and the denominator describes the production of turbulence due to shear. Neutral atmospheric conditions therefore have 𝑅𝑖 values of 0.  Substituting Equations A.19, A.20, and A.22, and rearranging Equation A.5 for 𝑇∗, the Richardson number can be simplified to:  𝑅𝑖 =𝑧𝐿1𝜓  (A.25) The Richardson number is negative in statically unstable (convective) atmospheres where the surface temperature is higher than air temperatures, and positive in stratified atmospheres. In these conditions, 𝐿 describes the height at which buoyancy and shear fluxes balance in the atmosphere; for a given convective flux, the lower the value of 𝐿, the more likely it is that the flow will be dominated by convection.  In the stable atmosphere, onset of turbulence occurs when 𝑅𝑖 drops below the critical Richardson number, with a value of approximately 0.25, and turbulent flow becomes laminar when 𝑅𝑖 rises above a value of 1. Convective Limit When wind speed drops to very low values, causing 𝑢∗ to approach zero, transfer coefficient calculations become asymptotic: 𝑐𝑑1/2 ∝ 1/𝑢∗   (Equations A.17, A.22, and A.23; Table A.1). Thus, although the heat and moisture fluxes appear at first glance to be directly proportional to 𝑢∗ in Equations A.4 and A.5, they instead approach infinity as 𝑢∗ approaches zero. The formulae for the transfer coefficients in an unstable atmosphere are therefore invalid when directly applied in very low winds. 97  In stratified conditions, this does not pose a problem, as the stability of the atmosphere prevents significant transfer in the absence of wind. However, in unstable or near-stable conditions, buoyancy fluxes become the dominant forces inducing mixing.  A weighting system to blend different calculation methods is therefore implemented:  𝝍 =𝟏𝟏+𝝃𝟐𝝍𝒙𝑲 +𝝃𝟐𝟏+𝝃𝟐𝝍𝒄 (A.26) Under this formula, where 𝜉2 <<1 (and so by definition, 𝐿<<𝑧), the profile function is effectively that at a marginally stable atmospheric condition. As 𝐿 approaches zero, the value of the profile function becomes similar to that at free convection.  Profile functions are assumed to be the same for temperature and humidity. Gustiness In addition to the problem of determining an appropriate value for the profile function in very low wind speeds, additional evaporation at low wind speeds is observed due to dry convective motions within the atmospheric boundary layer. These corrections for “gustiness” can lead to latent heat fluxes of 25 W∙m-2 in the tropics (Godfrey & Beljaars, 1991).  This phenomenon is analogous to penetrative convection in lakes, where cooling packets of water drop due to their increased density compared with the surrounding water. The velocity of these packets influences how much mixing they induce (Fischer, 1980). The mean wind speed used for calculating turbulent fluxes is therefore determined by augmenting the measured wind speed with a “gustiness” estimate:  𝑆2 = 𝑢2 + (𝛽𝑊∗)2 (A.27) where 𝛽 is an empirical constant with a value close to 1, and 𝑊∗ is the convection velocity, also known as the convective scaling velocity (Godfrey & Beljaars, 1991). In COARE, 𝛽 is assigned a value of 1.25.   98  The parameter 𝑊∗ describes the typical speed of the horizontal velocity fluctuations in the surface layer due to buoyant convection; this assumes dry convection, with neither cloud formation nor rainfall considered. 𝑊∗ is calculated as:  𝑊∗3 = −𝐹𝐵𝑧𝑖 =  𝑔𝑧𝑖𝑇𝑎[𝐻𝑆𝜌𝑎𝑐𝑝𝑎+0.61𝑇𝑎𝐻𝐿𝜌𝑎𝐿𝑒] (A.28) where 𝐹𝐵 is the buoyancy flux and 𝑧𝑖 is the depth of the convective boundary layer. 𝑧𝑖 generally ranges between 100 and 1000 m. The COARE algorithm assigns 𝑧𝑖 a value of 600 m based on a typical convective boundary layer height in the atmosphere. Refer to Stull (1988) for further details on the convective scaling velocity. Because 𝑊∗ is a function of 𝐻𝑆 and 𝐻𝐿, iteration to determine its value is necessary in unstable atmospheric conditions. When the atmosphere is stable, convection is limited and 𝑊∗ is assigned a value of zero. In COARE, this is increased to +0.5 in order to prevent numerical errors. Webb Correction When calculating latent heat flux, in addition to the corrections to the profile function at low wind speeds, and the increase in the effective wind speed due to convective motion, consideration must be made for the heat carried by the air from vertical air movements due to thermal instabilities. When the lake is losing heat, air particles close to the water are warmer, and therefore less dense, than the particles above. Convective exchange of these air masses must account for this density difference in order to meet the assumption that there is zero total vertical mass flux, with the lighter rising air particles having higher velocities than the falling cooler, denser air.  This additional velocity component, described in Webb et al. (1980), carries water vapour and heat. The new term for total latent heat flux then becomes:  𝐻𝐿 = 𝜌𝑎𝐿𝑒𝐶𝑒𝑢(𝑞𝑠 − 𝑞) + 𝜌𝑎𝐿𝑒?̅?𝑞 (A.29) That is to say, the original latent heat estimate is augmented with the flux due to the mean velocity, ?̅?.  Note that the Webb correction to the latent heat flux is evaluated separately from the latent heat in the Air-Sea toolbox; both terms must be considered when creating a heat budget. 99  A.4 Surface Temperature Estimation The input water temperature for the COARE model is the interfacial (surface) temperature of the water that is in direct contact with the atmosphere. However, as these measurements are difficult to obtain in practice, bulk sensors in the water column are typically used. Unlike radiative sensors that take measurements at the water surface, bulk sensors must be adjusted for diurnal and cool-skin effects. Fairall et al. (1996a) emphasizes the need to incorporate both cool-skin effects and warm-layer effects; a method for making these calculations is provided in Fairall et al. (1996b). Skin effects are included as an option in the Air-Sea toolbox; however, additional inputs (downwelling longwave radiation, net shortwave radiation, and measured insolation) are required for these calculations. Warm-layer effects are not included in the Air-Sea toolbox, and must be implemented separately. A.5 List of Assumptions The Monin-Obukhov scaling profile function for temperature and humidity is assumed to be the same as for wind. If calculating for precipitation effects, assume that rain is at the wet-bulb temperature. A non-comprehensive list of the values of parameters used by COARE are provided in Table A.3, below. Refer to Liu et al. (1979) and Fairall et al. (1996b) for a complete list. The values in the table are standard to COARE; however, in certain settings other values may be needed.    100  Table A.3 Parameter values as assigned in the Air-Sea toolbox. Modifications to values may be required in non-marine settings. Parameter  Component Value Units 𝑔 Gravity 9.81 m∙s-2 𝑧𝑖 Inversion height 600 m 𝛼 Charnock coefficient 0.011 [-] 𝛼 Mean albedo 0.055 [-] 𝛽 Convection scaling coefficient 1.25 [-] 𝛽𝑊∗ Minimum gustiness 0.5 m∙s-1 𝛾 Convective profile coefficient 12.87 [-] 𝜀 Mean emissivity 0.97 [-] 𝜅 von Karmán coefficient 0.4 [-]     101  Appendix B: Determining the Near-Surface Temperature The surface temperature of the water is a key input to the Air-Sea Toolbox, as the heat and moisture available for transfer to or from the atmosphere is a function of this variable (refer to Appendix A). This temperature can be determined directly, using infrared sensors (these would also account for the cool-skin effect on the water); however, in the absence of such data, estimation of the surface temperature must be done using models or empirical observations. Depending on the data available and the research goals, it may not be necessary or appropriate to estimate the near-surface temperature, and one can simply use the observations from the shallowest available thermistor. That is, in lakes with elevated mixing or incomplete data sets, or where precise heat flux estimates are not needed, attempting to estimate near-surface temperature may not yield more accurate results than what would be achieved otherwise. However, at Base Mine Lake (BML), the availability of high-quality data and need for accurate estimation of each heat flux demands that the warm layer of the surface be quantified for input into the model. An overview of the need and approach for estimating this parameter is presented below. During autumn overturn, when water temperatures are typically uniform in most lakes, near-surface stratification at BML can be observed during calm days (wind <5 m∙s-1) with net radiation above approximately 200 W∙m-2 (Figure B.1). During these periods, surface warming due to continued elevated solar heating and low mixing led to water temperatures at 0.15 m depth rising over 3 °C above the temperature at 0.5 m (𝑇0.5) well into fall turnover on 1 Oct 2015 (DOY 275).  For most of the study period, the near-surface sensors at P1, measuring temperatures at 0.05 m, 0.15 m and 0.3 m, are unavailable. The meteorological conditions at the lake and the temperature difference between the 0.5 m and 1 m sensors, 𝑇0.5 and 𝑇1, were therefore used to estimate the surface temperature immediately below the skin layer:  𝑇0𝑒𝑠𝑡 = [𝑇0.5 + 𝛼𝑡 (𝑇0.5 − 𝑇1)] + 𝛽  (B.1) where 𝛼𝑡 accounts for the overall level of stratification of the water column, and 𝛽 accounts for short-term effects of wind and solar radiation on the near-surface water. Values of 𝛼𝑡 were determined by visual examination of the 0.15 m sensor at P1, and the 0.5 m and 1 m sensors at  102   Figure B.1 Influence of wind and net radiation on stratification during fall overturn. From DOY 271.5 and DOY 273.4, water temperatures are from the P1 sensors, as the P2 and P3 sensors were undergoing servicing. The 0.15 m sensor is also located at P1. Grey bars in panels (a) and (c) correspond with low wind (<5 m∙s-1) and high net radiation (>200 W∙m-2), respectively. The thermistor record (panel b) is overlain with light grey bars present where only one of the two stratifying conditions occurs, and dark grey bars indicating periods with both high winds and low net radiation. The dashed line in panel (a) highlights where the wind speeds drop below 5 m∙s-1.103  P2 and P3, while 𝛽 was established by examining these data series as well as the meteorological conditions at the lake. The resultant values, a function of wind speed and net radiation, are described below and outlined in Table B.1.  The first parameter, 𝛼𝑡, functionally assumes that a similar level of stratification seen between the 0.5 m and 1 m sensors will also be present between 0 m and 0.5 m. Because the warm layer at the surface may persist after radiative heating ceases until strong winds or surface cooling induce mixing of the layer, the parameter 𝛼 was set at 1.1 except in windy conditions where net radiation was negative. In these conditions, 𝛼𝑡 = 0 (i.e. it was assumed that the surface layer is completely mixed down to a depth of 0.5 m).  The second parameter, 𝛽, addresses short-term stratification and mixing events. In windy conditions, 𝛽 was zero as it was assumed that formation of surface waves would reduce stratification in the top half-metre of the water column. In calm conditions, 𝛽 varied between 0 °C and 1.5 °C, with higher values set in periods with higher net radiation based on the assumption that higher solar radiation induced greater warming of the upper water column. Both parameters are needed, as 𝛽 alone would not address the persistence of a warm layer at night, while 𝛼 alone cannot take into consideration short-term wind and net radiation effects, particularly when the 0.5 m and 1 m sensors may be similar in temperature but strong radiation from above is stratifying the upper portion of the water column. Table B.1 Variation in 𝜶𝒕 and 𝜷 in Equation 3.5 based on wind and net radiation conditions.  Net Radiation < 0 W∙m-2 Net Radiation > 250 W∙m-2 Net Radiation > 450 W∙m-2 Wind  <5 m∙s-1 𝛼𝑡 = 1.1 𝛽 = 0 𝛼𝑡 = 1.1 𝛽 = 1 𝛼𝑡 = 1.1 𝛽 = 1.5 Wind  > 5 m∙s-1 𝛼𝑡 = 0 𝛽 = 0 𝛼𝑡 = 1.1 𝛽 = 0 𝛼𝑡 = 1.1 𝛽 = 0  104  Figure B.2 presents the effects of the above corrections. The temperatures recorded by the 0.15 m sensor at P1 are presented against 𝑇0.5 and 𝑇0𝑒𝑠𝑡, and show that 𝑇0.5 underestimates the near-surface temperatures, particularly at higher temperatures.   Figure B.2 Scatter plot of corrected and uncorrected surface temperature estimates compared against the 0.15 m sensor at P1. The black line is the 1:1 line. The temperatures recorded by the 0.15 m sensor at P1 are presented on the X-axis, while 𝑻𝟎.𝟓, the nearest-surface sensors at P2 and P3, and 𝑻𝟎𝒆𝒔𝒕, the estimated surface temperature after correction for warm-layer effects, are plotted on the Y-axis.       


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