EFFECT OF SPECTRAL AND TEMPORAL VARIATION OF SUBSURFACE IRRADIANCE ON THE HEATING OF LAKES by Yasmin Nassar B.Sc., Cairo University, 2003 M.Sc., Cairo University, 2005 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2013 © Yasmin Nassar, 2013 ii ABSTRACT This dissertation illustrates how the vertical distribution of temperature in lakes can be affected by the fact that the attenuation coefficient of light is often strongly dependent on wavelength. The potential importance of this spectral effect is first examined by considering the solar radiation in isolation, and then by including all non– penetrative heat fluxes using a modified version of the numerical model, DYRESM. Comparing the subsurface spectral irradiance of different lakes reveals that the spectral variability of the attenuation coefficient is more significant when calculating the light intensity in relatively clear lakes than in turbid lakes. Comparisons made between field measurements and theoretical predictions of hypolimnetic heating show the importance of accounting for the spectral irradiance for two relatively clear lakes: Pavilion Lake and Crater Lake. A new parameterization that better describes the spectral attenuation coefficient and the distribution of subsurface irradiance is added to DYRESM. The results obtained, when running the original and the modified DYRESM on Pavilion Lake, show a significant improvement in predicting the thermal structure of the lake with the modified version. The effects of the variation in solar angle and the seasonal variation in water quality on the attenuation coefficient are also examined for Pavilion Lake using DYRESM modified to accept a time varying attenuation coefficient. Simulations were performed for Pavilion Lake using the original and modified versions of DYRESM on diurnal and seasonal scales. Results show no significant improvement in the thermal evolution of the lake when considering the diurnal variations, while slight improvement was shown on a seasonal scale. iii TABLE OF CONTENTS ABSTRACT ....................................................................................................................... ii TABLE OF CONTENTS ................................................................................................ iii LIST OF TABLES ............................................................................................................ v LIST OF FIGURES ......................................................................................................... vi LIST OF SYMBOLS ....................................................................................................... xi ACKNOWLEDGMENTS .............................................................................................. xv DEDICATION................................................................................................................ xvi 1 INTRODUCTION ...................................................................................................... 1 1.1 1.2 2 Background ......................................................................................................... 1 Research objectives ............................................................................................ 2 RADIATION TRANSFER FROM SUN TO WATER BODIES AND ITS EFFECT ON WATER TEMPERATURE: A LITERATURE REVIEW ................... 4 2.1 Introduction ........................................................................................................ 4 2.2 Attenuation coefficient in physical limnology.................................................. 4 2.3 Factors affecting the subsurface radiative transfer ........................................ 9 2.4 Determination of attenuation coefficients ...................................................... 12 2.4.1 Non-spectral attenuation coefficient ........................................................... 12 2.4.2 Spectrally variable Attenuation coefficient ................................................ 14 2.4.3 Attenuation coefficient as a function of wavelength and solar angle ......... 16 2.5 Parameterization of the attenuation coefficient in physical limnology ....... 19 2.5.1 Vertical ........................................................................................................ 19 2.5.2 Temporal ..................................................................................................... 19 2.5.3 Physical Modeling ...................................................................................... 20 3 SITES AND METHODS ......................................................................................... 23 3.1 Introduction ...................................................................................................... 23 3.2 Study sites.......................................................................................................... 23 3.2.1 Pure Water .................................................................................................. 24 3.2.2 Crater Lake.................................................................................................. 24 3.2.3 Lake Superior and San Vicente Reservoir .................................................. 25 3.2.4 Tailings Lake .............................................................................................. 25 iv 3.2.5 Pavilion Lake .............................................................................................. 26 3.3 Methodology ..................................................................................................... 32 3.3.1 DYRESM .................................................................................................... 32 3.3.2 Calculation of attenuation coefficients ....................................................... 34 4 PARAMETERIZATION OF SHORTWAVE IRRADIANCE USING SPECTRAL ATTENUATION ...................................................................................... 39 4.1 Introduction ...................................................................................................... 39 4.2 Subsurface Spectral Irradiance ...................................................................... 41 4.3 Subsurface irradiance intensity ...................................................................... 46 4.4 Heating rates due to spectral and non-spectral irradiance .......................... 48 4.5 Comparison of predicted heating with field data .......................................... 50 4.5.1 Pavilion Lake .............................................................................................. 50 4.5.2 Crater Lake.................................................................................................. 56 5 THE EFFECT OF SPECTRAL ATTENUATION ON TEMPERATURE EVOLUTION FOR DIFFERENT WATER CLARITIES ......................................... 59 5.1 5.2 5.3 5.4 5.4 6 Introduction ...................................................................................................... 59 DYRESM........................................................................................................... 60 ����𝐧 .................................................................................................. 61 Calculating 𝐊 Model results ..................................................................................................... 66 Conclusions ....................................................................................................... 75 TEMPORAL VARIATION OF ATTENUATION ............................................... 77 6.1 Introduction ...................................................................................................... 77 6.2 Dependence of the attenuation coefficient on the solar zenith angle ........... 78 6.3 Field observations ............................................................................................. 81 6.3.1 Diurnal variation in the attenuation coefficient .......................................... 81 6.3.2 Seasonal variation in the attenuation coefficient ........................................ 88 6.4 Model results ..................................................................................................... 91 6.4.1 Diurnal Variation ........................................................................................ 91 6.4.2 Seasonal Variation ...................................................................................... 94 6.5 Conclusions ....................................................................................................... 96 7 CONCLUSIONS ...................................................................................................... 98 7.1 7.2 7.3 Thesis summary ................................................................................................ 98 Main findings .................................................................................................... 99 Future work .................................................................................................... 100 BIBLIOGRAPHY ......................................................................................................... 101 APPENDICES ............................................................................................................... 112 v LIST OF TABLES Table 2.1: Trophic class according to Secchi depth adapted from Preisendorfer (1986) . 13 Table 3.1: DYRESM default parameters. ......................................................................... 34 Table 4.1: Optical properties of the four water types studied. .......................................... 43 Table 4.2: Observed dT dt , calculated dT dt using spectral K(λ)|Pav; and calculated dT dt using � ls |Pav .................................................................................................... 54 the non-spectral K ����n values estimated using Prony’s method ...................................................... 63 Table 5.1: K Table A.1: χ(λ), e(λ), and Kw(λ) adapted from Morel (1988) ......................................... 112 Table B.2: Absorption (a) and Scattering (b) coefficients of pure water adapted from Kirk (1984) .............................................................................................................................. 113 vi LIST OF FIGURES Figure 2.1: Solar radiation as function of wavelength at the top of the atmosphere and at sea level, adapted from Mecherikunnel and Duncan (1982). ............................................. 6 Figure 2.2: Spectral distribution of a) Ediffuse, b) Edirect, and c) ET for different sun angles at sea level. Values plotted are relative to the maximum light intensity, I(λ = 0.48 μm) = 1566 W m-2 μm-1, for vertical sun (θ = 0), based on equations (2-4) and (2-5).................. 8 Figure 2.3: Average attenuation coefficient, K(λ), under clear sky and no wind, for pure water (black), Crater Lake (blue), Lake Superior (green), and San Vicente Reservoir (red). Data from Bukata (1995) and Tyler and Smith (1970). .......................................... 10 Figure 2.4: Spectral distribution of a) Kdiffuse(λ) and b) Kdirect(λ, θ) for pure water, under different sun angles, based on equations (2-9) and (2-10). ............................................... 18 Figure 3.1: Bathymetry of Pavilion Lake showing north, central, and south basins. (+)TP marks the location of the thermistor chain and (x) marks the location of the meteorological station. ...................................................................................................... 27 Figure 3.2: Line plot of water temperatures measured using temperature recorders located in the central basin of Pavilion Lake, 2006. Temperatures are measured at z = 9.5, 15.5, and 26.5 m for the period from June 1st 2006 to August 1st 2006. The thermistors were recovered to retrieve the data, and then redeployed on August 2nd 2006 at z = 9.5, 15.1, and 24.1 m, recording temperatures until November 20th 2006. ...................................... 28 Figure 3.3: Meteorological data measured at the meteorological-station located near the north sill in Pavilion Lake for the period June 1st 2006 to November 20th 2006: a) Wind Speed (m s-1), b) Relative Humidity (%), c) Air Temperature (oC), d) PAR (W m-2). ..... 29 Figure 3.4: Line plot for subsurface PAR irradiance measured in the central basin of Pavilion Lake at z = 5.8 m (black) and 17 m (grey) for the period from June 1st 2006 to November 20th 2006. The empty gap shows the period, August 1st 2006 to August 8th vii 2006, during which the Licor sensors were pulled out to retrieve the data and redeployed at the same depths. ............................................................................................................ 30 Figure 3.5: Typical spectral response of LI-COR quantum sensor versus wavelength, adapted from (Licor 2006). ............................................................................................... 31 Figure 3.6: Flow chart showing the major steps of the heater module in DYRESM, and the new module parameterizing the spectral variability in the attenuation coefficient of PAR, presented in the dashed boxes. ................................................................................ 35 Figure 4.1: a) Spectral distribution of shortwave radiation at the surface, I(0, λ), based on an incoming irradiance intensity of 1000 W m-2. b) Attenuation coefficient, K(λ), for pure water (black) from Bukata (1995); Tyler and Smith (1970), Crater Lake (blue) from Bukata (1995), Lake Superior (green) from Tyler and Smith (1970), and San Vicente Reservoir (red) from Tyler and Smith (1970). For Crater Lake, K(λ) is estimated in the UV domain (blue dashed line) using Equation (4-5) suggested by Morel and Antoine (1994). ............................................................................................................................... 42 Figure 4.2: a) Spectral distribution of irradiance at the surface I(0, λ) (solid line, left scale) and K(λ) for Crater Lake (black dashed line, right scale). Also shown is the best �𝑙𝑠 (dash-dotted line, right scale). I(0, λ) was calculated using the typical maximum fit, 𝐾 solar radiation on a sunny day, ET(0) = 1000 W m-2. b) I(z, λ) calculated using K(λ), for Crater Lake (W m-2 µm-1). Line contours show the 200, 400, 600, 800, 1000, and 1200 �𝑙𝑠 . ...... 44 Wm-2 µm-1 levels. The thick solid line indicates λmax. c) I(z, λ) estimated using 𝐾 Figure 4.3: a) Spectral distribution of irradiance at the surface I(0, λ) (solid line, left scale) and K(λ) for Lake Superior (black dashed line, right scale). Also shown is the best �𝑙𝑠 (dash-dotted line, right scale). I(0, λ) was calculated using the typical maximum fit, 𝐾 solar radiation on a sunny day, ET(0) = 1000 W m-2. b) I(z, λ) calculated using K(λ), for Lake Superior (W m-2 µm-1). Line contours show the 200, 400, 600, 800, and 1000 W m-2 �𝑙𝑠 . ............... 45 µm-1 levels. The thick solid line indicates λmax. c) I(z, λ) estimated using 𝐾 viii Figure 4.4: Irradiance intensity profiles based on the PAR irradiance, assuming a surface �𝑙𝑠 (dashed line), and 𝐾 �1% (dotted intensity of 1000 W m-2, calculated using K(λ) (black), 𝐾 line) for a) Pure Water, b) Crater Lake, c) Lake Superior, and d) San-Vicente Reservoir. The profiles are plotted down to the photic depth per each water body as indicated by the right scale. ......................................................................................................................... 47 Figure 4.5: Rate of temperature increase based on the PAR irradiance, assuming a surface �𝑙𝑠 (dashed line) for a) pure intensity of 1000 W m-2, calculated using K(λ) (black), and 𝐾 water, b) Crater Lake, c) Lake Superior, and d) San-Vicente Reservoir. The profiles are plotted down to the photic depth for each water type as indicated by the right scale. ..... 49 Figure 4.6: Measured a) incoming PAR (W m-2), and b) Wind Speed (m s-1). Running daily averages of the temperature at c) 15.5 m, and d) 26.5 m; observed (black), estimated � ls |Pav (dashed magenta), predicted using K(λ)|Pav (solid magenta). .................... 52 using K Figure 4.7: Calculated 𝐾(𝜆)|𝑃𝑎𝑣 (magenta), using Equation (4-5) by Morel and Antoine � ls |Pav (dashed magenta) for PAR domain, K � ls |Crat (dashed blue) for UV and (1994), K PAR domains, K(λ)|Crat (blue), and 𝐾𝑤 (𝜆) (black) adapted from Bukata (1995). Grey area indicates the upper and lower bounds for K(λ)|Pav for C between 2 and 0.5 mg m-3 53 Figure 4.8: a) The average rate of increase in temperature with depth for Pavilion Lake. Observed (black triangle), using K(λ)|Pav �𝑙𝑠 (dashed (solid magenta), and using 𝐾 magenta). b) Temperature profile measured in June in Pavilion Lake showing the depth of the thermocline. ............................................................................................................ 55 Figure 4.9: a) Averaged temperature profiles for July (triangle) and August (circle) between 1988 and 2000 for Crater Lake. The points represent the mean with the standard deviation for each depth (adapted from Larson et al. (2007)), b) 𝑑𝑇 𝑑𝑡 observed (circle), � ls |Crat (dashed black). ... 57 calculated using K(λ)|Crat (solid black), and calculated using K Figure 5.1: Observed (11th August 2006) and estimated a) Irradiance profiles for the PAR domain, and b) rate of temperature gain based on the incoming PAR irradiance only in Pavilion Lake. Observed distribution (black), and calculated using a one (red), two ix (green), three (cyan), and four (blue) bands. Note that there is no significant difference between the three and four division PAR estimates, and that the blue line lies on top of the cyan line. ..................................................................................................................... 63 Figure 5.2: Observed and estimated a) Irradiance profiles for the PAR domain, and b) rate of temperature gain based on the incoming PAR irradiance only in Tailings Lake. Observed distribution (black), and calculated using a one (red), two (green), three (cyan), and four (blue) bands. Note that there is no significant difference between the two, three, and four division ET or 𝑑𝑇 𝑑𝑡 estimates, and that the blue line lies on top of the cyan and green lines. ........................................................................................................................ 64 Figure 5.3: (a) PAR, (b) wind speed, and (c) air temperature observed at Pavilion Lake for the studied period in 2006. .......................................................................................... 67 Figure 5.4: The evolution of the thermal stratification for Pavilion Lake simulated using DYRESM with a) one, b) two, c) three, and d) four band exponential forms, from the top down, representing the subsurface PAR distributions. Contours plotted are for 6, 8, 10, 12, 14, 16, 18, 20, 22, and 24 oC.. ..................................................................................... 68 Figure 5.5: Observed (black) and simulated temperatures in Pavilion Lake at 9.5 m, 15.5 m, and 26.5 m depths, one (red), two (green), three (cyan), and four (blue) band exponential forms.............................................................................................................. 70 Figure 5.6: a) Modeled surface temperatures using DYRESM with, one (red), two (cyan), three (green), and four band (blue) exponential forms. b) Differences between Tw modeled using four band exponential form and one band exponential form. .................. 71 Figure 5.7: a) latent, b) sensible, c) longwave, and d) total non-penetrative radiation (Wm-2), using one (red), two (blue), three (green), and four band (black) exponential forms. ................................................................................................................................ 74 Figure 5.8: Differences between a) latent, b) sensible, c) longwave, and d) total heat fluxes (W m-2), calculated using one and four band exponential forms. .......................... 75 x Figure 6.1:a) Fraction of diffuse and direct radiation in the total incoming radiation, at the water surface, Ediffuse/ET and Edirect/ET ; b) Kdirect(λ, θ), Kdiffuse(λ), and K(λ, θ) for pure water under different sun angles for λ = 0.48 μm. ...................................................................... 79 Figure 6.2: Attenuation coefficients for pure water at noon and at zenith angle = 77o; Kdirect(λ, 77o) (red), Kdirect(λ, 0o) (black), and Kdiffuse(λ) (blue). ......................................... 80 Figure 6.3: a) Irradiance intensity at the surface, ET(0) (black), at 5.8 m depth ET(5.8) (blue), and at 17 m depth ET(17) (red) b) Calculated 𝐾(PAR, θ) at z = 5.8 m (blue), z = 17 m (red) and a depth averaged 𝐾(PAR, θ) (dashed black). Black arrows point to the diurnal spike and dip on June 28th 2006. ....................................................................................... 83 Figure 6.4: a) Kdirect(θ), and Kdiffuse for Pavilion Lake under different sun angles b) Depth � (PAR, θ) (dashed black) for sunny days in June. K(λ, θ) for Pavilion Lake averaged K under different sun angles calculated using empirical equations. Grey area indicates the upper and lower bounds for K(θ) for C between 2 and 0.5 mg m-3................................... 86 Figure 6.5: Progression of the distributions of PAR over Pavilion Lake on a typical summer day, June 28th a) 8:00 am, b) 8:30 am, c) 11:15 am, d) 2:15 pm, e) 2:45 pm, and f) 4 pm. Contour colors represent the intensity of the incoming solar radiation, ET in W m-2. The lake is outlined in black, (+) and ( ) are the locations of the thermistor chain and the weather station...................................................................................................... 87 Figure 6.6: Running weekly average of calculated attenuation coefficient at z = 5.8 m (blue), z = 17 m (red), and a depth averaged K(dashed black). Grey lines delineate the icecover period (Lim et al., 2009). ........................................................................................ 91 Figure 6.7: a) Diurnal PAR at the sea surface, b) Decomposition of PAR into Ediffuse (black dashed) and Edirect (black solid), c) Attenuation coefficient calculated based on the PAR decomposition without including any change in chlorophyll concentration. ......... 92 Figure 6.8: (a) Wind speed, and (b) air temperature observed at Pavilion Lake for the studied period in 2006....................................................................................................... 94 xi Figure 6.9: Temperature at 9.5 m, 15.5 m, and 26.5 m; observed (black), estimated using constant attenuation coefficient (blue), and variable attenuation coefficient (red). ......... 93 Figure 6.10: (a) Wind speed, (b) air temperature, and (c) Vapor Pressure observed at Pavilion Lake for the model period in 2006. .................................................................... 95 Figure 6.11: Running weekly averages of the temperature at 9.5 m, 15.1 m, and 24.1 m; observed (black), estimated using constant attenuation coefficient (blue), and temporally variable attenuation coefficient (red). ............................................................................... 96 xii LIST OF SYMBOLS Roman Characters: a Water absorption function apw Pure water absorption function b Water scattering function B(z) Damping term that depends on the particulate biogeochemical matter in the C Chlorophyll concentration c Speed of light Cturb The sum of the biogeochemical particulate matter concentrations cp Specific heat capacity of water Ediffuse Diffused solar radiation Edirect Direct solar radiation ET Total solar radiation ET(z) Irradiance intensity profiles F Fraction of Ediffuse in the incoming radiation Fr Froude number h Planck’s constant � 𝐾 Non-spectral attenuation coefficient �1% 𝐾 Spectrally independent attenuation coefficient calculated using the photic xiii �𝑙𝑠 𝐾 Spectrally independent attenuation coefficient defined as a least square fit � direct(θ) 𝐾 Non-spectral attenuation coefficient for direct PAR under different zenith � diffuse 𝐾 Non-spectral attenuation coefficient for diffuse PAR to a profile of light intensity angles � (PAR, θ) 𝐾 Bulk temporally variable attenuation coefficient for PAR Kc Self-shading constant Kdiffuse Attenuation coefficient for diffuse light Kdirect Attenuation coefficient for direct light Kw Attenuation coefficient for pure water m Ratio between the mass of the air column encountered by the solar Qlh Latent heat flux Qlw Longwave flux Qsh Sensible heat flux Rn Reynolds number T Absolute temperature t Time Ta Air temperature Tr Percentage transmission Tw Surface water temperature xiv Uw Wind speed z Depth below water surface zmax Depth where λ becomes λmax z1% Photic depth Greek Characters: α Water surface albedo β Mass attenuation coefficient ε Light backscattering probability in a cloudless atmosphere ηk Shear production efficiency ηs Wind stirring efficiency γ Latitude less the solar declination ρ Density σ Stephan Boltzmann constant τ Optical thickness of the atmosphere θ Solar zenith angle xv List of Abbreviations: AOP Apparent Optical Properties CAEDYM Computational Aquatic Ecosystem DYnamics Model CTD Conductivity–Temperature–Depth DYRESM DYnamic Reservoir Simulation Model ELCOM Estuary, Lake and Coastal Ocean Model GOTM General Ocean Turbulence Model IOP Inherent Optical Properties IR Infra-Red NIR Near Infra-Red PAR Photosynthetically Active Radiation RH Relative Humidity SD Secchi Depth TKE Turbulent Kinetic Energy UV Ultra-Violet UVA Ultra-Violet A UVB Ultra-Violet B xvi ACKNOWLEDGMENTS I would like to express my deepest thanks and gratitude to my advisor, Professor Gregory Lawrence. I am indebted to the time and effort he took upon himself to offer me advice. I am also thankful for his constant consideration and support. I will never forget his assistance and will always feel fortunate for having him as my advisor. I feel deep appreciation to Dr. Roger Pieters for providing me with continued exposure and advice. His perfection motivated me to achieve a higher level of understanding my work. His accessibility, and continued advice were very constructive. This work would not have been possible without Dr. Bernard Laval who offered data, and assistance. Special thanks go to Dr. Dan Moore for his helpful suggestions and positive feedback. I feel blessed to have close friends here in Vancouver. They have supported, encouraged me throughout the years. The laughter and joy they brought to me during my low moments are priceless. I will never forget their true friendship and support. I would like to thank Rio Tinto Alcan who has partially sponsored the research conducted in this dissertation. My dearest husband, friend, and study mate Yehya has always been by my side with endless support and encouragement, especially when I was in low spirits. Our valuable discussions, as colleagues, contributed significantly to the improvement of my research. Finally, special thanks are owed to my parents and brothers. They have encouraged and believed in me for many years. They are the ones who are far away but yet closest with their hearts. I hope I can always make them happy and proud. xvii DEDICATION To my parents Chapter 1. Introduction 1 1 INTRODUCTION 1.1 Background Knowledge of light levels below the water surface is crucial in both biological and physical limnology. From a biological viewpoint, the distribution in light is needed to assess the growth of phytoplankton, whereas from a physical perspective, it is of primary importance in determining water temperature. The distribution of light in lakes is determined by the solar radiation incident on the water surface, and its attenuation as it passes through the water column. Numerous factors play a role including cloud cover, the spatial distribution of the light field at the water surface, surface waves, air quality and water clarity, and each of these factors can vary temporally and/or spatially. Here we will focus on one facet of this very broad problem, namely: how is the temperature distribution in lakes affected by the fact that the attenuation coefficient is often strongly dependent on wavelength? While several oceanographic studies have investigated the effects of the spectral variability of light attenuation (Paulson and Simpson, 1977; Schneider and Zhu, 1998; Chang and Dickey, 2004; Simpson and Dickey, 1981), their primary concern was surface heating and they focused on infrared radiation (IR; wavelengths between 0.7 and 3 µm) rather than ultra-violet (UV; 0.02 0.4 µm) or Photosynthetically Active Radiation (PAR; 0.4 - 0.7 µm). In this study we focus on PAR irradiance since it can penetrate to great depths, particularly in clear waters. In limnological models the spectral variability of the attenuation coefficient is often simplified by assuming that UV and IR radiations are absorbed at the surface, while the PAR is distributed throughout the water column using a constant attenuation coefficient (Spigel and Imberger, 1980; Hornung, 2003). This parameterization is subject to error as a result of the high spectral variability of light attenuation within the PAR domain. Chapter 1. Introduction 2 In this study we will assume that each lake has its characteristic spectral variation in attenuation coefficient, but that this variation does not change with depth or position within the lake. We will also address diurnal variations due to changes in solar angle and seasonal variations in the attenuation coefficient due to changes in water clarity, but the major focus will be on the effects of the spectral variation in attenuation coefficient. 1.2 Research objectives The primary objective of this research is to assess the importance of accounting for the spectral variability in the attenuation coefficient, K(λ), for the PAR domain in particular, on the evolution of the thermal structure of inland water bodies with different water clarities. The secondary objective is to evaluate two aspects of the temporal change in the attenuation coefficient resulting from solar angle and from seasonal changes in water quality. To achieve these research objectives, the following steps were taken: 1. Summarize common approaches used in calculating the subsurface irradiance. 2. Understand how the subsurface irradiance calculated using the spectral variation in K(λ), differs from the subsurface irradiance estimated using a non-spectral attenuation coefficient, for different water clarities. 3. Test how the difference between these two profiles of subsurface irradiances influences the thermal structure of lakes: a. Using solar radiation in isolation and examining profiles of light and heat fluxes estimated from each, and b. Using a version of the Dynamic Reservoir Simulation Model (DYRESM) modified to account for the wavelength dependence. 4. Understand the effect of light attenuation on the development of heat fluxes at the water surface. 5. Examine the diurnal and seasonal variations in the attenuation coefficient. 6. Examine the significance of accounting for the temporal variations in light attenuation on thermal structure of a clear lake using a modified version of DYRESM. Chapter 1. Introduction 3 The structure of this thesis is as follows. Chapter 2 reviews the relevant literature. In Chapter 3, the study sites are described along with the available meteorological and optical data for the lakes investigated in this study. In addition, Chapter 3 describes the hydrodynamic model (DYRESM) that will be used and the modifications made to it. Chapter 4 illustrates the importance of accounting for the spectral variability in K(λ) on irradiance and thermal distributions for different lakes with different water clarities. In Chapter 5, the modified version of DYRESM is used to estimate the long-term effect of accounting for this spectral variability for a relatively clear lake. Chapter 6 illustrates the dependence of the attenuation coefficient on the solar zenith angle and how accounting for seasonal variations in the attenuation coefficient improves the accuracy of predictions of the thermal structure of a clear lake. Finally, Chapter 7 summarizes the results, conclusions and possible future work. Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 4 2 RADIATION TRANSFER FROM SUN TO WATER BODIES AND ITS EFFECT ON WATER TEMPERATURE: A LITERATURE REVIEW 2.1 Introduction A proper understanding of the fate of shortwave radiation as it enters, and propagates into natural waters is needed to accurately estimate the flux of heat into the water column. The intensity of shortwave radiation and its attenuation are important factors in the development and evolution of thermal stratification in lakes (Bloss and Harleman, 1979; Pimentel et al., 2008). The following literature review focuses on the physical meaning of the attenuation coefficient (K) and its parameterization in physical limnology. Section 2.2 discusses solar radiation and its propagation through air and water. In section 2.3 we introduce different factors that affect subsurface radiative transfer that play an important role in determining the value of the attenuation coefficient. Section 2.4 illustrates different approaches used for predicting the attenuation coefficient. Section 2.5 discusses the parameterization of the attenuation coefficient in several hydrodynamic lake models. 2.2 Attenuation coefficient in physical limnology Solar radiation is the primary source of energy input to lakes. The sun is a blackbody that emits at every wavelength as shown in Figure 2.1. The distribution of this intensity can be estimated using Planck’s law, Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 5 𝐼(λ) = 2ℎ𝑐 2 λ5 1 ℎ𝑐 𝑒 λσ𝑇 −1 (2-1) where h is the Planck constant, c is the speed of light, and σ is the Stephan Boltzmann constant (5.67*10-8 W m-2). T is the absolute temperature in Kelvin and the effective temperature of the sun is 5780 K. The peak in this distribution is given by λ= hc 1 4.96 σ T (2-2) which is often referred to as Wein’s law. For the sun, the irradiance peaks at λ = 0.502 μm (green light). Integrating Planck’s law over all wavelengths gives the total radiation from a black body as the Stephan Boltzmann law, ET = 𝜎𝑇 4 ; where ET (W m-2) is the total radiation per unit area from the surface of the black body. The black body radiation from the sun can be divided into bands, including Xrays (0.00001 – 0.02 μm), ultraviolet radiation (UV; 0.02 – 0.4 μm), photosynthetically active radiation (PAR; 0.4 – 0.7 μm), infrared (IR; 0.7 - 3 μm), and microwave (3 μm – 1 m). The distribution and intensity of shortwave radiation arriving at sea level is different than the extraterrestrial spectrum because much of this radiation is absorbed and scattered by gases and aerosols in the atmosphere. For X-rays and microwave radiation, not only the irradiance from the sun is small, but irradiance in these bands is further reduced by high absorption and scattering in the atmosphere. Therefore, in this study we are only concerned with shortwave radiation which includes IR, PAR, and UV. The percentage of each in the incoming shortwave radiation at sea level is 7% UV, 51% PAR, and 42% IR (Gast, 1960; Jerlov, 1968; Matthews, 1987). Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 6 2500 X UV PAR IR Microwave -2 -1 Irradiance (Wm µm ) 2000 Solar Irradiance at top of atmosphere 1500 Solar Irradiance at sea level 1000 500 0 0 0.5 1 1.5 2 2.5 3 3.5 4 Wavelength λ (µm) Figure 2.1: Solar radiation as function of wavelength at the top of the atmosphere and at sea level, adapted from Mecherikunnel and Duncan (1982). The total solar radiation reaching the surface of the earth is a function of the solar zenith angle, θ, ET = (E solar cosθ ) (1 + ετsecθ ) (2-3) where; Esolar is the extraterrestrial irradiance with an average value of 1367 W m-2; θ is the solar zenith angle measured from the vertical; ε is the light backscattering probability which is approximately 0.2 in a cloudless atmosphere (at λ = 0.55 μm), and τ is the Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 7 optical thickness of the atmosphere (at λ = 0.55 μm) which is approximately 0.3 (Bukata, 1995). For cloudless sky, the radiation reaches the earth as either direct Edirect(θ), or diffuse, Ediffuse(θ), radiation, ET (θ ) = Ediffuse (θ ) + E direct (θ ) (2-4) where Edirect can be calculated using, E direct (θ ) = E solar cosθexp(− τsecθ ) (2-5) The direct radiation is the radiation that travels on a straight line from the sun to the earth, and the diffuse radiation is the one that gets scattered by air molecules and other particles in the atmosphere. On a cloudless day with the sun vertically overhead (θ = 0), a typical value of ET(0) is approximately 1290 W m-2, of which Edirect(0) is 1010 W/m2 and Ediffuse(0) is 280 W m-2. The two components, Edirect(θ) and Ediffuse(θ), vary with time of day. The above equations are implicitly a function of beam absorption and scattering through the path length of the beam. As θ increases, the longer the straight path of light through the atmosphere gets, the higher is the percentage of Ediffuse due to scattering. Ediffuse, Edirect and ET are shown in Figure 2.2 over the PAR band which will be the focus here. Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 8 Figure 2.2: Spectral distribution of a) Ediffuse, b) Edirect, and c) ET for different sun angles at sea level. Values plotted are relative to the maximum light intensity, I(λ = 0.48 μm) = 1566 W m-2 μm-1, for vertical sun (θ = 0), based on equations (2-4) and (2-5). The previous is a brief explanation of the path of light from the top of the atmosphere to the water surface, showing how the radiation spectrum and composition reaching the water surface is considerably modified from the extraterrestrial radiation by the atmosphere. Further detail can be found in Jerlov (1968), Kirk (1994), and Preisendorfer (1965). Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 9 2.3 Factors affecting the subsurface radiative transfer Transfer of shortwave radiation in water is determined by the behavior of photons, which can be absorbed or scattered by water molecules, and which varies with wavelength (λ). The attenuation of solar radiation through water can be calculated using Beer’s law (also known as the Beer-Lambert law), 𝐼(𝑧, 𝜆) = 𝐼(0, 𝜆)exp−𝐾(𝜆)𝑧 (2-6) where, I(z, λ) is the downwelling intensity of shortwave radiation of wavelength λ at depth z (m), I(0, λ) is the intensity of the solar radiation of wavelength λ at the surface, and K(λ) (m-1) is the attenuation coefficient for a given λ. The spectral attenuation coefficient of light, K(λ), in pure water under clear sky conditions with no wind is presented in Figure 2.3. The attenuation coefficient can be influenced by many factors that can be classified as either inherent optical properties (IOP), properties that depend solely on the water and its impurities; or apparent optical properties (AOP), properties that depend on the geometrical distribution of light field such as the solar angle, clouds, and surface waves (Anthony et al., 2004). In this section the role of these parameters is briefly described. Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 10 10 10 -1 K(λ ) (m ) 10 10 10 10 10 10 5 4 UV PAR IR 3 2 1 0 San Vicente Superior -1 Crater Pure -2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 λ (µm) Figure 2.3: Average attenuation coefficient, K(λ), under clear sky and no wind, for pure water (black), Crater Lake (blue), Lake Superior (green), and San Vicente Reservoir (red). Data from Bukata (1995) and Tyler and Smith (1970). Different impurities present in natural waters that affect its clarity include: dissolved organic substances (yellow substances) resulting from the decomposition of plant tissues and industrial discharges; chlorophyll in phytoplankton; and organic and inorganic particulate matter (tripton). Yellow substances absorb the blue part of the PAR spectrum significantly. Chlorophyll strongly absorbs the red and the blue parts of the PAR spectrum. Particulate matter scatters the photons significantly. Any increase in the amount of these impurities increases the amount of light which is absorbed and scattered, through increasing the overall path of light (Jerlov, 1968; Hieronymi, 2011). Note that pure water has the lowest K(λ) at all wavelengths (Figure 2.3). Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 11 Another factor that plays a role in determining the subsurface irradiance is the fraction of Ediffuse in the incoming radiation. The fraction of Ediffuse, depends on the solar zenith angle and the cloud cover. Both Edirect and Ediffuse, have a different dependence on wavelength (λ) as will be further discussed in section 2.4.2. A study done in the Sargasso Sea observed a large variation in subsurface PAR from clear to cloudy sky (Stramska and Dickey, 1998). There is a strong relation between waves and wind, and the penetration of light into the water column. Wind and waves cannot be separated as the wave heights depend on the fetch and the duration the wind is active in. The slope and the height of the superposed waves act as converging lens altering the irradiance below the water surface (Dera and Gordon, 1968). The lensing effect of surface, capillary and gravity waves, are one of the main reasons behind irradiance fluctuations below the water surface (Hieronymi and Macke, 2010). However, different superimposed wave types affect the irradiance differently; for example waves of heights higher than 1.5 m can affect the irradiance distribution down to a depth of 100m in sea water (Hieronymi, 2011). The wave effect can be eliminated by analyzing a period with no winds (Stramska and Dickey, 1998). Wind velocities that significantly affect the spectral irradiance distribution depend on the given water body. In a study in the Mediterranean Sea, the maximum fluctuations observed in the irradiance intensity profiles, ET(z), were for periods with wind speeds between 1 and 5 m s-1, the wind generated waves of ~ 0.5 m height (Gernez and Antoine, 2009). In another study in the Atlantic Ocean, wind velocities that altered the irradiance intensely were between 2 and 7 m s-1 (Hieronymi, 2011). For the purposes of the present study it is assumed that water constituents are homogenously spread throughout the water column, and for the days that light profiles were measured, there were no clouds or wind. Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 12 2.4 Determination of attenuation coefficients In the following, some of the methods used to estimate the attenuation coefficient are described along with studies that have used these methods. 2.4.1 Non-spectral attenuation coefficient The accuracy of approximating the irradiance below the water surface using a temporally and spectrally constant attenuation coefficient depends on the water clarity. However, because it is easy to estimate, most physical limnological studies use a constant coefficient, e.g. Yeates and Imberger (2003). We shall adopt the convention that a nonspectral attenuation coefficient refers to an attenuation coefficient that is constant over the PAR domain; UV and IR are assumed to be absorbed near the water surface or, in the case of layered models, in the upper layer. We will now discuss the most common methods used to determine non-spectral attenuation coefficients. Quantum light meters Quanta meters measure the irradiance intensity for a specific band in the electromagnetic spectrum, usually PAR (0.4 - 0.7 μm). By measuring the intensity at two different depths, z1 and z2, a non-spectral attenuation coefficient can be derived from Beer’s Law: 1 𝐸(𝑧 ) 𝐾 = − ∆𝑧 ln 𝐸(𝑧1 ), 2 (2-7) where ∆𝑧 = z1-z2. Secchi Depth A non-spectral attenuation coefficient can be estimated from the Secchi depth. The Secchi depth is the depth at which a circular disk disappears in the water. For a valid Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 13 Secchi depth, the measurement should be taken at noon, under cloudless sky and low wind speed. Secchi depths are a rough measure of water clarity which are useful, for example, in determining the trophic status of lakes (Table 2.1). Table 2.1: Trophic class according to Secchi depth adapted from Preisendorfer (1986) Trophic Class Secchi Depth (m) Ultra-oligotrophic >8 Oligotrophic 8–4 Mesotrophic 4-2 Eutrophic 2 – 0.5 Hyper-eutrophic < 0.5 Several relationships have been suggested for calculating K from the Secchi depth (SD), such as K = SD/1.7 given by (Poole and Atkins, 1929). However, most relationships have been in the form of 𝐾 = 𝑎1 𝑆𝐷−𝑎2 . For example, 𝐾 = 2.14 𝑆𝐷−0.87; was used by Reinart and Herlevi (1999), who nevertheless doubted its validity for 𝐾 > 4.6 m−1. This type of relationship for K is inaccurate, and the value of K estimated from Secchi depth remains approximate at best; for further detail see Preisendorfer (1986). Temperature Profiles A non-spectral attenuation coefficient can also be estimated using a dense grid of temperature profiles as done in Lake Kinneret, Valle de Bravo, and Winam Gulf by Adiyanti and Imberger (2007). The heat was calculated in different layers and an averaged attenuation coefficient was estimated from the increase in temperature in Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 14 subsequent profiles. Nevertheless, the method requires fine temperature measurements in the photic zone and the surface layer needs to be stratified (e.g. under low wind speed) to isolate the solar heat flux. In this technique, it is assumed that the water is mixed in the horizontal direction and the results are not very accurate at the water surface (Adiyanti and Imberger, 2007). 2.4.2 Spectrally variable Attenuation coefficient The importance of estimating ET(z) using the full K(λ) divided into different bands versus a non-spectral K was first investigated by Dickey and Simpson (1983). Nevertheless, this study was done on clear oceanic water and the whole PAR region was presented using a single band. Calculations based on non-spectral attenuation coefficients were also proposed by Zaneveld (1989) and Anderson (1993) who suggested that the use of a spectral attenuation coefficient can reduce errors in the calculation of photosynthesis rates by up to 50% for oceanic waters. However, lake studies often lack the measurements required to determine the K(λ) spectra. Given below are different methods used to determine K(λ). Spectroradiometers Spectroradiometers are devices that can measure the spectral power distribution of light. They measure the power of each color separately. Attenuation coefficient spectra can be estimated by applying Beer’s law between surface and subsurface irradiance measurements, for each wavelength. The K(λ) spectrum for the Pacific Ocean near California has been evaluated using a single irradiance collector that measured the subsurface radiation for 13 different wavelengths in the PAR region (Stramski, Booth, and Greg Mitchell, 1992). Another example is the K spectrum for the Sargasso Sea which has been estimated by Stramska and Dickey (1998) with the aid of a spectroradiometer. This time Beer’s law was used between two subsurface spectral irradiances measured at 15 m and 35 m depths. Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 15 Biological models Several empirical relationships have been proposed in order to determine K(𝜆) from, for example, the scattering coefficient at 𝜆 = 0.55 μm, particulate organic concentrations, incoming radiation, and the irradiance intensity at several depths (Morel, 1988). These empirical relationships considered the water as two component systems (pure water and a biogenic component) and they are only applicable to waters in which chlorophyll is the only impurity. Most of these empirical relationships used chlorophyll concentrations in calculating the attenuation coefficient through calculating particle scattering and absorption for different wavelengths (Baker and Smith, 1982), and have been widely used despite their uncertainties. A linear relationship between the chlorophyll concentration and the attenuation coefficient was suggested by Vila et. al (1996); however, other non-linear relations have also been suggested, such as the relation by Morel and Maritorena (2001), 𝐾(𝜆, 𝐶) = 𝐾𝑤 (𝜆) + 𝜒(𝜆)𝐶 𝑒(𝜆) (2-8) where, 𝐶 is the chlorophyll concentration (mg m-3), 𝐾𝑤 (𝜆) is the spectral attenuation coefficient of pure water, and 𝜒(𝜆) is an empirical coefficient. Similar relations were suggested and tested by Reinart and Herlevi (1999) on oceanic waters, but this time using water properties at 𝜆 = 0.49 μm, the wavelength corresponding to the maximum incoming radiation. This model was further developed to include waters with impurities other than chlorophyll, and was tested on 43 Estonian lakes by Paavel et al. (2006). Note that the relationship suggested by Paavel et al. (2006) requires measurements at three different wavelengths. Kirk (1979) analyzed the dependence of the attenuation coefficient spectra on water impurities for some inland waters. The study gave a detailed explanation of the effect of different components such as the effect of yellow substances on different bands Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 16 in the K(𝜆) spectrum. Attenuation coefficient values were estimated using the irradiance spectra measured by the spectroradiometer at different depths, and the behaviors of the spectra were explained in terms of impurity concentrations. The resulting relationships and their applicability were further investigated by Morel (2009). Remote Sensing Remote sensing is one of the most popular methods to calculate the attenuation coefficient spectrum (Liu et al., 2006). Advances in accurately sensing water color by remote sensing have facilitated the determination of water optical properties. It is possible to estimate the spectral intensity of light based on reflection characteristics of the water surface using two methods. The first method uses the reflected light at λ = 0.49 μm from which K(0.49 μm) can be estimated and used as the non-spectral attenuation coefficient. However, this is only applicable for hyper-eutrophic waters K(0.49 μm) < 0.25 m-1 (Mueller, 2000). Estimating the attenuation coefficient spectrum using K(0.49 μm) is also possible based on the strong correlation between K(0.49 μm) and the rest of the spectrum. The second method in estimating the attenuation coefficient spectrum is through determining chlorophyll concentration from the absorption and scattering properties. Nevertheless, the attenuation coefficient spectrum estimated based on the sensed chlorophyll concentration is not as accurate as the one estimated using the absorption and scattering coefficients (Lee, Carder, and Arnone, 2002) and it is only reliable in waters where the chlorophyll is a dominating component in water clarity. Note that both methods are based on the assumption of vertical homogeneity (Lee et al., 2005). 2.4.3 Attenuation coefficient as a function of wavelength and solar angle Ideally, the attenuation coefficient spectra calculated in the previous section would vary with time. As mentioned before K(λ) is dependent on the solar zenith angle. The Monte Carlo simulation model first proposed by Kirk (1984) managed to estimate Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 17 attenuation coefficient spectra for diffuse and direct light, as a function of the zenith angle, as follows: K direct ( λ, θ ) = K diffuse ( λ ) = a (λ ) [1 + (0.425 cos (θ ) − 0.19b(λ ) / a (λ ))]0.5 cos (θ ) a (λ ) [1 + 0.162b(λ ) / a (λ )]0.5 0.856 (2-9) (2-10) where θ, is the solar zenith angle, and a and b are the absorption, and the scattering functions. Values of a and b depend on Edirect(λ) and Ediffuse(λ) at the surface, and on the amount of dissolved organic substances, and other impurities. All these factors affect Kdiffuse(λ) and Kdirect(λ, θ) (Jerlov, 1968). The previous equations give the role of θ in determining Kdirect(λ), while Kdiffuse(λ) is constant for all θ as shown in Figure 2.4. The spectral attenuation coefficient of the global radiation; K(λ) can be expressed as a function of F = Ediffuse/ET, where K (λ ) = FK diffuse (λ ) + (1 − F )K direct (λ ) (Bukata, 1995). For pure water, a and b values are known (Appendix B); nevertheless, they are hard to determine for water with impurities. Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 18 Figure 2.4: Spectral distribution of a) Kdiffuse(λ) and b) Kdirect(λ, θ) for pure water, under different sun angles, based on equations (2-9) and (2-10). Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 19 2.5 Parameterization of the attenuation coefficient in physical limnology 2.5.1 Vertical Many researchers have observed that the vertical distribution of solar radiation often does not satisfy Beer’s-law if a single (non-spectral) attenuation coefficient is used. This deviation has often been attributed to a vertical variation in K due to a vertical variation in water properties, such as chlorophyll concentration (Zielinski et al., 1998, Lee et al., 2005, Losordo and Piedrahita, 1991). However, the discrepancy may be due, at least in part, to the wavelength dependence of both the incoming solar radiation and the attenuation coefficient in the water column. 2.5.2 Temporal The effect of temporal variability in the attenuation coefficient, on the stratification of lakes has not been deeply investigated. It is known, for example, that heat transfer and thermocline depth are dependent on the biological composition of water (Pérez-Fuentetaja et al., 1999) and that this composition is temporally variable. Nevertheless, studies that focus on the temperature structure in lakes, and which do not consider biological processes, has often assumed a temporally constant attenuation coefficient. Dejak et al. (1992) studied the effect of the temporal variation of the incoming shortwave on the thermal distribution below the water surface, the heat budget and the timing of the spring and fall turnovers. The study was done without considering the temporal variation in K, however, deeper understanding of other heating factors was recommended. Tanentzap et al. 2008 have demonstrated how variations in the annual average attenuation coefficient play a major role in the thermal evolution of Clearwater Lake. Together with other research (Han et al., 2000; Tanentzap, Hamilton, and Yan, 2007; Perroud et al., 2009), they have suggested that a temporally variable attenuation Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 20 � (𝑡), should replace the constant coefficient in the one dimensional coefficient 𝐾 hydrodynamic model, DYnamic Reservoir Simulation Model (DYRESM). 2.5.3 Physical Modeling Surface heat fluxes and wind are the major drivers to the evolution of the thermal structure of lakes. Mixing in the surface layer is maintained through three mixing energies driven by those two forcings (Imboden and Wuest, 1995); wind, which triggers water stirring; buoyancy fluxes due to the cooling of the water surface; and finally shear, which can occur between water layers and is also driven by wind. These energies are usually redistributed in the surface layer. Under weak wind conditions, stirring and convection are the dominant sources of energy (Spigel and Imberger, 1980). Nevertheless, under strong wind conditions, shear dominates, until the interface between the surface layer and the layer below is set up, that is when the baroclinic circulation dominates (Imberger and Hamblin, 1982). If the density stratification is weak; those baroclinic waves can cause mixing in the hypolimnion. The role of the spectral attenuation coefficient was often omitted in the classical literature when dealing with mixing and the lake’s thermal structure. Processes including position of the thermocline, penetrative convection, and hypolimnetic mixing were often related to internal properties of the system (Imberger and Ivey, 1991) and not to the spectral attenuation of light. This being said, Jassby and Powell (1975) suggested adding a term in the hypolimnion heat budget that accounts for solar radiation. Following are some widely used models and how they distribute the shortwave irradiance below the water surface. ELCOM: The three dimensional model ‘Estuary, Lake and Coastal Ocean Model’ is used to predict the thermal evolution of reservoirs. It divides the incoming shortwave into four bands, Photosynthetically Active Radiation (PAR), Near InfraRed (NIR), Ultra Violet A (UVA), and Ultra Violet B (UVB). The contribution of each band Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 21 is 45%, 41%, 3.5%, and 0.5%, respectively, and each has a corresponding attenuation coefficient. Note that the summation is only 90%, underestimating the UV and IR by 10%, which might result in underestimating the predicted surface temperature. The K default values set by the model are 0.25 m-1, 1 m-1, 1 m-1, and 2.5 m-1, respectively, and can be altered by the user. The shortwave radiation is then distributed using Beer’s law throughout the depth (Hodges and Dallimore 2006). In case the water quality model ‘Computational Aquatic Ecosystem DYnamics Model’ (CAEDYM) is coupled with ELCOM, CAEDYM recalculates the attenuation coefficient for the PAR band based on the water quality which means that it becomes temporally variable; however, the shortwave radiation is still distributed using Beer’s law, with no account for spectral variation in any of the bands. CE-QUAL-W2: The two-dimensional model CE-QUAL-W2 is used in modeling the hydrodynamics and water quality of reservoirs, rivers, and river basin systems. In the hydrodynamic module of the model, the subsurface irradiance is calculated using a bulk extinction coefficient K, which is a single constant value entered by the user and used throughout the whole simulation period, to represent the full shortwave spectrum. On the other hand, in the biological module, a different approach is utilized; the user enters the fraction of the incoming shortwave that is absorbed at the surface, and a temporally variable attenuation coefficient, which depends on the concentration of the inorganic and organic suspended solids, is calculated and used to distribute the remaining shortwave radiation (Cole and Buchak, 1995). DYRESM: The one-dimensional hydrodynamic model, DYnamic REservoir Simulation Model (DYRESM) is used for predicting the evolution of temperature, salinity and density structure, in lakes and reservoirs. The heat module in this model divides the shortwave radiation into two fractions, 55% representing the IR and UV radiation. This fraction is transferred to the surface layer regardless of its thickness. The remaining 45% represents the PAR, and is distributed throughout the water column using Chapter 2. Radiation Transfer from sun to water bodies and its effect on water temperature: A literature review 22 Beer’s law. A single K value, which is not temporally or spectrally variable, is entered by the user (Antenucci and Imerito, 2001). Just like ELCOM, DYRESM can be coupled with CAEDYM. In that case CAEDYM recalculates K at every time step based on the concentration of chlorophyll and other parameters. GOTM: The one-dimensional model, General Ocean Turbulence Model, is used to model thermodynamic and hydrodynamic processes. Shortwave radiation in the temperature calculations is distributed using 𝐸𝑇 (𝑧) = 𝐸𝑇 (0)(𝐴𝑒 − 𝑧 𝜗1 + (1 − 𝐴)𝑒 − (2-11) 𝑧 𝜗2 )𝐵(𝑧) where 𝜗1 and 𝜗2 are absorption coeﬃcients in (m) suggested by Paulson and Simpson (1977) and depend on the water type (Jerlov, 1968). B(z) is a damping term that depends on 𝑒 the particulate 0 −𝑘𝑐 ∫𝑧 (∑ 𝐶𝑡𝑢𝑟𝑏(𝜀) 𝑑𝜀 ) biogeochemical matter in the water column, 𝐵(𝑧) = , where Kc is the self-shading constant, Cturb is the sum of the biogeochemical particulate matter concentrations (Umlauf et al., 2006). The model basically divides the total incoming shortwave radiation into two fractions; A and (1-A), and it distributes the total incoming irradiance using a bi-modal exponential form, taking into account the vertical variability in the distribution of biogeochemical matter using B(z). Chapter 3. Sites and methods 3 23 SITES AND METHODS 3.1 Introduction In this chapter, we describe the data and methods used throughout. In Section 3.2, we describe the physical characteristics of the studied lakes, and summarize the meteorological, optical and water temperature data for each. In Section 3.3, a commonly used one dimensional model is described; this model is used to predict the thermal structure of lakes in subsequent chapters. The modifications to this model, made to include the spectral and temporal variability of the attenuation, are also described. Section 3.3 also includes the methods used to calculate the attenuation coefficient for different bands in the photosynthetically active radiation (PAR; λ = 0.4 – 0.7 µm) domain. 3.2 Study sites This section describes the study sites and the various data that are used for each. The first four sections describe the spectrally-variable attenuation coefficient data, K(λ), that are available for pure water and three different lakes, which span a range of turbidity levels. Pure water is examined to test the effect of maximum clarity on the irradiance intensity below the surface. It is followed by Crater Lake which is one of the clearest lakes in the world, then Lake Superior and San Vicente Reservoir as examples of lesser clarity. The last two sections describe Tailings and Pavilion Lakes. For the ultraoligotrophic Pavilion Lake there are a wide range of temperature, light, weather, and topographic data available. These data are used to test the significance of accounting for Chapter 3. Sites and methods 24 the wavelength dependence and the time variation of K on the evolution of the thermal structure of the lake. 3.2.1 Pure Water For pure water, data were compiled by Buiteveld et al. (1994) who presented the attenuation coefficient spectrum at 20oC from λ = 0.02 – 3 μm, covering ultraviolet (UV; λ = 0.02 - 0.4 µm), PAR, and infrared (IR; λ = 0.7 – 3 µm) radiation bands. Pure water has the lowest K(λ) as shown in Figure 2.3, with a minimum K(λ) value of 0.01 m-1 at λ = 0.42 μm. Pure water has an estimated Secchi depth of 79 m based on its photic depth. 3.2.2 Crater Lake Crater Lake (42.9° N, 122.1° W) located in south-central Oregon, United States, is formed in the bowl shaped caldera basin of Mount Mazama. Its surface is almost circular in shape with a mean diameter of 8 km, and a surface area of 53 km2. It has a maximum depth of 595 m, and an average depth of 350 m. Crater Lake is the deepest lake in the United States, and is the seventh deepest in the world. The lake does not get ice covered during winter (Larson et al., 2007). The lake filled with the accumulation of rain and snow over approximately 250 years. Now the lake maintains its water level as incoming precipitation balances evaporation and seepage losses. There are no tributary inflows or outflows. Because Crater Lake is primarily fed by precipitation, it is ultra-oligotrophic, and one of the clearest lakes in the world with Secchi depth down to 41.5 m and a maximum photic depth of 79 m (Larson et al., 2007). There are few organic materials in the water, and the primary dissolved and particulate materials come from the caldera walls and natural algal blooms. The spectrum of the attenuation coefficient was collected at three different depths in August by Tyler and Smith (1970). No significant differences are observed between Chapter 3. Sites and methods 25 the spectra at the different depths and it is assumed that optically the lake is vertically homogenous. The average of the measured K(λ) at all three depths is presented in Figure 2.3. The minimum K(λ) value is 0.016 m-1 at λ = 0.425 μm (Figure 2.3). 3.2.3 Lake Superior and San Vicente Reservoir Lake Superior (48°N, 88°W) is located in North America, between Wisconsin and Michigan in the United States, and Ontario in Canada. The maximum depth of the lake is 405 m, making it the deepest of the Great Lakes. It has a surface area of 82,100 km2, which makes it the third largest freshwater lake in the world. Many tributaries flow into Lake Superior, and as a result, three types of dissolved and suspended materials are found in the Lake: yellow substances, phytoplankton, and organic and inorganic particulate matter. These materials do not vary much between summer and winter (Dobson, Gilbertson, and Sly, 1974). Lake Superior is oligotrophic, with a Secchi depth of 8 m, which is the maximum recorded for all the Great Lakes. The Secchi depth measurements remain almost constant throughout the year. The minimum K(λ) value is 0.35 m-1 at λ = 0.57 μm (Figure 2.3). San Vicente Reservoir (32.9° N, 116.9° W) is located 40 km northeast of San Diego, California, United States. The total surface area of the reservoir is 4.3 km2, with a maximum depth of 58 m. San Vicente Reservoir is eutrophic, with a very shallow photic depth. The minimum K(λ) value is 0.51 m-1 at λ = 0.58 μm (Figure 2.3). 3.2.4 Tailings Lake Tailings Lake (64.3°N, 115°W) is located in the Northwest Territories, 220 km north of Yellowknife, Canada. The lake is part of the Colomac gold mine. The length of the lake is 1.4 km, the width is 0.35 km, the surface area is 0.5 km2, and the maximum depth is 14 m, with a volume of 0.0034 km3. Chapter 3. Sites and methods 26 The lake filled from 1990 to 1997, with solid tailings and wastewater during mining operations (Pieters and Lawrence, 2009). Tailings Lake has a small drainage area of 3.16 km2 and precipitation contributes to the inflow to the lake. The lake is abundant with yellow substances, phytoplankton, and organic and inorganic particulate matter. Its water is very opaque with a Secchi depth of less than 0.5 m, which is the shallowest secchi depth for the cases described here. 3.2.5 Pavilion Lake Site description Pavilion Lake (50.9° N, 121.7° W) is located in Marble Canyon, a limestone valley in south central British Columbia, Canada. Pavilion Lake is long and narrow with two constrictions dividing the lake into three major basins: north, central, and south (Figure 3.1). The length of the lake is 5.7 km and the maximum depth is 65 m deep in the main basin. Inflow to Pavilion Lake is groundwater dominated, despite inflows from Big Bar Creek and Bonaparte River. Pavilion Lake is dimictic and ice forms annually. The lake is ultra-oligotrophic with very clear water and low suspended sediment; the Secchi depth is approximately 15 m. According to Reinart et al. (2003) classification, Pavilion Lake is considered a class C lake (clear). Chapter 3. Sites and methods 27 N ↑ 5637 20m 40m Met St. Northing (km) 5636 40m 60m TP 20m 5635 20m 40m 5634 5633 5632 587 588 589 590 591 Easting (km) Figure 3.1: Bathymetry of Pavilion Lake showing north, central, and south basins. (+)TP marks the location of the thermistor chain and (x) marks the location of the meteorological station. Data Collection Water Temperature Data. High-accuracy temperature (±0.005oC) and conductivity (±1 𝜇S cm ) profiles were collected using a Seabird SBE19 plus -1 conductivity–temperature–depth (CTD) profiler. The casts were measured on June 24th 2006, August 11th 2006, and October 17th 2006. Chapter 3. Sites and methods 28 20 18 16 Temp (°C) 14 9.5 m 12 10 15.5 m 15.1 m 8 6 26.5 m 24.1 m 4 1Jun 1Jul 1Aug 1Sep 1Oct 1Nov 2006 Figure 3.2: Line plot of water temperatures measured using temperature recorders located in the central basin of Pavilion Lake, 2006. Temperatures are measured at z = 9.5, 15.5, and 26.5 m for the period from June 1st 2006 to August 1st 2006. The thermistors were recovered to retrieve the data, and then redeployed on August 2nd 2006 at z = 9.5, 15.1, and 24.1 m, recording temperatures until November 20th 2006. Another method to measure the temperature was the deployment of four Onset Hobo Water Temp Pro temperature loggers, with an accuracy of ±0.2oC. They were deployed lying on the bottom of the lake near the deepest part of the central basin (Figure 3.1). The loggers were deployed at 9.5, 15.5, and 26.5 m depths. The data were retrieved from the loggers on August 2nd 2006 and then they were re-deployed at 9.5 m, 15.1 m, and 24.1 m depths (Figure 3.2). Chapter 3. Sites and methods 29 -1 Uw(m s ) 15 a) 10 5 0 100 RH (%) b) 50 0 40 Ta (°C) 30 c) 20 10 0 -2 PAR (W m ) -10 500 400 d) 300 200 100 0 1Jun 1Jul 1Aug 1Sep 1Oct 1Nov Date (2006) Figure 3.3: Meteorological data measured at the meteorological-station located near the north sill in Pavilion Lake for the period June 1st 2006 to November 20th 2006: a) Wind Speed (m s-1), b) Relative Humidity (%), c) Air Temperature (oC), d) PAR (W m-2). Chapter 3. Sites and methods 10 3 2 -2 PAR (W m ) 10 30 10 10 1 0 1Jun 1Jul 1Aug 1Sep 1Oct 1Nov Date (2006) Figure 3.4: Line plot for subsurface PAR irradiance measured in the central basin of Pavilion Lake at z = 5.8 m (black) and 17 m (grey) for the period from June 1st 2006 to November 20th 2006. The empty gap shows the period, August 1st 2006 to August 8th 2006, during which the Licor sensors were pulled out to retrieve the data and redeployed at the same depths. Meteorological Data. A meteorological station was positioned on a dock near the sill separating the north and central basins. The station was a Campbell Scientific CR1000 weather station which recorded wind speed (Uw) and direction at 15 minute intervals, and air temperature (Ta), relative humidity (RH) and the incoming PAR at 30 minute intervals (Figure 3.3). The incoming PAR was recorded using a PAR LITE Kipp & Zonen Quantum Sensor (above-water sensor) which provides spectral response for radiation within λ = 0.4 - 0.7 μm. Chapter 3. Sites and methods 31 Ideal Assumed response 100 Percent relative response LI-COR Quantum sensor 80 60 40 20 0 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 λ (µm) Figure 3.5: Typical spectral response of LI-COR quantum sensor versus wavelength, adapted from (Licor 2006). Optical Data. Optical data was collected in two ways. First, a LI-193 sensor was mounted on the Seabird SBE19 plus CTD (Conductivity- Temperature- Depth) profiler. During the study period, CTD profiles were collected at various locations in June, August, and October of 2006. Second, two LI-193 sensors were deployed in the central basin at 5.8 and 17 m below the water surface and recorded at 15 minutes intervals (Figure 3.4). The sensors were recovered on August 1st, 2006 to retrieve the data, and then they were redeployed on August 8th, 2006 at the same depths (Figure 3.4). The LI193 sensors are quantum sensors that measure photon flux. The output in micromoles s-1 m-2, has been multiplied by 1.3 to convert it to energy flux, in W m-2 as suggested by the manual. Chapter 3. Sites and methods 32 The response of the LI-COR sensors is shown in Figure 3.5. The response is assumed to be constant in the range λ = 0.4 - 0.7 μm, since the variation is within ±10%. The data For Pavilion Lake were collected as part of another project (Lim et al., 2009). Some of the data records extended beyond the period chosen here; however, the period chosen is the only period when all measurements were available simultaneously. 3.3 Methodology 3.3.1 DYRESM To test the importance of the spectral and temporal variation in the attenuation coefficient, K, to the evolution of the thermal structure of lakes, we use the one dimensional model, Dynamic Reservoir Simulation Model (DYRESM). DYRESM simulates the stratification with high temporal resolution (Weinberger and Vetter, 2012). DYRESM assumes that any longitudinal and transverse variations relax quickly by horizontal dispersion. The lake is divided into a set of horizontal layers that are variable in thickness (Lagrangian layer scheme). The thickness of the layers depends on the vertical gradients in the water column, or on volume changes due to inflows and outflows. To run DYRESM, the user is required to input the lake bathymetry, latitude, elevation, surface area, inflows and outflows. Temperature and salinity profiles are provided as the initial condition. The model requires meteorological data consisting of the wind speed, air temperature, rainfall, relative humidity, shortwave radiation and cloud cover. The model also requires a set of input parameters (Table 3.1) that includes a single attenuation coefficient. The core of DYRESM is divided into two main modules: heater and mixer. The mixer module calculates and applies the available turbulent kinetic energy to the surface layer. Turbulent kinetic energy (TKE) results from shear mixing, convective mixing, and Chapter 3. Sites and methods 33 wind stirring. These processes are parameterized using efficiency coefficients; default coefficients are given by DYRESM and no calibration is required. In this study we set these coefficients to their default values given in Table 3.1. The mixer module then stabilizes the water column and mixes layers if there is enough TKE to overcome the potential energy of the stratification. The other component of DYRESM is the heater module, which applies the heat fluxes at the surface and throughout the water column. A flow chart for the heater module in DYRESM is given in Figure 3.6. It deposits the non-penetrative energy exchanges (longwave radiation, and the latent and sensible heat fluxes) into the surface layer and distributes penetrative radiation (shortwave heat fluxes), throughout the water column. In order to do this, the module divides the incoming shortwave into two parts: the first part consists of UV and IR which account for 55% of the solar radiation, while the second part consists of the PAR, which accounts for the remaining 45%. The module delivers the IR and UV into the top bin at the beginning of each time step, which assumes that the attenuation coefficients for these components are high and that they cannot penetrate beyond the first bin. The remaining PAR component is distributed throughout the depth of the water column using Beer’s law. The model does not include vertical variation in optical properties which could result from variations in the concentration of chlorophyll or other dissolved or suspended materials. The model assumes that the water body is optically homogenous in the vertical direction. The model uses the single attenuation coefficient for the entire PAR component and does not account for either the spectral variation or the temporal variation of the attenuation during the simulation period. DYRESM also does not consider any changes in the vertical propagation of light under different solar zenith angles, i.e. it assumes that there are no angular distributions to the subsurface irradiance (Hieronymi, 2011). In this work, DYRESM is modified in two different ways. First, DYRESM is modified to parameterize the spectral attenuation coefficient for the PAR domain by dividing the PAR domain into two, three, and four bands. Each band is given a separate attenuation coefficient input by the user as shown in the flow chart presented in Figure Chapter 3. Sites and methods 34 3.6. Light is then distributed vertically in the water column using Beer’s law for each band separately. Table 3.1: DYRESM default parameters. Parameter Value Bulk aerodynamic momentum transport coefficient 1.3x10-3 Emissivity of water 0.96 Shear production efficiency (ηk) 0.06 Critical wind speed (ms-1) 3 BBL dissipation coefficient 1.4x10-5 Vertical mixing coefficient 200 Potential energy mixing efficiency 0.2 Wind stirring efficiency (ηs) 0.4 Second, DYRESM is modified to accept a time varying attenuation coefficient defined by the user. Note in this case the time varying attenuation coefficient is for the entire PAR band; an attenuation coefficient that is both spectrally and temporally variable is not introduced due to limitations in the data as discussed in Chapter 7. 3.3.2 Calculation of attenuation coefficients Spectral variation of the attenuation coefficient In order to estimate the attenuation coefficients in the different bands of PAR, measured profiles of PAR are used. To get good estimates, the profiles used are measured on cloudless days and during calm periods, to minimize changes due to surface waves and the movement of clouds. In the selected profiles, data below the LI-COR detection limit is not included. Chapter 3. Sites and methods 35 Start Calculate Longwave (LW), Sensible (S), and Latent (L) Heat fluxes Figure 3.6: Flow chart showing the major steps of the heater module in DYRESM, and the new module parameterizing the spectral variability in the attenuation coefficient of PAR, presented in the dashed boxes. Calculate total non-penetrative energy deposited on the surface = (LW + S + L) × surface area × time step Adjust surface layer mass, salinity, density and volume to allow for changes due to evaporation and rainfall If SW=0 Calculate the SW energy penetrating through the surface layer Deposit 55% of SW energy into the surface layer Distribute the rest of the SW heat energy through all the sub-surface layers of the water column using the Beer’s Law Determine the temperature change in each layer End Divide the SW into different bands Using the Beer’s law distribute each band separately with its K starting at the surface layer. Chapter 3. Sites and methods 36 A mathematical method developed for circuit design – Prony’s method (Potts and Tasche, 2010) is used to fit the sum of exponentials representing Beer’s law for each band, to the observed profile, in order to determine the K values for the set of bands. The following is an explanation of the application of Prony’s method to estimate a set of K values from a profile of underwater irradiance, ET(z). In the exact case, there is a separate attenuation coefficient at each wavelength, 𝐸𝑇 (z) = ∫ I(0, 𝜆)exp(−𝐾(𝜆). 𝑧)𝑑𝜆 (3-1) where ET(z) is the irradiance intensity at depth z, K(λ) is the wavelength dependent attenuation coefficient and I(0, λ) is the spectral solar radiation at the surface of the lake. However, since K(λ) is usually unknown, ET(z) can be expressed as a summation for a set of n different wavebands, 𝐸𝑇 (z) = ∑𝑛𝑗=1 I�0, 𝜆𝑗 �exp(−K�𝜆𝑗 �. 𝑧) (3-2) where n is the number of wavelength bands. In this study the number of bands varies from n = 1 (a single attenuation coefficient for PAR) to n = 4. I(0, 𝜆𝑗 ) is the irradiance for band j at the surface. For a given number of bands (n), depth is discretized as an integral number (i = 0 , 1, …, 2n-1) of constant depth intervals, ∆z. Thus we can write z = i∆z and replace the depth dependence by an index i, ET(z) = ET(i∆z) = Ei, so that, 𝐸𝑖 = ∑𝑛𝑗=1 𝐸0 �𝜆𝑗 �(𝑥𝑗 )𝑖 (3-3) 𝑥𝑗 = exp(−K�𝜆𝑗 �. ∆z) (3-4) where Assuming that xj satisfies the following polynomial 𝛽0 + 𝛽1 𝑥 + 𝛽2 𝑥 2 + ⋯ + 𝛽𝑛 𝑥 𝑛 = 0 (3-5) Chapter 3. Sites and methods 37 where βn = 1. By multiplying both sides of equation (3-3) by βi and summing over n, the non-linear system of equations presented by equation (3-3) can be expressed as 𝑛 𝑛 𝑛 � 𝐸𝑖+k β𝑖 = �{ � 𝐸0 �𝜆𝑗 �(𝑥𝑗 )𝑖+k }β𝑖 𝑖=0 𝑖=0 𝑗=1 𝑛 𝑛 𝑗=1 𝑖=0 = � 𝐸0 �𝜆𝑗 �𝑥𝑗 k {� β𝑖 (𝑥𝑗 )𝑖 } = 0 (3-6) where k = 0,1,…,n-1, and equation (3-5) has been used. Thus, the left hand side of equation (3-6) can be written as 𝑛−1 � 𝐸𝑖+𝑘 𝛽𝑖 = −𝐸𝑘+𝑛 (3-7) 𝑖=0 The roots for equation (3-5), can be found after solving equation (3-7) for 𝛽 i, and finally K can be found by rearranging equation (3-4) as follows: 1 𝐾�𝜆𝑗 � = − ∆z ln 𝑥𝑗 (3-8) As stated above, Prony’s method is used in this work to fit the sum of up to four exponential functions to observed PAR profiles. Prony’s method is not intended to be used for such small values of n. For example, for n = 1, Prony’s method simply fits a single exponential to PAR values at z = 0 and ∆z, making the fit extremely dependent on choice of ∆z if, as shown for Pavilion Lake, the data do not follow a single exponential. Chapter 3. Sites and methods 38 Thus, the method was adapted to minimize root mean square (RMS) error of fit at small n by adaptively selecting ∆z to minimize error. To enhance the performance of Prony’s method and reduce its sensitivity to the depth interval (∆z), the measured PAR profiles should be smoothed to eliminate any noise (Coluccio, Eisinberg, and Fedele, 2007). In the present study, no smoothing was required since all ∆z used were larger than length scales associated with noise. The PAR spectrum is divided into bands with equal surface irradiance intensity for n = 2, and 4. For n = 3, in order to expand wavelength resolution at low attenuation, the band corresponding to the high 𝐾 value, in the two band division, is left unchanged, and the band corresponding to the low 𝐾 value, is divided into two bands with equal surface irradiance intensity. Temporal variation in the attenuation coefficient As described in the previous section, temporally variable measurements for PAR irradiance at a fixed depth were available for Pavilion Lake. These data, along with PAR measurements from the lake shore meteorological station are used to calculate an attenuation coefficient which varies in time K(t), using, 1 𝐸 (𝑧,𝑡) � 𝑇 (0,𝑡) 𝑇 𝐾(𝑡) = 𝑧 ln �(1−α)𝐸 (3-9) where ET(0, t) is PAR intensity measured above water at time t, and ET(z, t) is PAR intensity measured at the depth z by the submerged sensor at time t. The water surface albedo, α, is calculated using α = 0.06/cos(γ / 1.2 ) , where γ is the latitude less the solar declination (Robock, 1980). Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 4 PARAMETERIZATION OF SHORTWAVE IRRADIANCE USING SPECTRAL ATTENUATION 4.1 Introduction Solar radiation reaching the surface of lakes and attenuation of this radiation as it travels through the water column, are both strong functions of wavelength, λ. In limnological models this spectral variability is often simplified by distinguishing between ultra-violet (UV; λ = 0.02 - 0.4 µm), photo-synthetically active (PAR; λ = 0.4 - 0.7 µm) and infrared (IR; λ = 0.7 - λ = 3 µm) radiations (Cole and Buchak, 1995; Antenucci and Imerito, 2001; Hodges and Dallimore, 2006). Spectral dependence of attenuation is typically modeled in the following way: UV and IR are assumed to be absorbed at the surface, while the PAR is distributed throughout the water column assuming a constant attenuation coefficient (Spigel and Imberger, 1980; Hornung, 2003). This approach is subject to error, particularly in clear water lakes where the attenuation coefficient can vary spectrally by over an order of magnitude within the PAR domain, and UV radiation can penetrate to considerable depths (Figure 4.1). While solar radiation reaching the surface of a lake generally varies with time, space and wavelength, in this chapter we limit our attention to variation with wavelength (λ). Figure 4.1a presents the spectral incoming solar radiation (Mecherikunnel and Duncan, 1982) within the UV, PAR, and part of the IR domain. The spectrum has a maximum at λ = 0.48 μm. While, the spectrum is relatively flat for λ > 0.48 μm, it decreases rapidly through the UV domain to zero at a wavelength of 0.3 µm, because UV radiation below that wavelength is absorbed and scattered in the atmosphere. Of the incoming energy, 7% is UV, 51% is PAR, and 42% is IR. Subsurface spectral irradiance (I) can be modeled as a function of depth (z) and λ using Beers’ law as follows: 39 Chapter 4. Parameterization of shortwave irradiance using spectral attenuation I(z, λ) = I(0, λ) 𝑒 −𝐾(𝜆)𝑧 40 (4-1) where K(λ) is the spectral attenuation coefficient. As with surface spectral irradiance (I(0, λ)) K generally varies with time, space and wavelength. In this chapter we limit our attention to its variation with wavelength. To highlight the impact of the spectral dependence of subsurface irradiance and heat transfer, we apply equation (4-1) to the following four cases: pure water, Crater Lake, Lake Superior, and San Vicente Reservoir. For each of these four examples, the minimum attenuation coefficient is in the PAR domain (Figure 4.1b), which, following equation (4-1) suggests that the deepest penetrating radiation is in the PAR bands. Overall, K(λ) is smallest for pure water indicating longest attenuation length (K-1). For pure water, the attenuation length is over 100 m at λ = 0.44 μm (violet light). In the IR domain, the attenuation length for pure water is small (< 5 m), while the attenuation length in the near UV remains large (Figure 4.1b), though it quickly decreases with decreasing wavelength. Crater Lake, known for its water clarity, follows pure water with K values that are slightly larger (i.e. smaller attenuation length). Nevertheless, the attenuation length in Crater Lake is still over 60 m at λ = 0.42 μm (similar to pure water). In contrast, for Lake Superior and San Vincente Reservoir, attenuation length is much smaller, with a maximum of 7 m or less at λ = 0.58 μm. In this chapter, we demonstrate the significance of accounting for the spectral variability of subsurface spectral irradiance on the heating of lakes with different turbidities. We start by comparing the subsurface spectral irradiance between Crater Lake and Lake Superior, followed by an illustration of the significance of the spectral variability of surface spectral irradiance and attenuation coefficient on the light intensity for pure water, Crater Lake, San Vicente Lake and Lake Superior. We then make comparisons between theoretical predictions and measurements of hypolimnetic heating due to spectral irradiance for two relatively clear lakes: Pavilion Lake and Crater Lake. Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 4.2 41 Subsurface Spectral Irradiance In order to highlight the effect of the spectral dependence of K we compare predictions of subsurface irradiance made using equation (4-1) using 1) a spectrally varying K with 2) a spectrally independent attenuation coefficient. As examples, we make this comparison for Crater Lake and Lake Superior. A spectrally independent attenuation �𝑙𝑠 ) can be defined as a least squares fit to a profile of light intensity. Where coefficient (𝐾 light intensity is defined as 𝐸𝑇 (z) = ∫ 𝐼(𝑧, 𝜆)𝑑𝜆 (4-2) Note that ET(z) is sometimes referred to as irradiance, but we use intensity to differentiate it from spectral irradiance, 𝐼(𝑧, λ). Note that we do not have measured light �𝑙𝑠 from intensity profiles for the four cases we are studying, therefore we determine 𝐾 ET(z) profiles generated using K(λ) within the PAR domain. For San-Vicente Lake and �𝑙𝑠 is an order of magnitude higher than for pure water and Crater Lake Lake Superior 𝐾 (Table 4.1). The combined effect of the spectral variation of attenuation coefficient and surface spectral irradiance leads to spectral distribution of irradiance varying with depth as shown in Figure 4.2b for Crater Lake. In contrast, when modeling attenuation using a �𝑙𝑠 ), the subsurface spectral distribution of irradiance constant attenuation coefficient (𝐾 reflects the surface distribution at all depths (Figure 4.2c). In this case, more heat is absorbed at shallower depths; 97% of the incoming intensity, ET(0), is absorbed in the top �𝑙𝑠 is used in the calculations, compared to 140 m when K(λ) is used. 78 m when 𝐾 42 2000 (a) UV PAR IR -2 -1 I(0,λ ) (W m µm ) Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 1000 0 0 10 10 0 -1 1 10 10 San Vicente Lake Superior Attenuation length (m) -1 K(λ ) (m ) (b) Crater Lake Pure Water 2 -2 10 10 0.3 0.35 0.4 0.45 0.5 0.55 λ (µm) 0.6 0.65 0.7 0.75 Figure 4.1: a) Spectral distribution of shortwave radiation at the surface, I(0, λ), based on an incoming irradiance intensity of 1000 W m-2. b) Attenuation coefficient, K(λ), for pure water (black) from Bukata (1995); Tyler and Smith (1970), Crater Lake (blue) from Bukata (1995), Lake Superior (green) from Tyler and Smith (1970), and San Vicente Reservoir (red) from Tyler and Smith (1970). For Crater Lake, K(λ) is estimated in the UV domain (blue dashed line) using Equation (4-5) suggested by Morel and Antoine (1994). Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 43 Table 4.1: Optical properties of the four water types studied. � 𝐥𝐬 (m-1) 𝐊 z1% � 𝟏% (m-1) 𝐊 Pure Crater Lake San-Vicente Water Lake Superior 0.04 0.05 0.52 0.74 257 179 9.8 6.6 0.02 0.03 0.47 0.7 As a comparison, subsurface spectral irradiance is also calculated for the relatively turbid Lake Superior (Figure 4.3). The higher values of K(λ) for Lake Superior result in a rapid attenuation of subsurface spectral irradiance; 97% of ET(0) is absorbed in �𝑙𝑠 is used. the top 7.3 m when K(λ) is used in the calculations, and in the top 6.7 m when 𝐾 In both Crater Lake and Lake Superior, the wavelength of maximum subsurface spectral irradiance (λmax) varies with depth. This variation can be directly observed from �𝑙𝑠 is used in the calculations, λmax remains Figures 4.2b and 4.3b. In contrast, when 𝐾 constant with depth. The lack of a continuous function for I(0, λ) and K(λ) precludes writing λmax explicitly as a function of depth from equation (4-1); however, the following implicit relation can be determined by differentiating equation (4-1) with respect to 𝜆 and equating to zero. 𝑧𝑚𝑎𝑥 = 1 𝑑𝐼0 (𝜆) � 𝑑𝜆 � 𝐼0 𝑑𝐾(𝜆) 𝑑𝜆 (4-3) For each wavelength, equation (4-3) gives the depth at which I(z, λ) is extreme. The right hand side of equation (4-3) can be multi-valued at a given depth; however, for I(0, λ) and K(λ) used there is only one solution corresponding to positive depths and this occurs when the derivatives of I(0, λ) and K(λ) have the same sign. 1000 0.6 a) 0.4 500 0.2 0 0 0 b) -1 1500 K(λ ) (m ) -1 44 -2 I(0,λ ) (Wm µm ) Chapter 4. Parameterization of shortwave irradiance using spectral attenuation z (m) 50 100 150 0 c) z (m) 50 100 150 0.3 0.35 0.4 0.45 0.5 0.55 λ (µm) 0.6 0.65 0.7 0.75 Figure 4.2: a) Spectral distribution of irradiance at the surface I(0, λ) (solid line, left scale) and K(λ) for Crater Lake (black dashed line, right scale). Also shown is the best fit, 𝐾�𝑙𝑠 (dash-dotted line, right scale). I(0, λ) was calculated using the typical maximum solar radiation on a sunny day, ET(0) = 1000 W m-2. b) I(z, λ) calculated using K(λ), for Crater Lake (W m-2 µm-1). Line contours show the 200, 400, 600, 800, 1000, and 1200 W m-2 µm-1 levels. The thick solid line indicates λmax. c) I(z, λ) estimated using 𝐾�𝑙𝑠 . Grey area indicates the accurate range of wavelengths corresponding to the maximum irradiance throughout the water column. 1000 1.5 a) 1 500 0.5 0 0 -1 1500 K(λ ) (m ) -1 45 -2 I(0,λ ) (W m µm ) Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 0 b) z (m) 2 4 6 8 0 c) z (m) 2 4 6 8 0.3 0.35 0.4 0.45 0.5 0.55 λ (µm) 0.6 0.65 0.7 0.75 Figure 4.3: a) Spectral distribution of irradiance at the surface I(0, λ) (solid line, left scale) and K(λ) for Lake Superior (black dashed line, right scale). Also shown is the best �𝑙𝑠 (dash-dotted line, right scale). I(0, λ) was calculated using the typical maximum fit, 𝐾 solar radiation on a sunny day, ET(0) = 1000 W m-2. b) I(z, λ) calculated using K(λ), for Lake Superior (W m-2 µm-1). Line contours show the 200, 400, 600, 800, and 1000 W m-2 µm-1 levels. The thick solid line indicates λmax. c) I(z, λ) estimated using 𝐾�𝑙𝑠 . Grey area indicates the accurate range of wavelengths corresponding to the maximum irradiance throughout the water column. For Crater Lake, λmax decreases steadily from 0.48 µm at the surface to nearly 0.42 μm at 100 m depth, then asymptotes towards 0.42 μm (Figure 4.2b). For Lake Chapter 4. Parameterization of shortwave irradiance using spectral attenuation Superior, λmax increases steadily from 0.48 µm at the surface to nearly 0.57 µm at 2 m depth, then asymptotes towards 0.57 µm (Figure 4.3b). In both cases, λmax at the surface corresponds to the wavelength of maximum I(0, λ), i.e. where 𝑑𝐼0,𝜆 𝑑𝜆 asymptotes towards the wavelength of minimum K(λ), i.e. where = 0. At depth, λmax 𝑑𝐾(𝜆) 𝑑𝜆 = 0. In contrast, �𝑙𝑠 is used, equation (4-1) when the spectral variability in K(λ) is ignored and 𝐾 becomes 𝐼(𝑧, λ) = 𝐼(0, λ) 𝑒 −𝐾�𝑙𝑠 𝑧 , which indicates that the spectral distribution of I(z, λ) at any depth is influenced only by the spectral distribution of I(0, λ), and λmax remains 0.48 μm at all depths, as shown in Figures 4.2c and 4.3c. 4.3 Subsurface irradiance intensity For pure water and the three lakes, subsurface profiles of ET(z) were calculated �𝑙𝑠 . In these examples, the profiles are using equation (4-2) (Figure 4.4) with K(λ) and 𝐾 computed by integrating over wavelengths in the PAR domain only, because this is the only domain that the K(λ) spectra are available for the four cases. Profiles are evaluated down to the depth at which intensity decreases to one percent of surface intensity (z1%), commonly referred to as the photic depth (Table 4.1). �𝑙𝑠 decrease exponentially while the profiles The PAR profiles calculated using 𝐾 calculated using K(λ) do not (Figure 4.4). Comparing subsurface intensity calculated �𝑙𝑠 and K(λ), the two profiles diverge at the surface and converge with depth. The using 𝐾 profiles cross at 49, 38, 4, and 4 m depth for pure water, Crater Lake, Lake Superior, and San Vincente Reservoir, respectively. Differences between the profiles are a function of the shape and magnitude of K(𝜆). For each water type, differences decrease as the minimum K(𝜆) for that water increases. In the examples of Figure 4.1, pure water has the lowest minimum K(λ), while San Vincente has the highest minimum K(λ). For comparison we also compute profiles using a commonly used non-spectral �1% , , defined as, attenuation coefficient𝐾 46 Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 0 0 0.1 0 a) b) 47 0 c) 20 1 0 d) 1 50 0.2 2 40 0.3 60 100 0.4 4 3 80 0.5 5 z (m) z / z1% 2 3 150 0.6 0.7 100 6 120 7 4 5 200 0.8 8 140 0.9 1 0 160 100 0 50 -2 ET(z) (W m ) 6 9 250 100 50 -2 ET(z) (W m ) 0 100 0 50 -2 ET(z) (W m ) 100 50 -2 ET(z) (W m ) Figure 4.4: Irradiance intensity profiles based on the PAR irradiance, assuming a surface �𝑙𝑠 (dashed line), and 𝐾 �1% (dotted intensity of 1000 Wm-2, calculated using K(λ) (black), 𝐾 line) for a) Pure Water, b) Crater Lake, c) Lake Superior, and d) San-Vicente Reservoir. The profiles are plotted down to the photic depth per each water body as indicated by the right scale. �1% = ln(100) = 𝐾 z1% 4.6 z1% (4-3) �𝑙𝑠 , 𝐾 �1% is an order of magnitude less for Crater Lake and pure water As with 𝐾 �1% is lower than than for San-Vicente Lake and Lake Superior (Table 4.1). In general 𝐾 �𝑙𝑠 (Table 4.1). The difference between 𝐾 �1% and 𝐾 �𝑙𝑠 is lowest for pure water and 𝐾 increases with increasing water turbidity. The difference between profiles of ET(z) using �1% and K(λ) is almost twice as large as the difference between profiles using 𝐾 �𝑙𝑠 and 𝐾 Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 48 K(λ). Assuming subsurface spectral irradiance modeled using K(λ) is the most accurate �𝑙𝑠 when using a constant attenuation coefficient. representation, we will hence use 𝐾 4.4 Heating rates due to spectral and non-spectral irradiance From conservation of heat applied to a point in the water column, the rate of temperature increase at a given depth, due to the absorption of solar radiation, can be defined as follows: 𝑑𝑇 𝑑𝑡 =− 1 𝑑𝐸𝑇 𝜌𝑐𝑝 𝑑𝑧 (4-4) where, cp is the specific heat capacity of water = 4186 J kg-1 ºC-1, and ρ is the water density (kg m-3). To illustrate the role of spectral variability of attenuation coefficient on heating of a water column, the rate of temperature increase is calculated for the same four water types described above. As before, we compare the effect of calculating ET(z) using �𝑙𝑠 (Figure 4.5). Heating in the epilimnion is affected by K(λ) with that calculated using 𝐾 surface mixing driven by wind and convection, and other heat fluxes such as the net longwave flux (blackbody radiation emitted at the water surface) sensible heat flux (heat lost or gained due to differences between the surface water and overlying atmosphere), and latent heat flux (heat lost due evaporation). Effects of these fluxes and processes diminish below the surface mixed layer. The following compares spectral and nonspectral heating for completeness, but the comparison is only valid above the thermocline in rare circumstances when penetrating solar irradiance is the only influence on local heating. Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 0 0 0 0.1 49 0 1 20 1 50 0.2 2 40 0.3 60 100 0.4 4 3 80 0.5 z(m) z / z1% 2 3 5 150 0.6 0.7 100 6 120 7 4 5 200 0.8 8 140 0.9 6 9 160 250 1 -8 10 10 -6 -4 10 -1 dT/dt(°C s ) -8 10 -6 10 10 -1 dT/dt (°C s ) -4 10 -5 10 -3 -5 10 -1 dT/dt (°C s ) -3 10 -1 dT/dt (°C s ) Figure 4.5: Rate of temperature increase based on the PAR irradiance, assuming a surface �𝑙𝑠 (dashed line) for a) pure intensity of 1000 W m-2, calculated using K(λ) (black), and 𝐾 water, b) Crater Lake, c) Lake Superior, and d) San-Vicente Reservoir. The profiles are plotted down to the photic depth for each water type as indicated by the right scale. For all four water types, at the surface, the temperature gain from subsurface �𝑙𝑠 is underestimated, compared to that modeled using K(λ) and intensity modeled using 𝐾 the residual is distributed with depth. For pure water, this difference in heating at the �𝑙𝑠 underestimates the rate of temperature increase in the top 7 m, surface is 78%; using 𝐾 �𝑙𝑠 underestimates and overestimates it between z = 7 and 95 m. For Crater Lake, using 𝐾 𝑑𝑇 𝑑𝑡 (0) by 72% and overestimates the heat between z = 6 and 71 m. The difference at the surface decreases for the more turbid cases; for Lake Superior, 𝑑𝑇 𝑑𝑡 (0) is underestimated Chapter 4. Parameterization of shortwave irradiance using spectral attenuation by 9% and overestimated between z = 1 and 6.5 m; and for San Vincente 50 𝑑𝑇 𝑑𝑡 (0) is underestimated by 15%, and is overestimated between z = 0.75 and 5.7 m. Note that the values in Figure 4.5 are of very small values, as they are based on assuming a surface intensity of 1000 W m-2 for a period of 1 second, same amount of PAR intensity used earlier. 4.5 Comparison of predicted heating with field data In this section, we show the significance of spectral variation of surface irradiance and attenuation coefficient on the evolution of temperature in Pavilion and Crater Lakes. The results are based on comparing the observed and predicted temperature below the thermocline, away from the effect of non-penetrative surface fluxes and mixing processes. 4.5.1 Pavilion Lake Pavilion Lake (50.9° N, 121.7° W) is located in Marble Canyon, a limestone valley, in south central British Columbia, Canada. The lake is long and narrow with two constrictions dividing it into north, central, and south basins (Figure 3.1). The maximum depth is 55 m, the length is 5.7 km, while the average width is 0.8 km. Pavilion Lake is ultra-oligotrophic with very clear water, low suspended sediment and Secchi depths of 15 m in summer. A meteorological station was located on a dock between the north and central basins, which recorded wind speed (Uw) and direction at 15 minute intervals, and the incoming PAR at 30 minute intervals (Figure 4.6a and 4.6b). The PAR was measured using a PAR LITE Kipp & Zonen Quantum Sensor which provides spectral response to λ = 0.4 - 0.7 μm. Water temperatures at 15.5 and 26.5 m depth were recorded every 30 minutes for the period starting on June 17th 2006 and ending on July 17th 2006, using Hobo Water Temp Pros, with an accuracy of ±0.2 ºC, deployed at the lakebed. At both Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 51 depths the temperature increases throughout the period of interest. The temperature at 15.5 m increases from 6.9 to over 8.0 °C during the study period (Figure 4.6c), while the temperature at 26.5 m increases from 5.0 to 5.3 °C (Figure 4.6d). The significance of the spectral variability in the attenuation coefficient on the rate of temperature increase is tested by comparing the observed temperature to calculated temperatures estimated from the incident light. Water temperatures in Pavilion �𝑙𝑠 . Lake are estimated using both K(λ) and 𝐾 To proceed, we first estimate the non-spectral attenuation coefficient for Pavilion � ls |Pav ) by a linear least squares fit of the natural logarithm of PAR profile, Lake (K measured on August 11th 2006 in Pavilion Lake, with the fit forced through the observed surface irradiance intensity. The profile was measured using a LiCor LI-193 Underwater Spherical Quantum Sensor, designed to measure light within the λ = 0.4 - 0.7 μm range only. The sensor was mounted on a Seabird SBE19 plus conductivity–temperature–depth (CTD) profiler. The profile was measured during a cloudless and calm period of the day. � ls |𝑃𝑎𝑣 for this profile is 0.14 m-1. Note that the sensor measures scalar irradiance and not K downwelling irradiance. However based on Davies-Colley et al. (1993), the attenuation coefficient for upwelling, downwelling, and scalar irradiance are invariable except very close to the water surface. Also, Spigel and Williams (1984) showed that for clear waters, the downwelling irradiance is almost equal to the scalar irradiance. The spectral attenuation coefficient for Pavilion Lake (K(λ)|Pav) is also estimated. Given that Pavilion Lake has a very low dissolved organic carbon (DOC) content (1.3 mg L-1) compared to DOC in other lakes which range between 0.5 to 102 mg L-1 (Lim et al., 2009), and given that the suspended sediment is negligible, K(λ) is estimated using the non-linear formula proposed by Morel and Antoine (1994): 𝐾(𝜆) = 𝐾𝑤 (λ) + 𝜒(λ)[𝐶]𝑒(λ) (4-5) 400 -1 a) 200 0 6 Uw (m s ) -2 PAR (W m ) Chapter 4. Parameterization of shortwave irradiance using spectral attenuation b) 4 2 0 8.4 c) Temp (°C) 8 7.6 7.2 Temp (°C) 6.8 5.5 d) 5.2 4.9 20Jun 25Jun 30Jun 05Jul 10Jul 15Jul Date (2006) Figure 4.6: Measured a) incoming PAR (W m-2), and b) Wind Speed (m s-1). Running daily averages of the temperature at c) 15.5 m, and d) 26.5 m; observed (black), estimated � ls |Pav (dashed magenta), predicted using K(λ)|Pav (solid magenta). using K where 𝐾𝑤 (λ) is the attenuation coefficient spectrum for pure water, 𝜒(λ) and 𝑒(λ) are coefficients (Appendix A), and 𝐶 is the chlorophyll concentration. While the chlorophyll concentration was not available for Pavilion Lake, we estimate it from the photic depth, 30.5 m, calculated from the measured PAR profile, to be [𝐶] = 1.1 mg m-3 using Morel (1988). The resulting K(λ)|Pav is shown in Figure 4.7 along with that for pure water and Crater Lake, shown for reference. 52 Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 10 53 0 -1 K(λ ) (m ) K(λ) Pure K(λ) Crater K(λ) Pav 10 10 -1 -2 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 λ (µm) Figure 4.7: Calculated K(λ)|Pav (magenta), using Equation (4-5) by Morel and Antoine � ls |𝑃𝑎𝑣 (dashed magenta) for PAR domain, K � |𝐶𝑟𝑎𝑡 (dashed blue) for UV and (1994), K (λ) PAR domains, K(λ)|Crat (blue), and 𝐾𝑤 (black) adapted from Bukata (1995). Grey area indicates the upper and lower bounds for K(λ)|Pav for 𝐶 between 2 and 0.5 mg m-3 During the studied period, (from June 17th 2006 to July 17th 2006) , and below the thermocline that extended down to about 12.5 m depth, the only time series of temperature that were available were at 15.5 and 26.5 m depths in the central basin. Since our study period is only one month during the summer, and the lake extends down to 65 m deep, we assume that the only source of heating is the shortwave radiation and that any heating due to groundwater inflow is negligible. Due to the high values of 𝐾(𝜆) in the UV and IR domains (Figure 4.7); 𝐾(𝜆) ≥ 0.5 m-1 which makes the photic depth for UV and IR less than 0.9 m, and since we are only comparing temperatures below 15.5 m, any temperature changes at these depths are attributable to the PAR domain only, and Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 54 therefore the analysis is done for PAR only. The comparisons are also done using running daily averages, to eliminate the effect of the internal waves in the observed temperatures. During the period of interest, the temperature calculated using equation (4-4) at 15.5 m depth using K(λ)|Pav increased from 6.8 ºC to over 8 ºC (Figure 4.6c). The grey shaded area in Figure 4.6c, indicates the upper and lower bounds for the temperature for 𝐶 between 0.5 and 2 mg m-3. The warming observed at 15.5 m is comparable to the solar � ls |𝑃𝑎𝑣 exceeds heating estimated using 𝐾(𝜆). In contrast, temperature estimated using K the observed (Figure 4.6c). Note that, even though the thermistors are lying on the bed of the lake, observed heating rates are similar to the estimated ones. The temperature calculated at 26.5 m depth using the spectral K(λ)|Pav, increased � ls |𝑃𝑎𝑣 is from 5ºC to 5.3ºC (Figure 4.6d). At this depth, the temperature predicted using K lower than observed and lower than the lower bound for the temperature calculated using K(λ)|Pav if 𝐶 was 2 mg m-3 (indicated by grey shaded area). Temperatures calculated using K(λ)|Pav matched the observed temperatures at 15.5 and 26.5 m depths. The grey area in Figure 4.6c indicates that proper parameterization of the spectral dependence of PAR is crucial. For both depths the � ls |𝑃𝑎𝑣 in observed temperatures were within grey area. Note that using the non-spectral K the calculations, resulted in overestimating 𝑑𝑇 underestimating it at z = 26.5 m. 𝑑𝑇 Table 4.2: Observed 𝑑𝑡 , calculated � ls |𝑃𝑎𝑣 the non-spectral K Observed 𝒅𝑻 (ºC day-1) 𝒅𝒕 z = 15.5 m 0.037 𝑑𝑇 𝑑𝑡 at z = 15.5 m (Table 4.2) and using spectral K(λ)|Pav; and calculated Calculated C=0.5 0.035 𝑑𝑡 𝒅𝑻 𝒅𝒕 (ºC day-1) using 𝐊(𝛌)|𝐏𝐚𝐯 C=1 0.037 C=2 0.041 𝑑𝑇 𝑑𝑡 using 𝒅𝑻 Calculated 𝒅𝒕 (ºC day-1) using � 𝐥𝐬 |𝑷𝒂𝒗 𝐊 0.045 Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 55 0 b) a) 5 z (m) 10 15 20 25 30 10 -2 10 -1 dT/dt (°C day-1) 0 10 0 10 20 T (°C) Figure 4.8: a) The average rate of increase in temperature with depth for Pavilion Lake. �𝑙𝑠 (dashed Observed (black triangle), using K(λ)|Pav (solid magenta), and using 𝐾 magenta). b) Temperature profile measured in June in Pavilion Lake showing the depth of the thermocline. The previous comparison was for two depths only. However, when comparing the 𝑑𝑇 𝑑𝑡 � ls |𝑃𝑎𝑣 with that calculated using K(λ)|Pav , it is profile (Figure 4.8) calculated using K found that at depths less than 2 m, 𝑑𝑇 𝑑𝑡 � ls |𝑃𝑎𝑣 is less than that calculated calculated using K using K(λ)|Pav , and greater from 2 m to 24.5 m. The black triangles in Figure 4.8, represent the average observed heating rates, which are better presented by the calculated using K(λ)|Pav than the 𝑑𝑇 𝑑𝑡 𝑑𝑇 𝑑𝑡 profile � ls |𝑃𝑎𝑣 . These results show profile estimated using K the improvement resulting from accounting for the variability in 𝐾(𝜆) in modeling the temperature evolution below the surface layer in Pavilion Lake. Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 4.5.2 Crater Lake Crater Lake located in Oregon, United States (42.9° N, 122.1° W) is one of the clearest lakes in the world, and the deepest lake in the United States. It is formed in the bowl shaped caldera basin of former Mount Mazama. Its surface is almost circular in shape with an area of 53 km2 and a maximum depth of 595 m. The lake filled from the accumulation of rain and snow over a period of approximately 250 years. There are no inflows or outflows to the lake and it is primarily fed by precipitation, which results in an ultra-oligotrophic basin, with Secchi depths down to 41.5 m in June (Larson et al., 2007). There are minimal chemicals or organic materials in the water, and the present contaminants are dominated by particulate materials from the calderal walls and natural algal blooms. Temperature profiles of July and August are given in Larson et al. (2007), and are reproduced in Figure 4.9a. The profiles were originally calculated to assess long term changes in the lake characteristics and they are based on averages of temperature measurements during each month for the years 1988 to 2000. Temperature data are available at 0, 5, 10, 20, 30, 40, 60, 80, 100, and 120 m depth and then every 20 m to a 260 m depth. Here, the comparison is carried out for depths from 10 m to 120 m. Data above 10 m are highly variable due to surface processes. Data below 120 m are also discarded because the changes are very small. In this section, the significance of accounting for the spectral variability in the attenuation coefficient on the thermal evolution of Crater Lake is tested. The increase in the heat content during the studied period, ѱ, is calculated by integrating the total 10 temperature change, between July and August over depth, ѱ= ∫120 𝜌𝑐𝑝 (𝑇𝐴𝑢𝑔 − 𝑇𝐽𝑢𝑙𝑦 )𝑑𝑧 = 371,000 KJ m-2. This increase in heat, is assumed to be solely due to incoming solar radiation. ѱ is then redistributed using both the spectral attenuation coefficient for Crater � ls |𝐶𝑟𝑎𝑡 to predict Lake, K(λ)|Crat , and a non-spectral attenuation coefficient, K temperature changes, which are then compared with observations. 56 Chapter 4. Parameterization of shortwave irradiance using spectral attenuation 57 0 b) Observed Predicted using K(λ) Predicted using Kls a) 20 z (m) 40 60 80 TJuly T 100 120 Aug 4 6 8 10 Temp (°C) 12 10-4 10 -3 -2 -1 10 10 -1 dT/dt (°C day ) 10 0 Figure 4.9: a) Averaged temperature profiles for July (triangle) and August (circle) between 1988 and 2000 for Crater Lake. The points represent the mean with the standard dT deviation for each depth (adapted from Larson et al. (2007)), b) dt observed (circle), � ls |Crat (dashed black). calculated using K(λ)|Crat (solid black), and calculated using K The attenuation coefficient spectrum for Crater Lake, K(λ)|Crat , was measured within the PAR domain during the month of August by Tyler and Smith (1970). UV radiation is expected to penetrate deep in the water column due to the low concentration of dissolved organic carbon (DOC) in Crater Lake, which is only 0.1 mg L-1 (Boss et al., 2007). Therefore, we extend K(λ)|Crat in the UV domain for 𝜆 from 0.3 μm to 0.4 μm Chapter 4. Parameterization of shortwave irradiance using spectral attenuation using equation (4-5) as shown in Figure 4.7. The estimated values of K(λ)|Crat in the UV domain are very close to that of pure water. The least square fit attenuation coefficient, � ls |𝐶𝑟𝑎𝑡 , is 0.051 m-1. Note that this value is slightly higher than the value in Table 4.1. K This value is determined by fitting an ET(z) profile generated using K(λ)|Crat within both the PAR and UV domains. The observed rate of temperature change per day between July and August is simply calculated by dividing the temperature difference TAug - TJuly by the number of days in July. The rate decreases with depth. At 10 m depth, while at 120 m depth, 𝑑𝑇 𝑑𝑡 𝑑𝑇 𝑑𝑡 is around 0.1 ºC day-1, decreased to 0.001 ºC day-1. The decrease is not exponential as indicated by the lack of linearity in the log plot (Figure 4.9b). The profile calculated � ls |𝐶𝑟𝑎𝑡 is. The profiles usingK(λ)|Crat is not exponential while that calculated using K � ls |𝐶𝑟𝑎𝑡 rather than follow the same trend previously observed in section 4.4. Using K K(λ)|Crat results in less change in temperature near the surface (above 8 m) and deep in the hypolimnion (below 45 m ), while, it estimates more heat transfer from 8 m to 45 m. The rates calculated using K(λ)|Crat at z = 60, 80, 100, and 120 m are closer to � ls |𝐶𝑟𝑎𝑡 , which suggests the importance of those observed than those calculated using K wavelength. Nevertheless, there are differences between the rates, which are not explained. A number of factors might be contributing to these differences such as the averaging of the profiles for multiple years. 58 Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 59 5 THE EFFECT OF SPECTRAL ATTENUATION ON TEMPERATURE EVOLUTION FOR DIFFERENT WATER CLARITIES 5.1 Introduction Attenuation controls the vertical distribution of water column heating by solar radiation. As seen in Chapter 4 attenuation depends on wavelength. However, to our knowledge, lake models, including popular models such as the DYnamic REservoir Simulation Model, DYRESM (Imerito, 2007), do not include the wavelength dependence of attenuation, K(λ), within the PAR (Photosynthetically Active Radiation) domain (λ = 0.4 - 0.7 µm). Instead these models distribute PAR vertically in the water column using a constant attenuation coefficient. One consequence of ignoring the wavelength dependence of attenuation in the PAR domain, is that heat is not distributed correctly with depth, and heating rates are often underestimated in the deep hypolimnion (De Stasio et al., 1996, Perroud et al., 2009, and Weinberger and Vetter, 2012). In the previous chapter, we examined the significance of accounting for the spectral variability in K(λ) on the heating in the hypolimnion of two lakes, Pavilion and Crater Lakes, by comparing temperatures predicted using K(λ) to those predicted using a � . Comparisons were carried out for the water non-spectral attenuation coefficient, 𝐾 column below the surface mixed layer where solar radiation is the dominant source of heating. Accounting for spectral variability improved predictions even though temperature below the surface layer can be altered by other processes (Losordo and Piedrahita, 1991; Dejak et al., 1992; Imberger and Hamblin, 1982). In this chapter we explore two questions raised by the results of the previous chapter. Is the spectral effect important when all heat fluxes present in a lake are considered and is it worth considering the spectral effect when K(λ) is rarely known? We account for the heat fluxes using a modified version of the numerical model, DYRESM, Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 60 and we introduce a method for estimating K(λ) from a measured PAR profile. In Section 5.2 we describe DYRESM and our modifications to it. In Section 5.3 we describe how we divide PAR into multiple bands and obtain a constant attenuation coefficient for each band, and in section 5.4 the original and modified model results are compared with the field data from Pavilion Lake for which we have PAR profiles. 5.2 DYRESM DYRESM simulates the vertical distribution of salinity and temperature using mass-conserving layers. Surface heat fluxes and internal mixing processes are parameterized in the model (Hamblin, and Imberger, 1984; Antenucci and Horn, 2002). Air temperature, cloud cover, rainfall, vapor pressure, wind speed and solar radiation are required inputs for the simulation period. The model also requires the total incoming shortwave radiation, which includes ultraviolet radiation (UV; 0.02 – 0.4 μm), PAR, and infrared (IR; 0.7 - 3 μm). DYRESM is divided into two main modules; heater and mixer. The heater module deposits non-penetrative heat fluxes (longwave, latent and sensible) into the top numerical layer, and distributes the penetrative heat flux, namely shortwave radiation, throughout the water column. It divides the input shortwave radiation into two parts. The first part, 55%, accounts for IR and UV, and is placed into the top layer at the beginning of each simulation time step, regardless of the thickness of the top layer and the extinction coefficient. The second part, 45%, accounting for PAR, is distributed through the water column using Beer's law with a single attenuation coefficient. To provide a more accurate parameterization of the spectral attenuation coefficient K(λ), we modified the heater module to divide the PAR domain into two, three and four different bands, with a separate non spectral attenuation coefficient ���𝑛�). assigned for each band (𝐾 Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 61 The other main module of DYRESM is the mixer. Turbulent kinetic energy is introduced through shear mixing, convective mixing, and wind stirring. These processes are parameterized using efficiency coefficients, set to the default efficiency coefficients used in DYRESM (Table 3.1). 5.3 Calculating ���� 𝑲𝒏 In lake modeling, users typically use a non-spectral attenuation coefficient as representative of the spectral K(λ) within the PAR domain. They rarely have the information necessary to determine K(λ), and therefore they estimate the non-spectral coefficient using conventional methods such as fitting irradiance intensity profiles, or using Secchi depth measurements (Adiyanti and Imberger, 2007). The direct method to determine K(λ) is to use a spectral radiometer (Austin and Petzold, 1986; Reinart and Herlevi, 1999), which is not often available. Alternatively K(λ) can be estimated using a Monte Carlo simulation for photon propagation in water with dissolved organic substances; however, this requires detailed knowledge of the water composition (Kirk, 1984). For both Pavilion and Tailings Lakes we estimate the attenuation coefficients for ���𝑛�, by fitting to a measured PAR profile (ET(z)). To estimate n different bands of PAR, 𝐾 ���� 𝐾𝑛 we chose Prony’s method (Potts and Tasche, 2010). This method was developed for circuit design to extract exponential signals from time series data (Osborne and Smyth, 1995). We apply the same method but to spatial data, ET(z) rather than temporal data, as previously discussed in Chapter 3. In Pavilion Lake, profiles of PAR were measured using a LI-COR LI-193SA spherical quantum sensor mounted on a Seabird SBE19 plus CTD (ConductivityTemperature-Depth) profiler. A profile was measured on 11th August 2006, a calm and cloudless day during the summer stratification period. The CTD was lowered from the sunny side to avoid the shadow of the boat. The profile was measured on a calm day to ensure that it was not affected by refraction from surface waves. In Tailings Lake a PAR Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 62 profile was measured every 0.25 m only down to 1.75 m using a LI-COR LI-193 sensor. The profile was measured on June 18th 2006, also on a sunny day with low wind. ���1 ) for both Pavilion and Tailings Lakes Non-spectral attenuation coefficients (𝐾 are given in Table 5.1. ��� 𝐾1 is smaller in Pavilion than in Tailings Lake indicating spectral irradiance can penetrate much deeper in Pavilion Lake. Also shown, are the attenuation coefficients when the incoming PAR spectrum is divided into two bands with equal surface irradiance intensity; note how one value is much smaller than the other in both cases. Next, PAR is divided into three bands by evenly splitting the surface irradiance intensity of the band with the lowest value and leaving the other band unchanged. Finally, the spectrum is divided into four bands of equal surface irradiance intensity. Note that when dividing the PAR for Pavilion Lake into two or more bands, the ���𝑛� in Table 5.1 are higher than those estimated in Figure 4.7. There are upper values of 𝐾 several possible reasons for this discrepancy, including the following. The fit for the high values of �𝐾��� 𝑛 may be especially sensitive to details of the PAR profile such as the zero depth and the depth variation in irradiance near the surface. On the other hand, the estimates of K(λ) in Figure 4.7 assume the only impurity is Chlorophyll, while the long wavelengths of PAR are especially sensitive to impurities, e.g., particulates (Bukata, 1995). Note that even though there is uncertainty in the upper values of the attenuation coefficient, this uncertainly will not affect the hypolimnion (i.e., the results in Chapter 4) where the distribution of light and heat is controlled by the smaller values of K(λ). Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 63 ���𝑛� values estimated using Prony’s method Table 5.1: 𝐾 Number of bands Pavilion Lake ���� 𝑲𝒏 (m-1) Tailings Lake ���� 𝑲𝒏 (m-1) 1 0.13 3.0 2 0.1 2.2 3 0.09 0.13 4 0.09 0.13 2.7 2.2 1.7 4.6 3.6 2.6 2.9 2.6 2.9 3.6 3.5 3.6 0 b) a) 5 z (m) 10 15 20 25 30 0 10 10 1 10 ET(z) (W m-2) 2 3 10 10 -8 10 -6 10 -4 10 -2 dT/dt (° C s-1) Figure 5.1: Observed (11th August 2006) and estimated a) Irradiance profiles for the PAR domain, and b) rate of temperature gain based on the incoming PAR irradiance only in Pavilion Lake. Observed distribution (black), and calculated using a one (red), two (green), three (cyan), and four (blue) bands. Note that there is no significant difference between the three and four division PAR estimates, and that the blue line lies on top of the cyan line. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 64 0 0.2 b) a) 0.4 0.6 z (m) 0.8 1 1.2 1.4 1.6 1.8 2 -1 10 10 0 10 1 ET (z) (W m-2) 10 2 10 -7 10 -6 10 -5 10 -4 dT/dt (° C s-1) Figure 5.2: Observed and estimated a) Irradiance profiles for the PAR domain, and b) rate of temperature gain based on the incoming PAR irradiance only in Tailings Lake. Observed distribution (black), and calculated using a one (red), two (green), three (cyan), and four (blue) bands. Note that there is no significant difference between the two, three, and four division ET or green lines. 𝑑𝑇 𝑑𝑡 estimates, and that the blue line lies on top of the cyan and ���� We used the estimated 𝐾 𝑛 values for the one, two, three, and four band divisions in Table 5.1 to calculate PAR profiles, for Pavilion Lake (Figure 5.1a) and Tailings Lake (Figure 5.2a). For Pavilion Lake using one band overestimated the calculated PAR profile over almost the whole depth (Figure 5.1a). This is an artifact of Prony’s method for n = 1 with ∆z = 30 m, which fits a straight line between the values at z = 0 and ∆z; however, in �𝑙𝑠 ) = 0.14 m-1 this case, ��� 𝐾1 = 0.13 m-1 is very similar to the linear least square fit (𝐾 calculated in section 4.5.1. This consistency between methods comes from the adaptation of Prony’s method to select ∆z such that the RMS fit error is minimized. Dividing the Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 65 spectrum into two bands with two distinct ���� 𝐾𝑛 values, improved the calculated light profile for Pavilion Lake significantly near the surface. Further division into three bands resulted in a slight improvement for Pavilion Lake, while division into four bands was almost identical to that of three bands (Figure 5.1a). This rapid improvement with increasing bands is driven mainly by ∆z being reduced to 2 m for n = 2 - 4. The adapted Prony’s method reduces RMS error by choosing a small ∆z to fit the large near-surface curvature in the profile. This preferential fitting of the near surface comes from the nonlinearity of the profile with much larger values of PAR being absorbed near the surface. A consequence of this is that increasing the number of bands from 1 to 2 actually worsens the fit at depth, while increasing the number of bands to 3 and 4 provides only marginal improvement. Local heating, however, is dependent on the vertical gradient of irradiance intensity, and this deterioration of the fit of irradiance intensity at depth does not significantly affect local heating (Figure 5.1b). For Tailings Lake, the observed light profile is almost exponential over the entire depth of the profile (Figure 5.2a), and the calculated profile using one attenuation coefficient does a reasonable job of representing the observed data at all depths even for only one band. Increasing the number of bands makes little difference. To identify at which depths heating errors might occur, the rate of heat gain was estimated from both the observed and calculated light profiles for Pavilion Lake (Figure 5.1b) and Tailings Lake (Figure 5.2b). For Pavilion Lake, using one band to estimate the heat gain profile, using equation (4-4), overestimates the heat transfer between 1.5 and 25 m depth (Figure 5.1b). The estimated 𝑑𝑇 𝑑𝑡 profile is significantly improved when PAR is divided into two bands, due to the high �𝐾��𝑛� value (2.2 μm) which captured the measured ���𝑛� is 2 m, and it PAR distribution at the surface. Again, note that the photic depth for this 𝐾 has no influence on the heating in the hypolimnion shown in chapter 4. Further division into three bands results in slight improvement, and is almost the same as that for four bands. In contrast, the observed 𝑑𝑇 𝑑𝑡 profile for Tailings Lake (Figure 5.2b), is almost Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 66 exponential with depth, and dividing the PAR into bands did not result in significant improvement in the predicted heat flux. 5.4 Model results In this section we will show the temperatures predicted by DYRESM using the different parameterizations of the attenuation, and compare to field data for Pavilion Lake. The first simulation is run using the original DYRESM, which treats PAR as monochromatic light; and the second, third, and fourth simulations are run using DYRESM modified to divide PAR into two, three, and four bands as given in Table 5.1. For Pavilion Lake, the model runs from June 24th to July 24th, 2006. Measured temperature and salinity profiles were used to initialize the model. A meteorological station located on the dock between the north and central basins (Figure 3.1) provided incoming shortwave radiation, wind speed (Ua), air temperature (Ta), and vapor pressure recorded at 15 minute intervals (Figure 5.3). Cloud cover was calculated using the no sky PAR and the measured incoming PAR. During sunny periods the air temperature peaks at 20 to 30 °C; the wind speeds shows a diurnal signal, with winds generally increasing in the afternoon (Figure 5.3). The model results are summarized in Figure 5.4. The four model runs predict epilimnia with different temperatures and depths. Using one band for PAR (Figure 5.4a) predicts an epilimnion which is cooler and deeper than the surface layer predicted using the multi band forms (Figure 5.4b, c, and d). As seen in Figure 5.1b, the major difference in heat transfer between the one and multiple bands occurs in the top 5 m; for multiple bands there is more heat content near the surface, resulting in increased stratification. For one band there is less stratification, wind mixed the epilimnion deeper, and the predicted thermocline has a steeper gradient. Note also that the 8oC contour, at 11 to 13 m depth, lies 1 m deeper for one band (Figure 5.4a) than for multiple bands (Figure 5.4b, c, and d). This extra heat at depth is a consequence of overestimating the heat flux at depth as shown in Figure 5.1b. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 67 600 -2 PAR (W m ) a) 400 200 0 35 b) Ta (°C) 30 25 20 15 10 c) -1 Uw (m s ) 6 4 2 0 24Jun 29Jun 04Jul 09Jul 14Jul 19Jul 24Jul Date (2006) Figure 5.3: a) PAR, b) wind speed, and c) air temperature observed at Pavilion Lake for the studied period in 2006. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 68 Figure 5.4: The evolution of the thermal stratification for Pavilion Lake simulated using DYRESM with a) one, b) two, c) three, and d) four band exponential forms, from the top down, representing the subsurface PAR distributions. Contours plotted are for 6, 8, 10, 12, 14, 16, 18, 20, 22, and 24 oC. To look at the effect on the deep water temperature in greater detail, line plots of observed and simulated temperature are given in Figure 5.5 for z = 9.5, 15.5 and 26.5 m. The average observed rate of heating at 9.5 m depth is 0.09 oC day-1, while that calculated Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 69 using the one band is much higher, 0.2 oC day-1. Dividing PAR into two, three, and four bands results in an average heating rate of 0.08 oC day-1 which is very similar to the observed rate. At z = 15.5 m, the observed heating rate is 0.04 oC day-1. While using one band over-estimated this rate significantly, the rates predicted using the multiple band divisions were again comparable to those observed. The rate predicted using the four band division is the closest to the observed rates, 0.04 oC day-1. Nevertheless, at z = 26.5 m all forms were accurate in predicting the observed heating rate. The high attenuation coefficient representing the one band division inaccurately concentrated most of the heat at the surface, while the low attenuation coefficient representing some of the bands in the multi-band divisions transported more heat deeper. At 26.5 m depth, most of the heat was already depleted using all band divisions. The high heating rates simulated using the one band division at 9.5 m and 15.5 m, support the conclusions from the previous chapter, and are consistent with Figure 5.1. More heat is expected to reach the water column from 2 m to 25 m depths if PAR is distributed using one band rather than multiple bands, which match the theoretical calculations presented in Figure 4.8. At 26.5 m depth, the heating rate estimated using �𝑙𝑠 is 0.33 oC/30 days, which again are K(λ), is 0.36 oC/30 days and that estimated using 𝐾 very close to the heating rates simulated using DYRESM. We now turn to examine the differences in the surface temperature, which resulted from the high ��� 𝐾𝑛� values estimated for some bands, in greater detail. As shown in Figure 5.4, the major difference between the different simulations is the distribution and the magnitude of heat content. This change in heat content is explicitly dependent on the surface temperature. The model outputs for Pavilion Lake show that, when PAR is dealt with as one band, a cooler surface temperature (Tw) is predicted. The low surface temperature predicted using one band, results in less heat radiated from the water surface. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 70 16 14 z = 9.5 m Temp(°C) 12 10 8 z = 15.5 m 6 z = 26.5 m 4 24Jun 29Jun 04Jul 09Jul 14Jul 19Jul 24Jul Date Figure 5.5: Observed (black) and simulated temperatures in Pavilion Lake at 9.5 m, 15.5 m, and 26.5 m depths, one (red), two (green), three (cyan), and four (blue) band exponential forms. This, along with reduced stratification and a deeper thermocline, leaves the lake with higher heat content in the surface layer and the upper hypolimnion, than that simulated using PAR divided into multiple bands. The non-penetrative latent (𝑄𝑙ℎ ), sensible (Qsh), and longwave (𝑄𝑙𝑤 ) heat fluxes are explicitly dependent on Tw. The effect of surface temperature on these fluxes is described below. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 71 28 a) Tw (°C) 26 24 22 20 δ Tw (°C) 4 2 0 b) -2 24Jun 29Jun 04Jul 09Jul 14Jul 19Jul 24Jul Date (2006) Figure 5.6: a) Modeled surface temperatures using DYRESM with, one (red), two (cyan), three (green), and four band (blue) exponential forms. b) Differences between Tw modeled using four band exponential form and one band exponential form. The surface temperature from the model runs is shown in Figure 5.6. The difference between the simulated surface temperatures decrease on sunny warm days, for example, during the periods starting on June 28th to July 1st, and July 18th to July 21st, the difference drops to 0.5 oC during nighttime. On the other hand, the maximum difference, 3.2 oC, is calculated on the calmest night with the lowest air temperature (Ta) recorded Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 72 during the simulated period is, the night of June 26th. The average Ta and wind speed (Ua) recorded for that night are, 10 oC and 0.3 m s-1. Latent heat flux ranged between 0 and -386 W m-2 (Figure 5.7a). During the simulated period, no water condenses on the surface; water only evaporates from the surface losing latent heat. Therefore, the latent heat flux is always negative. 𝑄𝑙ℎ is a function of the wind speed and Tw (Appendix D), the volume of evaporated water increases with the increase in wind speed and Tw. The sensible heat flux outgoing from the water surface is also proportional to Tw (Appendix D) and changes from one band division to another. During nighttime, the water retains heat more than air does, because the specific heat capacity of water (~ 4.2 Jg-1 K-1) is almost 4 times that of air (~ 1.0 J g-1 K-1). This in turn results in the water surface cooling down, losing sensible heat to the atmosphere. On sunny days, the rate of heating of the water surface is smaller than that of air, because of the difference in heat capacities. Therefore, water warms up as it gains heat from the atmosphere. On the other hand, cloud coverage reduces the available solar radiation, which in turn decreases the heating rate of the water surface, but causes even higher decrease in the heating rate of the atmosphere. The water loses heat to the atmosphere and cools down (Figure 5.7b). Qsh ranges between 50 W m-2 and -170 W m-2. Another non-penetrative flux is the net longwave radiative flux, the sum of emitted and absorbed longwave fluxes is shown in Figure 5.7c. The outward flux is the blackbody radiation emitted from the surface of the water. The amount of radiation depends primarily on the temperature of this body. The absorbed longwave flux is the blackbody radiation emitted from cloud, known as back radiation. This flux depends on the cloud coverage, Ta, and cloud temperature (Appendix D). Only during sunny periods does the atmosphere radiate more heat than water, resulting in net incoming longwave into the water surface. Maximum net longwave flux heating the water is 34 W m-2. During night time, the net longwave flux is leaving the water surface. Maximum net longwave going out from the water surface is 160 W m-2 (Figure 5.7c). Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 73 Differences in the simulated Tw are inter-linked with differences between the simulated heat fluxes (Figure 5.8). The low Tw predicted using the one exponential form, on cloudy days with cool 𝑇𝑎 , is associated with estimating less heat leaving the water surface than the outgoing heat estimated when using the four exponential form. However, the high Tw estimated using the one exponential form, for the period starting July 18th to July 21st, results in estimating more heat leaving the water surface than that simulated using the four exponential forms. Differences between the estimated latent heat fluxes contributes the most to the total difference between the heat fluxes estimated using one and four exponential forms (Figure 5.8). The low magnitude of heat flux leaving the surface, predicted when using the one band division, left the water column with more heat content than the observed and simulated using the multi band divisions. This difference in heat contents is concentrated near the surface, increasing the predicted heating rates at 9.5 and 15.5 m depths. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 74 -2 Qlh (W m ) 0 -100 -200 -300 a) -2 Qsh(W m ) -400 50 0 -50 -100 b) -2 Qlw(W m ) -150 100 50 c) 0 -50 -100 -2 Q(W m ) -150 500 0 d) -500 -1000 24Jun 29Jun 04Jul 09Jul 14Jul 19Jul 24Jul Date (2006) Figure 5.7: a) latent, b) sensible, c) longwave, and d) total non-penetrative radiation (W m-2), using one (red), two (blue), three (green), and four band (black) exponential forms. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 75 -2 δQlh (W m ) 80 60 40 20 0 a) -2 δQsh(W m ) -20 30 20 10 0 b) -2 δQlw(W m ) -10 20 10 0 -2 δQ(W m ) -10 150 c) 100 50 0 d) -50 24Jun 29Jun 04Jul 09Jul 14Jul 19Jul 24Jul Date (2006) Figure 5.8: Differences between a) latent, b) sensible, c) longwave, and d) total heat fluxes (W m-2), calculated using one and four band exponential forms. 5.4 Conclusions • The modified heater version in DYRESM improved the simulated heating rates in the epilimnion and hypolimnion for Pavilion Lake. Chapter 5. The effect of spectral attenuation on temperature evolution for different water clarities 76 • Heating rates of the hypolimnion simulated using DYRESM were often under or overestimated in previous studies. However, from the runs made for Pavilion, using the modified DYRESM, the heating rates are accurately predicted, which indicates that the attenuation of solar flux may be the dominant mechanism of heat transfer in the hypolimnion. • For the relatively clear Pavilion Lake, most of the improvement in the simulated temperature contours was obtained from dividing the PAR into two band divisions. Nevertheless, further division of the band with the high attenuation, (transition from three to four bands) slightly improved the calculated PAR profile at the surface which gave better approximations of the surface temperature and consequently the total heat content and heating rates simulated by the model. • The attenuation coefficient values representing different divisions of the K(λ) spectrum can be estimated by applying Prony’s method on a single PAR profile. However, ∆z must be chosen carefully for n = 1, and the profile has to be measured on a sunny day under low winds and all the noise has to be filtered. • More band divisions are required to accurately calculate PAR profiles for Pavilion Lake than for the turbid Tailings Lake. Chapter 6. Temporal variation of attenuation 6 TEMPORAL VARIATION OF ATTENUATION 6.1 Introduction In the previous two chapters we studied the wavelength dependence of the attenuation coefficient and its significance on the stratification of several lakes with different clarities. In these examples, the attenuation coefficients were assumed to be temporally constant, during the relatively short simulation periods. However, the attenuation of light in water depends on water properties and other external factors such as the solar altitude, cloud cover, and wave height. These factors can vary in time and therefore result in temporal variation in irradiance attenuation (Lee et al., 2005). A few studies have looked at the effect of the temporal variation in the attenuation coefficient as a result of changes in solar altitude (Stramska and Frye, 1997; Zheng et al., 2002; Lee et al., 2005), while others have been concerned with the temporal variability in the attenuation coefficient due to changes in chlorophyll concentrations (Paavel et al., 2008). Tanentzap et al. (2008) have shown how variations in the attenuation coefficient play a major role in the thermal evolution of Clearwater Lake. They, along with others (Han et al., 2000; Tanentzap et al., 2007; Perroud et al., 2009), have suggested that the constant coefficient in the one dimensional hydrodynamic model, DYnamic Reservoir Simulation Model (DYRESM), should be replaced by a time varying attenuation coefficient. In Section 6.2, we look at the dependence of the attenuation coefficient, K(λ, θ), on solar zenith angle (θ), and how this affects heat transfer in pure water. In Section 6.3 we calculate K(λ, θ) from data from Pavilion Lake, and illustrate how this K(λ, θ) varies diurnally due to solar altitude and topographic shading. In section 6.4, we modify DYRESM to accept a time varying K(λ, θ), and give the results for simulations of Pavilion Lake using both the original and modified version of DYRESM. We then compare the model results with observed field data to explore both diurnal and seasonal variation in K(λ, θ). 77 Chapter 6. Temporal variation of attenuation 6.2 Dependence of the attenuation coefficient on the solar zenith angle As solar radiation propagates through the atmosphere it is both absorbed, and scattered by air molecules. These two processes segregate the irradiance intensity ET into two components, direct, Edirect(θ), and diffuse, Ediffuse(θ), which are both functions of the solar zenith angle, θ (Figure 6.1a). Direct radiation travels in a straight line from the sun to the earth, while diffuse radiation is scattered within the atmosphere. The two components Edirect(θ) and Ediffuse(θ) vary with time of day and can be calculated using equations (2-4) and (2-5). The relative components Edirect/ET and Ediffuse/ET are shown in Figure 6.1a. When θ = 0 (sun directly overhead), almost 80% of total irradiance intensity is direct. As the solar zenith angle increases, the path length of the solar radiation through the atmosphere increases, a higher percentage of incoming solar radiation is absorbed and scattered, and the percentage of direct radiation decreases. Not only is it important to decompose the light arriving at the water surface into direct and diffuse light, but it is also important to consider separately the attenuation of the direct and diffuse components of light propagating in the water (Bukata et al., 1995). As described in previous chapters, the amount of irradiance in the water at a given wavelength depends on both the amount of light at the surface at that wavelength, and the attenuation coefficient at that wavelength. First, consider the spectral distribution of the direct and diffuse components at the surface; they are often considered to be the same as that of the total solar radiation ET at the water surface, and we make use of this here 1. To represent the variation of irradiance in the water, two attenuation coefficients are used Kdirect(λ, θ) and Kdiffuse(λ). While, both of the components depend on the wavelength of irradiance (λ), only Kdirect(λ, θ) depends on the solar zenith angle, θ. It is possible to estimate the spectral attenuation of both the direct light, Kdirect(λ, θ), and the diffuse light Kdiffuse(λ) for pure water based on the absorption and scattering coefficients (Appendix B) and using equations (2-9) and (2-10) (Kirk, 1984). In this chapter we will again only deal with photosynthetically active radiation (PAR; λ = 0.4 - 0.7µm), which we set to be the total, ET. 1 78 Chapter 6. Temporal variation of attenuation 79 1 0.8 a) Edirect/ET E/ET 0.6 0.4 Ediffuse/ET 0.2 0 0.2 b) Kdirect -1 K (m ) 0.15 0.1 K K 0.05 diffuse 0 0 30 θ (°) 60 90 Figure 6.1: a) Fraction of diffuse and direct radiation in the total incoming radiation, at the water surface, Ediffuse/ET and Edirect/ET ; b) Kdirect(λ, θ), Kdiffuse(λ), and K(λ, θ) for pure water under different sun angles for λ = 0.48 μm. Figure 6.1b shows the effect of solar zenith angle on Kdirect(λ, θ) at the wavelength of peak solar radiation, λ = 0.48 μm. At low solar angle, with the sun overhead, Kdirect(λ, θ) and Kdiffuse(λ) are small and very similar; as the solar angle (θ) increases, Kdirect(λ, θ) begins to increase, while Kdiffuse(λ) remains the same. Chapter 6. Temporal variation of attenuation 10 -1 K(m ) 10 10 10 10 80 1 0 -1 -2 -3 0.4 0.45 0.5 0.55 0.6 0.65 λ (µm) Figure 6.2: Attenuation coefficients for pure water at noon and at zenith angle = 77o; Kdirect(λ, 77o) (red), Kdirect(λ, 0o) (black), and Kdiffuse(λ) (blue). The wavelength dependence of Kdirect(λ, θ) and Kdiffuse(λ) for pure water is illustrated in Figure 6.2. At θ = 0, the path length is shortest and yields the lowest Kdirect(λ, θ). At small angles, Kdirect(λ, θ) is similar to Kdiffuse(λ) for almost all wavelengths λ = 0.44 μm to 0.7 μm. As the solar angle increases and the slant path length to achieve a given depth increases, Kdirect(λ, θ) increases at all wavelengths. The next step is to combine Kdirect(λ, θ) and Kdiffuse(λ) into a total attenuation coefficient. By assuming that the spectral variation in Edirect(θ) and Ediffuse(θ) are the same Chapter 6. Temporal variation of attenuation 81 as that of the total solar radiation ET (Bukata et al., 1995), the total attenuation coefficient at any wavelength λ, K(λ, θ), can be estimated as a linear function of Kdirect(λ, θ) and Kdiffuse(λ), 𝐾(λ, θ) = 𝐹(𝜃)𝐾𝑑𝑖𝑓𝑓𝑢𝑠𝑒 (λ) + (1 − 𝐹(𝜃))𝐾𝑑𝑖𝑟𝑒𝑐𝑡 (λ, 𝜃) (6-1) where, the fraction of diffuse light, F(θ) = Ediffuse(θ)/ET. The fraction of diffuse light is a minimum, F(θ) = 0.24, when the sun is vertical (θ = 0) and increases toward F(θ) = 1 with increasing zenith angle (Figure 6.1a). The effect of solar zenith angle on the total attenuation coefficient, K(λ, θ) is shown in Figure 6.1b for λ = 0.48 μm in pure water. As θ increases, Kdirect(λ, θ) increases nonlinearly while Kdiffuse(λ) remains constant, and K(0.48 μm, θ) increases from 0.022 m-1 at θ = 0o, to a maximum of 0.041 m-1 at θ = 77o. As θ increases further, almost all the light at the water surface becomes diffuse, and K(0.48 μm, θ) decreases again (Figure 6.1b). 6.3 Field observations 6.3.1 Diurnal variation in the attenuation coefficient In this section we compare the theoretical dependence of the attenuation in the � (PAR, θ), on the solar zenith angle θ, to observations from Pavilion entire PAR band, 𝐾 � (PAR, θ) Lake. The study period is from June 25th 2006 to July 2nd 2006. We calculate 𝐾 using measurements of irradiance intensity recorded using an above-water sensor at the weather station located in the north basin of Pavilion Lake (Figure 3.1), and irradiance intensity measured using submerged sensors at 5.8 m, ET(5.8), and 17 m ET(17) depths in the mid basin. The above-water sensor is a PAR LITE Kipp & Zonen Quantum Sensor which provides spectral response within λ = 0.4 - 0.7 μm. The submerged sensors are LICOR LI-193SA Underwater Spherical Quantum Sensors and recorded at 15 minute intervals. The irradiance intensity both at the water surface and at 5.8 and 17 m depths are shown in Figure 6.3a for a week of sunny days. Chapter 6. Temporal variation of attenuation 82 The attenuation coefficient is calculated at depth, z, using, � = −1 ln 𝐸𝑇 (𝑧). 𝐾 𝑧 𝐸 (0) 𝑇 (6-2) � (PAR, θ) used at both depths is the Note that the reference irradiance intensity for 𝐾 surface irradiance intensity, ET(0), incoming irradiance after accounting for surface � (PAR, θ) is an apparent attenuation as explained in Appendix C. albedo, and as a result 𝐾 The variation in the apparent attenuation coefficient with the solar zenith angle is shown in Figure 6.3b, and resembles the trend of that calculated for pure water (Figure 6.1b), � (PAR, θ) observed at noon, and as the sun gets lower in the sky, with a minimum 𝐾 � (PAR, θ) increases. 𝐾 The attenuation coefficient calculated at 5.8 m varies significantly with the solar � (PAR, θ) zenith angle for Pavilion Lake (Figure 6.3b). This significant variation in 𝐾 might be a result of the relatively clear water, giving a low scattering to absorption ratio, which indicates that cos(θ) dominates the path length of light photons (equation (2-9)) � (PAR, θ), as suggested by Kirk (1984). 𝐾 � (PAR, θ) calculated at 17 m and consequently, 𝐾 depth, shows less dependence on θ. At this depth, the percentage of short wavelengths (< 0.51 μm) comprising ET increases, because the longer wavelengths correspond to higher attenuation coefficients and are rapidly absorbed at shallower depths. These wavelengths are highly sensitive to scattering (Appendix B) and therefore scattering � (PAR, θ) (Zheng et al., 2002). Values of the dominates over cos(θ) when determining 𝐾 attenuation coefficient calculated at 5.8 m depth are on average, 0.12 m-1 higher than � (PAR, θ) at 17 m depth. 𝐾 � (PAR, θ) are observed. Clouds and the spatial Irregularities in the calculated 𝐾 separation between the above-water and the submerged sensors are possible reasons behind these irregularities. Moreover, these irregularities intensify during cloudy days. Nevertheless, daily consistent spikes and dips are observed. The reason behind these spikes and dips is that the submerged sensor and the meteorological station are at Chapter 6. Temporal variation of attenuation 83 500 a) z=0 -2 ET (W m ) 400 300 200 z = 5.8 m 100 z = 17 m 0 0.4 0.35 b) 0.25 -1 K (m ) 0.3 0.2 0.15 0.1 0.05 0 26Jun 28Jun 30Jun 02Jul Date (2006) Figure 6.3: a) Irradiance intensity at the surface, ET(0) (black), at 5.8 m depth ET(5.8) � (PAR, θ) at z = 5.8 m (blue), z = (blue), and at 17 m depth ET(17) (red) b) Calculated 𝐾 � (PAR, θ) (dashed black). Black arrows point to the 17 m (red) and a depth averaged 𝐾 th diurnal spike and dip on June 28 2006. different locations, receiving different magnitudes and compositions of light. Note that � (PAR, θ) values are not representative of during these spikes, the calculated 𝐾 � (PAR, θ) calculated at z = 5.8 m, a diurnal upward Kdirect(PAR, θ) and Kdiffuse(PAR). For 𝐾 spike ranging from 0.35 - 0.4 m-1 is observed at the beginning of each day. Another diurnal downward dip ~ 0.2 m-1 is observed daily closer to sunset (Figure 6.3b). Chapter 6. Temporal variation of attenuation -1 K (m ) a) 84 Kdirect 0.2 Kdiffuse 0.1 0.2 b) -1 K (m ) 0.15 0.1 0.05 0 04:00 06:00 08:00 10:00 12:00 14:00 Hour 16:00 18:00 20:00 22:00 Figure 6.4: a) Kdirect(θ), and Kdiffuse for Pavilion Lake under different sun angles, black � (PAR, θ) (dashed black) for sunny days lines are for C = 1.2 mg m-3 b) Depth averaged 𝐾 in June. K(λ, θ) for Pavilion Lake under different sun angles calculated using empirical equations. Grey area indicates the upper and lower bounds for K(θ) for C between 2 and 0.5 mg m-3. Topographic shading of one of the sensors can explain these consistent diurnal � (PAR, θ), as will be further discussed after comparing the calculated 𝐾 � (PAR, signals in 𝐾 θ) with those calculated using empirical equations. Chapter 6. Temporal variation of attenuation 85 The spectral attenuation coefficients, Kdirect(λ, θ) and Kdiffuse(λ), for Pavilion Lake are estimated using equations (2-9) and (2-10). The absorption coefficient in these equations is calculated, assuming that chlorophyll is the only component that affects the attenuation in the lake, following Bricaud et al. (1995), 𝑎 = 𝑎𝑝𝑤 + 𝐶𝐴(𝜆)𝐶𝐵(𝜆) (6-3) where, 𝑎𝐶ℎ𝑙 is the specific absorption coefficient for chlorophyll, 𝐴(𝜆) and 𝐵(𝜆) are empirical wavelength parameters. As with pure water Kdiffuse(λ) is constant during the day, while Kdirect(λ, θ) is variable with the zenith angle, (θ). Separate ET(z) profiles were calculated using Kdiffuse(λ), and Kdirect(λ, θ) under different zenith angles. By fitting single exponentials to these profiles, non-spectral � diffuse and 𝐾 � direct(θ) were estimated, for diffuse and direct attenuation coefficients, 𝐾 different zenith angles (Figure 6.4a). Finally, Equation (6-1) was used to estimate a bulk � (PAR, θ), using those values, 𝐾 � diffuse and diurnally variable attenuation coefficient, 𝐾 � direct(θ). Note that the depth averaged 𝐾 � (PAR, θ) values lie within the grey area (Figure 𝐾 6.4b). Because of the algal spring bloom and the high chlorophyll concentration in June, � (PAR, the values predicted using PAR measurements, are closer to the upper bound of 𝐾 θ) estimated using the empirical equations. Pavilion Lake is surrounded by mountains. At the beginning of the day, as the sun is low, the mountains on the east bank cast their shadow over the whole lake. As the sun rises, the shadow of the east mountains recedes and an increasing portion of the west bank becomes exposed to direct radiation. When the sun is high in the sky, the whole lake is exposed to direct radiation. As the sun starts to set, the mountains on the west bank begin to cast a shadow on the west shore, and gradually the shadow advances towards the east shore, until it covers the whole lake. To examine topographic shading in further detail, the spatial distribution of ET over the whole lake at six different times throughout a typical summer day; June 28th, is shown in Figure 6.5. The plotted data are predicted using the r.sun module of the Chapter 6. Temporal variation of attenuation Geographic Information System, GRASS GIS. Given solar parameters and the topography, the module computes Edirect(θ) and Ediffuse(θ). At the beginning of the day, while the sun is low, the east bank obstructs Edirect(θ) from reaching the submerged sensor, while nothing obstructs both components of ET(θ) from reaching the above-water sensor (Figure 6.3b). The spike observed on June 28th is associated with the reduction of ~160 W m-2 of Edirect(θ) from reaching the submerged sensor. This in turn results in an increased F(θ) observed at the submerged sensor (F(θ)sub); ET(θ) at the above-water and the submerged sensors are 215 W m-2 and 55 W m-2, while F(θ) at the above-water sensor (F(θ)dry) is 0.3 and F(θ)sub is 1, respectively (Figure 6.5b). Now we move on to the part of the day, when both sensors are exposed to the same magnitude and composition of light with F(θ)dry = F(θ)sub ranging between 0.18 and � (PAR, θ) depends exclusively on θ and 0.2 (Figure 6.5c and d). During this time, 𝐾 represents a composite of Kdirect(PAR, θ) and Kdiffuse(PAR). However, just before sunset, the shadow of the western mountains, reaches the above-water sensor, while nothing shades the submerged sensor. The incoming irradiance at the above-water and submerged sensors are 40 W m-2 and 120 W m-2, while F(θ)dry is 1 and F(θ)sub is 0.2. This divergence � (PAR, θ). The attenuation coefficient in F(θ) results in the sudden drop in the calculated 𝐾 calculated at the time when both sensors are only under Ediffuse(θ) right before sunset and right after sunrise (Figure 6.5a and f), can be estimated to be equal to Kdiffuse(PAR, θ) � (PAR, θ) are assumed to be ~0.32 m-1. Any other slight fluctuations in the calculated 𝐾 associated to the turbulence in the top layer caused by winds, internal waves, or scattered clouds. 86 Chapter 6. Temporal variation of attenuation Figure 6.5: Progression of the distributions of PAR over Pavilion Lake on a typical summer day, June 28th a) 7:30 am, b) 8:00 am, c) 1:00 pm, d) 2:15 pm, e) 6:00 pm, and f) 7:30 pm. Contour colors represent the intensity of the incoming solar radiation, ET in W m-2. The lake is outlined in black, (+) and ( ) are the locations of the thermistor chain and the weather station. 87 Chapter 6. Temporal variation of attenuation 6.3.2 Seasonal variation in the attenuation coefficient Development of algal blooms, surface water temperature, wind speed, and the sun’s altitude are all factors that affect the composition of the incoming radiation and � (PAR, θ) (Paavel et al., 2008). In the contribute to the seasonal variations observed in 𝐾 previous section we focused on the composition of ET and its effect on the diurnal � (PAR, θ). In this section we will illustrate the seasonal variation in the variability in 𝐾 � (PAR, θ) in Pavilion Lake. The study period extends from August 9th 2006, calculated 𝐾 to April 24th 2007. To filter out the diurnal variation, an energy weighted daily average is calculated � (PAR, θ) at both depths 5.8 m and 17 m. The resulting attenuation coefficient is for 𝐾 � (PAR, θ) values calculated from both depths follow similar shown in Figure 6.6. The 𝐾 seasonal trends. An abrupt increase in the attenuation occurs after the formation of ice. This is to be expected since ice-cover is not considered in the calculations, and while black ice has similar optical properties to clear water, snow and white ice have high albedo and light absorption. Consideration of ice cover is beyond the scope of this study. After ice off, the attenuation drops to 0.2 m-1 at z = 5.8 m and 0.06 m-1 at z = 17 m, � (PAR, θ) at 17 m depth, is clearer indicating high water clarity. The apparent value for 𝐾 � (PAR, θ) at 5.8 m depth. For the period with no ice, the and with less fluctuation than 𝐾 � (PAR, θ) with PAR at 5.8 m depth is significant, with a correlation relation between 𝐾 coefficient of 0.7, however at z = 17 m, the relation is even more significant, with a correlation coefficient of 0.84. These coefficients show the stronger dependence of � (PAR, θ) on the magnitude of ET at deep depths than at shallower depths. On seasonal 𝐾 � (PAR, θ) increase by 0.3 m-1 during the period starting mid-August to average both 𝐾 mid-January. 88 Chapter 6. Temporal variation of attenuation ↓ 300 ↓ -2 ET(0) (W m ) 89 200 100 a) 0 0.8 0.7 b) ICECOVER -1 K(PAR)(m ) 0.6 0.5 0.4 0.3 z = 5.8 m 0.2 0.1 z = 17 m 0 1Sept 01Oct 01Nov 01Dec 01Jan 01Feb 01Mar 01Apr Date (2006/2007) Figure 6.6: Running weekly average of calculated attenuation coefficient at z = 5.8 m � (dashed black). Grey lines delineate the (blue), z = 17 m (red), and a depth averaged 𝐾 ice-cover period (Lim et al., 2009). � almost Note that during the period of interest, September 1st to November 15th, 𝐾 doubled, from 0.1 m-1 to 0.2 m-1. Two factors contribute to this increase. First, the increase in chlorophyll concentration, which occurs as the epilimnion deepens, and nutrients are released to the rest of the water column (Pettersson et al., 2003). Second, the � . Figure 6.7, shows the increase in F, Ediffuse/ET, which consequently increases 𝐾 insignificant contribution of the second factor in the increase. Calculations for this figure Chapter 6. Temporal variation of attenuation 90 -2 PAR (Wm ) 200 150 100 50 a) 0 150 -2 E (Wm ) Edirect 100 50 b) Ediffuse 0 0.11 -1 K (m ) 0.1095 0.109 0.1085 0.108 01Sep c) 15Sep 30Sep 15Oct Date (2006) 30Oct 15Nov Figure 6.7: a) Diurnal PAR at the sea surface, b) Decomposition of PAR into Ediffuse (black dashed) and Edirect (black solid), c) Attenuation coefficient calculated based on the PAR decomposition without including any change in chlorophyll concentration. only consider the change in F without any consideration to the change in chlorophyll concentration. Chapter 6. Temporal variation of attenuation 6.4 Model results To provide a more accurate parameterization of the attenuation coefficient in DYRESM, the heater module was modified to use a different attenuation coefficient value for each time step. This was done to test the significance of accounting for the temporal variation in the attenuation coefficient on the evolution of the thermal structure of lakes. Two sets of simulations are run for Pavilion Lake: the first set uses the original model and the second set uses the modified version. Model results are compared with the observed field data. 6.4.1 Diurnal Variation First, we test for the significance of accounting for the diurnal variability in � (PAR, θ). The chosen study period is the same presented in Figure 6.3, June 25th 2006 𝐾 to July 2nd 2006. Two simulations are run for this study period. The first is using the modified DYRESM and a temporally variable attenuation coefficient. The value of � (PAR, θ) is the calculated as the depth averaged 𝐾 � (PAR, θ) calculated at z = 5.8 and 𝐾 17 m depths, shown in Figure 6.3. The second is using the original DYRESM version, with a temporally constant attenuation coefficient of value of 0.13 m-1. This value � (PAR, θ) used corresponds to the time average, of one week, of the depth-averaged 𝐾 above. The short studied period and the depth average resulted in a value nearly identical to that in chapter 5, though it was determined independently from PAR profiles presented in chapter 5. A meteorological station located on the dock between the north and central basins (Figure 3.1) provided incoming shortwave radiation, wind speed (Uw), air temperature (Ta), vapor pressure, and cloud cover recorded at 15 minute intervals. Wind and air temperature follow the same diurnal cycles as the solar radiation (Figure 6.8). Maximum wind speed and air temperature recorded were 5 m s-1 and 32 oC, while the minimum was 10 oC during night time. The minimum daily averaged wind speed, during this study period, was recorded on day June 28th while the maximum was on June 30th. 91 Chapter 6. Temporal variation of attenuation 92 5 (a) -1 Uw (m s ) 4 3 2 1 0 40 35 (b) Ta (°C) 30 25 20 15 10 26Jun 28Jun 30Jun 02Jul Date (2006) Figure 6.8: (a) Wind speed, and (b) air temperature observed at Pavilion Lake for the studied period in 2006. During the studied period, water temperatures at 9.5, 15.5 and 26.5 m depth were recorded every 30 minutes, using Hobo Water Temp Pros, with an accuracy of ± 0.2 ºC, deployed at the lakebed. Clear diurnal signals were observed at 9.5 m depth as shown in Figure 6.9. The average rates of increase at the three depths were 0.18 oC day-1, 0.055 oC day-1, and 0.02 oC day-1. Chapter 6. Temporal variation of attenuation 93 13 12 9.5 m 11 Temp(°C) 10 9 8 15.5 m 7 6 26.5 m 5 4 26Jun 28Jun 30Jun 02Ju Date (2006) Figure 6.9: Temperature at 9.5 m, 15.5 m, and 26.5 m; observed (black), estimated using constant attenuation coefficient (blue), and variable attenuation coefficient (red). The two models make almost identical predictions of temperature (Figure 6.9). Using a diurnally variable attenuation coefficient rather than the constant one did not improve the model results. Estimated heating rates are insignificantly different than the observed ones at the three different depths. In other words, at any depth the diurnal � ((PAR, θ), is similar to the diurnal integration of the heat flux calculated using 𝐾 integration of heat flux using a constant attenuation coefficient of 0.13 m-1, ∫ 𝐸𝑇 𝑒 −𝐾𝑇 (𝑃𝐴𝑅,θ)z 𝑑𝑡 = ∫ 𝐸𝑇 𝑒 −0.13𝑧 𝑑𝑡. Chapter 6. Temporal variation of attenuation 6.4.2 Seasonal Variation The significance of the seasonal variation of irradiance attenuation on the thermal structure of Pavilion Lake was also examined on a seasonal scale. Two simulations were performed, the first using the modified DYRESM with a temporally variable attenuation � (PAR, θ) used is the calculated depth averaged 𝐾 � (PAR, θ) shown in Figure coefficient. 𝐾 6.6. The second simulation was performed using the original DYRESM, and a temporally constant attenuation coefficient equal to 0.13 m-1. The study period starts on September 1st 2006 and ends on November 15th 2006. During this period, the daily averaged air temperature decreased relatively uniformly from 20oC to 4oC (Figure 6.10b), except for a cool period with an average below zero, starting on October 28th 2006 to October 30th 2006. Observed temperatures show the beginning of the fall turnover. Temperature recorded at z = 9.5 m cooled throughout the simulation period, and the rate of cooling increased on cloudy days (Figure 6.11b). Temperatures at z = 15.1 and 24.1 m are in constant warming throughout the studied period. The combination of high wind speed, low Ta and drop in the vapor pressure, after cloudy periods, result in mixing the surface layer down to 15.5 m on October 15th 2006 and down to 24.1 m on November 12th 2006. Both simulations captured the evolution of temperature and fall turnover, (Figure 6.11b). However, using a temporally constant attenuation coefficient underestimated the heat transferred from the surface starting on October 15th till the end of the simulation. It underestimated the cooling rate at z = 9.5 m and the heating rates at z = 15.1 and 24.1 m. Consequently, the mixing of the surface temperature down to 15.5 m, is delayed until October 28th 2006. The rates estimated using the temporally variable attenuation coefficients are similar to those observed. Nevertheless, the mixing of the surface layer down to 24.1 m is delayed by 2 days. Heat fluxes calculated from the two simulations are similar. Sensible, latent and longwave fluxes estimated using both parameterizations are equal, as a result of the agreement between the surface skin temperatures estimated. 94 Chapter 6. Temporal variation of attenuation 95 12 -1 Uw (m s ) 10 (a) 8 6 4 2 0 30 Ta (°C) 20 10 0 (b) Vapor Pressure(hPa) -10 15 10 5 0 01Sep (c) 15Sep 30Sep 15Oct 30Oct 15Nov Date(2006) Figure 6.10: (a) Wind speed, (b) air temperature, and (c) Vapor Pressure observed at Pavilion Lake for the model period in 2006. 400 -2 ET(0) (W m ) Chapter 6. Temporal variation of attenuation 96 a) 200 0 20 b) 9.5 m Temp (°C) 15 15.1 m 10 24.1 m 5 01Sep 15Sep 30Sep 15Oct 30Oct 15Nov Date (2006) Figure 6.11: Running weekly averages of the temperature at 9.5 m, 15.1 m, and 24.1 m; observed (black), estimated using constant attenuation coefficient (blue), and temporally variable attenuation coefficient (red). 6.5 Conclusions The total attenuation coefficient of water is sensitive to the solar zenith angle. While, Kdirect(λ, θ) is dependent on θ, Kdiffuse(λ) remains constant with θ. When the sun is vertical, Kdirect(λ, θ) and Kdiffuse(λ, θ) are equal. Chapter 6. Temporal variation of attenuation � (PAR, θ) for Pavilion Lake clearly illustrate Diurnal variations in the calculated 𝐾 the temporal dependence of the attenuation coefficient on the decomposition of ET, introduced through the variation of θ and the spatial distribution of light. For Pavilion Lake, the increase in chlorophyll concentration, which occurs as the � during the fall epilimnion deepens, is the dominant factor contributing in increasing 𝐾 �. period, while the increase in F, has a minor effect on 𝐾 The methodology used provided a reasonable approach to explore the significance of the diurnal and seasonal variation in K(PAR, θ) in improving the accuracy of predicting the thermal structure of the Lake. For this specific clear lake, seasonal variations in the attenuation coefficient proved to be important in determining surface and hypolimnion temperatures (Figure 6.11). However, diurnal variations showed no significant improvements in the predicted thermal structure of the lake (Figure 6.9). The attention given to the diurnal variability in K(PAR, θ) is not essential. 97 Chapter 7. Conclusions 7 98 CONCLUSIONS 7.1 Thesis summary In this study the significance of accounting for the spectral variability of subsurface spectral irradiance has been explored and the effect of addressing the temporal variation of light attenuation on the evolution of the thermal structure of lakes was tested. Under certain conditions, such as high water clarity, the spectral and temporal variability in the light attenuation coefficient can play a significant role in the evolution of the thermal structure of lakes. Chapter 4, illustrated the importance of accounting for the spectral variability of subsurface spectral irradiance on the heating of lakes with different turbidities. We started by comparing the subsurface spectral irradiance between Crater Lake and Lake Superior. Then the subsurface heating calculated assuming a monochromatic light, and using the spectral irradiance were compared, using both theoretical predictions and field measurements of hypolimnetic heating in Pavilion and Crater Lakes. In chapter 5, the heater module in DYRESM was modified to accept different attenuation coefficient values for different bands in the PAR domain. The modified version was used to consider other heat fluxes present in lakes. The original and modified model results were compared with the field data from Pavilion Lake. Chapter 6, examined how the attenuation coefficient varied diurnally due to solar altitude for pure water and calculated the temporal variation in the coefficient from real PAR measurements in Pavilion Lake. The significance of the variation in the attenuation coefficient, on both diurnal and seasonal time scales, on the thermal evolution of Pavilion Lake was tested using a modified version of DYRESM. Chapter 7. Conclusions 99 7.2 Main findings • The subsurface irradiance intensity is a function of the wavelength dependence of both the incoming shortwave flux and the attenuation coefficient of the water body. • Wavelengths corresponding to high spectral attenuation coefficient, K(λ), get absorbed at shallower depths in a water column. For any water body, the wavelength of the maximum irradiance at the surface corresponds to the wavelength of the maximum incoming shortwave. With depth, the wavelength of maximum irradiance gradually shifts to the wavelength corresponding to the minimum K(λ) value. • Significant errors can occur when estimating the subsurface irradiance intensity using a non-spectral attenuation coefficient instead of a spectral one. Inaccuracies increase with increasing water clarity and increasing variation in K(λ). • Low attenuation coefficients have greater effects on hypolimnion temperatures, as they imply irradiance propagating to deeper depths. Accounting for the spectral variation in K(λ) influences heating rates below the thermocline, away from the effect of non-penetrative surface fluxes and mixing processes in Crater and Pavilion Lakes. • High attenuation coefficients directly influence the near surface temperature in lakes, and thereby can have a significant affect on the non-penetrative surface heat fluxes. This, in turn, can affect the depth and temperature of the surface mixed layer as seen for Pavilion Lake. • Modeling of the hypolimnion of Pavilion Lake, using DYRESM, modified to distribute PAR using multiple spectral bands for attenuation coefficient that better describes the spectral variability in K(λ), accurately predicted observed heating rates. This indicates that the attenuation of solar flux may be an important mechanism of heat transfer in the hypolimnion. • Attenuation coefficient values representing different bands of the K(λ) spectrum can be estimated by fitting a sum of exponential functions with different attenuation coefficient values to a single PAR profile. However, the profile has to Chapter 7. Conclusions 100 be measured on a sunny day. More bands were required to accurately calculate PAR profiles for the relatively clear Pavilion Lake than for the turbid Tailings Lake. • Temporal variation in the decomposition of ET, to Edirect and Ediffuse, due to the variation of θ and the spatial distribution of light, resulted in diurnal variations in � (PAR, θ), for Pavilion Lake, with a the calculated total attenuation coefficient, 𝐾 • � (PAR, θ) value occurring during solar noon. minimum 𝐾 � (PAR) during fall The increase in F, Ediffuse/ET, has a minor effect in increasing 𝐾 in Pavilion Lake. • � (PAR, θ) showed no significant improvements in the The diurnal variations in 𝐾 predicted thermal structure of Pavilion Lake. However, seasonal variations in the attenuation coefficient proved to be important in determining surface and hypolimnetic temperatures. 7.3 Future work A better parameterization of the subsurface irradiance in hydrodynamic models is warranted, especially for modeling clear lakes. The subsurface irradiance parameterization should include variation both with wavelength and time. Further work on the likely associated errors estimated for different water turbidities, under different meteorological conditions should be explored when modeling the thermal structure of lakes, which would require the collection of wavelength dependent light data from a variety of lakes. 101 BIBLIOGRAPHY Adiyanti, S., Imberger, J., 2007. Derivation of total diffuse attenuation coefficient from water column temperature data and meteorological water surface fluxes: A simple management tool. International Journal of River Basin Management 5, 293–303. Anderson, Thomas R, 1993. A spectrally averaged model of light penetration and photosynthesis. Limnology and Oceanography 38(7), 1403-1419. Antenucci, J., Horn, D., 2002. DYRESM: Dynamic reservoir simulation model. Water 29. Antenucci, J.P., Imerito, A., 2001. The CWR Dynamic Reservoir Simulation Model (DYRESM) Science Manual. Report WP 1573 JA. Anthony, K.R.N., Ridd, P.V., Orpin, A.R., Larcombe, P., Lough, J., 2004. Temporal variation of light availability in coastal benthic habitats: Effects of clouds, turbidity, and tides. Limnology and Oceanography 49, 2201–2211. Austin, R.W., Petzold, T.J., 1986. Spectral dependence of the diffuse attenuation coefficient of light in ocean waters. Optical Engineering 25, 471–479. Baker, K.S., Smith, R.C., 1982. Bio-optical classification and model of natural waters. 2. Limnology and Oceanography 27(3), 500–509. Bloss, S., Harleman, D.R.F., 1979. Effect of wind-mixing on the thermocline formation in lakes and reservoirs. Massachusetts Institute of Technology, Department of Civil Engineering, Ralph M. Parsons Laboratory for Water Resources and Hydrodynamics. 102 Boss, E.S., Collier, R., Larson, G., Fennel, K., Pegau, W.S., 2007. Measurements of spectral optical properties and their relation to biogeochemical variables and processes in Crater Lake, Crater Lake National Park, OR. Hydrobiologia 574, 149–159. Bricaud, A., Babin, M., Morel, A., & Claustre, H., 1995. Variability in the chlorophyll-specific absorption coefficients of natural phytoplankton: Analysis and parameterization. Journal of geophysical Research, 100(C7), 13321-13. Buiteveld, H., Hakvoort, J.H.M., Donze, M., 1994. Optical properties of pure water. Ocean Optics XII, pp. 174-183. International Society for Optics and Photonics. Bukata, R.P., 1995. Optical properties and remote sensing of inland and coastal waters. CRC Pressl Llc. Chang, G.C., Dickey, T.D., 2004. Coastal ocean optical influences on solar transmission and radiant heating rate. Journal of Geophysical Research 109, C01020. Cole, T.M., Buchak, E.M., 1995. CE-QUAL-W2: A Two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model, Version 2.0. User Manual. Coluccio, L., Eisinberg, A., Fedele, G., 2007. A Prony-like method for nonuniform sampling. Signal Processing 87, 2484–2490. Davies-Colley, Robert J., William N. Vant, and David Glynn Smith, 2003. Colour and clarity of natural waters. Science and management of optical water quality. Blackburn Pr. De Stasio, B.T., Hill, D.K., Kleinhans, J.M., Nibbelink, N.P., Magnuson, J.J., 1996. Potential effects of global climate change on small north-temperate lakes: Physics, fish, and plankton. Limnology and Oceanography 41, 1136–1149. 103 Dejak, C., Franco, D., Pastres, R., Pecenik, G., Solidoro, C., 1992. Thermal exchanges at air-water interfacies and reproduction of temperature vertical profiles in water columns. Journal of Marine Systems 3, 465–476. Dera, J., Gordon, H.R., 1968. Light field fluctuations in the photic zone. Limnology and Oceanography 13(4), 697–699. Dobson, H.F.H., Gilbertson, M., Sly, P.G., 1974. A summary and comparison of nutrients and related water quality in Lakes Erie, Ontario, Huron, and Superior. Journal of the Fisheries Board of Canada 31, 731–738. Dickey, T. D., and J. J. Simpson, 1983. The influence of optical water type on the diurnal response of the upper ocean. Tellus B 35(2), 142-154. Gast, P. R., 1960. Solar radiation. Handbook of Geophysics, rev. ed. Macmillan. Hamblin, P. F., and J. Imberger, 1984. Classification and dynamic simulation of the vertical density structure of lakes. Limnology and Oceanography 29(4), 845-86. Han, B.P., Armengol, J., Carlos Garcia, J., Comerma, M., Roura, M., Dolz, J., Straskraba, M., 2000. The thermal structure of Sau Reservoir (NE: Spain): a simulation approach. Ecological Modelling 125, 109–122. Hargreaves, B.R., Girdner, S.F., Buktenica, M.W., Collier, R.W., Urbach, E., Larson, G.L., 2007. Ultraviolet radiation and bio-optics in Crater Lake, Oregon. Longterm Limnological Research and Monitoring at Crater Lake, Oregon 107–140. Hieronymi, M., 2011. Solar radiative transfer into the ocean: A study on underwater light fluctuations due to surface waves. PhD dissertation, ChristianAlbrechts-Universität, 2011. 104 Hieronymi, M., Macke, A., 2010. Spatiotemporal underwater light field fluctuations in the open ocean. Journal of the European Optical Society-Rapid Publications 5, 1-8. Hodges, B., Dallimore, C., 2006. Estuary, Lake and Coastal Ocean Model: ELCOM. Science Manual. Centre of Water Research. University of Western Australia. Hornung, R., 2003. Numerical Modelling of Stratification in Lake Constance with the 1-D hydrodynamic model DYRESM. M.S. thesis, University of Stuttgart, Germany. Imberger, J., Hamblin, P.F., 1982. Dynamics of lakes, reservoirs, and cooling ponds. Annual Review of Fluid Mechanics 14, 153–187. Imberger, J., Ivey, G.N., 1991. On the nature of turbulence in a stratified fluid. II, Application to lakes. Journal of Physical Oceanography 21, 659–680. Imboden, Dieter M., and Alfred Wüest, 1995. Mixing mechanisms in lakes. Physics and chemistry of lakes. Springer Berlin Heidelberg, 83-138. Imerito, A., 2007. DYRESM v4. 0 science manual. Perth, Australia: Centre for Water Research, University of Western Australia. Jassby, A., Powell, T., 1975. Vertical patterns of eddy diffusion during stratification in Castle Lake, California. Limnology and Oceanography 530–543. Jerlov, N.G., 1968. Optical Oceanography. Elsevier. Kirk, J.T.O., 1979. Spectral distribution of photosynthetically active radiation in some south-eastern Australian waters. Marine and Freshwater Research 30, 81–91. Kirk, J.T.O., 1984. Dependence of relationship between inherent and apparent optical properties of water on solar altitude. Limnology and Oceanography 350–356. 105 Kirk, J.T.O., 1994. Light and photosynthesis in aquatic ecosystems. Cambridge Univ Pr. Larson, G.L., Hoffman, R.L., McIntire, D.C., Buktenica, M.W., Girdner, S.F., 2007. Thermal, chemical, and optical properties of Crater Lake, Oregon. Hydrobiologia 574, 69–84. Lee, Z.P., Carder, K.L., Arnone, R.A., 2002. Deriving inherent optical properties from water color: a multiband quasi-analytical algorithm for optically deep waters. Applied Optics 41, 5755–5772. Lee, Z.P., Du, K.P., Arnone, R., Liew, S.C., Penta, B., 2005. Penetration of solar radiation in the upper ocean: A numerical model for oceanic and coastal waters. Journal of Geophysical Research: Oceans 110, C09019. LI-COR Underwater Radiation Sensors Instruction Manual, 2006. Licor Biosciences. Lim, D.S.S., Laval, B.E., Slater, G., Antoniades, D., Forrest, A.L., Pike, W., Pieters, R., Saffari, M., Reid, D., Schulze-Makuch, D., others, 2009. Limnology of Pavilion Lake, BC, Canada Characterization of a microbialite forming environment. Fundamental and Applied Limnology; Hydrobiologie 173, 329–351. Liu, C.C., Miller, R.L., Carder, K.L., Lee, Z., D’Sa, E.J., Ivey, J.E., 2006. Estimating the underwater light field from remote sensing of ocean color. Journal of Oceanography 62, 235–248. Losordo, T.M., Piedrahita, R.H., 1991. Modelling temperature variation and thermal stratification in shallow aquaculture ponds. Ecological Modelling 54, 189–226. Mecherikunnel, A., Duncan, C.H., 1982. Total and spectral solar irradiance measured at ground surface. Applied Optics 21, 554–556. 106 Michalski, J., Lemmin, U., 1995. Dynamics of vertical mixing in the hypolimnion of a deep lake: Lake Geneva. Limnology and oceanography 809–816. Matthews, P. C., and S. I. Heaney, 1987. Solar heating and its influence on mixing in ice covered lakes. Freshwater Biology 18(1), 135-149. Morel, A., 1988. Optical modeling of the upper ocean in relation to its biogenous matter content(case I waters). Journal of Geophysical Research 93, 749–10. Morel, A., 2009. Are the empirical relationships describing the bio-optical properties of case 1 waters consistent and internally compatible. Journal of Geophysical Research 114, C01016. Morel, A., Antoine, D., 1994. Heating rate within the upper ocean in relation to its bio-optical state. Journal of Physical Oceanography 24, 1652–1665. Morel, A., Maritorena, S., 2001. Bio-optical properties of oceanic waters: A reappraisal. Journal of Geophysical Research 106, 7163–7180. Mueller, J.L., 2000. SeaWiFS algorithm for the diffuse attenuation coefficient, K (490), using water-leaving radiances at 490 and 555 nm. SeaWiFS postlaunch calibration and validation analyses, part 3, 24–27. Murtugudde, R., Beauchamp, J., McClain, C.R., Lewis, M., Busalacchi, A.J., 2002. Effects of penetrative radiation on the upper tropical ocean circulation. Journal of Climate 15, 470–486. Nakamoto, S., Kumar, S.P., Oberhuber, J.M., Ishizaka, J., Muneyama, K., Frouin, R., 2001. Response of the equatorial Pacific to chlorophyll pigment in a mixed layer isopycnal ocean general circulation model. Geophysical Research Letters 28, 2021–2024. 107 Osborne, M.R., Smyth, G.K., 1995. A modified Prony algorithm for exponential function fitting. SIAM Journal on Scientific Computing 16, 119–138. Paavel, B., Arst, H., Reinart, A., 2008. Variability of bio-optical parameters in two North-European large lakes. European Large Lakes Ecosystem changes and their ecological and socioeconomic impacts 201–211. Paavel, B., Arst, H., Reinart, A., Herlevi, A., 2006. Model calculations of diffuse attenuation coefficient spectra in lake waters. Proceedings of the Estonian Academy of Sciences, Biology and Ecology 55, 61-81. Pan, X., Zimmerman, R.C., 2010. Modeling the vertical distributions of downwelling plane irradiance and diffuse attenuation coefficient in optically deep waters. Journal of Geophysical Research 115, C08016. Paulson, C.A., Simpson, J.J., 1977. Irradiance measurements in the upper ocean. Journal of Physical Oceanography 7, 952–956. Pérez-Fuentetaja, A., Dillon, P.J., Yan, N.D., McQueen, D.J., 1999. Significance of dissolved organic carbon in the prediction of thermocline depth in small Canadian Shield lakes. Aquatic Ecology 33, 127–133. Perroud, M., Goyette, S., Martynov, A., Beniston, M., Anneville, O., others, 2009. Simulation of multiannual thermal profiles in deep Lake Geneva: A comparison of onedimensional lake models. Limnology and Oceanography 54, 1574. Pettersson, K., Grust, K., Weyhenmeyer, G., & Blenckner, T., 2003. Seasonality of chlorophyll and nutrients in Lake Erken–effects of weather conditions. Hydrobiologia, 506(1-3), 75-81. Pieters, R., Lawrence, G.A., 2009. Effect of salt exclusion from lake ice on seasonal circulation. Limnology and Oceanography 54, 401–412. 108 Pimentel, S., Haines, K., Nichols, N.K., others, 2008. Modeling the diurnal variability of sea surface temperatures. Journal of Geophysical Research 113, C11004. Poole, H.H., Atkins, W.R.G., 1929. Photo-electric measurements of submarine illumination throughout the year. Journal of the Marine Biological Association of the United Kingdom (new series) 16, 297–324. Potts, D., Tasche, M., 2010. Parameter estimation for exponential sums by approximate Prony method. Signal Processing 90, 1631–1642. Preisendorfer, R.W., 1965. Radiative transfer on discrete spaces. Pergamon Press. Preisendorfer, R.W., 1986. Secchi disk science: Visual optics of natural waters. Limnology and Oceanography 31(5), 909–926. Reinart, A., Herlevi, A., 1999. Diffuse attenuation coefficient in some Estonian and Finnish lakes, in: Proceedings of the Estonian Academy of Sciences, Biology Ecology 48(4), pp. 267–283. Reinart, A., Herlevi, A., Arst, H., & Sipelgas, L., 2003. Preliminary optical classification of lakes and coastal waters in Estonia and south Finland. Journal of Sea Research, 49(4), 357-366. Robock, A., 1980. The seasonal cycle of snow cover, sea ice and surface albedo. Monthly Weather Review 108, 267–285. Schneider, E.K., Zhu, Z., 1998. Sensitivity of the simulated annual cycle of sea surface temperature in the equatorial Pacific to sunlight penetration. Journal of Climate 11, 1932–1950. 109 Simpson, J.J., Dickey, T.D., 1981. Alternative parameterizations of downward irradiance and their dynamical significance. Journal of Physical Oceanography 11, 876– 886. Spigel, R.H., Imberger, J., 1980. The classification of mixed-layer dynamics in lakes of small to medium size. Journal of Physical Oceanography 10(7), 1104–1121. Spigel, Robert H., and Clive Howard Williams, 1984. A method for comparing scalar irradiance measurements with upward and downward irradiance in lakes. Water resources research 20 (4), 507-509. Stramska, M., Dickey, T.D., 1998. Short-term variability of the underwater light field in the oligotrophic ocean in response to surface waves and clouds. Deep-Sea Research Part I 45, 1393–1410. Stramska, M., Frye, D., 1997. Dependence of apparent optical properties on solar altitude: Experimental results based on mooring data collected in the Sargasso Sea. Journal of Geophysical Research 102, 15–15. Stramski, D., Booth, C.R., Greg Mitchell, B., 1992. Estimation of downward irradiance attenuation from a single moored instrument. Deep Sea Research Part A. Oceanographic Research Papers 39, 567–584. Sweeney, C., Gnanadesikan, A., Griffies, S.M., Harrison, M.J., Rosati, A.J., Samuels, B.L., 2005. Impacts of shortwave penetration depth on large-scale ocean circulation and heat transport. Journal of Physical Oceanography 35, 1103–1119. Tanentzap, A.J., Hamilton, D.P., Yan, N.D., 2007. Calibrating the Dynamic Reservoir Simulation Model (DYRESM) and filling required data gaps for onedimensional thermal profile predictions in a boreal lake. Limnology and Oceanography: Methods 5, 484–494. 110 Tanentzap, A.J., Yan, N.D., Keller, B., Girard, R., Heneberry, J., Gunn, J.M., Hamilton, D.P., Taylor, P.A., 2008. Cooling lakes while the world warms: Effects of forest regrowth and increased dissolved organic matter on the thermal regime of a temperate, urban lake. Limnology and Oceanography 53(1), 404–410. Tyler, J.E., Smith, R.C., 1970. Measurements of spectral irradiance underwater. Gordon & Breach Publishing Group. Umlauf, L., Burchard, H., Bolding, K., 2006. GOTM: source code and test case documentation. Vila, X., Colomer, J., Garcia-Gil, L.J., 1996. Modelling spectral irradiance in freshwater in relation to phytoplankton and solar radiation. Ecological Modelling 87, 59– 68. Weinberger, S., Vetter, M., 2012. Using the hydrodynamic model DYRESM based on results of a regional climate model to estimate water temperature changes at Lake Ammersee. Ecological Modelling 244, 38–48. Yeates, P.S., Imberger, J., 2003. Pseudo two-dimensional simulations of internal and boundary fluxes in stratified lakes and reservoirs. International Journal of River Basin Management 1, 297–319. Zaneveld, J.R.., 1989. An asymptotic closure theory for irradiance in the sea and its inversion to obtain the inherent optical properties. Limnology and Oceanography 40(5), 1442–1452. Zheng, X., Dickey, T., Chang, G., 2002. Variability of the downwelling diffuse attenuation coefficient with consideration of inelastic scattering. Applied optics 41, 6477–6488. 111 Zielinski, O., Oschlies, A., Reuter, R., 1998. Comparison of underwater light field parameterizations and their effect on a 1-dimensional biogeochemical model at station estoc, north of the Canary islands. Ocean Optics XIV, Kailua-Kona, HI. Appendices 112 APPENDICES APPENDIX A: Empirical water functions Table A.1: 𝜒(λ), 𝑒(λ), and 𝐾𝑤 (λ) adapted from Morel (1988) 𝛌 (μm) 0.605 0.61 0.615 0.62 0.625 0.63 0.635 0.64 0.645 0.65 0.655 0.66 0.665 0.67 0.675 0.68 0.685 0.69 0.695 0.7 𝑲𝒘 (m-1) 0.279 0.296 0.303 0.31 0.315 0.32 0.325 0.33 0.34 0.35 0.37 0.405 0.418 0.43 0.44 0.45 0.47 0.5 0.55 0.65 𝝌(𝛌) 0.035 0.036 0.0375 0.0385 0.04 0.042 0.043 0.044 0.0445 0.045 0.046 0.0475 0.049 0.0515 0.052 0.0505 0.044 0.039 0.034 0.03 𝒆(𝛌) 0.63 0.634 0.638 0.642 0.647 0.653 0.658 0.663 0.667 0.672 0.677 0.682 0.687 0.695 0.697 0.693 0.665 0.64 0.62 0.6 Appendices 113 APPENDIX B: Absorption and scattering functions for pure water Table B.2: Absorption (a) and Scattering (b) coefficients of pure water adapted from Kirk (1984) λ (nm) 400 405 410 415 420 425 430 435 440 445 450 455 460 465 470 475 480 485 490 495 500 505 510 515 520 525 530 535 540 545 550 a(m-1) 0.007 0.006 0.0056 0.0052 0.0054 0.0061 0.0064 0.0069 0.0083 0.0095 0.011 0.012 0.0122 0.0125 0.013 0.0143 0.0157 0.0168 0.0185 0.0213 0.0242 0.03 0.0382 0.0462 0.0474 0.0485 0.0505 0.0527 0.0551 0.0594 0.0654 b(m-1) 0.0053 0.005 0.0048 0.0045 0.0043 0.0041 0.0039 0.0037 0.0036 0.0034 0.0033 0.0031 0.003 0.0028 0.0027 0.0026 0.0025 0.0024 0.0023 0.0022 0.0021 0.002 0.0019 0.0018 0.0018 0.0017 0.0017 0.0016 0.0015 0.0015 0.0014 λ (nm) 555 560 565 570 575 580 585 590 595 600 605 610 615 620 625 630 635 640 645 650 655 660 665 670 675 680 685 690 695 700 705 a(m-1) 0.069 0.0715 0.0743 0.0804 0.089 0.1016 0.1235 0.1487 0.1818 0.2417 0.2795 0.2876 0.2916 0.3047 0.3135 0.3184 0.3309 0.3382 0.3513 0.3594 0.3852 0.4212 0.4311 0.4346 0.439 0.4524 0.469 0.4929 0.5305 0.6229 0.7522 b(m-1) 0.0014 0.0013 0.0013 0.0012 0.0012 0.0011 0.0011 0.0011 0.001 0.001 0.001 0.0009 0.0009 0.0009 0.0008 0.0008 0.0008 0.0008 0.0007 0.0007 0.0007 0.0007 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0005 0.0005 0.0005 Appendices 114 APPENDIX C: Local and Apparent attenuation coefficients It is important to distinguish between different methods for estimating the irradiance profiles, to be able to use the relevant K value for each method. The difference in the methods lies behind the reference value chosen. The irradiance profile calculated using the incoming shortwave as a constant reference, is different than that using variable subsurface irradiances as a reference. We introduce equivalent attenuation profiles, K(z), which are depth variable K values. These profiles ease the calculation of accurate irradiance profiles, without going through integrating calculated I(z, λ) spectra each time ET(z) is required. Two different approaches can be used to calculate K(z). The first is calculating K(z) from the logarithmic slope of the tangent at the irradiance at depth z, K local (z) = − 𝑑𝑙𝑛(𝐸(𝑧)) 𝑑𝑧 . Where Klocal(z) is the instantaneous attenuation coefficient. The second approach is easier in calculating. K(z) is calculated from the logarithmic slope of the line connecting the total irradiance at the surface (ET(0)) to the irradiance at z (Iz), K app (z) = 1 𝐸 (𝑧) − 𝑧 𝑙𝑛 𝐸𝑇(0) . The reference irradiance value in this approach is ET(0). Both profiles are 𝑇 valid in ET(z) calculations, however, attention should be taken to use the same method K(z) was first calculated with. If Kapp(z) is used then ET(z) should be calculated using 𝐸𝑇 (z) = 𝐸𝑇 (0)𝑒𝑥𝑝−𝐾𝑎𝑝𝑝 (𝑧)𝑧 , while if Klocal(z) is used then ET(z) should be calculated 𝑧 using 𝐸𝑇 (z) = 𝐸𝑇 (0)𝑒𝑥𝑝− ∫0 𝐾𝑙𝑜𝑐𝑎𝑙 (𝑧)𝑑𝑧 . While, Klocal(z) is an exact function of K(λ) and I(z, λ), Kapp(z) is an average of Klocal(z) profile above z. Accurately calculated ET(z) profiles, calculated using K(λ) are used to calculate K(z) (Figure C.1a) which shows how the composition of the San-Vicente water makes it the most turbid of the four cases. The High content of impurities in the Reservoir absorb and scatter higher percentages of the irradiances at the surface than Lake Superior, followed by Crater Lake and finally pure water. Pure water depletes the least irradiance near the surface. The average slope of the pure water ET(z) profile is the steepest. The slope gets milder as we go from Crater Lake, to Lake Superior, and finally being the mildest for San Vicente Reservoir water. Appendices 115 0 (a) 5 (b) Klocal - Kapp -- 10 15 z (m) 20 0.85 0.05 25 30 35 40 45 50 -15 10 -10 -5 10 10 -2 Iz (Wm ) 0 10 0 0.5 1 K(z) Figure C.1: a) Depth profiles of downward irradiance (0.4 - 0.7 µm) assuming an incoming irradiance of 100 Wm-2 b) K profiles showing the increase in K as we go near the surface, local K (dashed), and apparent K(solid), for pure water (black), Crater Lake (blue), Lake Superior (green), and San Vicente Lake (red). Calculated Klocal(z) and Kapp(z) profiles abate with depth because of the gradual absorption of wavelengths with higher K(λ) at shallower depths. They are equal at the surface, z = 0, however, Kapp(z) is greater than the Klocal(z) elsewhere, because its reference irradiance ET(0), is always of a bigger value than the reference irradiance of Klocal(z) (Figure C.1b). The surface K(z) values are 0.88, 0.54, 0.16, and 0.15 m-1 for SanVicente Reservoir, Lake Superior, Crater Lake, and pure water, respectively. K(z) values decrease reaching an asymptotic minimum with depth. The previous observation should be added to Kirk (1977) conclusions, who mentioned that pure water is almost biphasic with significantly higher attenuation coefficients near the surface than at greater depths. The percentage decrease in Kapp(z) and Klocal(z) is less abrupt for San-Vicente and Lake Superior waters than in Crater Lake and Pure waters. However, values decrease more significantly in turbid waters. Minimum Klocal(z) values calculated are 0.51, 0.35, 0.016, Appendices 116 and 0.04 m-1. These values correspond to the minimum K(λ) associated with the four cases. Appendices 117 APPENDIX D: Non-penetrative heat fluxes Latent Heat 0 0.622 𝑄𝑙ℎ = min � 𝐶𝐿 𝜌𝐴 𝐿𝐸 𝑈𝑎 (𝑒𝑎 − 𝑒𝑠 (𝑇𝑤 ))∆𝑡 𝑃 where, P is the atmospheric pressure (hectopascals), 𝐶𝐿 is the coefficient for latent heat transfer under 𝑈𝑎 , 𝑒𝑎 is the air vapor pressure (hectopascals), and 𝑒𝑠 (𝑇𝑤 ) is the saturation vapor pressure at 𝑇𝑤 (hectopascals). Sensible Heat 𝑄𝑠ℎ = 𝐶𝑠 𝜌𝐴 𝐶𝑃 𝑈𝑎 (𝑇𝑎 − 𝑇𝑤 )∆𝑡 where, 𝐶𝑠 is a dimensionless coefficient for sensible heat transfer under wind speed referenced at 10 m; 𝑈𝑎 , above the water surface, 𝜌𝐴 is the air density (kgm-3), 𝐶𝑝 is the specific heat of air under constant pressure (Jkg-1K-1). Longwave The net longwave flux entering the water surface can be expressed by Appendices 118 𝑄𝑙𝑤 = (1 − 𝑟𝑎𝑙𝑤 )(1 + 0.17𝐶 2 )𝜀𝑎 (𝑇𝑎 )𝜎𝑇𝑎4 − 𝜀𝑤 𝑇𝑤4 (𝑙𝑤) where, 𝑄𝑙𝑤 is the net incoming long wave radiation, 𝑟𝑎 is the albedo for the long wave radiation and it is a dimensionless unit equal to 0.03, 𝑇𝑤 is the water surface temperature in Kelvins, 𝜀𝑤 is the water surface emissivity, 𝜎 is the Stefan Boltzmann constant, 𝑇𝑎 is the air temperature in Kelvin, and 𝜀𝑎 (𝑇𝑎 ) = 9.37𝑒 −6 (𝑇𝑎2 ) (Antenucci and Horn 2002).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Effect of spectral and temporal variation of subsurface...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Effect of spectral and temporal variation of subsurface irradiance on the heating of lakes Nassar, Yasmin 2013
pdf
Page Metadata
Item Metadata
Title | Effect of spectral and temporal variation of subsurface irradiance on the heating of lakes |
Creator |
Nassar, Yasmin |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | This dissertation illustrates how the vertical distribution of temperature in lakes can be affected by the fact that the attenuation coefficient of light is often strongly dependent on wavelength. The potential importance of this spectral effect is first examined by considering the solar radiation in isolation, and then by including all non–penetrative heat fluxes using a modified version of the numerical model, DYRESM. Comparing the subsurface spectral irradiance of different lakes reveals that the spectral variability of the attenuation coefficient is more significant when calculating the light intensity in relatively clear lakes than in turbid lakes. Comparisons made between field measurements and theoretical predictions of hypolimnetic heating show the importance of accounting for the spectral irradiance for two relatively clear lakes: Pavilion Lake and Crater Lake. A new parameterization that better describes the spectral attenuation coefficient and the distribution of subsurface irradiance is added to DYRESM. The results obtained, when running the original and the modified DYRESM on Pavilion Lake, show a significant improvement in predicting the thermal structure of the lake with the modified version. The effects of the variation in solar angle and the seasonal variation in water quality on the attenuation coefficient are also examined for Pavilion Lake using DYRESM modified to accept a time varying attenuation coefficient. Simulations were performed for Pavilion Lake using the original and modified versions of DYRESM on diurnal and seasonal scales. Results show no significant improvement in the thermal evolution of the lake when considering the diurnal variations, while slight improvement was shown on a seasonal scale. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-04-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0073753 |
URI | http://hdl.handle.net/2429/44319 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2013-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2013_spring_nassar_yasmin.pdf [ 1.28MB ]
- Metadata
- JSON: 24-1.0073753.json
- JSON-LD: 24-1.0073753-ld.json
- RDF/XML (Pretty): 24-1.0073753-rdf.xml
- RDF/JSON: 24-1.0073753-rdf.json
- Turtle: 24-1.0073753-turtle.txt
- N-Triples: 24-1.0073753-rdf-ntriples.txt
- Original Record: 24-1.0073753-source.json
- Full Text
- 24-1.0073753-fulltext.txt
- Citation
- 24-1.0073753.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0073753/manifest