UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A technique for sub-pixel temperature retrieval from multispectral imagery and its application to Quesnel… Sentlinger, Gabriel Isaac 2005

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

Item Metadata

Download

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

Full Text

A TECHNIQUE FOR SUB-PIXEL TEMPERATURE RETRIEVAL FROM MULTTSPECTRAL IMAGERY A N D ITS APPLICATION TO QUESNEL LAKE, BRITISH COLUMBIA by GABRIEL ISAAC SENTLINGER B.Sc, University of British Columbia, 2001 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA August 2005 © Gabriel Isaac Sentlinger, 2005 ABSTRACT Water temperature is a significant environmental indicator in lakes, rivers and streams. It is important to many resource management and environmental fields: forestry practices in adjacent lands modify runoff patterns which in-turn affects water temperatures; accelerated glacial melting caused by climate change increases the flux of cold water; and spawning salmon stocks are lethally affected by thermal shock from small (<5°C) temperatures changes. Quesnel Lake in the interior of British Columbia, Canada, is an important lake in all of these respects, and is the subject of this study due to recent mass mortality of spawning salmon believed to be caused by extreme temperature fluctuations. The conventional means of tracking temperature are in-situ temperature loggers installed in a small number of critical locations. This approach is costly and prone to data loss however, and does not provide a spatially comprehensive measure of lake temperature patterns. Coupled with satellite Thermal Infrared (TIR) imagery, the two measures can provide a contiguous dataset in space, depth, and time. Achieving an adequate balance between temporal and spatial coverage in TIR imagery can be difficult. Sensors with a sufficiently short revisit time often do not provide adequate spatial resolution to resolve narrow reaches of lakes and rivers. The MODIS sensor provides very high temporal coverage, up to four visits per day, but at a resolution of 1km, which is insufficient for Quesnel Lake's narrow reaches. This research presented in this thesis attempts to reconcile this deficiency in available TIR products, to produce an adequately accurate and fine spatial and temporal resolution water surface temperature dataset. By using a priori vectorized water boundary data, we estimate the sub-pixel surface temperature to within 1°C using a gradient descent solution to the mixed pixel at-sensor radiance equation. Ground-leaving radiance is estimated from standard MODIS temperature and emissivity data products for pure pixels and a simple regression technique to estimate atmospheric effects. The resulting algorithm is simple, effective, and useful for scientists without extensive remote sensing expertise, essentially unlocking a largely untapped resource for limnological and hydrological studies. i i TABLE OF CONTENTS ABSTRACT II TABLE OF CONTENTS Ill LIST OF TABLES V LIST OF FIGURES VI LIST OF SYMBOLS VIII GLOSSARY IX ACKNOWLEDGEMENT X CO-AUTHORSHIP STATEMENT XI 1 INTRODUCTION 1 1.1 Fundamentals of Thermal Infrared Imagery Analysis 6 •1.1.1 Remote Sensing by Satellite 6 1.1.2 Planck's Law of Blackbody Radiation 9 1.1.3 Kirchoff s Law 12 1.1.4 Radiative Transfer 13 1.1.5 Split Window Technique 15 1.2 Spectral Mixture Analysis 17 1.2.1 Temperature End-members 18 1.2.2 SMA in the VNIR 20 1.2.3 SMA in the TIR Domain 22 1.3 Structure and Intent of Thesis 24 1.4 References 25 2 SUB-PIXEL WATER TEMPERATURE ESTIMATION FROM THERMAL-INFRARED IMAGERY USING VECTORIZED LAKE FEATURES 27 2.1 Introduction 27 2.2 Theoretical Background 29 2.2.1 Fractional abundance from Vectorized feature data 33 2.2.2 Atmospheric Correction from pure pixels 35 2.2.3 Emissivity from Pure Pixels 37 2.2.4 Sub-Pixel Temperature Recovery 38 2.2.5 Georegistration and Autoregistration 40 2.3 Results 41 2.4 Discussion 45 2.5 Conclusions 49 2.6 References 50 i i i 3 CHARACTERISING LAKE DYNAMICS IN A DEEP FJORD LAKE USING SUB-PIXEL TEMPERATURES DERIVED FROM MODIS IMAGERY 52 3.1 Introduction 52 3.1.1 TIR Imagery 53 3.1.2 Seiching 54 3.2 Site Description and Instrumentation 55 3.2.1 Satellite Imagery 57 3.2.2 Field Measurements 60 3.3 Results and Discussion 61 3.3.1 Antecedent Conditions 61 3.3.2 SPTR validation and the skin effect 64 3.3.3 Day 214 upwelling event 69 3.3.4 Summer 2003 MODIS trends and variance 77 3.4 Conclusions 84 3.5 References 87 4 CONCLUDING REMARKS 89 4.1 Future Research 91 4.2 References 92 APPENDIX A: DAY 217 SURFACE FEATURE 93 iv LIST OF TABLES Table 1-1: Sensor specification for MODIS and ASTER sensors (TIR = Thermal Infrared, SWIR= Short Wave Infrared, VNIR = Visible and Near Infrared (NASA 2005) 8 Table 2-1: Summary of error estimates 43 Table 3-1: Characteristics of Quesnel Lake Basins. Total(North) is the entire lake going north, Total(East) is the entire lake going East 64 Table 3-2: Multiple Regression Statistics for all variables. P-value is the probability that there is not a relationship 69 Table 3-3: Multiple Regression statistics for significant variables ...69 Table 3-4: Multiple Regression Results 69 Table 3-5: Estimated interface deflection at the sill for the three main events, based on wind speeds at Goose Spit 75 v LIST OF FIGURES Figure 1-1 Landsat image of Quesnel Lake and area 2 Figure 1-2 Schematic of pixel mixing: two material model 4 Figure 1-3 Schematic of satellite in sun-synchronous orbit 8 Figure 1-4 Radiance Changes as a Function of Wavelength 11 Figure 1-5 The Increase in Radiance as a Function of Temperature 11 Figure 1-6. The interaction of light with matter 12 Figure 1-7 Atmospheric attenuation of a signal 14 Figure 1-8 Two Unmixing Situations, Three Possible Unmixing Products 20 Figure 2-1 Sub-pixel Temperature Retrieval flow diagram 32 Figure 2-2 Abundance from reverse geo-referencing vector data into original lo-res space 34 Figure 2-3 At-sensor Radiance vs Estimated Ground-leaving radiance, Ch32 37 Figure 2-4 5x5 roving net to determine local emissivity 37 Figure 2-5 Gradient Descent convergence 39 Figure 2-6 SPTR results before and after autoregistration 41 Figure 2-7 Sub-Pixel Temperature Retrieval Results 44 Figure 2-8 Sub-Pixel Temperature Retrieval results: scatter plot 45 Figure 2-9 SPTR toolset with a quicktime movie 47 Figure 2-10 Sub-pixel Temperature retrieval results with water smoothing 48 Figure 2-11 Top 5 metres of water column during wind event 48 vi Figure 3-1 Quesnel Lake 57 Figure 3-2 Quesnel Lake depth and temperature profiles, Day 212 63 Figure 3-3 Comparison of SPTR results at M2 and M8 with thermistor temperatures and wind data 67 Figure 3-4 Skin effects through the day 68 Figure 3-5 Quesnel Lake T-Chain data, M2 and M8 72 Figure 3-6 Quesnel Lake Skin temperature during 214 upwelling event 73 Figure 3-7 Quesnel Lake Skin temperature during 214 upwelling event ; 74 Figure 3-8 Simulated flows into the West Basin epilimnion 76 Figure 3-9 Number of cloud-free images along W-E and S-N thalweg 79 Figure 3-10 Average lake temperature difference from mean, summer 2003 80 Figure 3-11 Solar Illumination of Quesnel Lake 81 Figure 3-12 Map of lake temperature variance about the average lake temperature, summer 2003 83 Figure 3-13 Covariance matrix of all Quesnel Lake transects 84 Figure A- l Surface feature traveling westward 94 Figure A-2 217 Surface feature in M8 temperature data 95 vii L I S T O F S Y M B O L S X wavelength, urn u wave number, cm1 B(A,Ti) planck function at wavelength A. and temperature T, Wm^m'sr"1 U s at-sensor radiance, Wm"2um'1sr'1 ground-leaving radiance, Wm'2um W 1 upwelling radiance, Wm'2umW n grouped reflected radiance, adjacency effect and noise, Wm'!|xm'1sr Tx transmittance or transmissivity EX emittance or emissivity ax absorbtance or absorbtivity rx reflectance or reflectivity f. i fractional abundance of ith material T temperature, Kelvin or Celcius P density, kg/m3 kx • absorption coefficient,m2/kg ds path of radiation, m u optical depth vii i GLOSSARY VNIR Visible Near Infrared TIR Thermal Infrared SWIR Short Wave Infrared SMA Spectral Mixture Analysis TES Temperature Emissivity Seperation DEM Digital Elevation Model ENVI Environment for Visualizing Imagery (software published by Research Systems Inc.) GSD Ground Sampling Distance IDL Interactive Data Language (software published by Research Systems Inc.) ISAC In Scene Atmospheric Correction MODTRAN MODerate resolution TRANsmission model (software published by Spectral Science Inc.) PSF Point Spread Function SRF Spectral Response Function ix ACKNOWLEDGEMENT The author would like to thank his supervisor Bernard Laval for his insight, encouragement, and understanding. Support and guidance was provided from Phil Austin in the UBC EOS department, Mike Davenport at MacDonald Dettwiler and Associates, John Morrison and Eddie Carmack at IOS, Geoff Schladow and Todd Steissberg at UC Davis, and Simon Hook at JPL. Financial support was provided through an NSERC IPS grant in conjunction with my sponsor company, MacDonald Dettwiler and Associates, for which I am very grateful. Thanks to my mom, Marilyn Sentlinger, for help with the young ones and her unswerving support. Thanks Ryan North, Jeff Carpenter, Waukene Sentlinger, and Michael Sentlinger for help with field work. Dan Pott's thesis on Quesnel Lake was an invaluable source of information and direction. My appreciation for Cheryl's love and support, and the patience of my kids Mia, Eve, and Noah, is as deep as Quesnel Lake (deepest fjord lake in the world) and as strong as the winds that rock it. x CO-AUTHORSHIP STATEMENT SUB-PIXEL WATER TEMPERATURE ESTIMATION FROM THERMAL-INFRARED IMAGERY USING VECTORIZED LAKE FEATURES Gabe Sentlinger: 80% of the identification and design of research program. 90% of performance of research. 90% of data analysis. 90% of the manuscript preparation. Dr. Bernard Laval: 20% of the the identification and design of research program. 10% of performance of research. 10% of data analysis. 90% of the manuscript preparation. CHARACTERISING LAKE DYNAMICS IN A DEEP FJORD LAKE USING SUB-PIXEL TEMPERATURES DERIVED FROM MODIS IMAGERY Gabe Sentlinger: 80% of the identification and design of research program. 90% of performance of research. 90% of data analysis. 90% of the manuscript preparation. Dr. Bernard Laval: 20% of the the identification and design of research program. 10% of performance of research. 10% of data analysis. 90% of the manuscript preparation. NOTE: Neither of these manuscripts have been submitted yet, although we intend to. xi INTRODUCTION Water temperature is an important indicator of water quality and water body dynamics for fisheries, forestry, and climate change studies. Anadromous fish are cold-blooded which makes their metabolism directly regulated by external water temperature: increased predation and mortality can result from extremely warm water, or extreme temperature gradients (McCullough 1999) . Aggressive forestry harvesting regimes can modify the runoff patterns and thereby the temperature of the water entering lakes from streams and rivers, which also serve as spawning grounds for anadromous fish. Some of these watersheds are fed by glaciers in the late summer, which are believed to be subject to accelerated recession with global warming. These factors make accurate water temperature estimates, i.e. < 1C, extremely important from a short-term, and long-term management policy perspective. Quesnel Lake (see Figure 1-1) is a lake in the interior of BC where such temperature variations have occurred with devastating impact on spawning salmon. Quesnel Lake is 81 km long but less than 1km wide in the reach leading to the mouth of the Quesnel River. The lake is an interesting subject as it is the deepest fjord lake (formed by glaciers) in North America at 505m and home to 30% of the Fraser River's Sockeye salmon run. The western arm is much shallower (max depth 100m) and separated from the deep Main Basin by a sill that rises to 30m below the surface. 1 Figure 1-1 Landsat image of Quesnel Lake and area. The West Basin is defined as the narrow reach west of Cariboo Island and feeding the Quesnel River. The Main Basin extends from Cariboo Island to the Junction between the North Arm and the East Arm. The lake has an East-West extent of 81 km and a North-South extent of 36 km, although it is only 3km across at its widest point near the Junction. Extreme temperature fluctuations (11°C in the West Basin, 5m depth, on Augl9th, 2003 (Morrison 2004)) are believed to be the result of exchange flows over this sill after large wind events. How these events arise and propagate are still largely a mystery to researchers however, as the use of in-situ sensors, such as temperature loggers, to monitor these changes in temperature are insufficient for spatially distributed phenomena. In this case, the use of remote sensing in the thermal infrared domain (TIR) is proposed as a more effective tool for measurement. TIR imagery can only take snapshots in time however, either too infrequently to capture phenomena or too coarse to spatially resolve narrow reaches. The work presented in this thesis attempts to reconcile these issues to produce an adequately high spatial and temporal resolution image series of inland water bodies, such 2 as Quesnel Lake, in order to offer some insight into what may be causing such devastating temperature fluctuations. Thermal infra-red imagery from space-borne sensors (ASTER, MODIS) and air-borne sensors (MASTER, SEBASS) has been used in previous research to determine water skin (top lOOum) temperatures of both the ocean (Mcmillin et al. 1984; Emery 1997) and fresh water bodies (Szymanski et al. 1999; Kay et al. 2001; Steissberg et al. 2003). ASTER imagery is high resolution (hi-res) low revisit time (90m pixels, 16 days), while MODIS is low resolution (lo-res) high revisit time (1km pixels, 6 hours). Generally, determining fresh water temperature is a more difficult task than ocean temperatures due to spectral mixing of land-based materials and water within a pixel. This is the case for the narrow West Basin of Quesnel Lake, which happens to be the location of greatest temperature variance, salmon activity, and likewise interest for researchers. Kay et al. (2004) have shown that this mixing effect adds unacceptable error to temperature estimates if ignored, Szymanski et al. (1999) and Kay et al. (2004) have shown that improved subpixel temperature estimates are possible by applying spectral mixture analysis (SMA) techniques to solve the mixed pixel radiance equation as set out by Gillespie (1992): =Uf,eMB(A,Tl) + f2e,2B(A,T2)] + LAu+n/i (l-l) where X is the wavelength, L^s is the at-sensor radiance, T\ is the transmissivity of the atmosphere, is the emissivity of the jth material, B(X,T) is the blackbody radiation at the jth material's temperature, f is the fractional abundance of the jth material within a pixel, the [bracketed] term is the effective ground leaving radiance, L j ^ is the upwelling radiance, and n*. is the grouped reflected radiance, adjacency effect and noise of the pixel. A schematic illustration of this effect for a single pixel, two material case is shown in 3 Figure 1-2. It is important to make a distinction between mixed pixels and pure pixels at this time. For our purposes, a mixed pixel is one which contains a significant amount of water, specifically lake water; a pure pixel is one which is entirely land-based or water. For a complete discussion of choice of materials for the unmixing model, see section 1.2.1. Figure 1-2 Schematic of pixel mixing: two material model. Radiation is emitted (and reflected) by both materials on the surface as a function of their temperature and emissivity. It is then attenuated and re-emitted by the atmosphere. This modified composite signal is what the sensor measures, and adds it's own noise to. The at-sensor radiance equation is difficult to solve, as it is under-constrained given only one channel measurement of radiance, and an unknown material emissivity, temperature, and atmospheric effect. Recent advances in the above mentioned multispectral and hyperspectral instruments contribute several more channels to help solve this problem. To solve this equation for T, and T2 we need four components: • at-sensor radiance, L ^ ; • emissivity estimates for endmembers, ex, and ex,; • atmospheric components Tx, L x u , and n; and • fractional abundance of endmembers, f, and f2. 4 The at-sensor radiance is provided in the standard MODIS product MOD02. In order to find the emissivity of a material, it is necessary to separate the emissivity from the temperature terms in equation 1-1, this is commonly known as Temperature-Emissivity Separation (TES). Most TES algorithms assume a constant and homogenous emitting surface within a pixel, which is rarely the case. More likely the surface contained within a pixel is made of heterogeneous materials and at different temperatures. Determining the sub-pixel abundance of each material in a pixel is the goal of Spectral Mixture Analysis (SMA). Several validated (TES)/atmospheric correction algorithms exist already for pure pixels (Gillespie et al. 1998; Young et al. 2002), our approach is to use the results of the MODIS TES algorithm (Wan et al. 1996) as distributed in the standard MODIS product MOD11 Temperature Emissivity for Channels 31 and 32. The emissivity estimates are taken from the pure pixel values in this dataset. There is currently no atmospherically corrected MODIS product, but atmospheric parameters can be estimated from all pure pixel estimates of T and e in the MOD11 product. Several alternative approaches to determining pixel fraction were tried, and a vectorized water feature approach was used. This procedure has been dubbed Sub-Pixel Temperature Retrieval, or SPTR. Once the sub-pixel temperature estimates are generated, there are many issues still to be addressed in order to properly interpret the results from SPTR: upwelling, skin effects, wind cooling, clouds contamination. Work by Hook et al. (2002) suggest that the water skin (top 100(im) temperature varies by +-0.5 degrees on average from the bulk water temperature (5cm down), but can swing much further depending on heating and cooling effects. The final chapter of this thesis looks at a series of images from the summer of 2003, along with concurrent weather and thermistor chain data, in order to understand 5 significant upwelling events in the West Basin, as well as the general lake thermodynamic patterns. 1.1 Fundamentals of Thermal Infrared Imagery Analysis TIR imagery differs from Visible and Near Infrared (VNIR) imagery primarily in the origin of the radiation. VNIR imagery is made possible by reflected solar radiation, while TIR imagery measures primarily the amount of radiation generated from the material being sensed. TIR imagery is necessarily much lower intensity that VNIR, due to the shape of the Planck function, and in order to obtain adequate Signal to Noise Ratios (SNR), larger pixel sizes are used, which also allows for larger coverage. This section is an overview of TIR imagery and the issues involved. 1.1.1 Remote Sensing by Satellite Remote sensing is the science of deriving information about something without having a sensor actually at the location (or in-situ). In the case of satellite TIR imagery, we are counting the number of photons that are incident on a charged-coupled-device (ccd) array (the sensor or instrument) on a piece of metal orbiting earth (the satellite or platform). There can be many instruments on board a platform, this is the platforms payload. In the case NASA's TERRA platform, both ASTER and MODIS are on board, along with several other sensors, which allows concurrent imaging of locations with both sensors. Conversely, the MODIS sensor is onboard both the TERRA and AQUA platforms. This reduces the MODIS revisit time significantly. 6 The revisit time is the time between consecutive images of the same location on earth. For MODIS with both sensors the revisit time is approximately 6 hours, although this varies with latitude as the orbits converge as they pass over the poles. Both AQUA and TERRA are on sun-synchronous orbits, as illustrated in Figure 1-3, with AQUA following approximately 1 hour behind TERRA. This time scale is ideal for imaging thermal upwelling events in lakes which rarely move faster than lkm/hour. ASTER is a much higher resolution instrument, 90m in the TIR, than MODIS, 1km in the TIR. Sensor characteristics are summarized in Table 1-1. Data collection and acquisition issues must also be addressed in any study. ASTER imagery is much more difficult to obtain, as it does not always acquire images as it passes over a location. ASTER is a tasking instrument, meaning a request is put into a priority queue to image a location, and ASTER saves images as resources allow. This is not only a data storage and bandwidth issue, as ASTER must be pointed an a location to image it, which involves moving parts, and wear issues are always a concern for the ASTER team. MODIS, on the other hand is much more freely available. It always acquires and save images, and all of these images are downloadable off the internet at no charge. This offers a wealth of information, dating back to 1999, if the tools are accessible for accessing this information. The purpose of this research is to tap into that largely unused resource for limnological study purposes. 7 Figure 1-3 Schematic of satellite in sun-synchronous orbit. The TERRA and AQUA satellites' orbital plane transects the earths poles and the sun in a "sun-synchronous" orbit. They also are arranged to always cross the equator at the same local time. For TERRA (ASTER and MODIS), this is 10:30 AM (ascending orbit) and 10:30 PM (going down). AQUA follows approximately 1 hour behind TERRA. The revisit time is a function of the swath width, or the track of earth the sensor images as it orbits: the larger the swath, the shorter the revisit time. MODIS large swath allows it to image the same spot on the earth approximately every 12 hours, and more frequent at higher latitudes as the swaths converge. ASTER's smaller swath requires 16 days to return. Because MODIS is on AQUA as well, the revisit time is reduced to approximately 6 hours or more accurately 4 times/ day. The viewing angle on half of these images is not ideal and clouds affect many of the images. Table 1-1: Sensor specification for MODIS and ASTER sensors (TIR = Thermal Infrared, SWIR= Short Wave Infrared, VNIR = Visible and Near Infrared (NASA 2005). MODIS A S T E R Launched 1999 1999 Revisit 6 hours 16 days Swath 2330 km 60 km Bands 1-2 (VNIR) 1-3 (VNIR) resolution 250m 15m 3-7 (VNIR-SWIR) 4-9 (swiR) 500m 30m 8-36 (VNIR-TIR) 10-14 (TIR) 1000m 90m Error ±1.0 K ±1.5 K 8 1.1.2 Planck's Law of Blackbody Radiation Planck postulated that matter emits radiation based on the cumulative effect of a multitude of quantum oscillators. Given all possible states of the oscillators, a continuous emission spectra is observed, with the peak corresponding to the temperature of the emitter. The energy emitted over a range of wavelengths is (Liou 2002): where B is the energy flux density in Wm 2. The energy flux density per unit wavelength is given by: (1-2) o 2hc2 (1-3) Xs (e he IKXT -1) where c is the speed of light in a vacuum=2.998 x 108 ms"1, X is the wavelength in m, T is the temperature in Kelvin, h is the Planck constant=6.626 x 10 W J s, and K is the Boltzmarm constants 1.3806 x 1023 Jdeg1. Given that (1-4) then dBdX dX dv (1-5) 9 where C,= hc/K, and v is the wave number (Liou 2002). These equations are for a blackbody, i.e.: matter that absorbs and re-emits all incident radiation. Figure 1-4 shows radiance as a function of wavelength for a number of temperatures, and Figure 1-5 plots radiance as a function of temperature for several wavelengths. Theses figures illustrate a few basic physical principals: • all matter emits some amount of radiation at all wavelengths; • the amount of radiation emitted is roughly linear with respect to temperature; and • the intensity of the emitted thermal radiation is strictly a property of the temperature, that is, no lines cross each other. A distinct line exists for each temperature, the wavelength of maximum radiance moving towards shorter wavelengths as temperatures increase. This means that for a blackbody each radiance spectrum corresponds to a single temperature and it is possible to solve for that temperature given a radiance measure in only one channel. The challenge is introduced by two factors: the atmosphere which corrupts the signal and the fact that most objects are not blackbodies. The actual radiance at a given wavelength is: Lx{T) = exBx(T) (1-6) where e is the emissivity at a given wavelength. This makes it impossible to distinguish between a greybody radiance at a given wavelength, an atmospherically corrupt signal, or a blackbody emitter at a lower temperature. For this reason, atmospheric correction must be performed to remove atmospheric effects, several bands are required, as well as some preliminary knowledge of an object's emissivity spectrum. 10 25 ,,350 Wavelength fum) Figure 1-4 Radiance Changes as a Function of Wavelength The radiance increases, and the peak moves to shorter wavelengths, with increasing temperature. The blackbody curves never cross; at a given wavelength, the radiance is only a function of temperature. 0.7 r Ground Temperature (K) Figure 1-5 The Increase in Radiance as a Function of Temperature The local linearity of the Planck function is often used to remove wavelength dependence of radiance values at different temperatures and wavelengths. 11 1.1.3 Kirchoff's Law Based on the 2nd law of thermodynamics, Kirchoff noted that absorbtance equals emissivity: (1-7) where Ex is the emittance at wavelength X, and ax is the absorbtance at wavelength X. This means that if a body is at thermodynamic equilibrium, it is emitting as much energy as it is absorbing, and the temperature is constant. A corollary to this law is the conservation of energy which states that: (1-8) where T\ is the transmittance, and rxis the reflectance. Figure 1-6 shows the interaction of radiation with matter. In our case, the grey block would the atmosphere, or a layer in the atmosphere, and we are interested in removing the emittance and transmittance effects to recover the incident radiance. emitted=8 inc ident transmit ted=T ref lected = r Figure 1-6. The interaction of light with matter. The sum of the energy absorbed, reflected, and transmitted must equal the energy incident. 12 1.1.4 Radiative Transfer The purpose of a radiative transfer equation is to distinguish between the effects of atmospheric attenuation and emission in the remotely sensed signal from the ground leaving signal. The amount of attenuation of the signal, and the addition of atmospherically emitted radiation, is given by: dLA = -Lxdu + Bxdu du where, Lx is the radiance entering the layer, du is the differential optical depth of the layer, and Bxis the radiation emitted by the layer. The optical depth of the atmosphere is defined as: dux = -pkxds r i (1-10) ux(s\,s) = | pkxds where p is the density of the gas in kg/ m', k\ is the absorption coefficient in m2/kg, and ds is the path of radiation in metres. Equation 1-10 is referred to as the Schwarzchild's equation, and its solution is given by: Lx(Si) = Lx(0)e-"^0) + lBx[T(s)-\e-u^s)dux(sx,s) ( M l ) The radiance reaching the sensor is the sum of two components (Figure 1-7): • the ground leaving radiance, attenuated exponentially from the ground, s=0, to the sensor height at s,; and • the contribution of the atmosphere, which radiates as a blackbody at the temperature of the layer T(s), and is attenuated as it passes upward through the above layers from height s to s r 13 A + • + ' + ' + * + = radiance L a y e r m Ground . . . " • V V i .' -s i 1 1 : i L a y e r 1 Figure 1-7 Atmospheric attenuation of a signal. This is a schematic of radiation leaving the material of interest on the ground Iy being attenuated as it passes through the layers of the atmosphere, each with there own optical properties and temperatures, which in turn radiate in all directions. The radiation leaving layer m is the sum of the attenuated radiation from the ground and all lower layers, plus that emitted by layer m. The radiance measured at the sensor at the Top of the Atmosphere (TOA) is the sum of all of these contributions. The total amount of attenuation, e~ , is also referred to as the transmittance or transmissivity, \ \ . Neglecting the reflectivity of TIR radiation, the total absorption coefficient, ax, is (1-Tx). From Kirchoff's law, eq 1-6, the emissivity is then (1-Xx) also. For the special case where the atmospheric temperature is constant, equation 1-11 becomes: Z,( 5 , ) = Z / l ( 0 ) e ^ ( v " ( , ) + 5 / l ( l - ^ , v l ' 0 ) ) (1-12) in which case the second term is the emissivity of the entire atmospheric column. These equations only consider the radiation traveling nadir, or straight up; to consider radiation coming from a direction off-nadir these equations are multiplied by cos 9 (u) to yield: 14 Z , ( 5 1 ) = I,(0)e-"'(''-0) + ( , 5 , [ r ( S ) ] e ^ ( , , ' I ) ^ / l ( 5 1 , j ) (1-13) Integrating over all u will give the total radiation reaching the sensor. As seen from equation 1-10, determining uj.(s,sl) requires summing the k^  for all layers and all wavelengths. There are two main approaches to this problem: Line-By-Line integration and Correlated -K Distribution. Line-By-Line adds all of the contributions in a computationally intensive algorithm, while correlated-K takes advantage of the correlation between K values in different layers, dividing all K values into percentile weights to account for the majority of the information in a computationally relaxed routine. An alternative to either of these approaches is not to calculate the atmospheric affects, but simply exploit the local linearity of the Planck function seen in Figure 1-5, and solve for temperature using a split-window technique. .5 Split Window Technique The split window technique is a means of recovering the surface temperature from the at-sensor radiance using two measurements of the same atmospheric column/ ground sample in two different channels (Mcmillin et al. 1984). Each channel will have different k-coefficients, and therefore different optical depths, but the same atmospheric and surface temperatures. If two channels, centered on A,=l and X=2, are used to make two measurements of L \ the at-sensor radiance from equation 1-12 we have: LL=BI(TM) = B](Ts)Tt+BL(TA)(\-t[) L2=B2(TB2) = B2(TX)T2+B2(TA)(\-T2) 15 where Bx is the blackbody temperature at a given temperature, TBx is the effective temperature of the atmospheric column/ground sample, T s is the temperature of the surface, TA is the effective temperature of the atmosphere, and Tx is the transmissivity of the channel = e"u<0,s1). So we have two equations and three unknowns: T s , T A , and Tx We can make some justifiable assumptions to solve the equations, however. The first is that T s , TA, TB 1, and TB, are close enough to assume the curve SBx/ 8T is linear about TA, giving: B,(T) = BA(TA) + ^ (T-Ta) (1-15) We can then write each equation inl-14 in terms of this linear approximation, and eliminate (T-TA) between them to give: B2(T)^BI(TA) + rdB2/dT^ y dB] IdT j [BX{T)-BX{TA)\ (1-16) Which relates a change in B2 to a change in B r Substituting 1-16 back into 1-14 gives: B](TB2) = BI(TS)T2+B](TA)(\-T2) (1-17) and eliminating BjfTJ gives the split-window equation: Bl(Ts) = Bi(TBi) + r1[Bi(Tm)-Bi(TB2)] (1-18) where: 1 - T , 7 7 - 1 r i T2 (1.19) 16 Equation 1-19 can be used if the transmissivities are known for the atmospheric window in question, but usually the atmosphere is modeled using several different amounts of column water vapor to find an estimate of r\. Rearranging equation 1-18 results in: , B,{Ts)-B,(Tm) B\(TB\)-B\(TBl) ( 1 . 2 0 ) where T s is known for the simulation, and TB 1 and T B 2 are modeled. This estimate of r\ is then applied to the measured brightness in each of the channels to determine the surface brightness, which can then be used to find T5 if the emissivity is known. 1.2 Spectral Mixture Analysis Spectral Mixture Analysis (SMA) is based on the assumption that we can model the radiation leaving a surface as a linear combination of a finite number of radiating, or reflecting, materials. The goal of SMA is to determine the fractional abundance of each of those materials. This is a simpler problem in the VNIR, where terrestrial materials emit very little energy and act primarily as reflectors of solar energy. In the TIR, a material's radiation is a function of it's emissivity in a given wavelength and the material's temperature. There is no direct solution to this problem, so many techniques employing different assumptions have been employed to get at the sub-pixel temperature. 17 1.2.1 Temperature End-members The first step in SMA is to discriminate between the different elements in scene. The degree of discrimination required depends on the application: for determining water temperature spatial gradients, it is not necessary to determine the percent of each pixel in the scene that is fir or pine, or young pine and old, but we do require a high degree of accuracy in determining sub-pixel abundances of water at various temperatures. Alternatively, for a Mountain Pine Beetle assessment, water temperature would not be critical, whereas pine tree age would be. This simple example illustrates the challenged associated with choosing classes, or end-members, when unmixing a scene. The number of end-members that can be unmixed with an acceptable degree of accuracy is a function of the signal to noise ratio for the scene and how many bands are available. Strictly speaking, the number of bands is the maximum number of end-members possible, although generally much better results are obtained with fewer end-members. Before we proceed to thermal unmixing an image, we must first decide the desired output or unmixing end-members. Figure 1-8 illustrates two unmixing situations we are considering: 1. end-members with various emissivities, but roughly homogenous temperatures within the image, and 2. pixels filled with a homogenous emissivity surface, but varying temperature. The ultimate goal of this project is to unmix a pixel of varying emissivities and temperatures. The problem would still be to determine what an appropriate output of such an algorithm would be. The outputs for each possible case might be: 18 1. an abundance image with an image for each end-member's abundance as well as each end-members temperature within a pixel. Each pixel value in the temperature image would indicate the average temperature of the portion of the pixel filled by that material; 2. an image which is arbitrarily larger than the original in order to effectively increase the resolution, where each pixel value represents the average temperature of that hi-res pixel; or 3. a hybrid of 1 and 2 where there is an end-member for each temperature of the homogenous material (or the material of interest). This image would have the same dimensions as 1 with the extra temperature end-members tacked on. The choice of unmixing paradigm would depend on the application. Some possible scenarios are: 1. if an accurate measurement of a spatially contiguous shape is required, it may be appropriate to simply assign a temperature to each emissivity value, or unmix based on emissivity alone and output a single abundance image. 2. If the purpose to determine thermal anomalies within a homogenous emissivity surface, such as finding landmines buried beneath a field, an unmixing of the scene based on known temperatures of expected endmembers would potentially identify thermal anomalies. 3. If a higher resolution of a homogenous emissive surface is required, such as a lake or tarmac, the higher resolution temperature image may be appropriate. In this paper we use a combination of output 1 for shore pixels and output 2 for pure water pixels. 19 pixel composed of discrete materials at various temperatures pixel composed of same material at continuously varying temperatures WT1 WT2 WT3 WT4 %Water @ W T 4 % Water @ WT3 %Rock @ RT % Forest @ FT Figure 1-8 Two Unmixing Situations, Three Possible Unmixing Products The upper left pixel contains discrete materials, forest, rock, and water, each having a semi-discrete boundary; the bottom pixel contains only water. The top pixel would be better unmixed using configuration (1), an abundance image and a corresponding temperature image. The bottom pixel would be better suited to unmixing to a higher resolution, (2) where T l , T2, T3, and T4, are the average temperature of the sub-pixels. A hybrid solution (3) is also possible, shown in the bottom right. 1.2.2 S M A in the VNIR Thermal unmixing is a different problem than unmixing in the visible domain. For mixing in the VNIR domain, very little energy is emitted from the surface; most ground leaving radiation is reflected solar or sky light. The at-sensor radiance can be approximated (Adams et al. 1989) by the following equation: ^ = ^ [ ^ 1 ( 1 - 2 1 ) where T*, is the transmissivity of the atmosphere, r^ is the reflectivity, is the incident radiation, mostly solar, [fjJLjJ, is the effective ground leaving radiance, L^U is the 20 upwelling radiance, and n is the grouped emitted radiance, adjacency effect, and noise of the pixel. If we equate Equation 1-21 with a mixed pixel model (Adams et al. 1989), we obtain: firM+/2r2 =L>A where i. is the sub-pixel abundance of the jth end-member, and is the reflectivity of the jth end-member at a given wavelength. Generalizing equation 1-22 for many end-members leads to the standard mixing model matrix equation: r = y r f ' . p " ' (1-23) r = M f where r is a vector with each element representing the measured reflectivity in a single band, M is an array of mean end-member radiances. M has as many rows as there are bands, and each column represents one end-member, and f is a column vector of that pixel's abundances. The standard solution to this equation is : f = ( M T M ) - , M T r (1-24) where the pseudo inverse of M , (IV^M)"1 is used to account for a non-invertible matrix. 21 .3 SMA in the TIR Domain If we follow the VNIR SMA we can arrive at similar equation for determining the sub-pixel abundance of emitting surfaces. The single material at-sensor radiance is given by (Gillespie 1992): LAs =TA[£AB(A,T)]e +LAu+n (1-25) where T\ is the transmissivity of the atmosphere, ex is the emissivity, B(X,T) is the blackbody radiation at a given temperature, [exB(?i,T)](, is the effective ground leaving radiance, Lx,, is the upwelling radiance, and n is the grouped reflected radiance, adjacency effect and noise of the pixel. Al l variables in equation 1-25 are at a certain wavelength or integrated over a band for a given sensor. Equation 2-1 assumes e and T are homogenous for a pixel. The more accurate mixed pixel representation would be: L,s=r,[fleMB(A,Ti) + f2£,2B(A,T2)] + LMl +n 1-26 where f is the abundance of the jth end-member, and j=[l,2 ] for a simple two material case. If we were to equate the mixed pixel model with the homogenous pixel model, we would arrive at (Gillespie 1992): L/if„5(A,7;)+/ 2^ 25(A,r 2)] = [^5(A,r)] e where the effective ground-leaving radiance is the sum of the individual material ground leaving radiances. However, it is not correct to unmix a pixel using the homogenous pixel model of emissivity and temperature products because, as seen in Equation 2-5: 22 f\f\ ~*~ f'J-l ^ (1-28) We must take a different approach to the problem. One approach is to assume a characteristic temperature for each material and unmix either the ground leaving radiance or the at-sensor radiance. Analogous to equation 1-23, the matrix unmixing equation for a TIR image would be: 7^ (1-29) L = L f g m where L g is a column vector of ground leaving radiance with as many rows as bands, L m is a column vector of mean end-member effective ground-leaving radiances, (eiB(A,iT)), and f is an vector of each end-members abundance. A similar equation can be used to unmix the at-sensor radiance image as all of the L m vectors undergo the same gain and offset - transmission losses and additive atmospheric offset - leaving a linear system. If the materials are chosen from the at-sensor radiance image, the model takes the form: Lb =(Tx£AMfTl) + L^+n)fl+(t,enB(A,T2) + Llu+n)f2 (1-30) = TXEMB{A,TM +TA£Z2B(A,T2)f2+(LXM +n)(fl +/2) where f,+f2 sum to unity and therefore are in the same form as Equation 2-4, and can therefore be unmixed using a linear equation. The benefit to unmixing an at-sensor radiance image is that the values seen in each pixel are truly an additive mix of photons. Whereas unmixing an atmospherically corrected image will unavoidably contain artifacts from the correction and will introduce more error into the unmixing results. 23 In this paper TIR unmixing is used in an autoregistration routine, but not to solve the mixed-pixel ground-leaving radiance equation. This is due to the limited number of MODIS bands available, which restricts both the number of materials that can be solved for and increases the error in the unmixing results. Instead, an unmixing algorithm is presented which retrieves material fractional abundance from vectorized material features. This approach is only applicable for materials whose extent is known a priori, such as water and land. 1.3 Structure and Intent of Thesis The remainder of the thesis is divided into two chapters: one which addresses the issues outlined in this chapter and presents a solution, and the other which then applies this solution to a test case. The intention of the thesis is to prove that accurate sub-pixel temperature retrieval is possible using vectorized water features to determine fractional abundance of water and land within a pixel. We also hypothesize that standard temperature and emissivity products derived from accurate split-window algorithms can be reversed to extract atmospheric components, and solve the mixed-pixel ground-leaving radiance equation. The solution is simple, effective, and useful to a wide range of scientists with or without remote sensing background, and also applies to studies involving any TIR sensor, so long as there is vector data delineating the materials of interest, and an available temperature and emissivity dataset available. The solution is very economical, all datasets are free of charge and freely available on the internet, and has been packaged into a Graphical User Interface (GUI) in the popular remote sensing software platform, ENVI. The solution is presented in the next chapter. 24 1.4 References Adams, J. B., M. O. Smith and A. R. Gillespie, 1989. "Simple models for complex natural surfaces: a strategy for the hyperspectral era of remote sensing." IGARSS. Emery, W. J., and Yu Y., 1997. "Satellite sea surface temperature patterns." International Journal of Remote Sensing 18(2): 323-334. Gillespie, A., S. Rokugawa, T. Matsunaga, J. S. Cothern, S. Hook and A. B. Kahle, 1998. "A temperature and emissivity separation algorithm for Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images." Geoscience and Remote Sensing, IEEE Transactions on 36(4): 1113-1126. Gillespie, A. R., 1992. "Spectral mixture analysis of multispectral thermal infrared images." Remote Sensing of Environment 42(2): 137-145. Hook, S., F. J. Prata and S. G. Schladow, 2002. "Final Report: ASTER and MODIS Calibration on Lake Tahoe." http://shookweb.jpl.nasa.gov/validation. Kay, J., R. N. Handcock, A. Gillespie, C. Konrad, S. Burges, N. Naveh and D. Booth, 2001. "Stream-temperature estimation from thermal infrared images". Geoscience and Remote Sensing Symposium, 2001. IGARSS '01. IEEE 2001 International. Kay, J. E., R. N. Handcock, K. A. Cherkauer, S. K. Kampf, A. R. Gillespie and S.J. Burges, 2004. "Accuracy of lake and stream temperatures determined from atmospherically corrected thermal-infrared imagery,." Journal American Water Resource Association submitted. Liou, K.-N., 2002. An introduction to atmospheric radiation. Amsterdam ; Boston, Academic Press. McCullough, D. A., 1999. A review and synthesis of effects of alterations to the water temperature regime on freshwater life stages of salmonids, with special reference to chinook salmon. Seattle, Washington, U.S. Environmental Protection Agency. Mcmillin, L. M. and D. S. Crosby, 1984. "Theory and Validation of the Multiple Window Sea-Surface Temperature Technique." Journal of Geophysical Research-Oceans 89(NC3): 3655-3661. Morrison, J., 2004. Unpublished Data, Vynx Design, Sidney BC, V8L 4B2, Canada. NASA, accessed August 2005. Land Processes Digital Active Archive Center, lpdaac.usgs.gov/main.asp Steissberg, T. E., S. J. Hook and S. G. Schladow, 2003. "Characterizing Partial Upwellings and Surface Circulation at Lake Tahoe, CA/NV, USA with Thermal Infrared Images." submitted. 25 Szymanski, J. J., C. C. Borel, Q. O. Harberger, P. Smolarkiewicz and J. P. Theiler, 1999. "Subpixel temperature retrieval with multispectral sensors". Algorithms for Multispectral and Hyperspectral Imagery V, Orlando, FL, USA, SPIE. Wan, Z. and J. Dozier, 1996. "A generalized split-window algorithm for retrieving land-surface temperature from space." IEEE Trans, on Geoscience and Remote Sensing 34(4): 892-904. Young, S. J., B. R. Johnson and J. A. Hackwell, 2002. "An in-scene method for atmospheric compensation of thermal hyperspectral data." Journal of Geophysical Research (Atmospheres) 107(A14): ACH14-1. 26 2 SUB-PIXEL WATER TEMPERATURE ESTIMATION FROM THERMAL-INFRARED IMAGERY USING VECTORIZED L A K E FEATURES 1 2.1 Introduction Satellite imagery is valuable for recording and analyzing spatial-temporal patterns of water surface temperature, but is often limited by a sensor's ability to resolve narrow reaches of small lakes and rivers. Often a narrow reache leads to the outlet of the lake and is the most critical region for spawning salmon, which can be devastated by temperature fluctuations. The conventional approach to lake temperature measurement is spot measurements using thermistor chains. These can only report localized temperature patterns, which often do not offer insight into the dynamics of complex lake systems such as Quesnel Lake in central British Columbia, which is used as a test case in this paper. Coupled with satellite thermal infrared imagery, however, the two data sources are able to paint a more complete picture of the lake dynamics. Before this can happen, the narrow reaches of the lake must be resolved and with enough satellite images to create a 1 This paper is yet to be submitted for publishing. 27 meaningful time series. To date, no satellite sensor provides this product. The current research attempts to rectify this deficiency by extracting sub-pixel temperatures using vectorized water features and spectral unmixing techniques. This results in water temperature estimates which recover the high spatial resolution information and offer insight into the temporal and spatial distribution of water temperature. Thermal infra-red imagery from space-borne sensors (ASTER, MODIS) and air-borne sensors (MASTER, SEBASS) has been used in previous research to determine water skin (top 100(im) temperatures of both the ocean (Mcmillin et al. 1984; Emery 1997) and fresh water bodies (Szymanski et al. 1999; Kay et al. 2001; Steissberg et al. 2003). Generally, determining fresh water temperature is a more difficult task than ocean temperatures due to spectral mixing of land-based materials and water within a pixel. Kay et al. (2004) have shown that this mixing effect adds unacceptable error to temperature estimates if ignored, while Szymanski et al. (1999) and Kay et al. (2004) have shown that improved subpixel temperature estimates are possible by applying spectral mixture analysis techniques to solve the mixed pixel radiance equation as set out by Gillespie (1992). The at-sensor radiance equation is difficult to solve, as it is under-constrained and non-linear; 10 channel emissivities + 1 temperature are unknown with only 10 channels given, spectral mixing of different materials involves addition of the non-linear Planck function. Recent advances in the above mentioned multispectral and hyperspectral instruments and the accompanying algorithms contribute several channels to help solve this problem. The problem of temperature-emissivity separation (TES) has been addressed successfully in several algorithms (Gillespie et al. 1998; Young et al. 2002). Others avoid the task of TES (Wan et al. 1996) and assume the Planck function is locally linear in order to perform a 28 split-window technique together with emissivity lookup tables to extract temperature estimates. Both of these approaches only solve for the pure-pixel case. To solve the mixed pixel case, an estimate of the sub-pixel material fraction is necessary. Conventional approaches such as linear unmixing (Adams et al. 1989) break down due to non-linear effects, and other approaches involving classification of higher-resolution VNIR data (Hook 2004) only applies to daytime imagery. Our approach synthesizes the high-quality pure-pixel information publicly available with novel vector unmixing techniques to address these issues and unlock a wealth of archived low spatial resolution (lo-res) satellite imagery for use in limnological studies of narrow lakes and rivers. Given that several validated temperature emissivity separation (TES) and atmospheric correction algorithms already exist for pure pixels, we use as input pure-pixel emissivity and temperature values supplied by these algorithms as anchors for the derived sub-pixel temperature information. Using MODIS data, which has a short revisit time but low resolution —approximately 6 hours and 1km2 respectively —we generate a higher resolution temperature image series, both along the shore and for pure water pixels, for Quesnel Lake. The sub-pixel corrected imagery is compared with concurrent ASTER data and in situ temperature data for validation. 2.2 Theoretical Background The pixel value in any remotely sensed image is the radiance leaving the entire column beneath it. This includes the additive contribution from the atmosphere, and the attenuation of the radiance leaving the surface below. In order to recover subpixel 29 temperatures in lakes and rivers, we must solve the at-sensor radiance equation. The single material at-sensor radiance is given by: L^=ri[eiB(A,T)l+LAu + Ti(\-el)L^+nl . (2-1) where Xx is the transmissivity of the atmosphere, Ex is the emissivity, B(X,T) is the blackbody radiation at a given temperature, [ExB(?i,T)]i, is the effective ground leaving radiance, Lx^ is the upwelling radiance, and Lx d is the downwelling radiance, and nx is the grouped adjacency effect and noise in the pixel. All spectral variables in equation 2-1 represent mean values over a sensor's channel width. Although there are additional complicating factors in this equation, such as the sensor's point spread function (PSF) and spectral response function (SRF), we adopt this equation initially for simplicity. As we are generally dealing with near blackbody materials, water and vegetation, the reflected term can be neglected. We further simplify the problem by defining a single additive component: Lu=Liu+Tx{\-E,)LXd+ nx (2-2) to give: Lh=Tx[exB{Xj)\+LXA (2-3) The effective-ground leaving radiance term is set in brackets here because it is useful to define such a quantity, though physically unrealistic. Equation 2-3 assumes £ and T are homogenous for a pixel although the radiance leaving the surface is always a mixture of materials at different temperatures and with different emissivity spectra. The degree of generalization of the pixel components, or endmembers, depends on the application and spectral resources available. For MODIS imagery, there are 15 channels in the TIR, but only 5 that are in atmospheric windows. Of these 5, the MODIS product MOD11 only 30 provides two estimates of emissivity. Because we are only interested in two materials, water and not-water, this study uses MOD11 and a two material mixing model, which is a more accurate representation of at-sensor radiance than the single material case: Lh=tx[f^MB{fT,) + f28X2B{Xj2)] + LM (2-4) where f is the abundance of the jth end-member, j=[l,2 ] for a simple two material case. To solve this equation for T, and T, we need four components: • at-sensor radiance, Lx,; • fractional abundance of endmembers, f, and f2; • emissivity estimates for endmembers, Ex, and £x 2; and • atmospheric components Tx, Lx,,, and n. The at-sensor radiance is provided in the standard MODIS product MOD02. The fractional abundances can be generated by "reverse geo-referencing" vector data, that is, placing the vector into the original pixel space and determining sub-pixel fractions. The emissivity estimates are taken from the pure pixel training sets. The atmospheric parameters can be estimated from all pure pixels by fitting a line to equation 2-1. Because there is no analytical solution to equation 2-4 for T, and T,, we use a gradient descent method to minimize the error between the measured at-sensor radiance and an estimated at-sensor radiance. This procedure has been dubbed Sub-Pixel Temperature Retrieval, or SPTR. The SPTR flow diagram is shown in Figure 2-1. 31 At-Sensor Radiance Image Perform T E S T-e Image Estimate T,£ for End-members e.Tofpure pixels Retrieve Ground-Leaving Radiance (Atmo-correct) Ground-Leaving Radiance Image eofpurfy pixels Thermal Temperature Unmix Sub-pixel Temperature Image Figure 2-1 Sub-pixel Temperature Retrieval flow diagram This strategy uses validated pure-pixel temperature and emissivity products and vectorized water features as input, minimizing the load on the algorithm to produce accurate results. All inputs are readily and freely available on the internet. Dashed lines indicate optional user-supplied input. 32 .1 Fractional abundance from Vectorized feature data There are three options for determining fractional abundance: 1. perform spectral mixture analysis ( S M A ) on radiance data; 2. classify concurrent hi-res V N I R data as land/water and transform to lo-res space; or 3. determine pixel fraction from vectorized water features. S M A was tried and unsatisfactory results were obtained given the non-linear mix ing effects introduced when materials have variable temperatures. Concurrent V N I R data was used successfully, but only improved the resolution four times for M O D I S and was only available i n daytime images; half of the available images. Vectorized water feature data addresses both of these deficiencies but georeferencing needed to be addressed. G i v e n that water features rarely change position, it is reasonable to use vectorized water features to determine the fractional abundance of water. The degree of water movement , as w e l l as the positional accuracy of the vector data and raster georeferencing information, should be much smaller than the pixel we are trying to correct. In order to use this data, we must first transform the two data sets into a common space. Because m i x i n g occurs i n the original p ixel space, resampling dur ing georeferencing corrupts the mixed signal, so we cannot transform to a common geographic space. Instead, a technique had to be developed to transfer georeferenced vector data into the original low-resolution pixel space. The steps are: 1. georeference low-resolution (lo-res) radiance image into high-resolution (hi-res) space at approximately lOx the original resolution; 2. rasterize vector data into same hi-res space; 3. reverse geo-reference pixels by cross-referencing all hi-res pixels to original lo-res space, g iv ing approximately 100 hi-res pixels per lo-res pixel (Figure 2-2) 33 4. determine fractional abundance in each lo-res pixel by dividing the number of hi-res pixels of a certain endmember by the total number of hi-res pixels in a given lo-res pixel. The issues associated with this approach are georeferencing error and digitizing error in the vector data. This approach also does not take into account the sensor's PSF. In this research, we have aimed for error < 10% of pixel spacing. To accomplish this, we do the following: • use vector data with positional accuracy < 10% of pixel spacing. For 1km pixels, 1:50,000 vector data is used with a positional accuracy of 25m. • use an autoregistration routine to fine tune georeferencing of raster data. This approach uses a gradient descent algorithm to minimize the error between vector abundance data and spectral unmixing results to within 10% of a pixel spacing. The benefits of this approach are a reduction in the error of abundance estimates of stationary endmembers, especially those that vary in temperature such as lakes, rivers, tarmacs, buildings, and roads. Figure 2-2 Abundance from reverse geo-referencing vector data into original lo-res space. The rasterized lake vector data shown above is cross-referenced back into the original lo-res space using a geo-referenced radiance image. By counting the number of hi-res water pixels in the lo-res pixel, and dividing by the total number of hi-res pixels in the lo-res pixel, we find the lake abundance in pixel., to be approximately 40%. 34 2.2.2 Atmospheric Correction from pure pixels Atmospheric correction for SPTR is done using in-scene techniques and pure-pixels from abundance estimates. SPTR atmospheric correction takes as input a pre-generated temperature-emissivity (T-e) image and the at-sensor radiance image. This technique modularizes the temperature retrieval process, allowing the user to determine the best TES algorithm for his/her application or sensor. The accuracy of mixed pixel temperatures will be based on the accuracy of the input T-e image, but should retain the pure-pixel temperatures. For this study, we used the MODIS standard products: • MOD021KM MODIS/Terra Calibrated Radiances 5-Min LIB Swath 1km at-sensor radiance, radiometrically corrected; • MYD11_L2 and MODll_L2 MODIS/Aqua and MODIS/Terra Land Surface Temperature/Emissivity Daily 5-Min L2 Swath 1km for channels 31 (Hum) and 32 (12pm); and « MYD03 and MOD03 MODIS/Aqua and MODIS/Terra Geolocation Fields 5-Min L1A Swath 1km. All data products were swath data for 5 minutes of acquisition with no temporal averaging. They are available at no-charge at the EOS data gateway (NASA, 2005). A description of the products can be found at the NASA/ USGS MODIS website (NASA, 2005). Following the ISAC algorithm (Young et al. 2002), we set out to determine the atmospheric transmission and path-radiance using a line-fitting technique. Unlike Young et al. (2002) we don't try to determine blackbody materials and reference bands; we assume the T-e product input is correct (within error) for pure-pixels in a scene and solve equation 2-3 for T\ and L^A. This solution takes advantage of the very accurate temperature estimates possible from the split-window technique used in the MOD11 product which requires no 35 atmospheric correction in and of itself. The steps involved in atmospheric correction are as follows: 1. determine pure-pixels using the method outlined in section 2.2.1. 2. using the supplied T-e image, determine the ground-leaving radiance for all pure-pixels. 3. regress the ground leaving radiance against the measured at-sensor radiance; the slope of the line will be the transmissivity for a given channel, and the offset will be the combined path radiance and n from equation 2-1. 4. solve for ground-leaving radiance for all pixels using the derived transmissivity and path radiance/reflected radiance from equation 2-4. This technique assumes that these atmospheric parameters are constant for the entire scene, which has been proven to be reasonable in the ISAC algorithm (Young et al. 2002), and that noise, adjacency effect, and reflected radiation are negligible. Once X\ and LxA are found, we solve for ground leaving radiance: = [fl£uB(A,Tirlf2£2,B(A,T2)] 2-5 = L, where LxB is the ground-leaving radiance. Figure 2-3 shows a scatter plot of channel 32 estimated ground-leaving radiance and measured at-sensor radiance. 36 Figure 2-3 At-sensor Radiance vs Estimated Ground-leaving radiance, Ch32. By fitting all pure pixels in the scene to equation 2-4, we get a Tx of 0.85 and Tx and L x A 1.02 for this channel. Using all unclouded pixels, the correlation coefficient is 0.983 compared to 0.984 when using just pure unclouded pixels. The small difference in correlation is possibly because this scene is a forestry scene comprised of mostly near-black bodies. 2.2.3 Emissivity from Pure Pixels To estimate e we use a roaming net that takes the average of all adjacent pure pixels to the pixel being unmixed for water and non-water. This assumes the land coverage is homogenous within a 5x5 net. Figure 2-4 illustrates the technique. Figure 2-4 5x5 roving net to determine local emissivity. The red pixels in this images are pure shoreline pixels, blue are pure water pixels. To determine the local emissivity of the shore and water when unmixing the central pixel, we average pure pixel emissivity values that fall within this 5x5 net. 37 2.2.4 Sub-Pixel Temperature Recovery Once we have an estimate of the ground leaving radiance and fractional abundance of the water and land within each pixel, we can solve for sub-pixel temperature, T. We cannot solve for T explicitly as it is buried within the Planck Function. The Planck function is often assumed to be locally linear, which is the premise of the split-window technique. This technique, however, relies on there being a significant difference between emissivity spectra in order to reliably solve the sub-pixel temperature equation. Attempts at solving this with a linear approximation were made, but without success. A numerical method was used with much more success. The gradient descent method minimizes the error between the measured ground-leaving radiance and the sum of the estimated sub-pixel radiance. Minimizing the error between the two gives: Err2 ={LAg-[fleMB(A,Ti) + f2£,2B(A,T2)]}2 ^ = 2fcfe - \ffuB(X,Tx) + f2£A2B(A,T2)1 dB(A,Tx) dT, (2-6) and The gradient descent prescription is: dT. dt n = -v dErr2 K (2-7) 38 where r\ is the gradient descent gain, is the measured ground leaving radiance, and dTi/dt is the amount to adjust the nth temperature each iteration (t) until convergence. This typically required 5-100 iterations before convergence, as shown in Figure 2-5 . -0.OOO1 ^ f— -0.0002 - \ 1 \ t: -0.O0O3 - i w \ -D.0D04 - \ j -0.0005 • V Figure 2-5 Gradient Descent convergence Two example plots of the gradient descent algorithm converging for a two endmember case. Eq 2-7 solves for the average sub-pixel temperature for the fraction fn of the nth endmember. To further increase the resolution enhancement, we can assume the water temperature varies smoothly between lo-res pixels, and perform a modified convolution of the data. Using only the lake temperature pixels, we run a 20x20 (2kmx2km) filter over the data, then normalize by the number of lake pixels within each kernel. This preserves the average temperature for every kernel within the lake. The sub-pixel smoothed temperature is: j J (2-8) Where T n s is the hi-res sub-pixel temperature of the nth endmember after smoothing, Tn j is the is hi-res sub-pixel temperature before smoothing, K is the convolution kernel centered on T n s and j the jth element of the kernel, W is the pixel classification (lake =1, non-lake = 0). 39 2.2.5 Georegistration and Autoregistration One of the largest sources of error in this technique is georeferencing error in the MODIS image. For georeferencing, the georegistration product MOD03 was used, and initially the MRTSwath tool available from (NASA 2005). Artifacts from the nearest neighbour resampling technique in this algorithm rendered it unsuitable for this algorithm, however. In its stead, JPL's Win Vicar and SciExtract tools were used. These tools had to be modified to handle spectrally subset images, but eventually provided much better results. MOD03 inherently has georeferencing error, which adds unacceptable error to the results. To address this, an autoregistration algorithm was implemented. Autoregistration is accomplished using a gradient descent method which minimizes the error between a SMA fractional abundance image and a vector abundance image. The steps are as follows: 1. determine fractional abundance from initial georeferenced vector data as outlined in section 2.2.1; 2. use all pure pixels (f = 1) from vector abundance image to train unmixing algorithm; 3. measure the rms error between the abundance image and vector abundance image. 4. move the vector data l/10 l h of a pixel in one direction. 5. repeat steps 1-3. If error is less, continue in this direction; if error is greater, change direction. 6. repeat steps 1-5 until no movement or rotation reduces the error. Typically the initial georeferencing error was less than V i a pixel for the entire scene, but never more than 1 pixel. Error was worse for oblique angle images. In this technique, only scene-wide georeferencing error was considered, not relative error between pixels. Figure 2-6 shows that this approach can dramatically reduce the subsequent thermal unmixing error. Although autoregistration did not work in all cases, it worked > 50% of the time. This algorithm could be improved in subsequent research efforts. 40 Figure 2-6 SPTR results before and after autoregistration. The left image is the land and water image before SPTR. The middle image shows the retrieved water temperatures using the vector data without autoregistration. The right image shows the retrieved water temperatures with auto-registration; notice the offset between the vector data and adjusted temperature lake outline is 0.7 pixels down and 0.1 pixels left. 2.3 Results Using the SPTR technique outlined above, a test image was processed acquired from the multispectral sensor MODIS and compared to the hi-res ASTER T08 Temperature product. MODIS acquires in the TIR at 1km pixel spacing. Figure 2-7 shows SPTR results from an Aug 2001 image which was acquired concurrently with a hi-res TIR image from the ASTER instrument on-board the TERRA platform. The ASTER T08 temperature product is used here. In order to compare SPTR results with the ASTER image, the same reverse georeferencing lookup table was used to transform the corrected temperature image into a high resolution (100m pixels) georeferenced space. This image assigns the lo-res pixel value to all hi-res pixels that correspond to the lo-res pixel using a nearest neighbor resampling. This approach has the benefit of providing hi-res spatial information in a format that the analyst can readily interpret, discern temperature features, and identify potential 41 unmixing errors (usually caused by unflagged clouds). Another note to remember is that this format is not simply a georeferencing; the temperature value assigned to a cluster of high res pixels is the estimated average temperature of all of those hi-res pixels. Figure 2-7 shows the results. Based on Hook et al's. (2002) results, ASTER T08 temperature estimates have a cold bias of 0.4°C and MOD11 Temperature estimates have a hot bias of 0.3°C on average. In this scene, a difference of 0.9°C between average lake temperatures for MODIS and ASTER was removed by adding 0.9°C to the ASTER image, which is similar to the combined 0.7°C bias reported by Hook et al. (2002). Figure 2-8 shows a scatter plot of ASTER lake temperatures vs the original MOD11 temperature product for lake temperatures on the left, and ASTER vs SPTR on the right. This was a daytime image, so the majority of the MOD11 mixed pixels were hotter than the corresponding ASTER pixels. 99% of the MOD11 product lake temperatures are with 2.55°C of the ASTER temps, compared to 0.96°C for the SPTR results. Based on Hook's (2001) reported ASTER T error of 0.44°C (bias removed) and our 0.32°C error, the final RMS error estimate for the SPTR results is 0.54°C after the bias has been removed. This means that the error in the temperature difference between any two pixels will be less than 1.08°C , 95% of the time, or ±0.54°C, based on this proxy analysis. This implies that for a temperature difference to be significant, it must be >1.08°C, which represents a worst case-scenario given that we are referencing our error to another image product's error as opposed to a near-surface measure of radiation. Given the bias reported in Hook et al. (2002) of 0.3°C in the MODIS temperature product, MOD11 (which has probably been improved in more recent products) and practically no bias in the at-sensor radiance product, MOD02, which is also used in the SPTR algorithm, an absolute rms temperature error of ±0.62°C is expected. This is within the 1°C error originally aimed for. Lake surface radiometers would be 42 required to determine the actual relative error between pixels and the absolute error in skin temperature. Extracting the bulk temperature from the skin temperature is another matter. Sentlinger et al. (2005) describes the relationship between the 3m bulk temperature and derived SPTR results and finds a comparable 95% confidence of retrieving 3m depth bulk temperatures to within 1.0 °C after a multiple linear regression with wind speed, wind history, air temperature, cloud cover, and time of day parameters. The errors described here are shown in Table 2-1. In this table, the dataset is the data which has an error reported to the reference. The relative error is between pixels (T i x c l a - T ixclb) and represents the standard deviation about the mean error, the absolute error includes both the variance and any potential bias of the dataset against the reference. Table 2-1: Summary of error estimates. Dataset Reference Relative Abso lu t SPTR ASTER T8 0.32 ASTER T08 Skin 0.44 SPTR Skin 0.54 0.64 SPTR 3m depth 0.49 43 Figure 2-7 Sub-Pixel Temperature Retrieval Results. The upper right image shows the original MOD11 Temperature product. Notice the unacceptable contamination in the western reach and along the shore line. The SPTR results in the lower left image have greatly reduced this error, recovering the colder temperatures in the western arm which are apparent in the hi-res ASTER image in the upper left. The plot in the lower right compares transects for the lake from west to east (upper trace) and south to north (lower trace). Notice the MOD11 product is much too warm in the western arm while the corrected image tracks the ASTER temperature much better, resolving the periodic signal at locations 100. 44 Figure 2-8 Sub-Pixel Temperature Retrieval results: scatter plot. These plots show the ASTER lake temperature pixel values on the x-axis and MOD11 lake temperature pixel values on the y-axis, for corresponding pixels. The standard error here is 0.84°C. The SPTR results are shown on the right (ASTER x-axis, SPTR y-axis) and form a much tighter distribution around the ASTER values, here the standard error is 0.32°C. The solid lines show the equality relationship, the dashed lines show ±1°. 2.4 Discussion Results from this analysis appear promising, albeit there are still many issues to work out. Figure 2-10 shows the SPTR results after lake temperature smoothing. The smoothing algorithm averages over a 2km kernel and normalizes to average lake temperature, neglecting land influences. This is believed to be a more accurate representation of the smooth temperature gradients typically seen. This image is one of the more interesting images seen. Both the north and east basin are fed by cold rivers, and are flanked by tall ranges which block the sun in the late evening and early. The extension of the cold water into the Main Basin perhaps suggests these arms are dominated by the cold water fed into them, or that they receive significantly less solar radiation due to shadows. The spatially periodic signal in the West Basin is also an indicator of the very active temperature structure in this part of the lake. Whether this pattern is moving, or a standing wave, or 4 5 just localized heating or wind effect is uncertain without more images, which were unavailable due to clouds. One of the main challenges is interpreting the skin temperatures in terms of what is happening beneath the skin. Again, Hook et al. (2003) have done extensive work in monitoring surface skin response to external forcing such as wind, air temperature, and irradiance, and its correlation with bulk temperatures. On average, they found the skin temperature was 0.5°C colder than the bulk temperature measured 5cm below the surface. This varied throughout the day depending on air temperature wind, and solar radiation. Figure 2-11 again relates the top metre to temperatures 5m down, where thermistor chains often have their top thermistor. These profiles were taken using a Brancker TL1000, a watch, and a measuring tape at and around the junction of Quesnel Lake on August 13, 2004. The thermistor was lowered into the lake and held at each depth for 30seconds. The time was recorded and temperatures were averaged at each depth after removing the aclimatation period. The intention was not to record the effect of wind mixing on the top lm, which happened to result, but to get a spatial sampling of surface temperatures around the Junction of Quesnel Lake. Understanding this section of the water column is crucial if satellite TIR imagery is to be used to extrapolate to 5m and below. Modeling programs such as DYRESM do a good job of modeling 5m downwards, but require some calibration data. If a model could be built of the top 5m and the effects of diurnal heating and wind mixing, then TIR imagery could be used to calibrate such models and reduce the need for in situ instrumentation. Sentlinger et al. (2005) perform a simple multiple linear regression of this problem and achieve an estimate of the bulk temperature (3m depth) using SPTR temperature estimates and concurrent weather data to within ±0.5°C of the measured bulk temperature. 46 Much of this analysis was done using a Graphical User Interface (GUI) built in the IDL Language and using ENVI's custom remote sensing based libraries. One of the goals of the project was to make tools that a limnologist could use with limited knowledge of RS and atmospheric correction. This tool can process large datasets, 500 images have been done automatically, sorting them into clouded, oblique, or good quality, and performing SPTR with autoregistration on latter. It is hoped to translate these tools into MATLAB, which has a larger user-base in limnology than IDL. Figure 2-9 SPTR toolset with a quicktime movie. The two main tools are for performing batch SPTR on many images, and the second tool is for integrating wind and thermistor chain data into a multiband image which could be converted into a movie. In the future, the following improvements could be made to the algorithm: incorporate sensor PSF and SRF; use more MODIS channels; unmix for more than two endmembers; try alternative temperature and emissivity estimates. incorporate a skin-bulk temperature model to estimate the temperature 5m down based on time of day, cloud-cover, topography, wind data, and air-temperature. 4 7 Figure 2-10 Sub-pixel Temperature retrieval results with water smoothing This image shows the results of SPTR after a lake-normalized 1.5km x 1.5km kernel has been convolved with the initial SPTR results. This temperature transect is compared with a smoothed ASTER transect in the upper left for the west-east profile (upper) and south-north (lower). Agreement is within 0.32 °C except in the narrow passage by Cariboo island at location 200. This is likely because of registration error. This image was acquired August 9th, 2001 at 11:30 A M local time. Niagra Creek at the far end of the East Arm is glacier fed which may explain the cold water around its inlet. No local wind data is available for this time period, but prevalent wind directions have been overlaid from more recent field work along with 100m elevation contours. The relative arrow lengths indicate the relative frequency of a direction. - 1 r - — - — - 7 1 after wind during wind before wind 1 ) £ 3 a OI 5 18 18.5 19 19.5 20 20.5 21 Temperature (C) Figure 2-11 Top 5 metres of water column during wind event. This plot shows how the top 5m of the water column cools and mixes during a wind event. These profiles were taken around the central junction of Quesnel Lake on Aug 13 » 2004, just as a minor wind storm passed through. Note that these were not the same location, so heat is not conserved, but the structure of the column is representational of the heating effect at the surface approximately 2°, wind cooling causing instability, and the subsequent falling and mixing of the cooled water column. 48 2.5 Conclusions The SPTR algorithm is able to extract sub-pixel temperatures from MODIS data to +0.32°C of ASTER skin temperature, or ±0.54°C of the relative temperature between pixels. The absolute error in recovered skin temperature is expected to be ±.62°C, based on this proxy analysis and reported uncertainty in the MODIS products used in SPTR. This is comparable to the 1.0° degree accuracy required for the study of temperature effects on spawning salmon in lakes and rivers. By using a priori knowledge of lake water extent in the form of vectorized water features, SPTR iteratively solves the mixed pixel ground leaving radiance equation using a gradient descent method. Atmospheric correction is done using in-scene techniques and pure-pixel estimates of temperature and emissivity from standard MODIS products. A user-friendly GUI has been built in IDL to allow automatic processing of large numbers of images, sorting into clouded, oblique, or good quality, and allowing manual fine-tuning of the results for important images. This allows users outside of the remote sensing community to access a wealth of stored MODIS high temporal frequency datasets for limnological analysis. Sub-pixel temperatures from MODIS imagery are very important for many interior lakes which tend to be narrow and long, wedged between mountain ranges. These are often important salmon spawning grounds, and temperature estimates are crucial for resource management. 49 2.6 References Adams, J. B., M. O. Smith and A. R. Gillespie, 1989. "Simple models for complex natural surfaces: a strategy for the hyperspectral era of remote sensing." IGARSS. Emery, W. J., and Yu Y., 1997. "Satellite sea surface temperature patterns." International Journal of Remote Sensing 18(2): 323-334. Gillespie, A., S. Rokugawa, T. Matsunaga, J. S. Cothern, S. Hook and A. B. Kahle, 1998. "A temperature and emissivity separation algorithm for Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images." Geoscience and Remote Sensing, IEEE Transactions on 36(4): 1113-1126. Gillespie, A. R., 1992. "Spectral mixture analysis of multispectral thermal infrared images." Remote Sensing of Environment 42(2): 137-145. Hook, S., 2001. "Progress Report: ASTER and MODIS Calibration on Lake Tahoe." http:// shookweb.jpl.nasa.gov/validation. Hook, S., 2004. Discussion during JPL visit. Pasadena, CA. Hook, S., F. J. Prata and S. G. Schladow, 2002. "Final Report: ASTER and MODIS Calibration on Lake Tahoe." http://shookweb.jpl.nasa.gov/validation. Hook, S. J., F. J. Prata, R. E. Alley, A. Abtahi, R. C. Richards, S. G. Schladow and S. O. Palmarsson, 2003. "Retrieval of Lake Bulk and Skin Temperatures Using Along-Track Scanning Radiometer (ATSR-2) Data: A Case Study Using Lake Tahoe, California." Journal of Atmospheric and Oceanic Technology 20: 534-548. Kay, J., R. N. Handcock, A. Gillespie, C. Konrad, S. Burges, N. Naveh and D. Booth, 2001. "Stream-temperature estimation from thermal infrared images". Geoscience and Remote Sensing Symposium, 2001. IGARSS '01. IEEE 2001 International. Kay, J. E., R. N. Handcock, K. A. Cherkauer, S. K. Kampf, A. R. Gillespie and S.J. Burges, 2004. "Accuracy of lake and stream temperatures determined from atmospherically corrected thermal-infrared imagery,." Journal American Water Resource Association submitted. Mcmillin, L. M. and D. S. Crosby, 1984. "Theory and Validation of the Multiple Window Sea-Surface Temperature Technique." Journal of Geophysical Research-Oceans 89(NC3): 3655-3661. NASA, accessed August 2005. Land Processes Digital Active Archive Center, lpdaac.usgs.gov/main.asp Sentlinger, G. I., S. Hook, B. Laval and J. Morrison, 2005. "Characterising lake dynamics in a deep fjord lake using sub-pixel temperatures derived from modis imagery." To Be Submitted. 50 Steissberg, T. E., S. J. Hook and S. G. Schladow, 2003. "Characterizing Partial Upwellings and Surface Circulation at Lake Tahoe, CA/NV, USA with Thermal Infrared Images." submitted. Szymanski, J. J., C. C. Borel, Q. O. Harberger, P. Smolarkiewicz and J. P. Theiler, 1999. "Subpixel temperature retrieval with multispectral sensors". Algorithms for Multispectral and Hyperspectral Imagery V, Orlando, FL, USA, SPIE. Wan, Z. and J. Dozier, 1996. "A generalized split-window algorithm for retrieving land-surface temperature from space." IEEE Trans, on Geoscience and Remote Sensing 34(4): 892-904. Young, S. J., B. R. Johnson and J. A. Hackwell, 2002. "An in-scene method for atmospheric compensation of thermal hyperspectral data." Journal of Geophysical Research (Atmospheres) 107(A14): ACH14-1. 51 3 CHARACTERISING LAKE DYNAMICS IN A DEEP FJORD L A K E USING SUB-PIXEL TEMPERATURES DERIVED FROM MODIS IMAGERY 2 3.1 Introduction Quesnel Lake in the central interior of British Columbia is home to 30% of the Fraser River's sockeye salmon, making it of major importance to the health of the stock and the fisheries industry on the west coast of Canada. It is also the deepest Fjord Lake in the world, glacier fed, and in an active logging region. During the last few years, Fisheries and Oceans Canada (DFO) has noticed record numbers of pre-spawning salmon deaths near the mouth of the Quesnel River, which drains Quesnel Lake. Since this time, DFO scientists at the Institute of Ocean Sciences (IOS) have installed temperature loggers and weather stations in an attempt to understand what is causing the fatalities. The temperature loggers show drastic swings in temperature —10 degrees within 24 hours at the mouth of the Quesnel River — which is beyond the Sockeye salmon's coping range (McCullough 1999). But the cause of the temperature swings is still a mystery. 2 This paper is yet to be submitted for publishing. 52 Figure 3-1 shows a bathymetric map of Quesnel lake and relevant weather and temperature logger stations. It is a long and narrow lake with a E-W span of 81km which forks in the middle for a N-S span of 36km, although the widest transect is only 3km. A sill defining the eastern limit of the West Basin is only 30m below the surface, and is a suspected influence in what is responsible for the large temperature fluctuations in the West Basin and subsequently Quesnel River. Temperature data in the West Basin suggest large internal waves and upwelling events, but constructing the story of why these large temperature fluctuations occur is difficult given point measurements only. In an attempt to paint a synoptic portrait of these events, satellite Thermal Infrared Imagery was used to show the spatial structure of lake surface temperature during a series of major cold upwelling events in August of 2003. 3.1.1 TIR Imagery Thermal infra-red (TIR) imagery from space-borne sensors (ASTER, MODIS) and air-borne sensors (MASTER, SEBASS) has been used in previous research to determine water skin (top 100pm) temperatures of both the ocean (Mcmillin et al. 1984; Emery 1997) and fresh water bodies (Szymanski et al. 1999; Kay et al. 2001; Steissberg et al. 2003). Satellite TIR imagery of lakes is often unreliable evidence as the sensor cannot resolve the narrow reaches where the events took place, as is the case for Quesnel Lake, and the temperature values in these regions are contaminated by land-based signal mixing. Kay et al. (2004) have shown that this mixing effect adds unacceptable error to temperature estimates if ignored, while Szymanski et al. (1999) and Kay et al. (2003) have shown that improved sub-pixel temperature estimates are possible by applying spectral mixture analysis techniques 53 to solve the mixed pixel radiance equation as set out by Gillespie (1992). The current research applies the Sub-Pixel Temperature Retrieval (SPTR) algorithm (Sentlinger et al. 2005) to resolve the narrow reaches of Quesnel Lake in MODIS TIR imagery in an attempt to explain the lake dynamics leading up to the events of Aug 2003 — putting historical spectacles on the MODIS sensor. .1.2 S e i c h i n g Wind-driven upwelling occurs when a prolonged and sufficiently strong wind blows across the water's surface. In an unstratified lake, the free surface will pile up towards the down-wind end of the lake and oscillate about a steady state slope until the wind relaxes. For a stratified lake, the free surface again piles up at the down-wind end of the lake, but because of the stratification, the denser, colder water below moves to the upwind end of the lake, resulting in an interface slope opposite of the free-surface slope. The steady state slope of the interface is given by (Heaps et al. 1966): — = -4- (3-1) dx g\\ Where is the vertical deflection of the interface from it's mean, g' is the reduced gravity gAp/p, h, is the thickness of the epilimnion, and u.2 is the wind stress at the water surface given by: ul = C°P°U" (3-2) Po 54 where CDis the wind drag coefficient, pais the air density, p() is the water density. The wind speed at 10m above ground, u^ is given by u10= u,/(2/10)36 (Tennessee Valley Authority 1972). Integrating along x, the vertical deflection of the interface at the boundaries is: ul L A ^ > - T 7 ( 3 " 3 ) Where L is the length of the lake in the wind direction. The interface oscillates about this steady state slope with a period given by: 2L (3-4) T = c. where: hh, <3"5) H. = '1"2 H c is the wave phase speed, H, is the effective height, h, is the height of the epilimnion, is the height of the hypolimnion. The wind must blow for a time T/4 to push the interface to it's steady state slope in order for a basin scale seiche to be established. Using a 12 hour filter will remove transient wind signals not applicable to this basin, as discussed in Stevens and Lawrence (1997). 3.2 Site Description and Instrumentation Quesnel Lake sits in the interior western hemlock biogeoclimatic zone at 121°00' W and 52°30'N, while the upper reaches of its catchments are in the subalpine Engelmann spruce zone. The average residence time of the lake's 41.8 km3 of water is 10.1 years, based on a 55 mean annual outflow of 131 mV 1 (Potts, 2004). The lake elevation is 730m above MSL. Quesnel Lake is long and thin, <lkm at the mouth of the Quesnel River, which is typical of other deep fjord lakes carved from glaciers. The lake can be divided into 4 Main Basins, with their August 2003 characteristics as follows: • the West Basin has a maximum depth of 100m and is 16km long. It is separated from the Main Basin by a sill at Cariboo Island which comes up to 30m below the surface. The West Basin is the outlet to the Quesnel River and a major resting and feeding area for migrating salmon, especially at the sill (Morrison 2005). Its calculated internal seiche period is 17 hours. • the Main Basin is defined as the moderately deep, 250m maximum depth, basin running east of Cariboo Island to the junction of the East and North Arms. It is 26km long with an estimated seiche period of 1 day. The Horsefly River enters into Horsefly Bay just south of Cariboo Island. • the North Arm's southern extent is marked by a sill at 100m depth, and it's northern extent by the mouth of the Mitchell. The basin's length is 27km with an estimated seiche period of 27 hours. • the East Arm is the deepest basin of the lake, dropping to a maximum depth of 505m at the bottom of sheer lake walls. The bottom water is not anoxic, surprisingly, but how and when the lake turns over is not completely understood. The Niagra Creek feeds into the far end of the East Arm and is generally quite cold due to it's high elevation catchment and glacier. 56 1 2 t ° 3 0 ' w I21°15"W 121°w 120°45W 120°30W 120°15'w Figure 3-1 Quesnel Lake. Bathymetric contours show 100m intervals. Quesnel Lake can be characterized by four bathymetric areas: West Basin, Main Basin, North Arm, and East Arm, with characteristics shown in Table 3-1. The West Basin, west of the 30m sill, has a maximum depth of 100m, while the East Arm reaches depths of 500m. Weather stations are shown along with arrows indicating the relative prevalence of a wind direction—the largest arrow indicates the wind direction with the highest frequency, not necessarily the strongest winds. The major inflows and outflows are also shown with relevant characteristics. 3.2.1 Satellite Imagery MODIS imagery was collected from the EOS Data Gateway (NASA 2003) for the period June 1, 2003 to Aug 31, 2003. There were a total of 520 images available which imaged Quesnel Lake for this period, or approximately 5 images/ day, or a revisit time of approximately 5 hours. Of these 520 images, 145 contained Quesnel Lake at a viewing angle which was useful for analysis and in which the lake was not completely covered by clouds. Clouds presence was determined from the MOD11 Data set Quality Control, 57 which flags clouds with a 99% confidence over land and 66% confidence over lakes (Hook et al. 2002). Of these 145 images, 75 images were considered very good quality, with an elevation angle (between sensor and horizon) that was 70° or larger, and with a large portion of the lake cloud-free. Although the revisit time is very short for MODIS imagery, the corresponding large pixel size of 1km is too coarse to resolve the crucial temperatures from the West Basin. The current research employs a new vector/pure pixel based TIR unmixing algorithm termed Sub-Pixel Temperature Retrieval (SPTR, Sentiinger et al. 2005). This algorithm requires only standard MODIS products as input without the need for ancillary atmospheric parameters. For this study, we used the MODIS standard products: • MOD021KM MODIS/Terra Calibrated Radiances 5-Min LIB Swath 1km at-sensor radiance, radiometrically corrected; • MYD11_L2 and MODll_L2 MODIS/Aqua and MODIS/Terra Land Surface Temperature/Emissivity Daily 5-Min L2 Swath 1km for channels 31 (Hum) and 32 (12um); and • MYD03 and MOD03 MODIS/Aqua and MODIS/Terra Geolocation Fields 5-Min L1A Swath 1km. All data products were swath data for 5 minutes of acquisition with no temporal averaging. They are available at no-charge at the EOS data gateway (NASA 2003). A description of the products can be found at on the MODIS website (NASA 2005). The goal of this research was to essentially unlock the high-resolution data from this extremely valuable dataset. The short revisit time is ideal for monitoring upwelling events, and the acquisition policy of storing all images is very conducive to experimental research such as this. SPTR is a spectral mixture analysis technique which uses a combination of pure-pixel temperature-emissivity values and fractional abundance from vectorized water features. 58 The algorithm solves the two material at-sensor radiance mixing model as described by Gillespie (1992): ^=T,[f]eA]B(A,T]) + f2£,2B(A,T2)) + Ll4 (3-6) where Xx is the transmissivity of the atmosphere, A. is the wavelength, f is the fractional abundance of the jth end-member where j=[l,2 ] for a simple two material case, £x„ is the emissivity, B(A,T) is the blackbody radiation for the jth materials temperature T, L x A is the combined upwelling radiance, reflected the downwelling radiance, adjacency effect and noise in the pixel. All spectral variables in equation 3-6 represent mean values over a sensor's channel width. To solve this equation for T, and T 2 we need four components: • at-sensor radiance, L x s ; • fractional abundance of endmembers, f, and f2; • emissivity estimates for endmembers, Ex, and Ex,; and • atmospheric components Xx, L x A The at-sensor radiance is provided in the standard MODIS product MOD02. The fractional abundances can be generated by "reverse geo-referencing" vector data, that is, placing the vector into the original pixel space and determining sub-pixel fractions. The emissivity estimates are taken from the pure pixel training sets. The atmospheric parameters can be estimated from all pure pixels by fitting a line to equation 3-6, following the IS AC algorithm (Young et al. 2002). Because there is no analytical solution to equation 3-6 for T l and T2, we use a gradient descent method to minimize the error between the measured at-sensor radiance and an estimated at-sensor radiance. This procedure has been dubbed Sub-Pixel Temperature Retrieval, or SPTR. The final SPTR results are geo-referenced back into a hi-res image, the average lake temperature is subtracted in order to partially remove the skin effect. This is based on the assumption that the skin effect is 59 generally homogenous across the lake. A full description of the SPTR algorithm can be found in Sentlinger et. al. (2005). .2.2 Field Measurements Field measurements used in this study consists of thermistor chain (t-chain), weather station, and CTD data. The location of the t-chain moorings and weather stations are shown in Figure 3-1, while CTD stations are shown in Figure 3-2. Eleven temperature moorings have been installed in Quesnel Lake to date, although only two, M8 and M2 are used in this study. M2 was installed by the Institute of Ocean Sciences (IOS) scientists in the West Basin and has been acquiring data since 2002 (Morrison 2005). This is a combination of Stowaway Hobo thermistors and Brancker TLlOOOs moored nominally at 5m ±3m below the surface. M8 is installed on a spar buoy jointly by IOS and UBC (Potts 2004). It is located at the junction of the North and East Arms at a depth of 293 m, with 10 Brancker thermistors beginning just below the surface at 3m and going to depth of 283m. This t-chain was installed July 31s1, 2003, and is the primary temperature measurement for validation purposes of the SPTR results owing to its near-surface temperature logger. Thermistor data in this report has not been filtered in order to capture the instantaneous upwelling signals. Shore-based weather stations were mounted at Elysia fishing resort south of the Main Basin, Goose Spit in the North Arm, and Neilson's Fishing resort at the mouth of the Quesnel River. All stations were Davis Instruments Vantage Pro mounted approximately 2m above the land surface, measuring wind speed and direction, precipitation, humidity, air pressure, and air temperature. CTD casts were made on Julian day 212 (July 29lh) using a Seabird 911 CTD . 60 3.3 Results and Discussion In order to construct a synoptic picture of the upwelling events in the West Basin, we look at four major data sources: temperature profiles from CTD casts, in-situ temperature loggers, wind and weather data, and satellite imagery. We begin by establishing the antecedent conditions from CTD temperature data and establishing the validity of the satellite imagery; we then examine a large upwelling event and a subsequent surface feature; finally we look at the average surface temperature and variance over the summer of 2003. 3.3.1 Antecedent Conditions There are three significant upwelling events in the West Basin identified in the thermistor chain record at M2 during the summer of 2003 on Julian dates 206, 214, and 231. Fortunately, CTD casts, weather stations, and the mooring M8 were installed just days before the second major event on Julian day 214 (August 2nd). The event on day 206 had no Goose Spit wind data, which we demonstrate to be critical for forcing the large internal motions that lead to upwelling in the West Basin. The major event on day 231 was cloudy so MODIS imagery is. sparse for this event. For this reason, we focus on the day 214 event. Figure 3-2 shows the depth profile and CTD cast profiles for the day 212, 2003 casts. Casts within the same basins as defined in Figure 3-1 are very similar. CTD stations QL000-003 are in the Main Basin, North and East Arms. QL004-009 are in the West Basin at various depths. It is interesting that the temperature below 40m is much warmer in the West Basin than in the Main Basin; the sill apparently acts as a cold water filter allowing only the 61 upper layers to mix horizontally. The west basin below 30 metres is warmer probably it concentrates the solar energy absorbed in a smaller volume. Also, the Main Basin, North and East Arms appear warmer than the West Basin in the upper 20m. Figure 3-2 also shows the longitudinal depth profile of the basin from east to west (solid) and from the Junction northward (dotted). This plot shows the two sills separating the West Basin and the North Arm from the Main Basin, as well as the extreme depths of the East Arm. Table 3-1 shows the characteristics of the four Quesnel Lake basins and significant combinations. The lake is complex and the purpose of this table is to establish some ballpark estimates for internal wave velocities and seiche periods. Notice that the Main+North basin has a seiche period of 2 days, which requires only 12 hours of strong winds to establish an internal seiche. This suggests that a 12 hour wind pattern could establish a harmonic standing wave in this basin system. 62 Figure 3-2 Quesnel Lake depth and temperature profiles, Day 212 Shown on this map are the CTD stations for this date and water depth. The upper left panel shows thalweg depth for the west to east (solid) and south to north (dotted) transects. To the right of this is the time of each CTD cast. Profiles within a basin are similar, so representative profiles are shown. Al l stations have established a thermocline at approximately 30m. The lower left plot shows a close-up of the elbow at 30m. Notice that the West Basin profile (QL008) is much different than the Main Basin (QL000&QL002) below the sill level of 30m, and they are similar above the sill. QL005 is near the sill and shows intermediate temperatures. The lower right plot shows the epilimnion, or water at approximately the same temperature of 17-18°, extending from the surface down to 10m in these CTD casts. 63 Table 3-1: Characteristics of Quesnel Lake Basins. Total(North) is the entire lake going north, Total(East) is the entire lake going East. Total(North) Total(Easf) Main+North West Main Norlh East Range (km) 0-69 0-93 16-69 0-16 16-42 42-69 42-93 L (km) 69 93 53 16 26 27 51 H(m) 130 230 150 80 200 100 400 h1(m) 30 30 30 30 30 30 30 g' (m/sA2) 0.015 0.015 0.015 0.015 0.015 0.015 0.015 c (m/s) 0.59 0.63 0.60 0.53 0.62 0.56 0.65 T (days) 1/4T (hours) ^.71 3.44 2.04 0.70 0.97 1.11 1.83 16 29 20 65 12 27 4.19 5 84 6 68 10 98 3.3.2 SPTR validation and the skin effect SPTR was used to estimate the water skin temperature for pure and mixed pixels. The SPTR algorithm has been validated against a concurrent ASTER image by (Sentlinger et al. 2005), and found to be within 0.64° of the ASTER water surface temperature 95% of the time, after the bias between the two products is removed. Given that ASTER temperature error is 0.44° (Hook 2001) after biasing effects are removed, we expect a SPTR RMS error of 1.08° with a 95% confidence for the relative error between pixels. This error estimate does not take into account large areas of homogenous temperatures, which would reduce the uncertainty associated with the predicted temperatures within each area. Determining the absolute accuracy of the SPTR results is more difficult as it requires very accurate radiometers to be mounted on buoys, then radiative transfer algorithms, such as MODTRAN be employed to construct a forward model of what radiance should be reaching the sensor. Hook et al. (2002) reports a bias of 0.3°C too hot for the MOD11 temperature product, but no significant bias for the at-sensor radiance product. Because the SPTR algorithm uses both of these products, and considering all of these factors, the SPTR algorithm absolute error in skin temperature is approximately 1.24° at a 95% confidence interval, or ±0.64°C. For this study, we are primarily interested in the relative 64 error between pixels after subtracting the average lake temperature as this indicates upwelling by partially removing the skin effect (Steissberg et al. 2003). This first order removal of the skin effect is based on the assumption that the skin effect is nearly homogenous across the lake surface. Figure 3-3 shows the retrieved SPTR surface temperatures at the M2 and M8 mooring locations. Unfortunately, M2 is a subsurface mooring so near-surface temperature data is not available. M8, however, has a surface buoy with a thermistor mounted just below the surface at 3m depth. Notice the SPTR temperatures at M8 follow the expected skin effect: hotter midday after being warmed by solar radiation and air temperatures; cooler that the bulk temperature at night due to radiative losses and conductive losses to the air. The fact that the SPTR results oscillate around the in-situ temperature indicates an accurate skin temperature retrieval, albeit M8 is always within a pure-pixel. This pattern is in agreement with Hook et al. (2002) work on Lake Tahoe, which also showed a diurnal variance of ± 1° about the bulk temperature. In Hook et al.'s (2002) work, the bulk temperature refers to the temperature 5cm below the surface; in this paper, the bulk temperature is the shallowest available temperature data at 3m depth. The measured skin effect is compared with time of day for this time period in Figure 3-4, which shows significant variance in both day and night skin effects. In order to determine the influence of the weather variables on the skin effect, and the localized variation in the skin effect, a multiple regression analysis was performed. A total of 65% of this variance is explained by cloud cover and air temperature alone, but surprisingly only 25% by wind speed. It is speculated that this is because the bulk temperature is measured at 3m depth, while the variance due to diurnal wind is primarily 65 above lm (Imberger 1985). A multiple regression analysis was performed on 9 weather variables to determine their influence on the skin effect, with results shown in Table 3-2. Of these 9, only Goose wind, Goose wind past 4 hours, air temperature, cloud cover, and time of day were significant influences on the skin effect. In this analysis, these 5 variables explained 76% of the variance seen in the skin effect, which is summarized in Table 3-3 and Table 3-4. This results lends credibility to the assumption that the skin effect is nearly homogenous across the lake, as cloud-cover, time of day, and air temperature are nearly homogenous. 66 22 21 20 19 18 17 | M 8 S F T R I a • I • | • o n • 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 213 214 215 216 217 218 219 220 - 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 •5T 8 | I 6 I * s j ! d) I Neilson's 12 hr Average Wind Speed I 2 0 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 ^ _ | uysia nr Average wina speea | / j £ 6 ^ 4 2 IS o 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 Julian Day 2003 Figure 3-3 Comparison of SPTR results at M 2 and M 8 with thermistor temperatures and wind data. Plot a) shows retrieved skin temperatures at the M8 location compared to M8 thermistor chain data at 3m below the surface. Notice how the SPTR results track the M8 thermistor temperature, with a skin effect of -1°C during the night, and +1°C during the day. The same comparison is made for M2 in b), however the upper thermistor at M2 is at 5m below the surface. M8_3m is also shown in this plot to compare the skin effect at M2 with M8. Goose Spit wind speed averaged over 12 hours is shown in c), Neilsons in d) and Elysia in e). The large upwelling events at M2 on day 214 and 231 correlate with the two largest wind events at Goose Spit; Elysia and Neilson's have equally large winds on day 214, 225, and 231, although there is no upwelling on day 225. 67 2.00 1.50 1.00 _ 0.50 o m 0.00 « -0.50 -1.00 -1.50 -2.00 1 2 3 4 5 7 8 9 10 11 1213 14 15 16 17 18 19 20 2)1 22 23 24 local time sunrise sunset Figure 3-4 Skin effects through the day. The skin effect is usually positive during the day and negative at night. There are several factors involved however, and many non-linear relationships. Imaging of the lake occurred at three times throughout the summer, approximately 12noon, 10pm and 2am. SPTR temperatures at M2 also appear to exhibit the skin effect, although swing farther either way. This is likely due to the large variation in the bulk temperature over this period. Notice the continual cooling of the M2 SPTR night temperatures at day 214-216 as the cool water that was at 5m rises to the surface at day 216. The lower plots shows the wind magnitude at the three weather stations averaged over 12 hours. Goose spit at the end of the North Arm appears correlated to the large upwelling events at day 214 and day 231. Goose spit 12hour average peaks occur on day 214 and day 231. Neilsons and Elysia peak on day 214, day 225, and day 231, although a major upwelling did not occur on day 224. Lower isotherms rose approximately 5m on day 224, but did no exhibit the prolonged cooling (>44 hours) shown on day 214 and 213, as seen in Figure 3-5. The next section discusses the day 214 upwelling event in more detail. 68 Table 3-2: Multiple Regression Statistics for all variables. P-value is the probability that there is not a relationship. Coefficients Standard Error tStat P-value Intercept -1.33 0.75 -1.77 0.09 hour since sunrise -0.05 0.02 -2.17 0.04 %Cloud -0.97 0.36 -2.70 0.01 ElysiaAirT 0.11 0.03 3.37 0.00 GooseWind 0.11 0.07 1.44 0.16 GooseWindPast -0.14 0.07 -2.10 0.05 ElysiaWind 0.01 0.13 0.09 0.93 ElysiaWindPast 0.06 0.07 0.75 0.46 NeilsonsWind 0.11 0.08 1.39 0.18 NeilsonsWindPast -0.06 0.11 -0.54 0.60 3-3: Multiple Regression statistics for significant variables. Coefficients Standard Error t Stat P-value Intercept -0.79 0.65 -1.21 0.24 hour since sunrise -0.07 0.02 -4.71 0.00 %Cloud -0.99 0.32 -3.12 0.00 ElysiaAirT 0.10 0.03 3.61 0.00 GooseWind 0.18 0.06 2.98 0.01 GooseWindPast -0.15 0.07 -2.27 0.03 Table 3-4: Multiple Regression Results Regression Statistics Multiple R R Square Adjusted R Square Standard Error Observations 0.87 0.76 0.72 0.49 35.00 .3.3 Day 214 upwel l ing event. The day 214 upwelling event is documented in Figure 3-6 and Figure 3-7. The event appears to be a result of prolonged winds at Goose Spit. Based on the error analysis above and in Sentlinger et al. (2005), only temperature differences >1°C are considered as we can be 95% confident that there is a significant difference in surface temperatures. Also the color ramp has been stretched to highlight changes that are greater than this 1.08°C difference. The winds begin midday on day 213, and the images from this time are: 69 a) 213.93: This is the first clear MODIS image, the night of the strong winds. It can be seen from the image that the North Arm surface is cooler than the Main Basin by approximately 1°. This suggests a basin scale internal seiche is relaxing, which is also suggested by the warmer temperature at 5m in M2 and M8. b) 214.89: The next image is immediately after the second day of strong South West winds. In this image, the North Arm surface water is warmer than the West Basin. The thermistor plots shown in the lower left corner of this image show an internal bore traveling through the West Basin, while something is also seen the West Basin in this image. The M2 5m temperatures have dropped 6° at this point. c) 215.45: This image is taken midday on the third day of strong SW winds at Goose Spit. As it is a daytime image, the skin effect is expected to be more variable due to preferential heating and cooling of the surface. d) 215.53: Taken 1:40 later, this image again shows no significant features on the surface. The West Basin surface temperature is approximately the same temperature as the Main Basin. This image pair shows no moving features, however, it is valuable in that the two images are consistent across time, lending credibility to images that do show a moving feature. Note the cloud contamination in the North Arm. e) 215.92: Clouds covered much of this image, although enough of the North Arm and West Basin were exposed to suggest the West Basin was beginning to overturn. f) 216.09: This is the most interesting image of the West Basin upwelling, showing 16° water, 2° colder than the average lake temperature of 18°, fill the surface of the West Basin and travel over the sill into the Main Basin. Although the M2 5m records do not show much variation after day 214, the contour plots in Figure 3-5 show an upwelling on day 215 which continues through day 216. Figure 3-5 shows the contour lines from the M2 and M8 t-chains. Examining this plot we can make several observations: • there are several approximately 2day features in the M2 data; • the upwelling at 214 and at 231 are more extreme than that at 222 and 225. • there is an increase of 6° water and then 9° water in the 214 event, and an almost complete evacuation of the epilimnion. • almost all activity takes place above 30m, the sill level, which is also the bottom of the thermocline. • the 6° isotherm rises approximately 10m during the 214 event. 70 • there is one day of strong winds and a minor upwelling at 214 before the epilimnion is completed evacuated on 215-216, and again there is a minor upwelling on 230 before a complete evacuation of the epilimnion on 231. This does not occur on the 222 event. It appears that a large upwelling in the West Basin do not occur every large wind event; only when the average wind is >5m/s at Goose Spit for more than two consecutive 12 hour periods. Based on the calculated seiche period in Table 3-1, a prolonged wind of 12 hours is long enough to establish an internal seiche in the main + north basin. Although the surface extends into the West Basin and the East Arm is a complicating factor, and indeed there does appear to be several internal waves of different periods, and different modes, bouncing around the entire basin. Our hypothesis is that the sill at 30m and the hydraulic constriction at the sill is enough of a boundary condition to establish the 2day seiche calculated for the main+north basin. The fact that the wind blows in the opposite direction over night may help to force this basin to establish the seiche. 71 o £ 1 3 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 1 1 1 1 T Wind Energy 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 i f f ^f^m'^m 20 Figure 3-5 Quesnel Lake T-Chain data, M2 and M8. This is a plot of t-chain contours from M2 in the West Basin and M8 in the Main Basin. The black lines show a 44 hour period. The top plot show wind energy, u, 3in(m/s) 3, at Goose Spit. A l l winds during the day blow from the South West, while night winds blow from the Northeast, as seen in the panels in Figure 3-6 and Figure 3-7. Although the wind begins to blow on day 213, the isotherms are all pushed up approximately 10m on day 215 where they stay for approximately 44 hours. The largest volume increase appears to be of 6° water. On 217, there appears to be a flood of mostly 9° water, although there is more warm water in general. Isotherms rise on midnight of 218-219, then contours settle back to their pre-event levels. There is a minor wind event on day 222, which does not appear to cross a critical threshold which would cause a complete evacuation of the epilimnion, which does again happen on day 231. In both 214 and 231 event, the 6° water appears to rise 10m, while on the 222 wind event, it only rises 5m. 72 f& 55 B « m | & [ f W 76 £6 915 TOO r Don / —|—s—r 1—5—S 5 1" Figure 3-6 Quesnel Lake Skin temperature during 214 upwelling event. Time series moves from left to right then top to bottom, (a) Temperatures at M8 5m and M2 5m are nearly the same before the upwelling event. Wind has been blowing at Goose Spit during the day previous to this night image, (b) This is the second day of heavy winds at Goose Spit. The North Arm is approximately 1°C warmer than the West Basin, which is starting to show some of the turbulence below. The water at 5m has cooled 6°C from 16°C before the event to 10°C after. M8 cools shortly after, (c) during the night, the surface is nearly homogenous. 73 rj—r;—<r—j—r-Figure 3-7 Quesnel Lake Skin temperature during 214 upwelling event. (d) 1 3/4 hours later is nearly identical, (e) Although clouded, the third day of heavy winds at GooseSpit have now pushed the temperature at M2_5m down to 5.6°, and the surface water has cooled to 16° compared to the 18° of the previous night and the M8 skin temperature, (f) Upwelling at M2 is now very large, 16°C and 3° cooler than the North Arm, and extends into the Main Basin, while M8_5m is warming. 74 The question is begged: what is the difference between the 222 & 225 events and the 214 & 231 events? There is apparently something critical occurring on 214 & 231, whose threshold is not breached in the 222 & 225 events. We put forward two possible explanations: one based on the temperature difference between the two basins below the sill, and the second based on hydraulic control at the sill. Using equation 3-3, we can calculate the steady state deflection of the 6° isotherm with a u2 of 7m/s (Goose Spit 12hr average wind speed on day 212). Table 3-5 shows the parameters used and steady state deflections for the three main events during this period. Table 3-5: Estimated interface deflection at the sill for the three main events, based on wind speeds at Goose Spit. Event 214 222 231 u2 (m/s) 9 7 8 u10(m/s) 16.1 12.5 14.3 L(m) 53000 53000 53000 Cd 0.0013 0.0013 0.0013 pa (kg mA-3) 1.2 1.2 1.2 po (1000kg m A 1000 1000 1000 u*2 (m/s)A2 0.0004026 0.000243545 0.0003181 hi (m) 30 30 30 g1 (m/sA2) 0.015 0.015 0.015 d Interface (m) 23.7 14.3 18.7 Note that all of these values are estimates with many assumptions but they are meant to show the relative difference between the events, and the order of magnitude of the interface displacement. We hypothesize that something critical happens with a 20m displacement at the sill, and looking at the CTD profile in Figure 3-2, we suggest that the 5-6° water that exists on the Main Basin side of the sill is forced up over the sill in the two larger events, causing a rush of 6° water into the West Basin which in turn evacuates the top 10m of epilimnion on the 2nd period of the internal seiche set up in the main +north basin from the previous day's blowing. This is perhaps why the top 10m is evacuated causing the upwelling of the 16 degree water from 10m below. We suggest this is not 75 apparent in the satellite imagery until day 216, even though there is upwelling below beginning day 214, because the top 10m of the epilimnion is approximately the same temperature. Perhaps the 222 event is not large enough to reach this critical point where the West Basin water is warmer than the Main Basin water at the same depth. The second explanation is based on the exchange flows over the sill from work by Potts (2004). His modeled exchange flows, shown in Figure 3-8, across the sill are compared against a theoretical hydraulic limit. The heat content method, uses a conservation of energy approach, predicts flows to exceed the theoretical hydraulic limit on both the day 214 and 231 events, but not the 222 and 225. There have been suggestions that this hydraulic control causes the thermocline to rear up and cause upwelling in the West Basin. 20000 213 215 217 219 221 223 225 227 229 231 233 235 Julian Date in 2003 Figure 3-8 Simulated flows into the West Basin epilimnion. A negative flow rate indicate water leaving the epilimnion, and entering the metalimnion/hypolimnion. Notice that the only two instances that the flows as predicted by the heat content method occur on day 214 and 231; the events on 221 and day 223 do not cross this threshold. This is perhaps an over-simplification of a difficult and non-linear problem, but it is one way of incorporating all of the important suspects into the deduction, and gives scientists something to watch for. If the wind at Goose Spit exceeds 5m/s two days in a row, the top 76 5m, which is the most crucial layer as it drains into the Quesnel River, will evacuate and colder water will be entering the river. .4 Summer 2003 MODIS trends and variance. Using SPTR results for all 145 good quality MODIS images from the summer of 2003, average temperature and variance images were generated. This was done by co-registering all images, masking clouds and 2 lo-res pixels away from the clouds to cover any possible contamination, then calculating the statistical moments in the temporal direction of all good quality pixels: where <AT>, is the average delta temperature for the ith pixel, N is the number of images averaged, T u is the temperature of the i"1 pixel in the tlh image, TLlis the spatially averaged lake temperature of the ta image. The standard deviation is given by: StdDev^jjjj^fa-Tj-iATtf (3-8) These are the moments of the delta images, meaning the value of each of these delta image pixels is the difference of a given pixel from the average lake temperature for the scene in which the pixel exists. This is not the same as the difference of a pixel's temperature from the temporal average of that location. Equation 3-7 and 3-8 measures the i1" pixel's difference from the lake average, intending to remove diurnal temperature patterns to get at the spatial temperature patterns. This may illuminate local prevalent wind patterns as 77 opposed to bulk water patterns, but the averaging over time should remove random lake skin temperature patterns. The variance image, calculated using equation 3-20, depicts 1 standard deviation from the mean delta temperature, or the variance around a given pixel's time averaged difference from the spatial mean. Figure 3-9 shows the number of images in which the pixel location along the E-W thalweg and S-N thalweg are cloud-free. This is an indication of the confidence interval of the measurement; all pixels have > 50 observations which indicates a fairly good statistical measure. This image also shows that the Main Basin is the most frequently cloud free, while the east and North Arm, and the West Basin, are only clear 60% of this time. Many interesting features are found in the average temperature delta image Figure 3-10. Most noticeably are the relatively cold East Arm, and slightly colder North Arm. The East Arm is fed by Niagra creek, which in turn is glacier fed and mostly snowmelt in the summer. Niagra Creek is significandy colder than East Arm surface water during the summer (summer average T = 10°C); however, this cold water will quickly plunge below the surface and only the water surface near the mouth of Niagra Creek will be directiy influenced. It is unknown to what depth Niagra Creek plunges, but it may be partly responsible for the cooling the large region observed. Also, the relatively cooler East Arm may be due to shading from the relatively steep sides flanking the East Arm and the relatively little solar radiation reaching the surface due to cloud cover. Figure 3-11 shows the relative solar illumination hillshade images for three times of the year, December 21, March 15, and July 15, with maximum solar elevation angles of 17°, 40°, and 60° respectively. The only significant times that the east and North Arm are in shadow are in 78 the mornings in the winter, and in the evenings in the summer, as the sun sets in the northwest. Also of note is the colder West Basin, which is believed to be due to these upwelling events documented during day 214 and day 231. ^ . . . i . . . i 0 20 40 60 s 80 Distance along Thalweg (km) Figure 3-9 Number of cloud-free images along W-E and S-N thalweg.. Of the 132 images used in this analysis, approximately 100 of them were cloud free in over the Main Basin, while only approximately 60 were cloud-free in the farthest reach of the west arm and North Arm, and peculiarly, over the 70km mark in the East Arm. The variance image Figure 3-12 shows some unusual areas of low variance. These images show the temporal deviation of the pixels temperature from it's particular scene's spatial average temperature. This is to remove lake-scale skin effect and seasonal heating and cooling. Notice the high degree of variance in the west end of the West Basin and lack thereof at location 8km, suggesting a 1st mode horizontal basin scale seiche, also noted by other researchers from t-chain data in this basin (Morrison 2004). Similarly, the variance of the west half of the Main Basin suggests a node midway across the Main Basin, between the two sills shown in Figure 3-2. Strangely, the north basin has the least amount of variance over time. This doesn't mean that it was always the same temperature, it means it didn't vary in it's difference from the average lake temperature; it was also 0.2° cooler 79 0 20 . 4 0 60 , 80 distance along thalweg (km) Figure 3-10 Average lake temperature difference from mean, summer 2003 This is the difference between a pixel and the lake average in any particular scene, averaged over June-August 2003. For example, the East Arm is usually approximately 1 degree colder than average lake temperature of any given scene. than the average lake temperature (see Figure 3-10). Nodes indicated in this image were taken from local minima in the variance transects shown in Figure 3-12 b). SO Figure 3-11 Solar Illumination of Quesnel Lake These hillshade images model the relative lake illumination as a function of time of day and year. The morning images on the left are 1 hour after sunrise, and the evening images are 1 hour before sunset. The sun only reaches a maximum elevation angle of 17° on December 21 and the East and North Arms are shaded in the morning but not for the remainder of the day. In March, the morning shadows are reduced as the sun rises almost due east. In the summer, there are little shadows in the morning as the sun rises in the Northeast but the North and East Arm are shaded in the evening as the sun sets to the Northwest. Figure 3-13 show the covariance matrix of all transects extracted for the East - West thalweg. The covariance is given by: C o v * r = 4l(A7V, -<A7^,)lATYl -(ATy)) (3-9) Where ATX. is the temperature delta of the ith transect at location X along the transect, <ATX> is the average delta temperature at that location, and similarly for location Y. If the temperature at position Y is hotter than the average delta temperature at the Y position in the ith image, and it is also hotter at position X, and if this happens in many images, X and Y will have a positive covariance. Conversely, if it is often cooler at X when it is hotter at SI Y, these locations will have a negative covariance. When X = Y, as is the case along the diagonal of the covariance matrix, this is simply the equation for variance; the upper right covariance values, Covxy, is a mirror image of the lower left, Covyx. Positive values (red) show positive correlation, while negative values (blue-black) show negative correlation. There is apparent correlation between the two sides of the West Basin node as well as near Niagra Creek, which suggests a possible 2nJ horizontal mode basin scale seiche within these basins. The apparent anti-correlation between the Main Basin temperature and the first leg of the East Arm suggests a possible l s l horizontal mode basin scale seiche contained within this basin system. The TIR imagery analysis in this report would not have been possible with data from a sensor with an infrequent revisit time such as ASTER, or with standard MODIS data products. The statistical analysis requires many images for significance and the standard lkm pixel size of MODIS is too coarse to give any credibility to the patterns seen in the West Basin during the day 214. SPTR unlocks this massive data resource for limnological analysis of narrow lakes, and rivers, maximizing the resource with a limited amount of user expertise and ancillary products, such as atmospheric radiosondes. As shown in this analysis, this resource helps to create a synoptic portrait of major upwelling events, which in turn helps researchers understand complex, multi-variable systems such as Quesnel Lake 82 0 10 20 30 40 50 60 s 70 distance along thalweg (km) Figure 3-12 Map of lake temperature variance about the average lake temperature, summer 2003. This is a map of the standard deviations from the average after removing the lake average for each scene. For example, the west end of the West Basin is generally 0.5°C colder than the rest of the lake (see Figure 3-10), but also varies by 1.3 °C about that -0.5°C. Conversely, the skin temperature at M4 just west of Cariboo Island varies little. Arrows indicate apparent nodes taken from the minima of the variance plot in the upper left of the image. 83 distance along thalweg distance along thalweg Figure 3-13 Covariance matrix of all Quesnel Lake transects The covariance matrix measures degree of correlation, or covariance. The diagonal is the variance at a given location along the west-east transect, off-diagonal is the correlation, or anti-correlation of two locations; the upper right is a mirror image of lower left. Apparent nodes are marked on this plot with points and the corresponding location along the thalweg. Correlations appear between the peaks of variance on either side of the West Basin middle node at 5km. There is also a correlation on either side of the 81.5 km node in the East Arm. The strongest anti-correlation is between the Main Basin and first leg of the East Arm to 63km, suggesting that when the Main Basin is warm, the first leg of the East Arm is cold. A thumbnail of the entire covariance matrix is shown in the upper right. 3.4 Conclus ions MODIS TIR imagery can be used to help understand the internal lake dynamics of narrow lakes. Using the SPTR algorithm, we are able to extract the most limnologically important information from the MODIS dataset, the water temperature, by exploiting several key properties of inland water bodies: they don't move very much, they are near blackbodies, and temperature gradients are very small over large areas. This technique has been applied successfully to a reduce a large dataset — 500 images, 30 channels each — to a relatively small dataset of critical temperature information: 75 clear, near-nadir images of 84 water temperature. The absolute error on the sub-pixel extraction technique is ±0.64° , and ±0.54°C for relative temperature differences within a single image. This latter measure has been employed to remove the skin effect on a lake-wide scale and examine the relative surface temperature patterns in Quesnel Lake, B.C. Sub-pixel retrieved satellite temperature images are useful for determining the spatial extent of upwelling events, transient internal disturbances, and statistical measurements of the lake surface. The upwelling event in the West Basin on Julian day 214 (Aug 2nd) appears linked to strong winds at the Goose Spit weather station that lasted for at least two days. The sub-pixel temperature image showed upwelling had reached the surface 2 days later on day 216 (Aug 4lh; this gives an indication of how deep the isothermal epilimnion was, and the spatial extent of the West Basin colder water upwelling respectively. This has important implications regarding transport from the West Basin into the Main Basin over the natural cold-water filter of the 30 m sill. West Basin upwelling may be the only mechanism by which deep (beneath the sill) West Basin water is able to pass over the sill and into the Main Basin. This mechanism is an important consideration for fisheries management, and the sub-pixel temperature images generated by SPTR give defensible proof of this mechanism. The average lake image shows the peculiar trend of the East Arm to be consistently colder than the rest of the lake. By tracking this trend over many consecutive summers, it may be possible to determine long-term trends of glacial melt which feed into the East Arm, or the impact of forestry practices within the watershed on the temperature of water flowing into the East Arm. The nodes apparent in the variance image of both the West and Main Basin suggest possible 2nd horizontal mode basin seiches within these basins, and that these patterns can exist independently of one another. This has implications for predicting when, and the 85 magnitude of, the upwelling events which devastated the salmon stocks in recent years. Based on these findings, we suggest TIR imagery corrected for sub-pixel mixing with the SPTR algorithm can be a valuable source of information, in conjunction with in-situ sensors, for the analysis of narrow lakes and rivers. The MODIS TIR dataset is free of charge, always acquiring and archiving, and easily accessible, which makes it an ideal data source for future limnological studies of Quesnel Lake and other inland water bodies. 86 3.5 References Emery, W. J., and Yu Y., 1997. "Satellite sea surface temperature patterns." International Journal of Remote Sensing 18(2): 323-334. Gillespie, A. R., 1992. "Spectral mixture analysis of multispectral thermal infrared images." Remote Sensing of Environment 42(2): 137-145. Heaps, N. S. and A. E. Ramsbottom, 1966. "Wind effects on the water in a narrow two-layered lake." Philos. Trans. R. Soc. London A259: 391-430. Hook, S., 2001. "Progress Report: ASTER and MODIS Calibration on Lake Tahoe." http://shookweb.jpl.nasa.gov/ validation. Hook, S., F. J. Prata and S. G. Schladow, 2002. "Final Report: ASTER and MODIS Calibration on Lake Tahoe." http://shookweb.jpl.nasa.gov/validation. Imberger, J., 1985. "The diurnal mixed layer." Limnol. Oceanog. 30(4): 737-770. Kay, J., R. N. Handcock, A. Gillespie, C. Konrad, S. Burges, N. Naveh and D. Booth, 2001. "Stream-temperature estimation from thermal infrared images". Geoscience and Remote Sensing Symposium, 2001. IGARSS '01. IEEE 2001 International. Kay, J. E., R. N. Handcock, K. A. Cherkauer, S. K. Kampf, A. R. Gillespie and S.J. Burges, 2004. "Accuracy of lake and stream temperatures determined from atmospherically corrected thermal-infrared imagery,." Journal American Water Resource Association submitted. Kay, J. E., S. K. Kampf, R. N. Handcock, K. A. Cherkauer, A. Gillespie and S. Burges, 2003. "Lake and stream temperatures from thermal-infrared imagery: Estimation strategies and practical accuracies." Journal American Water Resource Association submitted. McCullough, D. A., 1999. A review and synthesis of effects of alterations to the water temperature regime on freshwater life stages of salmonids, with special reference to chinook salmon. Seattle, Washington, U.S. Environmental Protection Agency. Mcmillin, L. M. and D. S. Crosby, 1984. "Theory and Validation of the Multiple Window Sea-Surface Temperature Technique." Journal of Geophysical Research-Oceans 89(NC3): 3655-3661. Morrison, J., 2004. Unpublished Data, Vynx Design, Sidney BC, V8L 4B2, Canada. Morrison, J., 2005. Unpublished Data, Vynx Design, Sidney BC, V8L 4B2, Canada. NASA, accessed July 2003. "EOS Data Gateway", edcimswww.cr.usgs.gov/pub/imswelcome/, NASA and USGS. 87 NASA, accessed Aug 2005. "MODIS Data Product Description." eosdatainfo.gsfc.nasa.gov/eosdata/aqua/modis/modis_dataprod.html, NASA Potts, 2004. The Heat Budget of Quesnel Lake, British Columbia. Civil Engineering. Vancouver, UBC. Sentlinger, G. I., S. Hook, B. Laval and J. Morrison, 2005. "Sub-pixel water temperature estimation from thermal-infrared imagery using vectorized lake features." Yet to be submitted. Steissberg, T. E., S. J. Hook and S. G. Schladow, 2003. "Characterizing Partial Upwellings and Surface Circulation at Lake Tahoe, CA/NV, USA with Thermal Infrared Images." submitted. Stevens, C. L. and G. A. Lawrence, 1997. "Estimation of wind-forced internal seiche amplitudes in lakes and reservoirs, with data from British Columbia, Canada." Aquatic Sciences 59: 115-134. Szymanski, J. J., C. C. Borel, Q. O. Harberger, P. Smolarkiewicz and J. P. Theiler, 1999. "Subpixel temperature retrieval with multispectral sensors". Algorithms for Multispectral and Hyperspectral Imagery V, Orlando, FL, USA, SPIE. Tennessee_Valley_Authority, 1972. Heat and mass transfer between a water surface and the atmosphere. Laboratory Report No. 14. T. V. Authority. Norris, Tennessee, Tennessee Valley Authority. Young, S. J., B. R. Johnson and J. A. Hackwell, 2002. "An in-scene method for atmospheric compensation of thermal hyperspectral data." Journal of Geophysical Research (Atmospheres) 107(A14): ACH14-1. 88 CONCLUDING REMARKS The work presented in this thesis began with my supervisor Bernard Laval, saying "There's no use for satellite imagery in this project as it can't resolve the western arm". I hope we've addressed this need of limnologists, and provided some insight into how a very important BC Lake works. Many remote sensing scientists and limnonlogists have decried the lack of sufficiently high spatial and temporal TIR imagery for limnological analysis (Kay et al. 2001; Steissberg et al. 2003). MODIS imagery is a prime candidate for the high temporal resolution with up to 4 images/ day, and nobody had used vector data to estimate abundances up until this research. The purpose of having two manuscript chapters was to appeal to both remote sensing scientists interested in the algorithm, and limnologists, interested in the results. The algorithm's main measure of success is what it can help us understand about the lake dynamics, not the absolute error of measurement. The overlap between the manuscript chapters was intended to provide enough background of the algorithm for limnologists, but to really entice them with the results to want to find out more from reading the first manuscript. And the reason the first chapter focused on a lake was to suggest to the 89 remote sensing community that this is not just a theoretical paper, that it has immediate and useful results. Some papers have presented similar approaches (Szymanski et al. 1999) to what is presented here but did not go on to apply the algorithms to real data. Others have done very thorough analysis of the error associated with atmospheric contamination and land mixing (Kay et al. 2004), but did not present a comprehensive approach to handling the problem. What is novel in this thesis' approach to the limnological analysis is the attempt to present all relevant data in an easily digestible manner, together with the satellite image to aid the researcher to get a synoptic picture of the lake dynamics, as opposed to presenting the data in a piecemeal fashion, which relies on the reader to piece together the data. With regards to the analysis of Quesnel Lake, much of this data had not yet been analyzed, so no defensible theories as to what could be causing the upwelling had been put forth. One of the great hurdles in limnological analysis is the intimidating amount of data that must be processed before an understanding can be achieved. What is presented here is a very cursory understanding of a very complex system, but I believe provides a foundation for further inquiries. One of strengths of the research is that is provides a relatively straightforward though effective solution to a complicated problem using a few easily accessible products and tools, and achieves useful results. Many researchers have focused on atmospheric correction techniques (Kay et al. 2003) and not taken advantage of existing corrected products, such as the MOD11 temperature emissivity product. Perhaps one of the weaknesses, or at dependencies, of this approach is on the MOD11 product. Scientists have remarked on the error in this product (Hook et al. 2002; Steissberg et al. 2003), but 90 agree the estimates are within the published error estimates. A further dependency of the approach is the correct georeferencing of the MODIS data. I've tried to address this with an autoregistration component, however, this algorithm works only half of the time. Because of this a manual registration GUI was developed that provided the user with instant feedback of the results and suggestions on how to improve them. The working hypothesis that sub-pixel water temperatures estimates within 1° of the absolute temperature could be achieved through the use of vector/pure-pixel unmixing algorithm was successful. The degree to which it is successful is not fully explored, due to lack of concurrent field measurements and/or hi-res TIR imagery from satellites such as ASTER. The autoregistration algorithm should work much better than it does, and may have a bug buried in the code somewhere. 4.1 Future Research Elaborating on the items in Chapter 1. the future research work would be to: • incorporate sensor point spread function and spectral response function; • use more MODIS channels; • unmix for more than two endmembers; • try alternative temperature and emissivity estimates. Hook et al. has worked out a split-window solution for MODIS channels 31 and 32. • incorporate a skin-bulk temperature model to estimate the temperature 5m down based on time of day, cloud-cover, topography, wind data, and air-temperature. • improve autoregistration results. • convert IDL/ENVI code into MATLAB. The success of the vector based unmixing is applicable to any sensor, not just MODIS, and not only TIR but all spectral ranges. More complex models need to be built to remove the 9 1 ambiguity around the relationship between water skin and bulk temperature, and especially down to 5m. Perhaps incorporating research (Hook et al. 2003) involving the effects of weather variables on the skin-bulk relationship with a model such as DYRESM would allow researchers make inferences about the entire water column profile for every pixel over time, and providing one more source of information to help understand complex lake systems such as Quesnel Lake. 4.2 References Hook, S., F. J. Prata and S. G. Schladow, 2002. "Final Report: ASTER and MODIS Calibration on Lake Tahoe." http://shookweb.jpl.nasa.gov/vahdation. Hook, S. J., F. J. Prata, R. E. Alley, A. Abtahi, R. C. Richards, S. G. Schladow and S. O. Palmarsson, 2003. "Retrieval of Lake Bulk and Skin Temperatures Using Along-Track Scanning Radiometer (ATSR-2) Data: A Case Study Using Lake Tahoe, California." Journal of Atmospheric and Oceanic Technology 20: 534-548. Kay, J., R. N. Handcock, A. Gillespie, C. Konrad, S. Burges, N. Naveh and D. Booth, 2001. "Stream-temperature estimation from thermal infrared images". Geoscience and Remote Sensing Symposium, 2001. IGARSS '01. IEEE 2001 International. Kay, J. E., R. N. Handcock, K. A. Cherkauer, S. K. Kampf, A. R. Gillespie and S.J. Burges, 2004. "Accuracy of lake and stream temperatures determined from atmospherically corrected thermal-infrared imagery,." Journal American Water Resource Association submitted. Kay, J. E., S. K. Kampf, R. N. Handcock, K. A. Cherkauer, A. Gillespie and S. Burges, 2003. "Lake and stream temperatures from thermal-infrared imagery: Estimation strategies and practical accuracies." Journal American Water Resource Association submitted. Steissberg, T. E., S. J. Hook and S. G. Schladow, 2003. "Characterizing Partial Upwellings and Surface Circulation at Lake Tahoe, CA/NV, USA with Thermal Infrared Images." submitted. Szymanski, J. J., C. C. Borel, Q. O. Harberger, P. Smolarkiewicz and J. P. Theiler, 1999. "Subpixel temperature retrieval with multispectral sensors". Algorithms for Multispectral and Hyperspectral Imagery V, Orlando, FL, USA, SPIE. 92 APPENDIX A: DAY 217 SURFACE FEATURE A cold surface feature was noticed in the image at 217.52, although no feature is present in the 217.45 image, 1:40 earlier. This sequence is shown in Figure A- l and the corresponding M8 temperature trace in Figure A-2. There are clouds in the area at this time, and are more prevalent in the 218.08 image. It is difficult to say with certainty that this is not an error, however the disturbance in the M8 top 5m thermistor data suggests that the feature in 217.52 and again at 218.08 passed through the M8 location at approximately 217.7. This works out to a propagation speed of 0.32m/s from 217.52-217.7, and a wave speed of 0.34 m/s from 217.7 to 218.08. The period shown in the 5m trace is 3.6 hours, and using a wave speed of 0.32 m/s gives a wavelength of 3.8 km. The width of the feature from the satellite data is 3.6km. This does not prove its existence, but all evidence points to a legitimate internal wave appearing at the North Arm sill and traveling westward at a velocity of approximately 0.32m/s. As seen in Figure A-2, the feature does not affect the 3m temperature significantly. Examining the wind plots in Figure A- l , however shows that the two satellite images were acquired during a low wind period, while the M8_3m feature occurred during strong winds, suggesting the feature may not have been apparent at the surface due to mixing. Based on the SPTR results the feature is also circular as opposed to basin-wide. 93 Figure A-l Surface feature traveling westward. A cold feature appears at 217.52 although it is not present at 217.45,1:40 earlier. The feature appears at the sill of the north basin and travels west, seen in the M8_5m thermistor trace, and again at 218.08 midway across the Main Basin. It travels at a consistent speed of 0.3m/s, and based on this speed, has a wavelength of 3km from the thermistor trace; approximately the width of the feature in the images 94 Figure A-2 217 Surface feature in M8 temperature data. The width of the disturbance is 3.6 hours. If traveling at .32m/s, this gives a wavelength of 3.8 km. There has been speculation that this is perhaps a Rossby wave. The Rossby radius has been worked out as: f where Rr is the Rossby radius, sqrt(g'h,) is the internal wave speed approximately 0.6m/s, and f is the coriolis parameter = lx l0 4 s '. This gives a Rossby radius of 6km, which is wider than the channel at this point, indicating a Rossby wave is not possible. However, if the internal wave speed is 0.32, which is the speed of this feature, the Rossby Radius is 3.2 km, the width of the channel at the feature's location. This latter wave speed suggests the rotation of the earth may play a part in Quesnel Lake internal waves, although this is a very cursory study. Further studies may look at other MODIS SPTR images in attempt to find this pattern again. 9 5 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items