The Estimation of the North American Great Lakes Turbulent Fluxes Using Satellite Remote Sensing and MERRA Reanalysis Data

: This study provides the ﬁrst technique to investigate the turbulent ﬂuxes over the Great Lakes from July 2001 to December 2014 using a combination of data from satellite remote sensing, reanalysis data sets, and direct measurements. Turbulent ﬂuxes including latent heat ﬂux (Q E ) and sensible heat ﬂux (Q H ) were estimated using the bulk aerodynamic approach, then compared with the direct eddy covariance measurements from the rooftop of three lighthouses—Stannard Rock Lighthouse (SR) in Lake Superior, White Shoal Lighthouse (WS) in Lake Michigan, and Spectacle Reef Lighthouse (SP) in Lake Huron. The relationship between modeled and measured Q E and Q H were in a good statistical agreement, for Q E , R 2 varied from 0.41 (WS), 0.74 (SR), and 0.87 (SP) with RMSE of 5.68, 6.93, and 4.67 W · m − 2 , respectively, while Q H , R 2 ranged from 0.002 (WS), 0.8030 (SP) and 0.94 (SR) with RMSE of 6.97, 4.39 and 4.90 W · m − 2 respectively. Both monthly mean Q E and Q H were highest in January for all lakes except Lake Ontario, which was highest in early December. The turbulent ﬂuxes then sharply drop in March and are negligible during June and July. The evaporation processes continue again in August.


