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

Abstract

This study provides the first technique to investigate the turbulent fluxes 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 fluxes including latent heat flux (Q_E) and sensible heat flux (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² varied from 0.41 (WS), 0.74 (SR), and 0.87 (SP) with RMSE of 5.68, 6.93, and 4.67 W·m⁻², respectively, while Q_H, R² ranged from 0.002 (WS), 0.8030 (SP) and 0.94 (SR) with RMSE of 6.97, 4.39 and 4.90 W·m⁻² 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 fluxes then sharply drop in March and are negligible during June and July. The evaporation processes continue again in August.

1. 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.

2. Materials and Methods

2.1. Datasets

2.1.1. 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.

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.

2.1.2. 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.

2.1.3. 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. Summarized all the data used in this study.

Data	Quantities Estimated	Spatial Resolution	Temporal Resolution
MOD11L2	Water surface temperature	1 km	Daily
MOD06L2	Water surface temperature	5 km	Daily
MERRA IAU T10M	Air Temperature	1/2° × 2/3°	Hourly
MERRA IAU U10M	Wind Speed	1/2° × 2/3°	Hourly

).

2.2. Methods

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:

$${\mathrm{e}}_{\mathrm{a}}=6.11\exp [\frac{Lv}{Rv}\left(\frac{1}{273.15}-\frac{1}{{\mathrm{T}}_{\mathrm{d}}}\right)]$$

(1)

