Satellite-Derived Light Extinction Coefficient and its Impact on Thermal Structure Simulations in a 1-D Lake Model

One essential optical parameter to specify in lake models is water clarity, which is parameterized based on the light extinction coefficient (Kd). A global constant value of Kd is usually specified in lake models. One-dimensional (1-D) lake models are most often used as lake parameterization schemes in numerical weather prediction and regional climate models. This study aimed to improve the performance of the 1-D Freshwater Lake (FLake) model using satellite-derived Kd for Lake 10 Erie. The CoastColour algorithm is applied to MERIS satellite imagery to estimate Kd and evaluated against Kd derived from Secchi disk depth (SDD) field-based measurements collected during Lake Erie cruises. A good agreement is found between field and satellite-derived Kd (RMSE = 0.63 m, MBE = -0.09 m, I_a = 0.65) (in situ data was collected in 2004, 2005, 2008, 2011, 2012). The constant (0.2 m) and satellite-derived Kd values as well as radiation fluxes and meteorological station observations are then used to run FLake at the location of a buoy where lake surface water temperature (LSWT) was measured 15 in 2008. Results improved compared to using a constant Kd value (0.2 m) (lake-specific yearly average Kd value: RMSE=1.54 oC, MBE= -0.08 oC; constant Kd value: RMSE=1.76 oC, MBE= -1.26 oC). No significant improvement is found in FLake simulated LSWT when Kd variations in time are considered using a monthly average. Therefore, results suggest that a timeindependent, lake-specific, and constant satellite-derived Kd value can reproduce LSWT with sufficient accuracy for Lake Erie NDBC station. 20 A sensitivity analysis is also performed to assess the impact of various Kd values on the simulation of LSWT, mean water column temperature (MWCT), lake bottom water temperature (LBWT), mixed layer depth (MLD), water temperature isotherms as well as ice dates and thickness. Results show that FLake is sensitive to variations in Kd to estimate the thermal structure of Lake Erie. Dark waters result in warmer spring and colder fall temperatures compare to clear waters. Dark waters always produce warmer MWCT, shallower MLD, longer ice cover duration, and thicker ice. The sensitivity of FLake to Kd 25 variations is more pronounced in the simulation of MWCT and MLD. The model is particularly sensitive to Kd values below 0.5 m. This is the first study to assess the value of integrating Kd from the satellite-based CoastColour algorithm into the FLake model. Satellite-derived Kd is found to be a useful input parameter for simulations with FLake and possibly other lake models, and with potential for applicability to other lakes where Kd is not commonly measured.


Introduction
There has been significant progress made in recent years in the representation of lakes in regional climate models (RCMs) and numerical weather prediction (NWP) models.Lakes are known to be an important continental surface component affecting weather and climate, especially in lake-rich regions of the Northern Hemisphere (Eerola et al., 2010;Martynov et al., 2012;Samuelsson et al., 2010).They can influence the atmospheric boundary layer by modifying the air temperature, wind, and precipitation.Therefore, consideration of lakes in NWP and RCMs is essential (Kheyrollah Pour et al., 2012, 2014b;Martynov et al., 2010).In order to account for lakes in NWP and RCMs, a description of energy exchanges between lakes and the atmosphere is required (Eerola et al., 2010).Lake surface water temperature (LSWT) is one of the key variables when investigating lakeatmosphere energy exchanges (Kheyrollah Pour et al., 2012).There are various approaches to obtaining LSWT and integrating it in NWP models, such as through climatic observations, assimilation and/or lake parameterization schemes (Eerola et al., 2010;Kheyrollah Pour et al., 2014a).Currently, LSWT is broadly modeled in NWP models using one-dimensional (1-D) lake models as lake parameterization schemes (Martynov et al., 2012).For instance, the 1-D freshwater lake (FLake) model performs adequately for various lake sizes, from shallow to relatively deep (artificially lim-Published by Copernicus Publications on behalf of the European Geosciences Union. ited to 40-60 m depth; Kourzeneva et al., 2012a), located in both temperate and warm climate regions (Kourzeneva, 2010;Martynov et al., 2010Martynov et al., , 2012;;Mironov, 2008;Mironov et al., 2010Mironov et al., , 2012;;Samuelsson et al., 2010;Kourzeneva et al., 2012a, b).
One of the optical parameters required as input in the FLake model is water clarity.This variable is considered as an apparent optical property and is parameterized using the light extinction coefficient (K d ) to describe the absorption of shortwave radiation within the water body as a function of depth (Heiskanen et al., 2015).A global constant value of K d is usually used to run lake models, including FLake.For example, Martynov et al. (2012) coupled FLake with the Canadian Regional Climate Model (CRCM) by specifying a K d value equal to 0.2 m −1 (A.Martynov, personal communication, 2015) for all North American lakes, including Lake Erie, for the years 2005-2007.Heiskanen et al. (2015) evaluated the sensitivity of two 1-D lake models, LAKE and FLake, to seasonal variations and the general level of K d for simulating water temperature profiles and turbulent fluxes of heat and momentum in a small boreal Finnish lake.Modeled values were compared to those measured for the lake during the ice-free period of 2013.The study found a critical threshold for K d (0.5 m −1 ) in 1-D lake models.Heiskanen et al. (2015) concluded that for too-clear waters (K d < 0.5 m −1 ), the model is much more sensitive to K d .The study recommends a global mapping of K d to run the FLake model for regions with clear waters (K d < 0.5 m −1 ) for future use in NWP models.The authors also suggest that this global mapping can be time-independent (i.e., with a constant value per lake).
The global mapping of K d can be derived from satellite imagery.Potes et al. (2012) used empirically derived water clarity from space-borne medium resolution imaging spectrometer (MERIS) measurements to test the sensitivity of FLake to this parameter.The sensitivity analysis was conducted using two K d values, representing the expected extreme water clarity cases for their study (1.0 m −1 for clear water and 6.1 m −1 for dark water).The importance of lake optical properties was evaluated based on the evolution of LSWT and heat fluxes.Their results showed that water clarity is an essential parameter affecting the simulated LSWT.The daily mean LSWT increased from 1.2 • C in clear water to 2.4 • C in dark water (Potes et al., 2012).Water clarity measurements are included in water quality monitoring programs; however, global measurements of clarity are not yet available.Satellite remote sensing can provide water clarity observations to the modeling communities at higher spatial and temporal resolutions, to fill the gap of field measurements.
In recent years, a number of algorithms have been devised to retrieve different water optical parameters, including water clarity, from satellite observations for coastal (ocean) and lake waters (Attila et al., 2013;Binding et al., 2007;Binding et al., 2015;Olmanson et al., 2013;Potes et al., 2012;Wu et al., 2009;Zhao et al., 2011;Zolfaghari and Duguay, 2016).Turbid inland and coastal waters are optically more complex compared to open ocean, and large optical gradients exist.There is more than only one component (phytoplankton species, various dissolved and suspended matters with non-covarying concentrations) in coastal waters and lakes that determines the variability of water-leaving reflectance.Considering this complexity, the development of algorithms for coastal waters and lakes is more challenging.MERIS, which operated from March 2002 to April 2012, collected data from the European Space Agency's (ESA) Envisat satellite.The spatial resolution and spectral bands settings were carefully selected in order to meet the primary objectives of the mission; addressing coastal monitoring from space.The best possible signal-to-noise ratio, additional channels to measure optical signatures, and the relatively high spatial resolution of 300 m are some of the specific instrument characteristics (Ruescas et al., 2014).In 2010, ESA launched the CoastColour (CC) project to fully exploit the potential of MERIS instrument for remote sensing of coastal zone waters.CoastColour is providing a global data set of MERIS full resolution data of coastal zones that are processed with the best possible regional algorithms to produce water-leaving reflectance and optical properties (Ruescas et al., 2014).
The objectives of this study were to: (1) evaluate satellitederived K d values for a large lake in the Great Lakes region; (2) apply the evaluated satellite-derived K d in FLake model to investigate the improvement of model performance to reproduce LSWTs, compared to previous studies using a constant K d value of 0.2 m −1 .Therefore, three different values of K d were used in the simulations: yearly average, monthly average, and a constant value of 0.2 m −1 to evaluate the impact of a time-independent, lake-specific K d value in simulating LSWT; and (3) understand the sensitivity of the FLake model to variations in K d , based on the analysis of simulated LSWT, mean water column temperature (MWCT), lake bottom water temperature (LBWT), mixed layer depth (MLD), and water temperature isotherms during the ice-free season on Lake Erie (from April to November).The impact of K d variations on ice dates (freeze-up, break-up, and duration) and ice thickness was also evaluated.

