Remote Sensing of the North American Laurentian Great Lakes’ Surface Temperature

: The Great Lakes Surface Temperature (GLST) is the key to understanding the effects of climate change on the Great Lakes (GL). This study provides the ﬁrst techniques to retrieve pixel-based GLST under all sky conditions by merging skin temperature derived from the MODIS Land Surface Temperature (MOD11L2) and the MODIS Cloud product (MOD06L2) from 6 July 2001 to 31 December 2014, resulting in 18,807 scenes in total 9373 (9434) scenes for MOD11L2 (MOD06L2). The pixel-based GLST under all sky conditions was well-correlated with the in situ observations ( R 2 = 0.9102) with a cool bias of ´ 1.10 ˝ C and a root mean square error (RMSE) of 1.39 ˝ C. The study also presents the long-term trends of GLST. Contrary to expectations, it decreased slightly due to the impact of an anomalously cold winter in 2013–2014. this study provides the ﬁrst techniques to retrieve pixel-based GLST under all sky conditions by merging skin temperature derived from MODIS Land Surface Temperature (MOD11L2) and MODIS Cloud product (MOD06L2). The GLST was derived from the daily swath ﬁles ranging from 6 July 2001 to 31 December 2014, resulting in 18,807 scenes in total, and 9373 (9434) scenes for MOD11L2 (MOD06L2). MODIS-derived skin temperature was validated with water surface temperature observed from buoys. A cool bias (~ ´ 1.1 ˝ C) was observed for MODIS-derived GLST when compared to the in situ temperature measurements. The monthly mean GLSTs were interpolated from daily MODIS-derived surface temperatures. The lake-wide monthly mean water temperatures were observed to peak in July for the lower Great Lakes (Erie and Ontario) and in August for the upper Great Lakes (Superior, Michigan and Huron) followed by a rapid decline in the autumn in all lakes. The months of January (Erie and Ontario) and February (Superior, Michigan and Huron) had the lowest water temperatures then steadily increased through March to July. Mean monthly maps of GLST reveal that surface temperatures decreased from the lateral zone toward the lake center and from the upper Great Lakes (Lake Superior, Lake Michigan and Huron) to the Lower Great (Lake Erie and Lake Ontario) during the fall throughout the winter. Spatially, in early spring, surface temperatures started to rise from the shorelines toward the lake center and from the lower Great Lakes toward the upper Great Lakes. The peak temperatures were observed in July for the lower Great Lakes (Erie and Ontario) and in August for the upper Great Lakes (Superior, Michigan and Huron). Nonetheless, for future study, several limitations of the method proposed have to be considered. The accuracy of