where Lv is the latent heat of vaporization (2.5 × 10⁶ 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:

$${\mathrm{Q}}_{\mathrm{H}}=-\mathsf{\rho }{c}_p{\mathrm{k}}^2\frac{\Delta \overline{\mathrm{u}}·\Delta \overline{\mathrm{T}}}{{\left[\ln \left(\raisebox{1ex}{${\mathrm{z}}_2$}\!\left/ \!\raisebox{-1ex}{${\mathrm{z}}_1$}\right.\right)\right]}^2}{({\mathrm{\Phi }}_{\mathrm{M}}{\mathrm{\Phi }}_{\mathrm{H}})}^{-1}$$

(2)

$${\mathrm{Q}}_{\mathrm{E}}=-Lv{\mathrm{k}}^2\frac{\Delta \overline{\mathrm{u}}·\Delta \overline{{\mathsf{\rho }}_{\mathrm{v}}}}{{\left[\ln \left(\raisebox{1ex}{${\mathrm{z}}_2$}\!\left/ \!\raisebox{-1ex}{${\mathrm{z}}_1$}\right.\right)\right]}^2}{({\mathrm{\Phi }}_{\mathrm{M}}{\mathrm{\Phi }}_{\mathrm{V}})}^{-1}$$

(3)

where ρ is air density (kg·m⁻³); $c_p$ is specific heat of air at constant pressure (J/kg/K); Lv is the latent heat of vaporization (2.5 × 10⁶ J/kg); k is the von Kármán constant (0.41); $\Delta \overline{\mathrm{u}}$ is wind speed gradient derived from MERRA Reanalysis (m·s⁻¹); Δ$\overline{\mathrm{T}}$ is temperature gradient (°C) between Z₁ and Z₂; Δ$\overline{{\mathsf{\rho }}_{\mathrm{v}}}$ is vapor pressure gradient (hPa); $\ln \left(\raisebox{1ex}{${\mathrm{z}}_2$}\!\left/ \!\raisebox{-1ex}{${\mathrm{z}}_1$}\right.\right)$ is natural logarithm of two different heights (approximately 30 m); Z₁ is lower height (at the water surface); Z₂ is upper height (at approximately 30 m).

${({\mathrm{\Phi }}_{\mathrm{M}}{\mathrm{\Phi }}_{\mathrm{H}})}^{-1}\mathrm{and}{({\mathrm{\Phi }}_{\mathrm{M}}{\mathrm{\Phi }}_{\mathrm{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:

$$\mathrm{Ri}=\frac{\mathrm{g}}{\overline{\mathrm{T}}}\frac{(\Delta \overline{\mathrm{T}}/\Delta \mathrm{z})}{{(\overline{\mathrm{u}}/\Delta \mathrm{z})}^2}$$

(4)

where g is acceleration due to gravity (m·s⁻²); $\overline{\mathrm{T}}$ is mean temperature in layer $\Delta \mathrm{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)

$${({\mathrm{\Phi }}_{\mathrm{M}}{\mathrm{\Phi }}_{\mathrm{x}})}^{-1}={(1-5\mathrm{Ri})}^2$$

(5)

Unstable case (Ri negative)

$${({\mathrm{\Phi }}_{\mathrm{M}}{\mathrm{\Phi }}_{\mathrm{x}})}^{-1}={(1-16\mathrm{Ri})}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$$

(6)

We then converted latent heat fluxes (W·m⁻²) to evaporation rate (mm per day) by dividing the latent heat of vaporization (assumed to be 2.5 × 10⁶ J/kg) with 84,000 (seconds in a day). The density of fresh water is assumed to be 1000 kg/m³.

3. Results and Discussion

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). The modelled air temperature from MERRA showed a significant agreement with measured air temperature [22]. The correlation coefficients (R²) 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. Root Mean Square Error (RMSE) and R² 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.

Station Name	R2	RMSE (°C)	Coefficients
			Intercept (°C)	Slope (°C)
SR (Lake Superior)	0.96	1.6	0.91	1.18
WS (Lake Michigan)	0.88	2.1	−5.96	1.00
SP (Lake Huron)	0.88	2.1	−6.34	0.92

).

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.

The statistical relationship between the measured dewpoint temperatures and air temperature from MERRA from all stations also showed a strong correlation. The R² 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 bias, RMSE, and R² between MERRA air temperatures and dew point temperatures from in situ measurements (p-value < 0.0001).

Station Name	R2	RMSE (°C)	Coefficients
			Intercept (°C)	Slope (°C)
WS	0.88	2.0	−3.93	0.85
SP	0.88	2.0	−4.46	0.79
SR	0.88	1.8	−3.67	0.75

).

3.2. 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² varied from 0.41 (WS), 0.74 (SR) and 0.87 (SP) with RMSE of 5.68, 6.92, and 4.67 W·m⁻² respectively. The correlation between Q_H and EC were also significant with the exception of WS R² ranged from WS (R² = 0.002), SP (R² = 0.80) and SR (R² = 0.94) with a RMSE of 6.96, 4.39 and 4.90 W·m⁻², respectively (

Table 4. Summary of R², 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⁻²).

Stations	QE				QH
	R2	RMSE	Intercept	Slope	R2	RMSE	Intercept	Slope
SR	0.87	6.92	−67.91	2.84	0.94	4.90	−56.65	2.8
SP	0.74	4.67	14.47	0.7	0.80	4.39	−28.2	1.10
WS	0.41	5.68	5.54	0.72	0.002	6.96	8.803	0.05

). 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 Q_E and Q_H estimated from bulk aerodynamic approaches with atmospheric stability correction for each of the Great Lakes (

Table 5. Summary monthly mean Q_E and Q_H from Bulk Aerodynamic approach (W·m⁻²) for each of the Great Lakes from July 2001 to December 2014.

Month	Superior		Michigan		Huron		Erie		Ontario		All Lakes
	QE	QH	QE	QH	QE	QH	QE	QH	QE	QH	QE	QH
January	199.7	297.0	133.1	138.1	121.5	141.3	46.3	−3.2	86.8	51.1	314.0	198.6
February	153.4	189.8	98.4	77.3	89.7	79.2	28.9	−10.1	66.6	15.5	54.1	129.2
March	40.5	3.7	26.0	−7.3	28.9	−6.0	20.3	−14.4	23.1	−12.0	7.6	−5.3
April	17.4	−16.3	17.4	−14.6	17.4	−13.9	17.4	−14.5	17.4	−11.2	−0.7	−17.2
May	5.8	−13.4	11.6	−13.5	14.5	−13.2	17.4	−15.6	11.6	−12.2	−5.0	−15.3
June	0.0	−12.2	14.5	−12.8	17.4	−12.5	23.1	−13.6	23.1	−13.9	−5.6	−15.6
July	14.5	−12.6	23.1	−16.3	23.1	−14.5	28.9	−13.0	26.0	−15.3	1.0	−16.4
August	23.1	−15.3	26.0	−11.8	28.9	−12.1	28.9	−9.6	28.9	−11.3	4.4	−13.6
September	28.9	−12.1	37.6	−7.8	37.6	−8.4	40.5	−7.3	37.6	−9.7	10.0	−10.9
October	43.4	−4.6	55.0	−4.1	52.1	−3.7	57.9	−5.7	55.0	−6.9	20.3	−6.3
November	81.0	32.0	78.1	23.6	69.4	19.9	60.8	2.4	57.9	2.0	33.7	26.5
December	156.3	109.1	112.8	104.1	112.8	120.3	60.8	13.0	75.2	37.1	63.4	136.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.0 Wm²) then gradually decreased through the spring and reached its minimum in June (−29.5 Wm²). 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²) and steadily dropped to the minimum in July (−9.0 Wm²). 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].)