Study site and station observations
Lake Erie (42 • 11 N, 81 • 15 W; Fig. 1) is a large shallow temperate freshwater lake covering a surface area of 25 700 km 2 .The lake is characterized by three basins: shallow western, central, and deep eastern basins with maximum depths of 19, 25, and 64 m, respectively.Lake Erie is monomictic with occasional dimictic years (Bootsma and Hecky, 2003).It is the shallowest and smallest by volume of the Laurentian Great Lakes (Daher, 2000).These characteristics make Lake Erie unique from the other Great Lakes.
The meteorological forcing variables required for FLake model runs include solar (shortwave) and longwave irradiance, air temperature, air humidity, wind speed, and cloudiness.These data were collected from different stations shown in Fig. 1.Mean daily air temperature, wind speed and water temperature measurements were obtained for the years 2003-2012, from the National Data Buoy Center (NDBC) of NOAA, station 45005 (41 • 40 N, 82 • 23 W, and depth: 12.6 m).Air temperature is measured 4 m above the water surface and anemometer height is 5 m above the water surface to measure the wind speed, whereas the water surface is at 173.9 m above mean sea level.Water temperature is measured at 0.6 m below the water surface.The NDBC station was selected to perform simulations with FLake, since water temperature observations collected at the buoy station can be used to evaluate the model output.The other meteorological forcing variables required for model simulations at the NDBC station were obtained from nearby stations.Air humidity and cloudiness were available in a daily format from the Ontario Climate Center (OCC), Environment and Climate Change Canada (ECCC) for the Windsor station (climate ID: 6139525) (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).This station is a near-shore station close to the NDBC station.The distance between the OCC and NDBC stations is less than 81 km.Incoming radiation flux data were supplied by the National Water Research Institute (NWRI), ECCC, from a station located in the western basin of Lake Erie.Daily shortwave irradiance measurements were available only for 2004 and 2008.Therefore, a daily time series of solar irradiance for the entire study period (2003-2012) was completed for the NDBC station using solar irradiance model data (see Sect. 2.2).Longwave irradiance was measured only in 2008 at the NWRI station.An empirical equation (see Sect. 2.2) was therefore employed to obtain longwave irradiance for the full period of study (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).
FLake requires information on water transparency (downward light K d ) as input for model runs.MERIS satellite imagery was used to derive K d for the NDBC station during the study period.The method is described in details in Sect.2.3.Available Secchi disk depth (SDD) field measurements were collected by ECCC research cruises on board the Canadian Coast Guard ship Limnos and utilized in this study to evaluate the satellite-derived water clarity.The cruise visited Lake Erie at a total of 89 distributed stations in five different years (September 2004;May, July, and September 2005;May and June 2008;July and September 2011;and February 2012).