Introduction
Lake surface temperature is one of the most important variables in the study of the surface energy budgets of large lakes [1].The temperature difference between the water and the overlying air play a large role in the magnitude and the direction of turbulent fluxes of heat and water vapor in the air-water boundary layer.Subsurface temperatures also control the exchange of heat within the water column.This entire process affects the stability of the water column, which is important in expanding our understanding of the characteristics of the thermal regime of the Great Lakes [2].There has been an overall warming trend of the surface waters in the past decade in the Great Lakes [3] and consequently a decline of Great Lakes ice cover, both in duration and extent [4].It is also clear that the presence or absence of ice has a significant effect on thermal moderation of the Great Lakes [5].Most of the water surface temperature measurements rely on meteorological data from the buoys, rafts, and near-shore weather stations.The water temperature sensors measure approximately one meter below the surface of the water.In many lakes, buoys can only be deployed during the open-water season and are decommissioned in wintertime due to the ice cover.Previous studies have shown that the lake-atmosphere temperature gradients are largest during winter months when over-lake measurements are limited [6].The Great Lakes Surface Temperature (GLST) is typically defined as the temperature from the uppermost layer to a depth of one meter [7].The factors that control the GLST are solar radiation, cloud cover, wind-driven mixing, water clarity, ice cover, and lake bathymetry.Clouds are also considered one of the most important factors that dictate water surface temperature.In generally, the presence of clouds is associated with a cooler surface due to the cloud's high albedo [8].
Over the past three decades, remote sensing has been a crucial tool for evaluating the surface temperature of large lakes, reservoirs, and saline lakes.Several studies successfully applied the thermal infrared channel from a variety of satellite remote sensing platforms to retrieve water surface temperature.For example, Bussieres, Verseghy and MacPherson [9] used the thermal infrared channel onboard AVHRR to derive the water surface temperature of 132 bodies of water in western Canada.In Lake Malawi, the surface temperature distributions were derived from MODIS and AVHRR [10].In Brazil, over 1200 images of clear sky MODIS land surface temperature were used to study water temperature and surface heat fluxes of Itumbiara hydroelectric reservoir [11].Similar remote sensing methods were also employed in Lakes Vänern and Vättern using MODIS Sea Surface Temperature product [12].In Lake Tahoe, California-Nevada, the average brightness temperatures from the second Along-Track Scanning Radiometer (ATSR-2) nadir 11-µm and 12-µm channels were correlated with in situ observations of radiometric skin temperature and bulk temperature from four monitoring stations to retrieve skin temperature and bulk temperature [13].Here, skin temperature is defined as the temperature within the first few millimeters of the surface, and the bulk temperature is defined as the temperature below the surface to a few meters' depth [14].The same methods were applied to determine the skin effect over a large, high-altitude freshwater lake [15].
The study using thermal infrared (TIR) to retrieve water surface temperature has been applied not only to freshwater bodies but also saline lakes, for example in the Great Salt Lake [16,17], the Dead Sea [18], and Urmia Lake [19].In the North American Great Lakes, Schwab, Leshkevich and Muhr [20] proposed a technique for generating daily cloud-free maps of surface water temperature based on AVHRR imagery from NOAA's Coast Watch program.The maps have a 2.6 km spatial resolution and daily coverage of the Great Lakes.Surface water temperature estimates from this method were well correlated with water temperatures from NOAA weather buoys in the lakes (ice-free conditions only).The mean difference between the buoy temperature and the satellite-derived surface temperature was less than 0.5 ˝C for all buoys.Clear sky conditions are usually required to obtain these surface temperatures using the TIR methods [21].However, in many areas, cloud free images are rarely available, as is the case in the winter over the Great Lakes.The cloud properties together with their spatiotemporal variations create uncertainty in retrieving radiometric temperatures from the thermal infrared sensors.Therefore, GLST data availability is sparse, particularly during the autumn and winter when the sky conditions are mostly cloudy.Also, data validation is almost impossible since the data from buoys are also not available during that time.
This study provides the first technique to retrieve pixel-based data under all sky conditions over the Great Lakes by merging skin temperature derived from MODIS Land Surface Temperature (MOD11) and MODIS Cloud product (MOD06).After regressing MODIS-derived skin temperature against water temperature one meter below the surface observed from buoys, GLST is then estimated from regression analysis for each of the Great Lakes.The results are validated with in situ data from multi-platform meteorological stations (NOAA buoys and lighthouses).Similar approaches for retrieving GLST were used to observed seasonal lake turnover date.The main objective of this study is to present a framework for the GLST estimations under all sky conditions and to understand the characteristics of its spatiotemporal variability.

Study Area
The surface area for all of the Great Lakes covers a total area of 244,000 km 2 , with 16,000 km of coastline.Lake Superior is the largest, the deepest and has the highest elevation of the lakes in the Great Lakes drainage basin.Lake Superior drains through St. Mary's River into Lake Huron.Lakes Michigan and Huron are considered one system due to having the same surface elevation, as they are connected through the Straits of Mackinac.Water then flows through St. Clair and Detroit River into Lake Erie, the shallowest lake in the system.The largest volume of flow occurs between Lake Erie and Lake Ontario through the Niagara River.The Great Lakes waters flow through the St. Lawrence River, and end their journey at the Atlantic Ocean in the Gulf of St. Lawrence (Figure 1).Due to the large spatial extent and the sheer volume of water in the Great Lakes, the lakes have been shown to play a large role in the surface climate of the region [22].In the autumn, the heat stored in the water bodies is released in the form of latent heat of evaporation.Thus, the highest evaporation rates usually occur during this time of year (late fall and early winter).Shorter day lengths together with an increase in cloudiness gradually decrease the air temperatures.The temperature of the lakes continues to decline in the winter, when the frigid and dry Arctic air from the northwest travels over the moist warmer lakes.In early spring, increasing solar radiation begins to heat the surface and the overlying atmosphere, resulting in the thaw of snow and lake ice (if any is present).During this period, the lakes are usually warmer than the land and tend to keep the surrounding area cool, thus occasionally calm conditions continue well into April.In the summer, the longer sunshine duration warms the surface layer of water in the lakes, making surface water warmer and less dense than the heavier and colder water below.Due to the large spatial extent and the sheer volume of water in the Great Lakes, the lakes have been shown to play a large role in the surface climate of the region [22].In the autumn, the heat stored in the water bodies is released in the form of latent heat of evaporation.Thus, the highest evaporation rates usually occur during this time of year (late fall and early winter).Shorter day lengths together with an increase in cloudiness gradually decrease the air temperatures.The temperature of the lakes continues to decline in the winter, when the frigid and dry Arctic air from the northwest travels over the moist warmer lakes.In early spring, increasing solar radiation begins to heat the surface and the overlying atmosphere, resulting in the thaw of snow and lake ice (if any is present).During this period, the lakes are usually warmer than the land and tend to keep the surrounding area cool, thus occasionally calm conditions continue well into April.In the summer, the longer sunshine duration warms the surface layer of water in the lakes, making surface water warmer and less dense than the heavier and colder water below.