Introduction
The North American Laurentian Great Lakes are among the most valuable freshwater resource in the world.In the last fifty years, the water levels have been fluctuating [1], especially in the lower Great Lakes (Erie and Ontario) [2].While water surface temperatures and evaporation rates of the Great Lakes have been increasing, winter ice coverage has also been decreasing rapidly in response to regional warming [3].The study of the surface energy and water balances of the Great Lakes system will expand our limited understanding of how the warming climate affects the hydro-climate of the Great Lakes and their surrounding terrain [4].The greatest challenge in understanding the energy and water balances of the Great Lakes is how to accurately quantify the surface heat fluxes over the water surface.Therefore, acquiring the dynamics of the surface energy and water balances of large lake systems requires, not only direct over-lake measurements, but also further assessment of the spatial variability across the entire lake surface.Nonetheless, measuring the surface fluxes of the Great Lakes is extremely difficult due to their massive size.Lenter et al. [4] pointed out that the annual maximum turbulent fluxes rates on the Great Lakes occurs in late fall or early winter when air temperature is much lower than the water temperature.However, during this time, the "direct measurements" from moving platforms (e.g., The National Oceanic and Atmospheric Administration (NOAA) buoys) are decommissioned due to ice extent.
One of the most accurate ways to measure turbulent fluxes is the eddy covariance technique.This technique measures high-frequency wind and atmospheric data to provide accurate estimates of sensible and latent heat fluxes from the water surface.In order to get year-round direct measurements of the Great Lakes turbulent fluxes, eddy covariance instrumentation must be installed on tall, stable platforms such as lighthouses or small islands [5][6][7][8].Most of these ideal platforms are located in remote areas and are difficult to access.In the Great Lakes, the eddy covariance method has been employed for all-year direct measuring of surface energy balance on Lakes Superior, Huron, and Michigan since 2007, 2009, and 2012, respectively.The main findings from year-round measurements are that the losses of latent and sensible heat are highest in the winter and insignificant in the summer, with a six-month lag between energy inputs and outputs [7].However, the measurements sample a relatively small area (approximately 8 km upwind) relative to the large size of the lakes and, therefore, the spatial variation across the lakes should be determined.For this reason, Earth observation systems including remote sensing and reanalysis are the most convenient means to adequately capture the spatiotemporal variability of the surface energy balance for large lakes like the North American Great Lakes.
Over the past three decades, the study of surface fluxes using remote sensing techniques has advanced greatly with the development of more complex retrieval algorithms.Many studies have presented sophisticated methods on the estimation turbulent heat flux over the open water.For example, Boisvert et al. [9] combined remotely sensed data and reanalysis to calculate moisture flux from the North Water polynya using bulk aerodynamic formulas.Bentamy et al. [10] used the European Remote Sensing (ERS) satellite scatterometer on ERS-2, NASA scatterometer (NSCAT) for wind at 10 m and several Defense Meteorological Satellite Program (DMSP) radiometers (Special Sensor Microwave Imager (SSM/I)) on board the satellites F10-F14 for brightness temperature and later for surface layer air specific humidity.Alcântara et al. [11] used MODIS (Moderate Resolution Imaging Spectroradiometer)-derived land-surface temperature (LST) level 2, 1-km nominal resolution data (MOD11L2, version 5) all available clear-sky from 2003 to 2008, resulting in a total of 786 daytime and 473 nighttime images to derive water surface temperature and heat fluxes over a hydroelectric reservoir in Brazil.In the Great Lakes, Schwab, Leshkevich and Muhr [12] proposed a procedure for producing daily cloud-free maps of surface water temperature based on satellite derived AVHRR (Advanced Very High Resolution Radiometer) imagery.The maps have a nominal resolution of 2.6 km and provide the Great Lakes daily coverage by using previous imagery to estimate temperatures in cloud covered areas.Surface water temperature derived from this procedure compared well with water temperatures measured at the eight NOAA buoys in the lakes.The average difference between the buoy temperature and the satellite-derived temperature estimates was less than 0.5 • C for all buoys.Lofgren and Zhu [13] applied large near shore and offshore gradients in surface temperature and associated heat fluxes.However, these methods were poorly correlated with the network of Great Lakes surface buoys.Although a variety of sources for remotely sensed data are available to calculate evaporation and surface heat fluxes, those data need to be validated with direct measurements.
Retrospective analysis-or reanalysis-provides continuous and consistent information by combining numerical modeling of atmospheric processes with conventional and satellite observations through data assimilation.One of the significant advantages of using reanalysis is that the data are available both spatially and temporally for most atmospheric variables.Nonetheless, reanalysis products contain uncertainty in several variables.Such data from reanalysis products also need to be compared with in situ measurements.Yi et al. [14] evaluated land surface variables from the Modern-Era Retrospective analysis for Research and Applications (MERRA) product against in situ measurements and similar products derived from satellite remote sensing.The results demonstrate that surface air temperature derived from MERRA and two remotely sensed datasets were in significant agreement; however, moderately they were correlated with middle latitude regions and relatively varied in large spatial variation.In this paper, we investigate the turbulent fluxes over the Great Lakes from July 2001 to December 2014 based on the bulk aerodynamic approach by using a combination of data from satellite remote sensing, reanalysis datasets, and direct measurements.

Over Water Direct Measurements
Data are continuously collected from SR in Lake Superior (47.184 • N, 87.225 • W) from 2008-present, SP in Lake Huron (45.773 • N, 84.137 • W) from 2009-present, and WS in Lake Michigan (45.842 • N, 85.136 • W) from 2012-present.Each station was equipped with a 10-Hz eddy covariance system (sonic anemometer, krypton hygrometer, and infrared gas analyzer) providing direct measurements of the sensible and latent heat fluxes and, therefore, evaporation rates.All equipment was mounted at the top of the lighthouse approximately 30 m above the lake surface to measure the sensible heat and latent heat fluxes (evaporation), air temperature, humidity, wind speed and direction, barometric pressure, rain rate, incoming shortwave and long-wave radiation as well as lake surface temperature year-round (Figure 1).All hydro-meteorological data was collected by a data logger connected to a telecommunication system.The percentage of the 30-min eddy covariance and meteorological data that was available during calibration/validation period 93% (7 June 2010 to 31 December 2011; SR), 93% (1 January 2010 to 31 December 2010; SP), and 99% (1 January 2014 to 31 December 2014; WS).Given the completeness of the data during this time period, and to avoid any potential errors, we did not perform any gap-filling.We also validated modeled net radiation with the CNR 1 net radiometers (Campbell Scientific Inc., Logan, Utah) and the modeled heat fluxes with eddy covariance measurements over the lakes surface.Our spatial-temporal assessment and model are used to understand and predict the surface energy balance and hydro-climatology of the Great Lakes.