Shortwave and longwave irradiance
The SUNY model, a satellite solar irradiance model, has been developed to exploit Geostationary Operational Environmental Satellites (GOES) for deriving solar irradiance using cloud, albedo, elevation, temperature, and wind speed observations (Kleissl et al., 2013).The basic principles of solar-irradiance modeling based on inputs from geo- stationary satellites and atmospheric models are described in Kleissl et al. (2013).Data from these sources are used to generate site-and time-specific high-resolution maps of solar irradiance with the SUNY model.The daily mean solar irradiance data for the present study was obtained from the second version of the SUNY model (Version 2.4), available in SolarAnywhere ® (https://www.solaranywhere.com).The model provides a gridded data set with a spatial resolution of 1/10 of a degree (ca. 10 km).The solar irradiance data was extracted from a tile corresponding to the NWRI station for 2004 and 2008, when observations were available for evaluation, and also for FLake model run on Lake Erie for the full study period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).There is a strong agreement (R 2 = 0.93) between model-derived and measured solar irradiance at the NWRI station.The SUNY model slightly underestimates observations by 2.18 W m −2 (N = 362, RMSE = 21.58W m −2 , MBE = −2.18W m −2 , I a = 0.88; see Sect.2.5 for details).
Longwave irradiance was computed on a daily basis using the equation of Maykut and Church (1973), as implemented in the Canadian Lake Ice Model (CLIMo) (Duguay et al., 2003): where T is the air temperature at screen height ( • K) and G is the cloudiness in tenth from meteorological stations.Longwave irradiance calculated from Eq. ( 1) was evaluated against observations from the NWRI station, only available in 2008.The two data sets are highly correlated (R 2 = 0.74) with the equation underestimating measured irradiance by 0.86 W m −2 (N = 194, RMSE = 17.74W m −2 , MBE = −0.86W m −2 , I a = 0.76).Model-derived incoming shortwave and longwave fluxes were used as input in FLake model simulations for subsequent analyses in the NDBC station over the 2003-2012 period.

Satellite-derived extinction coefficient
MERIS operated on-board the ESA Envisat polar-orbiting satellite until April 2012.The sensor was a push-broom imaging spectrometer which measured solar radiation reflected from the Earth's surface at high spectral and radiometric resolutions with a dual spatial resolution (300 and 1200 m).Measurements were obtained in the visible and near-infrared part of the electromagnetic spectrum (across the 390-1040 nm range) in 15 spectral bands during daytime, whenever illumination conditions were suitable, and with a full spatial resolution of 300 m at nadir, with a 68.5 • field of view.MERIS scanned the Earth with a global coverage of every 2-3 days.
In this study, a total of 326 full-resolution archived MERIS images encompassing the NDBC station in Lake Erie were acquired from CC (Version 2) products through the Calvalus on-demand processing service for the period of 2003-2012.The level 2 products are generally geolocated geophysical products and CC Level2W products are the result of in-water processing algorithms to derive optical parameters from the water-leaving reflectance.These parameters include inherent optical properties (IOPs), concentrations of some water constituents, and other optical water parameters such as spectral vertical K d .The IOP parameters are first derived applying two different inversion algorithms: neural network (NN) and quasi-analytical algorithm (QAA).The derived IOPs are then converted to estimate constituents' concentrations and apparent optical properties (AOPs), including diffuse K d for different spectral bands applying Hydrolight simulations (Ruescas et al., 2014).
The diffuse K d product (the average value between visible spectral bands) in CC Level2W data was extracted for the pixel at the geographic location of the NDBC station.The satellite-derived K d values were also extracted for pixels on the same day and location as the Limnos cruise stations, to evaluate the CC-derived diffuse K d values against SDD in situ data collected during Limnos cruises.A valid pixel expression was defined in all pixel extraction steps that excluded pixels with properties listed in Table 1.

FLake model and configuration
The FLake model is a self-similar parametric representation (assumed shape) of the temperature structure in the four media of the lake including water column, bottom sediments, and the ice and snow.The water column temperature profile is assumed to have two layers: a mixed layer with constant temperature and a thermocline that extends from the base of mixed layer to the lake depth.The shape of thermocline temperature is parameterized using a fourth-order polynomial function of depth that also depends on a shape coefficient C T .The value of C T lies between 0.5 and 0.8 so that the thermocline can neither be very concave nor very convex.FLake has an optional scheme for the representation of bottom sediments layer, which is based on the same parametric concept (De Bruijn et al., 2014;Martynov et al., 2012).The system of prognostic equations for parameters is described in Mironov (2008).
The prognostic ordinary differential equations are solved to estimate the thermocline shape coefficient, the mixed layer depth, the bottom, mean and surface water column temperatures, and also parameters related to the bottom sediment layers (Martynov et al., 2012;Mironov, 2008;Mironov et al., 2010).The same parametric concept is applied for the ice and snow layers, using linear shape functions (Martynov et al., 2012).The mixed layer depth is calculated considering the effects of both convective and mechanical mixing, also accounting for volumetric heating which is through the absorption of net shortwave radiation (Thiery et al., 2014).The non-reflected shortwave radiation is absorbed after penetrating the water column in accordance with the Beer-Lambert law (Gordon, 1989).
Stand-alone FLake simulations were conducted for the NDBC station.The setup condition of the NDBC buoy station, such as height of wind measurement (5 m), height of air temperature sensor (4 m), and the geographic location and depth of this site were used to configure the model.The measured meteorological parameters and model-derived irradiance were also used to force the FLake model.A fetch value of 100 km was used to run all simulations.It was found that there is only a little sensitivity to modifications in this parameter for Lake Erie.The same result was found for Lake Kivu in Thiery et al. (2014).The bottom sediments module was switched off in all simulations and the zero bottom heat flux condition was adopted.The initial temperature value for the upper mixed layer and the lake bottom were 4 • C. Mixed layer thickness had the initial value of 3 m.The simulations were run in a daily time step (using daily forcing data) for 2003-2012.
The ability of FLake to reproduce the observed temperature variations using different K d values was tested by comparing the simulated LSWT to the corresponding in situ observations in the NDBC station.Also, the model sensitivity to variations in water clarity was assessed by studying the LSWT, MWCT, LBWT, MLD, isotherms, ice phenology, and ice thickness.

Accuracy assessment
To assess the model outputs, three statistical indices were calculated: the root mean square error (RMSE), the mean bias error (MBE), and the index of agreement (I a ).RMSE is a comprehensive metric that combines the mean and variance of model errors into a single statistic (Moore et al., 2014).
The MBE is calculated as the mean of modeled values minus the in situ observations.Therefore, a positive (negative) value of this error shows an overestimation (underestimation) of the parameter of interest.I a is a descriptive measure of model performance.It is used to compare different models and is also modeled against observed parameters.I a was originally developed by Willmott in the 1980s (Willmott, 1981) and a refined version of it was presented by Willmott et al. (2012).The refined version, which was adopted in this study, is dimensionless and bounded by −1.0 (worst performance) and 1.0 (the best possible performance).These statistical indices are considered to be robust measures of model performance (e.g., Hinzman et al., 1998;Kheyrollah Pour et al., 2012;Willmott and Wicks, 1980).Hence, these values, identified as the average, the lower, and the upper limits of clarity at the NDBC station, were used to carry out a sensitivity analysis with FLake (see Sect. 3.2.2).

Evaluation of CoastColour K d
The validation of satellite observations against in situ data is important, because the in situ data are still considered as the most accurate measurement of water clarity.The assessment of the satellite-derived K d retrieval reliability highly depends on the comparison with independent in situ SDD measurements.The general form of the relationship between K d and SDD was established by the pioneer study of Poole and Atkins (1929): where K is a constant value of 1.7 (Poole and Atkins, 1929).Following this important work, there were other studies that found a high variability of the constant value (K) depending on the type of the lake considered (Koenings and Edmundson, 1991).Armengol et al. (2003) showed that K d and SDD are negatively correlated and they developed an empirical power relation between these two parameters.In this study, applying a cross-validation approach, an empirical relation was developed between in-situ-measured SDD and CC-derived K d .SDD measurements were conducted 117 times during cruises on Lake Erie from 2004 to 2012.These spatially distributed measurements had minimum, maximum, mean, and standard deviation values of 0.2, 11, 3.69, and 2.68 m, respectively.CC Level2W satellite products were acquired on the same day as the in situ measurements.Applying defined flags produced 49 data pairs (matchup data set) of CC observations of K d and SDD in situ data that were collected on the same day and location.
The matchup data set was divided into training and testing data in 100 iterations.In each iteration, the data used for the equation's training and evaluation were kept independent, where 70 % of the sample was used for equation calibration and 30 % for evaluation.Ordinary least-square regression was used in the calibration step of each iteration to relate the in situ measurements of SDD to the CC-derived K d .Locally tuned equations were derived from this step and applied on SDD observations to predict K d in testing matchup data.The statistical parameters of the model performance were derived between the estimated K d from SDD observations and the paired CC-derived values.These steps were repeated for 100 iterations, and the final statistical indices, slope, and power of the locally tuned equation was reported as the average of those derived over all iterations.
Results from the above procedure show that K d can be derived from SDD, using the equation K d = 1.64 × SDD −0.76 , with a strong determination of coefficient value (R 2 = 0.78).Arst et al. (2008) obtained a similar regression formula between SDD and K d for the boreal lakes in Finland and Estonia representing different types of water, expanding from oligotrophic to hypertrophic.Because there is a good agreement between K d and the corresponding ones estimated from in-situ-measured SDD (N = 49, RMSE = 0.63 m −1 , MBE = −0.09m −1 , I a = 0.65; Fig. 3), the satellite-derived water clarity measurements were considered to be representative of K d and were used in the modeling for this study.
However, SDD is not always describing K d values.SDD is a suitable characteristic to describe water transparency for small values of K d .For high values of K d (ranging above 4 m −1 ), Arst et al. (2008) and Heiskanen et al. (2015) suggested that SDD is unable to describe any changes in K d .Figure 3 also shows that SDD cannot describe the scatter of K d for values above 4 m −1 .Therefore, the estimation of K d from in situ measurements of SDD should be used with caution.Direct measurements of K d in the field are not widely available.These limitations motivate the investigation into the potential of integrating satellite-based estimations of K d into lake models (Arst et al., 2008).

Improvement of LSWT simulations with
satellite-derived K d Martynov et al. (2012) focused on 2005-2007 to run FLake at the NDBC station using a constant value of 0.2 m −1 for K d .They simulated the lake properties using both realistic and excessive depths of 20 and 60 m, respectively, for a grid tile corresponding to the NDBC station.They showed that apply-      The overall performance of each simulation is summarized in Table 3 during the period of data availability.For all years, the average and merged simulations perform better than simulations using K d (0.2 m −1 ) as in Martynov et al. (2012), with improvement in RMSE and MBE for both real depth and tile depth.In all 3 years, LSWT simulated from the K d value employed in Martynov et al. (2012) resulted in an underestimation (CRCM-12.6:MBE = −1.52,−0.98, −1.08 • C; CRCM-20: MBE = −1.54,−1.09, −1.35 • C; during years 2005, 2006, and 2007, respectively).In 2005, the average of K d for the year demonstrates a better performance compared to the merged results, contrary to the results of 2007.However, for the merged results in 2006, the MBE was improved compared to the simulation using the average K d , whereas its performance decreased in terms of RMSE.The extent of K d variations in each month might not be captured by the available MERIS images due to cloud coverage in MERIS images or the absence of any satellite overpass.Therefore, a yearly-average K d can be potentially closer to the actual value of K d .For this reason, the merged results cannot always perform better than average simulations.
Figure 5 illustrates the scatterplots of simulated LSWT for all four different runs including 3 years of data (2005)(2006)(2007), in comparison with the corresponding in situ observations.All simulated results were in a high agreement with in situ measurements.Figure 5a  ).Under-prediction of LSWT decreased when the yearly-average CC-derived K d values were used, rather than a generic constant value (0.2 m −1 ).Heiskanen et al. (2015) suggested that the effect of K d seasonal variations on LSWT simulations are not significant for lakes with K d values higher than 0.5 m −1 (e.g., Lake Erie).Therefore, in the absence of reliable values of the temporal evolution of K d , a lake-specific, time-independent, and constant value of K d can be used in 1-D lake models when the K d values are higher than 0.5 m −1 .Martynov et al. (2012) concluded that applying a more realistic lake depth parameterization improves the FLake model performance.Using the realistic lake depth (12.6 m) at the NDBC station slightly improves the model performance in reproducing LSWT compared to a simulation employing the corresponding tile depth (20 m) (CRCM-12.6:5c, d).