Methods
Wan and Dozier [23] explained a fundamental theory for the land surface temperature (LST) algorithm based on the Planck function as: where L is emitted spectral radiance; B is brightness temperature; λ is wavelength; Ts is surface temperature (K); and ε is spectral emissivity Planck function.However, Equation (1) can only be applied during clear-sky conditions when the spectral radiance (L) emitted from the surface is detectable.In cloudy conditions, L cannot be retrieved from The Great Lakes basin boundaries and geographic locations for water surface temperature measurements.

Methods
Wan and Dozier [23] explained a fundamental theory for the land surface temperature (LST) algorithm based on the Planck function as: Lpλ, Tq " εpλq, Bpλ, T s q (1) where L is emitted spectral radiance; B is brightness temperature; λ is wavelength; T s is surface temperature (K); and ε is spectral emissivity Planck function.
However, Equation (1) can only be applied during clear-sky conditions when the spectral radiance (L) emitted from the surface is detectable.In cloudy conditions, L cannot be retrieved from the surface since the cloud tops, not the surface, is what is detected by the satellites.Therefore, in this study, the 5 km skin temperature from the MOD06L2 product was used when clouds were present in the scene by merging the LST from MOD11L2 with the 5 km skin temperature from MOD06L2 using the equation: where MODIS-LST is the daily skin temperature derived from MODIS; MOD06 is the skin temperature derived from MOD06L2; MOD11 is the LST from MOD11L2; and MASK is the Cloud mask from MOD06L2.
To distinguish between cloudy and clear skies, the MOD06L2 cloud mask product was applied.Images with less than 20% cloud coverage over the study area were recognized as cloud-free.MODIS-LST was then regressed with one-meter below surface temperature from NOAA buoys to create the GLST as a linear function.The GLST was simply formulated as: where GLST is the daily average based on hourly means Great Lakes surface temperature at one-meter below the surface a (1.1 ˝C) and b (0.921 ˝C/ ˝C) are empirical coefficients determined by the least square fit (R 2 = 0.9102).

MODIS Products
In this study, two MODIS products, MOD11L2 and MOD06L2, were used in GLST estimations.MOD11 level 2 collection 5 (MOD11L2) nominal 1 km resolution data was generated with MODIS bands 31 (11 µm) and 32 (12 µm).Though designed for retrieving the radiometric temperature related to the thermal infrared radiation emitted from a broad range of land cover types, LST can only be retrieved for the clear-sky pixels with high confidence [24].Consequently, there are large LST spatial-temporal gaps over the Great Lakes especially during the winter due to the high percentage of cloud cover.Therefore, in our analysis, 49.84% of MOD11L2 and 50.16% of MOD06L2 data were used.Figure 2 indicates that both clear-sky MOD11L2 data were mostly available in summer, and peaked in June.MOD06L2 data were available nearly twice that of the clear sky data.MOD06 level 2 collection 051 (MOD06L2) contained many different cloud properties derived from 14 bands in total.The MOD06L2 surface temperature was obtained from various sources including the MOD11L2 product, the National Centers for Environmental Prediction (NCEP) gridded analysis, and the Data Assimilation Office data of the Global Data Assimilation System (GDAS) as input to the cloud mask algorithm available in 5 km resolution as detailed in [25].
Remote Sens. 2016, 8, 286 4 of 15 the surface since the cloud tops, not the surface, is what is detected by the satellites.Therefore, in this study, the 5 km skin temperature from the MOD06L2 product was used when clouds were present in the scene by merging the LST from MOD11L2 with the 5 km skin temperature from MOD06L2 using the equation: where MODIS-LST is the daily skin temperature derived from MODIS; MOD06 is the skin temperature derived from MOD06L2; MOD11 is the LST from MOD11L2; and MASK is the Cloud mask from MOD06L2.
To distinguish between cloudy and clear skies, the MOD06L2 cloud mask product was applied.Images with less than 20% cloud coverage over the study area were recognized as cloud-free.MODIS-LST was then regressed with one-meter below surface temperature from NOAA buoys to create the GLST as a linear function.The GLST was simply formulated as: where GLST is the daily average based on hourly means Great Lakes surface temperature at one-meter below the surface a (1.1 °C) and b (0.921 °C/°C) are empirical coefficients determined by the least square fit (R 2 = 0.9102).