Over Water Direct Measurements
Data are continuously collected from SR in Lake Superior (47.184°N, 87.225°W) from 2008-present, SP in Lake Huron (45.773°N, 84.137°W) from 2009-present, and WS in Lake Michigan (45.842°N, 85.136°W) from 2012-present.Each station was equipped with a 10-Hz eddy covariance system (sonic anemometer, krypton hygrometer, and infrared gas analyzer) providing direct measurements of the sensible and latent heat fluxes and, therefore, evaporation rates.All equipment was mounted at the top of the lighthouse approximately 30 m above the lake surface to measure the sensible heat and latent heat fluxes (evaporation), air temperature, humidity, wind speed and direction, barometric pressure, rain rate, incoming shortwave and long-wave radiation as well as lake surface temperature year-round (Figure 1).All hydro-meteorological data was collected by a data logger connected to a telecommunication system.The percentage of the 30-min eddy covariance and meteorological data that was available during calibration/validation period 93% (7 June 2010 to 31 December 2011; SR), 93% (1 January 2010 to 31 December 2010; SP), and 99% (1 January 2014 to 31 December 2014; WS).Given the completeness of the data during this time period, and to avoid any potential errors, we did not perform any gap-filling.We also validated modeled net radiation with the CNR 1 net radiometers (Campbell Scientific Inc., Logan, Utah) and the modeled heat fluxes with eddy covariance measurements over the lakes surface.Our spatial-temporal assessment and model are used to understand and predict the surface energy balance and hydro-climatology of the Great Lakes.Locations of over-lake meteorological stations at SR (Lake Superior), SP (Lake Huron) and WS (Lake Michigan).
The turbulent flux footprint, calculated using the Schuepp model [15], stated that during unstable atmospheric stability conditions, 80% of the turbulent fluxes occurred over an upwind area Locations of over-lake meteorological stations at SR (Lake Superior), SP (Lake Huron) and WS (Lake Michigan).
The turbulent flux footprint, calculated using the Schuepp model [15], stated that during unstable atmospheric stability conditions, 80% of the turbulent fluxes occurred over an upwind area of 7-8 km where no land surface existed (given the offshore locations and absence of any land surface at the lighthouse sites).Concerning energy balance closure, we are limited in our ability to calculate this in the traditional manner because the heat storage term of the lakes is so large and unrelated to surface fluxes at the 30-min time scale.During the summer, the net radiation roughly equals the heat storage term, and during the winter, the heat storage term equals the sum of the turbulent fluxes.However, the work of Blanken et al. [7] showed that at the Lake Superior site, the residual of the energy balance (i.e., the heat storage term) was very comparable to independent thermistor-based profile calculated heat storage, and based on those calculations we are likely within 10%-20% of closure.

Satellite Remote Sensing Data
Aboard the Terra (EOS AM) and Aqua (EOS PM) satellites, MODIS provides high temporal resolution viewed over the entire surface of Earth every day or two [16,17].Here, we used MODIS data to estimate the Great Lakes Surface Temperature (GLST) following the scheme provided by Moukomla and Blanken [18].Pixel-based GLST under all sky conditions was created by merging skin temperature derived from MODIS Land Surface Temperature (MOD11L2) and MODIS Cloud product (MOD06L2) from 6 July 2001 to 31 December 2014.All MODIS data are converted into a geospatial database, then processed for each layer in the geographic information system using a geospatial overlay technique.Pixel-based surface heat fluxes can then be calculated through a process of pixel-by-pixel of satellite data and interpolated climatic data.In the validation process, selected scenes were compared with in situ data from meteorological stations that measured shortwave and longwave radiation, for example, meteorological station observations from lighthouses located at SR, SP and WS.