Sensitivity of FLake to K d variations
The sensitivity of FLake to different values of K d to reproduce LSWT, MWCT, LBWT, MLD, isotherm, and ice phenology and thickness was investigated in this section for the year 2008.As indicated previously (Sect.2.1), shortwave irradiance measurements were available in that year and longwave irradiance was also measured from May to October 2008.Therefore, longwave irradiance for the other months of 2008 was modeled as described in Sect.2.2 to fill the temporal gaps.Figure 6 presents simulation results for LSWT, MWCT, and LBWT using the real lake depth at the NDBC station, and the lowest, average, and highest values of K d observed in the study period (minimum K d = 0.58 m −1 , average K d = 0.90 m −1 , maximum K d = 3.54 m −1 ).The water temperature simulation from CRCM-12.6 (using K d = 0.2 and realistic depth at station) simulation was also plotted.
In the case of extreme clear water (CRCM-12.6),LSWT showed smoother variations during the open water season in 2008 as opposed to the darkest water simulation (maximum, or Max) which displayed more abrupt LSWT variations (Fig. 6).This is because solar radiation is absorbed more in waters with low clarity due to existing particles in water.It penetrates less deeply and warms up only the shallow surface layer (which shows in lower LBWT; see Fig. 6c) causing thinner mixing depth (Fig. 6d).The high temperature of this shallow layer causes an increase in latent and sensible heat fluxes.Therefore, the shallow mixed layer exchanges heat faster with the atmosphere, resulting in sudden surface water temperature variations as opposed to clear waters.The fast heat exchange with atmosphere resulted in warmer LSWT during spring (start of heating season) and colder LSWT in fall for dark water as opposed to clear water.On average, the darkest water simulation (Max) resulted in 0.09 • C higher LSWT compared to the average (Avg) sim-  sults showed that FLake-simulated LSWT was not significantly sensitive to K d values when this value varied in the range of our Min to Max K d .However, the sensitivity increased rapidly for K d values less than our Min (0.58 m −1 ).This result supported the study of Rinke et al. (2010) that the thermal structure of lakes is particularly sensitive to changes in K d when its value is below 0.5 m −1 .More recently, Heiskanen et al. (2015) confirmed the critical threshold of K d (ca.0.5 m −1 ).They suggested that the response of 1-D lake models to K d variations is nonlinear.The models are much more sensitive if the water is estimated to be too clear.Heiskanen et al. (2015) recommended to use a value of K d that is too high rather than too low in lake simulations, if the clarity of lake is not known exactly.
The MWCT and LBWT in the darkest condition (Max) were less than for all other clear-water simulations.This is because the lower layers in dark waters accumulate less heat during the heating season as opposed to clear waters, which results in less heat storage and lower water column temperature in dark waters (Heiskanen et al., 2015;Potes et al., 2012).The MWCT decreased by 0.94 • C (increased by 0.63 • C) when maximum (minimum) K d value was used instead of its average value during the study period.The MWCT increased by 2.25 • C when using a K d value of 0.2 m −1 rather than the average value.Changes in K d value from its maximum (minimum) to its average value also caused a decrease (increase) of −0.67 • C (0.67 • C) in the LBWT.The increase in LWBT was even larger when the K d value of 0.2 m −1 was used instead of its average value (6.96 • C).Therefore, K d variations had a larger impact on MWCT and LBWT than on LSWT, and the largest difference was when K d was estimated to be extremely clear.
Figure 7 displays the simulated isotherms derived from using different K d values.It shows that the mixed layer in dark waters was warmer in spring and summer and colder in fall.There are a number of factors determining the mixed-layer temperature in lakes, including the radiation fluxes (sensible heat, latent heat, and longwave radiation), and cooling effects from the water below.Persson and Jones (2008) concluded that for dark waters, the combination of these heating and cooling effects leads to a warmer epilimnion initially.The radiation is used to warm up a thinner layer in dark waters leading to higher (lower) temperatures in spring and summer (fall).However, a lower temperature in the mixed layer is followed due to the gradual decrease in radiative forcing and increased effect of cooling from the layers below.Figure 7 also supports observations by Persson and Jones (2008) and Heiskanen et al. (2015) that the depth of the thermocline layer is always deeper in clear waters due to the faster heat distribution between different underneath layers.The deepening of the thermocline layer in clear waters is faster compared to dark waters.The reason is related to heat transfer through convection, wind-induced mixing, and internal waves.The heat transfer in dark waters is slower due to the sharp density gradient between layers which forms an effective barrier for the mixing to deepen the thermocline.
Figure 6d is focusing on the variations of the MLD in 2008, using different values of K d (Min, Avg, and Max K d , and CRCM-12.6) in simulations.All simulations showed two turnover (complete mixing) events, spring and fall.Full mixing in spring was at the same time for all simulations; however, fall full mixing occurred at different dates for each simulation.Fall turnover in CRCM-12.6 was at the end of summer (28 August), while the other three runs show that the fall turnover took place in late fall, before ice forms.Full mixing in the Min simulation was in early November (3 November), earlier than the Avg and Max simulations (21 November).
In the darkest water simulation (Max), the MLD was shallower than the other simulations (an average difference of 4.94 m in 2008 between two simulations of Max and CRCM-12.6, with extreme K d values).Clear waters have a deeper mixed layer when the solar radiation can penetrate further and distribute to a larger volume in the water column.Also, due to the weak density gradient in clear waters, wind-induced turbulent kinetic energy can destroy the density stratification to a deeper layer and form the mixed layer.This layer is shallower in dark waters, even with the same wind forcing.CRCM-12.6 produced a MLD of 3.47 m deeper compared to Avg simulation, whereas the Min (Max) simulations resulted in MLD 1.15 m (1.47 m) deeper (shallower) compared to the Avg simulation.Hence, clear water simulated deeper MLD; and the effect of K d on the MLD was larger when the K d value was estimated to be too clear.
Figure 8 shows the impact of K d variations on lake ice phenology and thickness in winter 2008 (January-March).Freeze-up corresponds to the earliest date that the NDBC station is completely covered by ice, and the earliest date the station is completely free of floating ice is defined as break-up.The Avg simulation reproduced similar ice phenol- ogy as the Max simulation, whereas Min and CRCM-12.6 resulted in the similar break-up and freeze-up dates.The breakups in CRCM-12.6 and Min simulations were on 23 March, 1 day earlier than Max and Avg simulations, and freeze-up occurred on 24 January, 2 days after Max and Avg simulations.CRCM-12.6 and Min simulations reproduced 1.28 and 1.27 cm thinner ice than Avg simulation in 2008, respectively.The darkest water (Max) reproduced 0.21 cm thicker ice in 2008 compared to the Avg simulation.The ice sheet formed later in clear waters (CRCM-12.6and Min) and disappeared earlier compared to dark waters (Max and Avg), resulting in a shorter ice-cover duration (3 days) and hence thinner ice in clear-water simulations.
Lake morphological properties determine ice cover as well as climatic factors.Among morphological aspects, lake depth is the most important factor that can impact the ice cover by influencing the amount of heat storage in the water and hence the time needed for the lake to cool and ultimately freeze (Brown and Duguay, 2010).For a given depth and climatic condition, however, the amount of heat storage is determined by water clarity.Dark waters store more heat in a shallower layer.Therefore, the heat can be transferred faster to the atmosphere through the lake surface, resulting in an earlier freeze-up as mentioned in Heiskanen et al. (2015), who reported that freeze-up occurs earlier in darker waters.However, as shown by simulations with 12.6 m, ice phenology in the NDBC station was minimally affected by K d value in FLake.It must be noted that these results could not be verified due to the lack of ice phenology observations.For a larger depth or in a different model, the impact of K d values in ice onset should be investigated.