MODIS Products
In this study, two MODIS products, MOD11L2 and MOD06L2, were used in GLST estimations.MOD11 level 2 collection 5 (MOD11L2) nominal 1 km resolution data was generated with MODIS bands 31 (11 µm) and 32 (12 µm).Though designed for retrieving the radiometric temperature related to the thermal infrared radiation emitted from a broad range of land cover types, LST can only be retrieved for the clear-sky pixels with high confidence [24].Consequently, there are large LST spatial-temporal gaps over the Great Lakes especially during the winter due to the high percentage of cloud cover.Therefore, in our analysis, 49.84% of MOD11L2 and 50.16% of MOD06L2 data were used.Figure 2 indicates that both clear-sky MOD11L2 data were mostly available in summer, and peaked in June.MOD06L2 data were available nearly twice that of the clear sky data.MOD06 level 2 collection 051 (MOD06L2) contained many different cloud properties derived from 14 bands in total.The MOD06L2 surface temperature was obtained from various sources including the MOD11L2 product, the National Centers for Environmental Prediction (NCEP) gridded analysis, and the Data Assimilation Office data of the Global Data Assimilation System (GDAS) as input to the cloud mask algorithm available in 5 km resolution as detailed in [25].

Surface Measurements
The NOAA National Data Buoy Center (NDBC) has continuously operated buoys in the Great Lakes since 1979.All of the buoys are deployed during the ice-free season, generally from April to December.The buoys measure hourly "bulk" water temperatures, at a depth of 1 m below the water surface.The thermistor can measure a temperature range from ´5 to 40 ˝C with the accuracy of ˘0.1 ˝C.For the purposes of this study, we used the available data from NDCB buoys during the satellite overpass time (approximately 10:00 a.m.UTC) from 2010 to 2013 as a representative dataset for regression analysis plotted against radiometric "skin" temperature from MODIS to create GLST then validated with surface temperature from in situ measurements in 2014.

Processing of MODIS Data
Many studies used only cloud-free surface temperature (e.g., [10][11][12]16,19,26]) because during the cloudy season, the lack of clear-sky images complicates satellite-derived lake surface temperature, particularly regarding the accuracy of interrannual variability statistics.The MODIS land surface temperature product (MOD11L2) and MODIS cloud product (MOD06L2) were downloaded from the NASA website [27] in the hierarchical data format (HDF).The daily swath files range from 6 July 2001 to 31 December 2014, resulting in 18,807 scenes in total (9373 scenes for MOD11L2 and 9434 scenes for MOD06L2).In general, the MODIS swath width was large enough to cover all the Great Lakes.However, in case the satellite data did not cover the whole study area in one scene, the adjacent scene(s) were downloaded to cover a study area.When the coverage required more than two scenes in one day, scenes were mosaicked into one scene.All data were reprojected into grid format with the MODIS Conversion Tool (available for downloading at [28]) then clipped with the Great Lakes boundary vector [29].It should be noted that skin temperature extracted from MOD11L2 data were available at 1-km resolution only for clear sky conditions.Nonetheless, due to high cloud cover of the study area, especially during the winter months, a large amount of spatial and temporal data were missing.We therefore applied the skin temperatures derived from MOD06L2 for filling those data gaps.However, data from MOD06L2 were available only at 5-km spatial resolution; thus, it was essential to resize all data from MOD06L2 to 1-km spatial resolution.We used a gdal_translate function from Geospatial Data Abstraction Library (GDAL) open source software for resizing MOD06L2 from 5-km spatial resolution to 1-km spatial resolution.The clear sky-pixels were separated from the cloudy pixels by using a cloud marking technique.Then skin temperature extracted from both MOD11L2 and resized-MOD06L2 were merged into one scene using Equation (2) for the spatial gap filling.