Reanalysis Data: MERRA
The Modern Era Retrospective-analysis for Research and Applications (MERRA) was generated with version 5.2.0 of the Goddard Earth Observing System Version 5 (GEOS-5) and data assimilation system (DAS) [19].The project focused on historical analyses of the hydrological cycle.MERRA products have been available since 1979.The GEOS-5 data assimilation system used for MERRA implements Incremental Analysis Updates (IAU) to adjust the model with the observations.The temporal resolution of MERRA is one hour for two-dimensional diagnostics (surface fluxes, single level meteorology, vertical integrals and land states) and six hours for three-dimensional atmospheric analyses.The spatial resolution is 1/2 • latitude by 2/3 • longitude.However, extensive three-dimensional three-hour atmospheric diagnostics on forty-two pressure levels are also available, but at coarser (1.25 • ) resolution.Among the other modern reanalysis (e.g., ERA-interim, CFSR, NCEP/NCAR), MERRA provides the best spatial and temporal resolution suitable for the study area.Temperature at 10 m above the displacement height (T10M) and Eastward wind at 10 m above displacement height (U10M) were downloaded from the NASA website (http://disc.sci.gsfc.nasa.gov/).It should be noted that the spatial resolution of MERRA is coarser than MODIS so a decision was made to resize MERRA into a one-kilometer resolution.MODIS data products and MERRA data were parameterized for the estimate of the Great Lakes turbulent fluxes based on bulk aerodynamic scheme (Table 1).

Surface Turbulent Fluxes
The bulk aerodynamic method was used to estimate Q H and Q E using wind speed, air temperature and humidity at the lower height (water surface) and upper height (approximately 30 m above the surface).
The near-surface vapor pressure was calculated from the dew point temperature (T d ) using the Clausius-Clapeyron equation as follows: e a = 6.11 exp [ Lv Rv where Lv is the latent heat of vaporization (2.5 × 10 6 J/kg); Rv is the gas constant for water vapor (461 J/kg/K); T d is dew point temperature (K) estimated from MERRA reanalysis data (see Section 3.1).
The basic aerodynamic method applied during neutral stability when the wind gradient is found to increase inversely with height from the surface.However, in the case of non-neutral stability, the stability correction terms were taken into account in the equations.Sensible and latent heat fluxes from bulk aerodynamic method can be expressed as: where ρ is air density (kg•m −3 ); c p is specific heat of air at constant pressure (J/kg/K); Lv is the latent heat of vaporization (2.5 × 10 6 J/kg); k is the von Kármán constant (0.41); ∆u is wind speed gradient derived from MERRA Reanalysis (ms −1 ); ∆T is temperature gradient ( • C) between Z 1 and Z 2 ; ∆ρ v is vapor pressure gradient (hPa); ln( z 2/z 1 ) is natural logarithm of two different heights (approximately 30 m); Z 1 is lower height (at the water surface); Z 2 is upper height (at approximately 30 m).
(Φ M Φ H ) −1 and (Φ M Φ V ) −1 are the dimensionless stability correction terms estimated based on the Richardson number.The relationship between stability correction terms and the Richardson number depends on the temperature gradient and wind speed gradient.The Richardson number is given by: where g is acceleration due to gravity (m•s −2 ); T is mean temperature in layer ∆z (K).
The Richardson number also distinguishes the transition between forced and free convection [20].For instance, when the Richardson (Ri) number is between ±0.01, the atmospheric stability is neutral.The greater instability in the buoyancy term grows in mixed convection where Ri between −0.01 to −1.Free convection presents at the value of Ri larger than −0.1.On the contrary, negative buoyancy reduces turbulent motion so that the Ri value is roughly +0.25 when the flow is laminar and vertical mixing is absent [21].Atmospheric stability based on Ri was calculated through Equations ( 5) and ( 6) as given below: Stable case (Ri positive) We then converted latent heat fluxes (W•m −2 ) to evaporation rate (mm per day) by dividing the latent heat of vaporization (assumed to be 2.5 × 10 6 J/kg) with 84,000 (seconds in a day).The density of fresh water is assumed to be 1000 kg/m 3 .The modelled air temperature from MERRA showed a significant agreement with measured air temperature [22].The correlation coefficients (R 2 ) were 0.96, 0.88, and 0.88 from Lake Superior, Lake Michigan, and Lake Huron respectively, while Root Mean Square Error (RMSE) in MERRA air temperature varied from 1.6 • C to 2.1 • C (Table 2).

Results and Discussion
Here, we used the Clausius-Clapeyron equation (Equation ( 1) for the vapor pressure estimation at 30 m above the water surface.The Clausius-Clapeyron equation requires air temperatures to compute vapor pressure.To achieve this, we regressed the daily average hourly air temperature derived from MERRA reanalysis data with the daily average hourly dew point temperature from the lighthouse-based observations.The scatterplots of air temperature obtained from MERRA versus dew point temperature from the lighthouses is presented in Figure 2.

Estimating Air Temperature and Dew Point Temperature
To accurately estimate dew point temperature, we first validated hourly mean air temperature derived from MERRA reanalysis (T10M) with 30-mins mean air temperatures obtained from yearround measurements from the lighthouses.For the purpose of comparison, only air temperature data during TERRA satellite overpass (~17:00 UTC) was matched with the observed from each station.Also, it should be noted that the observed period was limited by data availability from 25 September 2013-31 December 2014 (WS), 24 September 2009-31 December 2014 (SP), and 11 June 2008-1 April 2012 (SR).The modelled air temperature from MERRA showed a significant agreement with measured air temperature [22].The correlation coefficients (R 2 ) were 0.96, 0.88, and 0.88 from Lake Superior, Lake Michigan, and Lake Huron respectively, while Root Mean Square Error (RMSE) in MERRA air temperature varied from 1.6 °C to 2.1 °C (Table 2).
Here, we used the Clausius-Clapeyron equation (Equation ( 1) for the vapor pressure estimation at 30 meters above the water surface.The Clausius-Clapeyron equation requires air temperatures to compute vapor pressure.To achieve this, we regressed the daily average hourly air temperature derived from MERRA reanalysis data with the daily average hourly dew point temperature from the lighthouse-based observations.The scatterplots of air temperature obtained from MERRA versus dew point temperature from the lighthouses is presented in Figure 2. The statistical relationship between the measured dewpoint temperatures and air temperature from MERRA from all stations also showed a strong correlation.The R 2 ranged from 0.88 (SR), 0.88 (SP), and 0.88 (WS) while RMSE were observed from 2.0 °C (SP), 2.0 °C (WS), and 1.8 °C (SR) (Table 3).The statistical relationship between the measured dewpoint temperatures and air temperature from MERRA from all stations also showed a strong correlation.The R 2 ranged from 0.88 (SR), 0.88 (SP), and 0.88 (WS) while RMSE were observed from 2.0 • C (SP), 2.0 • C (WS), and 1.8 • C (SR) (Table 3).

The Great Lakes Turbulent Fluxes
The bulk aerodynamic method requires measurements of wind speed, air temperature and humidity at the two different heights and atmospheric stability correction term to estimate Q H and Q E .Thus, the MERRA-based air temperature, dew-point temperature, and surface temperature were used to calculate the lake-wide turbulent fluxes.We first calculated the atmospheric stability correction using the Richardson number (Ri).We categorized the Ri number into three atmospheric stability classes: (1) unstable; (2) neutral; and (3) stable.Given the average monthly atmospheric stability as the difference of air and water temperature, wind speed, and vapor pressure, the atmospheric conditions over the Great Lakes were mostly unstable during the fall and late winter, but stay neutral in mid-spring (March) and mid-early fall (October) and were stable from April through August.During the winter, the water temperatures were greater than the air temperatures as well as higher wind speed than during the summer.Monthly mean spatial-temporal variation of the atmospheric stability over the Great Lakes estimated from the Richardson number is illustrated in Figure 3.
The monthly average turbulent heat fluxes estimated from bulk aerodynamic scheme were compared with direct measurement using eddy covariance (EC) method from three different stations.Due to the broad range of data variability especially during the summer, we aggregated daily EC data to monthly means as the representative dataset.Q E and Q H from bulk aerodynamic scheme were correlated well with EC.For Q E, R 2 varied from 0.41 (WS), 0.74 (SR) and 0.87 (SP) with RMSE of 5.68, 6.92, and 4.67 W•m −2 respectively.The correlation between Q H and EC were also significant with the exception of WS R 2 ranged from WS (R 2 = 0.002), SP (R 2 = 0.80) and SR (R 2 = 0.94) with a RMSE of 6.96, 4.39 and 4.90 W•m −2 , respectively (Table 4).In the bulk aerodynamic approach, the atmospheric stability correction and wind speed were taken into account in the equation as it may increase the accuracy of the estimations.Interestingly, we find that there was no significant relationship between modeled and measured Q H over White Shoal Lighthouse.One possible reason may be because of the modeled over White Shoal Lighthouse might influenced by the ice cover.Additionally, the atmospheric stability correction calculated based on the temperature gradient and wind speed, hence, Q H over the ice surface may introduce a non-linear relation between heat flux and air temperature [23,24].We also find that surface temperature from MODIS-derived GLST usually overestimate when ice is present in the pixel.Monthly mean QE and QH estimated from bulk aerodynamic approaches with atmospheric stability correction for each of the Great Lakes (Table 5) exhibited that turbulent fluxes were high in the winter but negligible during the summer.The input energy in form of net-all wave radiation, also peaked during summer [22].The lake-wide maximum mean QE estimated in January (308.0Wm 2 ) then gradually decreased through the spring and reached its minimum in June (-29.5Wm 2 ).In August, QE began to increase again and continued through to December.Whereas monthly mean QH were highest in January (198.6 Wm 2 ) and steadily dropped to the minimum in July (-9.0Wm 2 ).In the winter, QE and QH in Lake Superior are more pronounced than the other Lakes because of greater heat storage capacity of Lake Superior.During the summer, when the water temperature was lower than air temperature, the Great Lakes began to gain heat therefore the negative QE and QH fluxes are corresponded to a heat gain for all months during late spring through late fall.Monthly mean Q E and Q H estimated from bulk aerodynamic approaches with atmospheric stability correction for each of the Great Lakes (Table 5) exhibited that turbulent fluxes were high in the winter but negligible during the summer.The input energy in form of net-all wave radiation, also peaked during summer [22].The lake-wide maximum mean Q E estimated in January (308.0Wm 2 ) then gradually decreased through the spring and reached its minimum in June (−29.5Wm 2 ).In August, Q E began to increase again and continued through to December.Whereas monthly mean Q H were highest in January (198.6 Wm 2 ) and steadily dropped to the minimum in July (−9.0Wm 2 ).In the winter, Q E and Q H in Lake Superior are more pronounced than the other Lakes because of greater heat storage capacity of Lake Superior.During the summer, when the water temperature was lower than air temperature, the Great Lakes began to gain heat therefore the negative Q E and Q H fluxes are corresponded to a heat gain for all months during late spring through late fall.The Great Lakes receive energy from incoming solar radiation from April through July.Therefore, during this time the latent heat flux usually flow downward directly into the water surface (negative flux representing condensation events) so that evaporation during these months is negligible.Great Lakes evaporation occurs during late fall through early winter when the water temperature is much higher than the overlying atmosphere.The cold, dry air flowing down from the Canadian continental climate meets the relatively warm lake and produces strong vapor pressure gradients in addition to high winds, causing the very high lake evaporation rates.The intense convection and the strong local atmospheric respond through lake-effect precipitation [25].The climate teleconnection patterns, such as the North Atlantic Oscillation (NAO), the Arctic Oscillation (AO), and the El Niño-Southern Oscillation (ENSO) might also impact on the Great Lakes ice cover [26][27][28][29][30][31][32].Many studies show the significant impacts of the ENSO including the increasing gustiness of winds (e.g., Li et al. [31]; Enloe et al. [30]) and ice cover (e.g., Bix et al. [27], Assel and Rodionov [26].)

Overlake Evaporation
During the winter, lake-effect snow is common, which influences the local terrestrial water balance.An excellent example of the typical synoptic atmospheric conditions during cold-air outbreaks is lake-effect clouds.The highest evaporation rate was found in Lake Superior during January with a monthly average of 6.9 mm per day, followed by Lake Michigan (4.6 mm per day), Lake Huron (4.2 mm per day), Lake Ontario (3.0 mm per day), and in Lake Erie during December (2.1 mm per day), respectively.The evaporation rate then decreased sharply in March, and the processes continued again in August.During the study period, we observed that Lake Ontario evaporation (not shown) extended through the month of March in a high-ice-extent year (2002 and 2013).Lake bathometry might also affect the amount of energy, which goes into warming the lakes in the summer (the heat storage term).For example, shallow, warm, southern Lake Erie's evaporation occurs earlier in the season, and evaporates more than deep, northern Lake Superior.In warm years, evaporation from Lake Superior is greater than during cold years, but there is a seasonal shift with increased evaporation during the spring and fall during the warm years (more mid-winter evaporation during the cold years).In contrast, Lake Erie's evaporation, typically the largest of all the Lakes, shows a similar pattern in exceptionally warm or cold years; weaker vapor pressure gradients during warm years mimic high ice     The over-lake cumulative evaporation was calculated based on monthly mean over-lake evaporation from July 2001 to December 2014.However, the calculations of the cumulative evaporation were based on the water year period beginning in October and ending in September of the following year (Table 7).Therefore, data from July to September 2001 and from October to December 2014 were not included.The highest annual cumulative evaporation was in 2014 in Lake Michigan (1162.1 mm) follow by Lake Huron (1103.9mm).During the winter of 2013-2014, the Great Lakes experienced consistent lowest temperatures and highest ice cover in recent history [33].Additionally, during this time Lake Superior had approximately 88% ice cover.With extent ice cover, therefore, the monthly cumulative evaporation of Lake Superior was found slightly lower than the Lake Michigan-Huron system.Evaporation of the Great Lakes is not spatially uniform, as evaporation rates vary with the movement of synoptic-scale air masses over the lake [34].The evaporation begins in Lake Erie, Lake Ontario, and south of Lake Michigan as early as July.The highest evaporation rates tend to occur in the lateral regions of Lake Superior in January, particularly in the Thunder Bay area and then switches to offshore regions by January and February when ice cover begins to limit evaporation in near shore regions.The evaporation rate drops close to zero in April, May, and June before Lakes Erie, Ontario and South Michigan begin to evaporate again in July.It has also shown that Lake Superior does not start to evaporate again until August.Recent extreme ice extent in Lake Superior during winter 2013-2014 slightly hold the overall evaporation trend down.While other Lakes also experienced extreme ice extent, they were only partial covered, which in turn intensified evaporation.The evaporation rate in Lake Erie tends to increase nearly 0.9 mm per day while the Michigan-Huron system tendency was to roughly increase 0.4 mm per day.Lake Ontario had the lowest positive trend (~0.25 mm per day).

Conclusions
We provide the first technique to investigate the turbulent fluxes including latent heat flux (Q E ) and sensible heat flux (Q H ) based on bulk aerodynamic scheme over the Great Lakes during July 2001 to December 2014.The bulk aerodynamic approach was calculated from wind speed, temperature gradient and vapor pressure gradient the atmospheric correction term calculated based on the Richardson number.The proposed method used a combination of data from satellite remote sensing, reanalysis data sets and direct measurements.
Estimated Q E and Q H were compared with the direct eddy covariance measurements from the rooftop of three lighthouses-SR in Lake Superior, WS in Lake Michigan, and SP in Lake Huron.The relationship between modeled and measured Q E and Q H were in good statistical agreement, for Q E , R 2 varied from 0.41 (WS), 0.74 (SR), and 0.87 (SP) with a RMSE of 5.68, 6.93, and 4.67 W•m −2 , respectively, while Q H , R 2 ranged from 0.002 (WS), 0.80 (SP) and 0.94 (SR) with a RMSE of 6.97, 4.39 and 4.90 W•m −2 , respectively.Both monthly mean Q E and Q H were highest in January for all Lakes with the exception of Lake Ontario where they were highest in early December then sharply dropped in March before the evaporation processes continued again in August.Recent extreme ice extent in Lake Superior during winter 2013-2014 slightly hold the overall evaporation trend down.While other Lakes also experienced extreme ice extent, they were only partially covered, which in turn intensified evaporation.The evaporation rate in Lake Erie tends to increase nearly 0.9 mm per day while the Michigan-Huron system tendency roughly increased 0.4 mm per day.Lake Ontario had the lowest positive trend (~0.25 mm per day).

Figure 1 .
Figure 1.Locations of over-lake meteorological stations at SR (Lake Superior), SP (Lake Huron) and WS (Lake Michigan).

Figure 1 .
Figure 1.Locations of over-lake meteorological stations at SR (Lake Superior), SP (Lake Huron) and WS (Lake Michigan).

3. 1 .
Estimating Air Temperature and Dew Point Temperature To accurately estimate dew point temperature, we first validated hourly mean air temperature derived from MERRA reanalysis (T10M) with 30-min mean air temperatures obtained from year-round measurements from the lighthouses.For the purpose of comparison, only air temperature data during TERRA satellite overpass (~17:00 UTC) was matched with the observed from each station.Also, it should be noted that the observed period was limited by data availability from 25 September 2013-31 December 2014 (WS), 24 September 2009-31 December 2014 (SP), and 11 June 2008-1 April 2012 (SR).

Figure 3 .
Figure 3. Spatial variation of monthly average of atmospheric stability based on the Richardson number from July 2001 to December 2014.

Figure 3 .
Figure 3. Spatial variation of monthly average of atmospheric stability based on the Richardson number from July 2001 to December 2014.

Figure 4 .
Figure 4. Long-term trend of monthly mean over-lake evaporation rate each of the Great Lakes observed from July 2001 to December 2014.

Figure 4 .
Figure 4. Long-term trend of monthly mean over-lake evaporation rate each of the Great Lakes observed from July 2001 to December 2014.

Figure 4 .
Figure 4. Long-term trend of monthly mean over-lake evaporation rate each of the Great Lakes observed from July 2001 to December 2014.

Figure 5 .
Figure 5. Spatial variation of the Great Lakes monthly average evaporation (mm) from July 2001 to December 2014.

Figure 5 .
Figure 5. Spatial variation of the Great Lakes monthly average evaporation (mm) from July 2001 to 2014.

Table 1 .
Summarized all the data used in this study.

Table 2 .
Root Mean Square Error (RMSE) and R 2 between Modern-Era Retrospective analysis for Research and Applications (MERRA) air temperature and measured air temperature at a height of ~30 m above the water surface from in-situ measurements during 2014.

Table 3 .
The bias, RMSE, and R 2 between MERRA air temperatures and dew point temperatures from in situ measurements (p-value < 0.0001).

Table 4 .
Summary of R 2 , RMSEs and coefficients of monthly average QE and QH measured by eddycovariance method versus modeled from bulk aerodynamic using data from Stannard Rock lighthouse, Spectacle Reef lighthouse, and White Shoal lighthouse (Wm −2 ).

Table 4 .
Summary of R 2 , RMSEs and coefficients of monthly average Q E and Q H measured by eddy-covariance method versus modeled from bulk aerodynamic using data from Stannard Rock lighthouse, Spectacle Reef lighthouse, and White Shoal lighthouse (W•m −2 ).

Table 5 .
Summary monthly mean Q E and Q H from Bulk Aerodynamic approach (W•m −2 ) for each of the Great Lakes from July 2001 to December 2014.

Table 7 .
Over-lake monthly average cumulative evaporation (mm) based on the water year.