Spatial and temporal variations in K d
As described in the previous section, variations in water clarity play an important role in defining lake heat budget and thermal stratification and thus is a significant parameter for processes in the air-water interface.However, the long-term spatial and temporal trends of water clarity cannot be achieved through discontinuous conventional pointwise in situ sampling.These observations can be provided from satellite measurements.This section demonstrates the strength of satellite observations to detect the spatial and temporal variations of K d in Lake Erie.Spatial variations of K d derived from the CC algorithm are shown in Fig. 9 for a selected day (3 September 2011).This particular day of 2011 was selected as the lake experienced its largest algal bloom in its recorded history in that year, before the new recent record of 2015 (Michalak et al., 2013;NOAA, 2015).The bloom was expanding from the western basin into the central basin.Algal bloom is one of the factors affecting the water clarity of Lake Erie (NOAA, 2015).Other parameters include the concentrations of suspended and dissolved matters in the lake.The western basin is the shallowest region of the lake, and therefore is the most vulnerable to sediment re-suspension that also results in reducing water clarity.The map shows that Lake Erie experienced different levels of clarity in various locations, with an average K d value of 0.90 m −1 (with standard deviation of 0.80 m −1 , shown as 0.90 ± 0.80 m −1 hereinafter) over the entire lake on this particular day.The NDBC station is also shown on the satellite-derived map as a reference (with K d = 0.87 m −1 on 3 September 2011).
Since fully cloud-free MERIS satellite images for consecutive months were only available in 2010, 4 months (May-August 2010) were selected to illustrate temporal variations in K d on a monthly basis for 1 selected year (Fig. 10).The spatial averages of K d over the full lake for the specific days in May, June, July, and August were 0.82 ± 0.85 m −1 , 0.72 ± 1.10 m −1 , 0.73 ± 1.20 m −1 , 0.78 ± 0.55 m −1 , respectively.The western basin was always experiencing the lowest levels of water clarity in comparison to other regions of the lake, with a maximum K d in May.This can be the result  of a spring algal bloom, and also wind-driven re-suspension of sediments.K d at the NDBC station for these selected days varied between 0.68, 0.62, 0.66, and 0.85 m −1 from the months of May to August 2010, respectively.
Two MERIS images with full coverage of Lake Erie were only available in the month of May for 2 selected consecutive years (2008 and 2009) to show the inter-annual changes in K d value.Hence, the MERIS images of May 2008 and May 2009 were selected to show variations in K d between the two years.Although the images are for the same month of the year, K d still varied across the lake (Fig. 11).In the selected day of May 2008, a spatial average value of 0.77 ± 0.49 m −1 was estimated for the entire lake, while on 17 May 2009 the spatial average value was 0.90 ± 0.93 m −1 .Comparing the estimated maps for the two years suggested that the spring bloom in 2009 was stronger than the one in 2008 for the western basin.However, algal bloom in all basins of Lake Erie for the complete year of 2008 was recorded as the third-largest that the lake experienced before the occurrence of the breaking record blooms in 2011 and 2015 (Michalak et al., 2013;NOAA, 2015).K d values estimated for the NDBC station were 0.69 and 0.62 m −1 on 29 May 2008 and 17 May 2009, respectively.
Spatial variability of K d in Lake Erie shows that the simulated thermal structure of the eastern basin would potentially differ significantly from the one simulated for the western basin.The spatial variations of K d have to be considered in Lake Erie simulations, specifically for the eastern basin, which has K d values in the critical threshold range (less than 0.5 m −1 ).Therefore, in 3-D lake models, the spatial variations in K d need to be taken into account.As well as this, a lake-specific constant value cannot be used for simulating the thermal structure of the full lake.Finally, the temporal variations of K d did not significantly change the simulation results for the NDBC station.However, this needs to be confirmed for other locations of the lake, due to the importance of depth on the simulation results.