Validation
At the water surface, there is a thin layer a few millimeters thick that exchanges heat and moisture with the atmosphere and emits infrared radiation; the so-called "cool skin layer" [13].This layer is generally observed to be approximately 0.1 ˝C-0.5 ˝C lower than the bulk temperature of the water beneath [14].Also, below the cool layer, there is a layer a few centimeters thick which is warmer than the "cool skin layer" above, caused by the absorption of sunlight; the so-called warm layer [30].Both observations of bulk surface temperature from NDBCs are generally made a few centimeters below the cool skin and warm layer [15].These surface temperatures are called "bulk surface temperatures" [7].The bulk surface temperature can be significantly different from skin temperature, especially during the summer under weak winds and high amounts of insolation [31].To create the MODIS-Derived one-meter below surface temperature (bulk surface temperature), it is essential to validate the MODIS-derived skin temperature against the one-meter below surface temperature from the in situ measurements.Satellite infrared sensors during completely cloud-free conditions only observed the temperature from the immediate surface or "skin" of the water rather than the bulk surface temperature as measured from the buoy.The results of the validation between daily MODIS-derived skin temperature and the one-meter below surface temperature measure from nine NDBC buoys during satellite overpass in 2010-2013 are presented here.The daily MODIS-derived skin temperature was extracted from a pixel centered at the buoy's geographic coordinates.Next, the data were matched up with one-meter surface temperature measured from nine NDBC buoys at the time that satellite overpassed.It should be noted that the data from buoy NDBC 45006 (West Superior) were excluded from the calculation due to the low correlation of skin temperature and bulk temperature (R 2 = 0.8067).
The linear regression from synchronized data of MODIS-derived skin temperature versus in situ measurements are shown in Table 1.The relationships between MODIS-derived skin temperature and in situ measurements from NDBCs were in excellent agreement (R 2 = 0.9102).The excellent relationship between MODIS-derived skin temperature and bulk surface temperature from NDBC buoys was due to the physical properties of the water.As water allows shortwave radiation to penetrate through deeper layers, this allows energy absorption through the water column.Also, the fluid motion permits heat transfer from uppermost layer through the deeper layer.As a result, validated skin temperature from satellite remote sensing can be used as an estimate of the bulk surface temperature.However, the results showed a cool bias of -1.10 ˝C and a root mean square error (RMSE) of 1.39 ˝C between the satellite-derived and the in situ observations.Furthermore, MOD11L2 can only be used as acceptable data if the error is better than 1 ˝C.However, this analysis combined MOD06L2, in which the error is assumed to be 5 ˝C higher or lower than the estimated temperature.

GLST Spatial and Temporal Variability Monthly Variation
Daily GLSTs were interpolated into monthly mean temperatures based on the period July 2001 to December 2014 to observe the annual variation (Figure 3).In general, the monthly average GLSTs peaked in July for the lower Great Lakes (Erie and Ontario) and in August for the upper Great Lakes (Superior, Michigan, and Huron), followed by a rapid decline in the autumn in all lakes.The lowest GLST was in the month of January (Erie and Ontario) and February (Superior, Michigan and Huron), and then the temperatures steadily increased during March to July (Table 2).Lake Superior had the lowest minimum monthly temperature in February (−4.7 °C) whereas the monthly mean GLSTs for February and March were slightly different (0.1 °C).From February to August, the temperature increased approximately 12 °C.From September to January, the temperatures decrease from approximately 12 °C (September) to 3 °C (January).Lake Michigan and Huron showed the same pattern; the maximum temperatures peaked in August at approximately 20 °C, and afterward, the water temperature rapidly dropped by around 15 °C from September to January.Like Lake Superior, the minimum temperature from Lake Michigan and Huron were in the month of February (−2.7 °C and −0.4 °C respectively).Although GLST should represent surface temperature below the water surface, the minimum surface temperature below 0 °C could be caused by the ice cover.Unlike the upper Great Lakes, however, the mean water temperatures of Lake Erie and Ontario were highest in July (20.8 °C ± 1.6 °C and 21.2 °C ± 1.4 °C respectively), then the mean temperature decreased quickly to 2.1 °C ± 1.7 °C (Erie) and 1.8 °C ± 2.4 °C (Ontario) in January.In Lake Ontario, the lack of ice cover might contribute to high standard deviation (SD) of monthly mean temperature, especially in January.Since Lake Ontario rarely freezes over, the temperature difference between warm water surface and overlying cold air create unstable atmospheric conditions [32].From late February or early March, the water temperatures slowly increased and reached their peak in July.A small summer temperature plateau was present in the maximum temperature in July and August.In general, the monthly average GLSTs peaked in July for the lower Great Lakes (Erie and Ontario) and in August for the upper Great Lakes (Superior, Michigan, and Huron), followed by a rapid decline in the autumn in all lakes.The lowest GLST was in the month of January (Erie and Ontario) and February (Superior, Michigan and Huron), and then the temperatures steadily increased during March to July (Table 2).Lake Superior had the lowest minimum monthly temperature in February (´4.7 ˝C) whereas the monthly mean GLSTs for February and March were slightly different (0.1 ˝C).From February to August, the temperature increased approximately 12 ˝C.From September to January, the temperatures decrease from approximately 12 ˝C (September) to 3 ˝C (January).Lake Michigan and Huron showed the same pattern; the maximum temperatures peaked in August at approximately 20 ˝C, and afterward, the water temperature rapidly dropped by around 15 ˝C from September to January.Like Lake Superior, the minimum temperature from Lake Michigan and Huron were in the month of February (´2.7 ˝C and ´0.4 ˝C respectively).Although GLST should represent surface temperature below the water surface, the minimum surface temperature below 0 ˝C could be caused by the ice cover.Unlike the upper Great Lakes, however, the mean water temperatures of Lake Erie and Ontario were highest in July (20.8 ˝C ˘1.6 ˝C and 21.2 ˝C ˘1.4 ˝C respectively), then the mean temperature decreased quickly to 2.1 ˝C ˘1.7 ˝C (Erie) and 1.8 ˝C ˘2.4 ˝C (Ontario) in January.In Lake Ontario, the lack of ice cover might contribute to high standard deviation (SD) of monthly mean temperature, especially in January.Since Lake Ontario rarely freezes over, the temperature difference between warm water surface and overlying cold air create unstable atmospheric conditions [32].From late February or early March, the water temperatures slowly increased and reached their peak in July.A small summer temperature plateau was present in the maximum temperature in July and August.