3.3. 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 cover during cold years. Monthly mean lake-wide evaporation rate (mm per day) are summarized in

Table 6. Summary of monthly mean lake-wide evaporation rate (mm per day).

Lake		Jan.	Feb.	Mar.	Apr.	May	Jun.	Jul.	Aug.	Sept.	Oct.	Nov.	Dec.
Superior	Mean	6.9	5.3	1.4	0.6	0.2	0.0	0.5	0.8	1.0	1.5	2.8	5.4
	SD	2.1	1.8	0.5	0.1	0.5	0.9	0.3	0.1	0.1	0.5	0.7	1.5
Michigan	Mean	4.6	3.4	0.9	0.6	0.4	0.5	0.8	0.9	1.3	1.9	2.7	3.9
	SD	1.3	1.9	0.2	0.1	0.2	0.2	0.1	0.1	0.3	0.5	0.7	1.0
Huron	Mean	4.2	3.1	1.0	0.6	0.5	0.6	0.8	1.0	1.3	1.8	2.4	3.9
	SD	1.5	1.2	0.2	0.1	0.3	0.2	0.1	0.1	0.2	0.5	0.7	1.2
Ontario	Mean	3.0	2.3	0.8	0.6	0.4	0.8	0.9	1.0	1.3	1.9	2.0	2.6
	SD	1.3	1.0	0.1	0.1	0.5	0.2	0.2	0.2	0.2	0.5	0.6	1.0
Erie	Mean	1.6	1.0	0.7	0.6	0.6	0.8	1.0	1.0	1.4	2.0	2.1	2.1
	SD	0.6	0.2	0.1	0.2	0.2	0.2	0.2	0.2	0.2	0.6	0.7	0.7

n/a is evaporation rate during the period of negative downward fluxes that can be ignored.

.

The difference of monthly median evaporation rate between the upper Great Lakes (Superior, Michigan and Huron) and the lower Great Lakes (Erie and Ontario) might be due to their location. The upper Great Lakes are located in the interior part of the continent that might be affected by continental climate, while the lower Great Lakes situated near the North Atlantic Ocean might be influenced by mild temperature and humidity from a maritime climate. The long-term trend of the monthly mean over-lake evaporation rate of each of the Great Lakes observed from July 2001 to December 2014 is shown in Figure 4, whereas the monthly median spatial variation of the Great Lakes evaporation (mm) from July 2001 to December 2014 is shown in Figure 5.

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. Over-lake monthly average cumulative evaporation (mm) based on the water year.

Water year	Superior	Michigan	Huron	Erie	Ontario
2002	732.4	749.6	692.5	425.6	368.7
2003	867.4	889.5	989.1	594.2	496.8
2004	635.3	768.3	776.4	488.3	515.3
2005	871.7	835.1	701.1	526.4	430.4
2006	718.0	734.9	699.1	532.9	505.3
2007	745.3	710.1	630.1	535.2	363.8
2008	769.1	737.4	710.4	552.6	416.8
2009	825.8	789.5	783.3	501.2	435.8
2010	674.6	693.6	740.7	455.4	475.7
2011	857.9	873.4	829.5	626.6	487.3
2012	619.2	605.5	582.5	395.6	416.0
2013	1047.0	990.2	890.2	861.4	579.5
2014	1093.5	1162.1	1103.9	1017.8	539.1

). 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.9 mm). 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).

4. 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 with 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² 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⁻², respectively, while Q_H, R² 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⁻², 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).