Summary and conclusion
Spatial and temporal variations of K d in Lake Erie were derived from the globally available satellite-based CC product during open water seasons 2003-2012.The CC product was evaluated against SDD in situ measurements.CC-derived K d values and modeled incoming radiation flux data, in addition to complementary meteorological observations during the study period, were used to force the 1-D FLake model.The model was run for a selected site (NDBC buoy station) on Lake Erie, a large shallow temperate freshwater lake.
FLake was run with the range of clarity values acquired from satellite observations.Results were compared to a pre-Hydrol.Earth Syst.Sci., 21, 377-391, 2017 www.hydrol-earth-syst-sci.net/21/377/2017/ vious study which assumed a constant K d value due to the lack of data.Results clearly showed that applying satellitederived K d values improves FLake model simulations using a derived yearly average value as well as monthly averaged values of K d .Although K d varies in time, a time-invariant (constant) annual value is sufficient for obtaining reliable estimates of lake surface water temperature (LSWT) with FLake for Lake Erie NDBC station.It was also shown that the model is very sensitive to variations in K d when the value is less than 0.5 m −1 .This finding is in agreement with the study of Rinke et al. (2010) and the recent study of Heiskanen et al. (2015) who determined that the impact of seasonal variations of K d on the simulated thermal structure is small, for a lake with K d values larger than 0.5 m −1 .The studies suggested that the response of 1-D lake models to K d variations is nonlinear.The models are much more sensitive if the water is estimated to be too clear.The results of our study showed that the sensitivity to K d variations was more pronounced in simulation results for mean water column temperature (MWCT), lake bottom water temperature (LBWT), and mixed layer depth (MLD) compared to LSWT.
Results of this study have important implications for the lake modeling community, demonstrating that integrating satellite-derived lake-specific K d values can improve the performance of 1-D lake models compared to using a "generic" constant K d value.Although field measurements of K d are not widely available, this study evaluated the strength of satellite observations and introduces them as a reliable data source to provide lake models with global estimates of K d with high spatial and temporal resolutions.However, the weakness of this method is that the availability of satellitederived K d product can be limited due to cloud coverage or satellite overpass.Also, the in situ measurements are still required for validating satellite observations, because the in situ data collection remains the most accurate solution for water clarity measurement.The accuracy of the satellitederived K d product has to be verified for the water body of interest, especially for the ones with complex optical properties.After validation, the on-demand globally available CC product can be simply used for the water body of interest, as a source to fill the gaps in K d in situ observations, and improve the performance of parameterization schemes and, as a result, further improve the NWP and climate models.Although MERIS is no longer active, the Ocean and Land Colour Instrument (OLCI), to be operated on the ESA Sentinel-3 satellite (launched on 16 February 2016), will provide continuity of MERIS-like data.OLCI has MERIS heritages and improves upon it with an additional six spectral bands.Therefore, investigation of the Sentinel-3 potential to provide lake modeling community with the water clarity information is the next step of the current study.Also, the possible improvement in FLake output, when forcing the model with air humidity data collected directly at the station, can be examined in the future studies.