Stratification Onset Date
One of the most important properties of water is the thermal expansion.Pure water reaches its maximum density at a temperature of 3.98 ˝C and becomes less dense at temperatures above or below this value.From heating or cooling processes, lake stratification is created as the result of the stable vertical density gradient.As the estimation of lake-wide mean daily surface temperatures were made available, lake stratification onset date can be determined.The lake stratification start date was defined as when the water temperature showed a consistent temperature of between 3.8 ˝C and 4.2 ˝C.The recent warming trends are responding to the decrease in Great Lakes ice cover.Lake Superior alone has experienced a 79% decrease in ice cover over the past few decades.The reduction in ice cover is related to the increases in summer water temperature and the timing and duration of the summer stratification [22].However, our analysis showed the opposite trends due to the anomalously cold conditions in winter 2013-2014.During the winter of 2013-2014, ice covered approximately 80% of the Great Lakes surface and extended through late spring 2014.At the same time, input of heat into the lakes decreased considerably, because the ice albedo reflected most incoming solar radiation back to the atmosphere, delaying the onset date of stratification in late summer of 2014.With the impact of extreme winter 2013-2014 conditions, trends of surface temperature may have shifted to pre 1997-1998 El Niño conditions [33,34].The first day of lake turnover and stratification duration is shown in Table 3.We assume that there is no lake turn over in this year.
The analysis shows lake stratification begins in mid-April to late June (DOY 100-160) and ends around the last week of December to early January of the following year (DOY 350-20).Summer GLST were sensitive to the timing of the start of stratification, and therefore summer GLST were strongly related to winter ice conditions [35].To delineate summer stratification, pixels that temperature equal or greater than 4 ˝C in May and June were selected.Figure 4 presents the analysis of lake-wide summer stratification maps using surface temperatures from July 2001 to December 2014.(2002, 2005, 2006, 2007, 2010, and 2011) on the contrary, in the colder years, stratification was delayed until June (2003, 2004, 2008, 2009, 2013, and 2014.Note that in 2012 surface temperature from our analysis showed no stratification (GLST > 4°C).

Interannual Variation
The pattern of surface temperature was well defined in the time series (Figure 5).From the period of analysis, the time series revealed a warm winter in 2011-2012.Monthly anomalous surface temperatures were calculated by subtracting from the monthly mean temperature during the study period (not shown here).The anomalously high surface temperature can be seen from spring to summer 2011, around 5-10 °C in all Lakes while anomalously cold years were also observed in winter 2011-2012 and winter 2013-2014 (~6 °C).In the warmer years, stratification began in May (2002, 2005, 2006, 2007, 2010, and 2011) on the contrary, in the colder years, stratification was delayed until June (2003, 2004, 2008, 2009, 2013, and 2014.Note that in 2012 surface temperature from our analysis showed no stratification (GLST > 4 ˝C).

Interannual Variation
The pattern of surface temperature was well defined in the time series (Figure 5).From the period of analysis, the time series revealed a warm winter in 2011-2012.Monthly anomalous surface temperatures were calculated by subtracting from the monthly mean temperature during the study period (not shown here).The anomalously high surface temperature can be seen from spring to summer 2011, around 5-10 ˝C in all Lakes while anomalously cold years were also observed in winter 2011-2012 and winter 2013-2014 (~6 ˝C).During an abnormally warm year, the lakes will experience an earlier ice-melt, allowing more solar radiation to be absorbed into the water and resulting in higher surface temperatures, and therefore a delay in the following winter's freeze-up period.An early spring thaw leads to a late winter freeze-up, creating a positive feedback loop [1,36].Studies have shown that changing climate During an abnormally warm year, the lakes will experience an earlier ice-melt, allowing more solar radiation to be absorbed into the water and resulting in higher surface temperatures, and therefore a delay in the following winter's freeze-up period.An early spring thaw leads to a late winter freeze-up, creating a positive feedback loop [1,36].Studies have shown that changing climate has resulted in decreased ice cover over the Great Lakes [4,[37][38][39][40].Now, this study is also finding that the changes in the ice cover themselves have an impact on the Great Lakes' hydro-climate and thermal stratification as a whole.The consequences of North Atlantic Oscillation (NAO) and El Niño-Southern Oscillation [41] on Great Lakes ice cover were investigated using lake ice observations from winters 1963-2010 and National Centers for Environmental Prediction (NCEP) reanalysis data [5].The Arctic Oscillation (AO) is defined as the shifts in multiple features of the air pressure, temperature, and the strength of the jet stream.Both AOs (most of the time associated with NAO) had impacts on Great Lakes ice cover, with a tendency for lower ice cover during the positive phase of the NAO.On the contrary, in the negative phase of AO and NAO, the jet stream shifts southward and can develop troughs that drive the frigid, polar air southward.Frigid polar air overlying the Great Lakes area create abnormal ice cover and extent.While ENSO events often influence the Great Lakes' ice cover, strong La Niña events are often associated with lower ice cover.Ice cover also affects surface temperature.The cold, dry Arctic air encroached into the Great Lakes region in both December 2008 and January 2009 led to the extensive ice cover in winter 2008-2009, and again in the winters of 2013-2014 and 2014-2015.This behavior is more likely due to the dynamics of a low-pressure system surrounding Iceland also known as the Icelandic low.Strangely, the Icelandic low was persistent in December 2008, with the anomaly centered on Greenland and extending over Hudson Bay.In January 2009, the anomalous Icelandic low developed dual centers, an event that rarely happens in winter.These dual centers were displaced westward-one stayed over Iceland, and the other stayed over the Labrador Sea.Also, both low centers in January 2009 were shifted southward compared with December 2008.Therefore, there was a deep trough of low pressure from the Great Lakes to the Southeastern United States [39].In winter 2013-2014, anomalous cold weather reached across the eastern United States.Frigid temperatures stayed most of the month of January.As a consequence, many areas experienced record low temperatures while the Great Lakes ice coverage reached up to 75% of the surface, the highest observed since 1996 [42].The significant loss of Arctic sea ice might weaken the temperature gradient and consequently weaken the jet stream over North America [43].Anomalous freezing air flows vertically as the process of heat transfer from the Arctic Ocean to the atmosphere during autumn and early winter.Ballinger, Allen and Rohli [44] pointed out that the significant correlations between the AO, PDO, and PNA teleconnection might be influencing the spatial variability of the cold weather outbreak.

Spatial Variations
The spatial variations for GLST were computed from daily MODIS-Derived surface temperatures.The lake-wide monthly mean water temperatures peaked in July for the lower Great Lakes (Erie and Ontario) and in August for the upper Great Lakes (Superior, Michigan and Huron), followed by a rapid decline in the autumn in all lakes.The months of January (Erie and Ontario) and February (Superior, Michigan and Huron) had the lowest water temperatures and then a steady increase during March to July.The peak temperature was observed in July and August due to the thermal inertia property of the Lakes.Spatial GLST maps (Figure 6) revealed a decrease in surface temperatures from the lateral zone toward the lake center and from the upper Great Lake (Lake Superior, Lake Michigan and Huron) to the lower Great Lakes (Lake Erie and Lake Ontario) throughout the fall and the winter.In early spring, temperatures started to rise from the shorelines toward the lake center and from the lower Great Lakes toward the upper Great Lakes.The temperature peak was observed in July for the lower lakes and in August for the upper Great Lakes.Heat storage in the deeper lakes (e.g., Lake Superior) was much larger than the shallower lakes (e.g., Lake Erie), thus the smaller and shallower lakes tended to lose heat faster than the deeper lakes.
throughout the fall and the winter.In early spring, temperatures started to rise from the shorelines toward the lake center and from the lower Great Lakes toward the upper Great Lakes.The temperature peak was observed in July for the lower lakes and in August for the upper Great Lakes.Heat storage in the deeper lakes (e.g., Lake Superior) was much larger than the shallower lakes (e.g., Lake Erie), thus the smaller and shallower lakes tended to lose heat faster than the deeper lakes.

Conclusions
Remote sensing offers the potential for monitoring the surface temperature of large areas where in situ measurements are difficult, such as the Great Lakes.To deal with the cloudy conditions over the Great Lakes region, this study provides the first techniques to retrieve pixel-based GLST under

Conclusions
Remote sensing offers the potential for monitoring the surface temperature of large areas where in situ measurements are difficult, such as the Great Lakes.To deal with the cloudy conditions over the Great Lakes region, this study provides the first techniques to retrieve pixel-based GLST under all sky conditions by merging skin temperature derived from MODIS Land Surface Temperature (MOD11L2) and MODIS Cloud product (MOD06L2).The GLST was derived from the daily swath files ranging from 6 July 2001 to 31 December 2014, resulting in 18,807 scenes in total, and 9373 (9434) scenes for MOD11L2 (MOD06L2).MODIS-derived skin temperature was validated with water surface temperature observed from buoys.A cool bias (~´1.1 ˝C) was observed for MODIS-derived GLST when compared to the in situ temperature measurements.The monthly mean GLSTs were interpolated from daily MODIS-derived surface temperatures.The lake-wide monthly mean water temperatures were observed to peak in July for the lower Great Lakes (Erie and Ontario) and in August for the upper Great Lakes (Superior, Michigan and Huron) followed by a rapid decline in the autumn in all lakes.The months of January (Erie and Ontario) and February (Superior, Michigan and Huron) had the lowest water temperatures then steadily increased through March to July.Mean monthly maps of GLST reveal that surface temperatures decreased from the lateral zone toward the lake center and from the upper Great Lakes (Lake Superior, Lake Michigan and Huron) to the Lower Great (Lake Erie and Lake Ontario) during the fall throughout the winter.Spatially, in early spring, surface temperatures started to rise from the shorelines toward the lake center and from the lower Great Lakes toward the upper Great Lakes.The peak temperatures were observed in July for the lower Great Lakes (Erie and Ontario) and in August for the upper Great Lakes (Superior, Michigan and Huron).Nonetheless, for future study, several limitations of the method proposed have to be considered.The accuracy of GLST during the winter needs to be improved due to the lack of bulk surface temperature measured from NDBCs.We also noticed that the surface temperature derived from MODIS cloud products have limited ability to estimate radiative skin temperature (accuracy ˘5 ˝C) especially when the surface is too cold, for example when snow and ice surfaces are present in the pixel.There are also possible effects of GLST accuracy from the resizing and merging processes.

Figure 1 .
Figure 1.The Great Lakes basin boundaries and geographic locations for water surface temperature measurements.

Figure 1 .
Figure 1.The Great Lakes basin boundaries and geographic locations for water surface temperature measurements.

Figure 2 .
Figure 2. Comparison of percent data coverage between MOD11L2 and MOD06L2 over the study area with an example from Terra MODIS acquired on 1 January 2014 data (a); Percent cloud cover

Figure 2 .
Figure 2. Comparison of percent data coverage between MOD11L2 and MOD06L2 over the study area with an example from Terra MODIS acquired on 1 January 2014 data (a); Percent cloud cover over each of the Great Lakes estimated using MOD06L2 MODIS Cloud product from July 2001 to December 2014 (b).

Figure 3 .
Figure 3. Monthly mean GLST computed from July 2001 to December 2014; vertical bars represent one standard deviation.

Figure 3 .
Figure 3. Monthly mean GLST computed from July 2001 to December 2014; vertical bars represent one standard deviation.

15 Figure 5 .
Figure 5. Lake-wide time series evaluated from GLST monthly means of daily averages during July 2001 to December 2014 (°C).

Figure 5 .
Figure 5. Lake-wide time series evaluated from GLST monthly means of daily averages during July 2001 to December 2014 ( ˝C).

Figure 6 .
Figure 6.Spatial variation for monthly mean GLST ( ˝C) from July 2001 to December 2014.

Table 1 .
Station name and location for water surface temperature measurements and linear regression between MODIS derived skin temperature and surface temperature observed from buoys.

Table 2 .
Monthly mean GLST of each lake ( ˝C) observed from July 2001 to December 2014.

Table 3 .
Start date (day of year) and duration (day) of the lake-wide average stratification.
S = Start, D = Duration.During 2012 water temperature remained above 4 ˝C.