Figure 1 .
Figure 1.Maps showing Lake Erie in Laurentian Great Lakes and the location of stations where different parameters were measured.NDBC: National Data Buoy Center.NWRI: National Water Research Institute.OCC: Ontario Climate Center.CCGS: Canadian Coast Guard Ship.Vertical dashed lines separate different basins in the lake.

3
Figure2shows the variations of CC-derived K d for the NDBC station during the full study period(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).Lake Erie (specifically its shallow regions) is more susceptible to re-suspension of bottom sediments compared to the other Great Lakes, which leads to lower water clarity(Binding et al., 2010).The results from applying the CC algorithm on MERIS satellite imagery shows that the highest K d values in the NDBC station occur in the turn-over times in spring and fall.The maximum value of K d was 3.54 m −1 , estimated in April 2003.A minimum value of 0.58 m −1 was estimated in June 2007.The average value of K d during the study period was 0.90 m −1 with a standard deviation of 0.38 m −1 .Hence, these values, identified as the average, the lower, and the upper limits of clarity at the NDBC station, were used to carry out a sensitivity analysis with FLake (see Sect. 3.2.2).

Figure 2 .
Figure 2. Variations of CoastColour-derived K d for the selected location during the study period (2003-2012).

Figure 3 .
Figure 3. Relation between satellite-derived K d and in situ SDD matchups.

Figure 4
compares the results of different LSWT FLake simulations with observations at the NDBC station.LSWT observations had maximum values of 27.53, 26.48, and 25.46 • C in August during 2005, 2006 and 2007, respectively.The minimum values of 2.71, 7.3, and 3.42 • C were observed in December 2005, and April in 2006 and 2007, respectively.The average LSWT observations in 2005, 2006, and 2007 had values of 18.45, 17.12, and 17.75 • C, respectively.Four different simulation schemes were made which were then compared to the observed LSWT.The simulated LSWT values in Fig. 4 were produced by first applying K d = 0.2 m −1 from Martynov et al. (2012) using both the real lake depth at the station (12.6 m: CRCM-12.6) and also a tile depth corresponding to the station in their study (20 m: CRCM-20).Then, simulations using the yearly average CCderived K d for each year of study were plotted (Avg).The K d values derived from the monthly average of each year were used to simulate the surface water temperature and produce a merged LSWT product (Merged).Both Avg and Merged simulations used the real lake depth at the NDBC station (12.6 m).Comparing LSWT in situ observations (Obs) with the modeled values in Fig. 4 demonstrated that in Avg and Merged simulations for 2005-2007, surface temperature was modeled warmer in spring (April-June) and colder in summer (July-September) and fall (October-November) than in situ observations (spring: MBE Avg = 1.31 • C, MBE Merged = 1.25 • C; summer: MBE Avg = −0.72 • C; Hydrol.Earth Syst.Sci., 21, 377-391, 2017 www.hydrol-earth-syst-sci.net/21/377/2017/

Figure 4 .
Figure 4. Daily LSWT simulation results in 2005 (a), 2006 (b), 2007 (c).Avg.simulation is the CoastColour-derived average value for K d during selected months of each year (0.81, 0.71, and 0.73 m −1 , respectively).Merged simulation is based on merging simulation results for monthly average values of K d .CRCM-12.6 and CRCM-20 used a constant value of K d (0.2 m −1 ) with depth values of 12.6 and 20 m, respectively.The corresponding observations for LSWT are also plotted.Missing lines indicate no data.

Figure 5 .
Figure 5. Modeled (y axis) versus observed (x axis) LSWT for yearly average, merged, CRCM-12.6, and CRCM-20 simulations during the ice-free seasons in 2005-2007.A linear fit (dashed line) and its coefficients are shown on the plot.The statistics related to the regression of parameters, and a 1 : 1 relationship (solid line) are also shown.The average LSWT values of Obs, Avg, Merged 18.56, 18.50, 17.38, and 17.27  • C, respectively.
and b show that the resulting LSWT from yearly average (Avg) and monthly average (Merged) K d were not significantly different, whereas simulations with yearly average K d reproduced LSWT with improved RMSE and MBE values compared to monthly averages (Avg: RMSE = 1.54 • C, MBE = −0.08 • C; Merged: RMSE = 1.57• C, MBE = −0.14• C).It is possible that the actual K d value is best represented by the yearly average value.Therefore, using a constant annual open water season value for K d could be potentially sufficient to simulate LSWT in 1-D lake models with relatively high accuracy (the range of K d variations that brings the most sensitivity for the modeling is discussed in Sect.3.2.2).Both CRCM simulations (Fig. 5c: depths of 12.6 and Fig. 5d: depth of 20 m) under-predicted LSWT (for LSWT values larger than ca.7 • C), with MBE values of −1.26 and −1.37 • C, respectively.The under-prediction of these model runs was stronger, particularly for LSWT above 12 • C, which can be explained by the K d value used.This is because, no matter what depth is used in simulations (either actual or tile depth), both CRCM runs have larger MBE compared to Avg and Merged simulations.However, the CRCM-20 simulation tended to produce the coldest LSWT (the most underpredicted; MBE = −1.37 • C).This is due to the lake depth value considered for the model run, which corresponds to the tile depth as opposed to the other simulations that were based on using the actual depth at the station.The time-dependent (monthly average) K d did not improve simulation results for Lake Erie (K d ranging from 0.58 to 3.54 m −1 with an average value of 0.90 m −1 during open water seasons of 2003-2012).However, comparing results from Fig. 5a and c showed improvement in LSWT simulations when a lake-specific value of K d is used Hydrol.Earth Syst.Sci., 21, 377-391, 2017 www.hydrol-earth-syst-sci.net/21/377/2017/ (Avg: RMSE = 1.54 • C, MBE = −0.08 • C; CRCM-12.6:RMSE = 1.76 • C, MBE = −1.26• C

Figure 9 .
Figure 9. Spatial variation of satellite-derived K d in Lake Erie, on 3 September 2011.Location of NDBC station is shown on the map as a solid dot.

Figure 10 .
Figure 10.Temporal and spatial variation of satellite-derived K d in Lake Erie for different months of a year: May-August 2010.Location of the NDBC station is shown on the map as a solid dot.

Figure 11 .
Figure 11.Temporal and spatial variation of K d in Lake Erie during May of 2 consecutive years: 2008 and 2009.Location of the NDBC station is shown on the map as a solid dot.

Table 1 .
Flags of excluded pixels.

Table 2 .
CC-derived average values of K d for each month(2005)(2006)(2007).The values correspond to the time of year when water LSWT observations, as well as the CC derived K d values, are available.

Table 3 .
Simulated LSWT compared to in situ observations (2005- 2007).Period corresponds to the time of year when LSWT and K d values were available.