Optimizing Lake Surface Water Temperature Simulations Over Large Lakes in China With FLake Model

Abstract

Lake surface water temperature (LSWT) is sensitive to climate change; however, simulated LSWT and its response to climate change remain uncertain. In this study, FLake, a one‐dimensional freshwater lake model, is optimized to simulate the LSWTs of 94 large lakes with surface areas greater than 100 km² in China. While most of these lakes are seasonally ice‐covered over the Tibetan Plateau, FLake with default parameters significantly underestimated LSWT in spring and winter and slightly overestimated LSWT in summer and autumn in seasonally ice‐covered lakes. We performed sensitivity experiments and calibration in the trial lake (Qinghai Lake). Then, parameter calibrations of three lake‐specific properties (albedo, lake mean depth and light extinction coefficient) were performed in all the studied lakes. The optimized FLake substantially improved the simulations of seasonal and interannual variations in LSWT. The root mean square error decreased from 3.64 ± 1.54°C to 1.97 ± 0.72°C, and the mean bias of 96% of the lakes decreased to less than 1°C. Our study showed that the optimized FLake can reproduce the temporal variations in LSWT across China with optimized parameters, providing the possibility to simulate and project the response of LSWT to rapid climate change.

Introduction

Lakes, as important natural resources, hold enormous amounts of liquid surface freshwater (Gleick, 1993), support global biodiversity and provide ecosystem services to humans (Xu et al., 2018). However, lakes are facing numerous threats due to climate change (Woolway et al., 2020). Evidences have emerged that lake temperature (O'Reilly et al., 2015; Verburg & Hecky, 2009), lake area (Yang & Lu, 2014; Zhang et al., 2017), lake water level (Chen et al., 2017; Lei et al., 2019), lake mixing regimes (Butcher et al., 2015; Woolway et al., 2017) and lake ice phenology (Cai et al., 2019; Sharma et al., 2019) are responsive to climate change. Among these phenomena, the most direct impact of climate change on lakes is changes in lake surface water temperature (LSWT).

Lake surface warming affects the internal processes of lake ecosystems by changing the physical and biogeochemical properties of lakes, such as more stable stratification, decreased lake mixing, increased algal growth, more methane emissions and decreased primary productivity (Bastviken et al., 2004; Carter & Schindler, 2012; O'Reilly et al., 2003; Paerl & Paul, 2012), which may affect lake ecosystem services (Jeppesen et al., 2014). Moreover, increases in LSWT will affect the hydrological process in a basin (Chen et al., 2017; Gronewold & Stow, 2014) and even affect the regional climate, leading to an increasing intensity of lake-effect precipitation (Steiner et al., 2013).

There are ∼2,600 lakes (>1 km²) in China, accounting for ∼1% of the total land area (Ma et al., 2011). The largest lake (Qinghai Lake) has a surface area of 4,267 km². The altitudes of lakes span a large range from 0 to 5,198 m. Seventy-five percent of lakes are seasonally ice-covered and concentrated in high-latitude or high-altitude regions (Ma et al., 2011). Millions of people live around lakes and depend on the ecosystem services of lakes. From the perspective of both lake ecosystem services and regional climate effects, it is necessary to analyze and predict the temporal and spatial variations in LSWT. Although the diversity and importance of Chinese lakes make them ideal regions to study LSWT, previous studies were mostly on a single lake or over parts of China (Su et al., 2019; Wu et al., 2019; Zhang et al., 2014), rendering the LSWT dynamics and its response to climate change uncertain.

The lake model is the most common and effective tool to study LSWT and lake-air interactions. During recent decades, research institutions have developed many one-dimensional lake models and three-dimensional lake models based on different key physical processes and parameterization schemes. Although three-dimensional lake models contain more comprehensive and complicated physical processes within lakes, both their computational cost and time efficiency constrain their application. Moreover, accurate parameters of three-dimensional lake models are difficult to obtain. Considering these limitations, one-dimensional lake models have been more widely used in previous studies (Lv et al., 2019; Turuncoglu et al., 2013; Wu et al., 2020). One-dimensional lake models assume zero lateral heat/water exchange. Studies have proven that one-dimensional lake models can meet the needs of a range of lake simulations (Fenocchi et al., 2018; Thiery et al., 2015). These common one-dimensional lake models include the freshwater lake model (FLake) based on self-similarity theory, the eddy-diffusive Hostetler model (Hostetler et al., 1993; Subin et al., 2012), the Mixed-Layer model (Goyette et al., 2000), and the $k-\epsilon$ turbulence closure model (Perroud et al., 2009). The Lake Model Intercomparison Project (LakeMIP) compared the performances of different lake models and revealed that no single model outperformed the others globally, with performances varying among lakes due to differences in lake parameters (Stepanenko et al., 2010, 2013).

FLake is a model capable of reproducing the surface thermal behavior and vertical temperature structure in lakes of various depths (Mironov, 2008; Mironov et al., 2010). In recent decades, many studies have evaluated the performances of FLake in reproducing lake thermodynamics (Lv et al., 2019; Rooney & Bornemann, 2013; Thiery et al., 2015). Despite the difference between simulations from FLake and observations (Bernhardt et al., 2012; Layden et al., 2016; Pour et al., 2012), FLake is thought to be a reliable model for global-scale studies after tuning trials (Layden et al., 2016; Woolway & Merchant, 2019).

In this study, we improved the performance of FLake in large lakes over China to accurately simulate the thermal behaviors of lake surfaces. First, we conducted preliminary trial work in Qinghai Lake, the largest lake in China, and set up a series of sensitivity experiments to calibrate key lake-specific parameters with an improved observation-driven method integrating both satellite observations and the model-tuning algorithm recommended by Layden et al. (2016). Based on the trial in Qinghai Lake, we further optimized FLake model performances in all 94 large lakes during 1996–2011. We organized the paper as follows: Section 2 describes the studied lakes, data sets, an overview of the FLake model, model modifications, and methods for model optimization. Section 3 compares and improves the model performance against satellite-derived data sets. Section 4 summarizes and discusses the results.

Materials and Methods Study Area

Ninety-four large lakes (>100 km²) with LSWT observations of the Along Track Scanning Radiometer Reprocessing for Climate: Lake Surface Water Temperature and Ice Cover (ARC-Lake) data set in China were chosen as study areas to calibrate and improve the FLake model performance. Lake boundaries were extracted from the HydroLAKES data set (Messager et al., 2016), as well as lake areas, lake depths and altitudes. These lakes account for nearly 50% of lake areas in China and vary in surface area (100.8–4,266.6 km²), depth (1–120 m), altitude (0–5,198 m), and geographic location (Figure 1). The studied lakes were divided into six subregions to demonstrate the simulations of LSWT according to geographical and climatic features (Ma et al., 2011; Zhang et al., 2019). These regions included the Tibetan Plateau (TP), Xinjiang (XJ), Inner-Mongolian Plateau (MP), Northeast Plain and Mountain (NE), Eastern Plain (EP), and Yunnan-Guizhou Plateau (YG). The lake numbers in the subregions are 58, 9, 3, 1, 20, and 3. Non-ice-covered lakes are mostly located in the EP region and YG region. Seasonally ice-covered lakes are concentrated in the other four subregions. The general characteristics of the studied lakes are shown in Table S1.

Figure 1. Spatial distribution of the studied lakes. The lakes are divided into six subregions: Tibetan Plateau (TP), Xinjiang (XJ), Inner-Mongolia Plateau (MP), Northeast Plain and Mountain (NE), Eastern Plain (EP), and Yunnan-Guizhou Plateau (YG). The underlying background displays the elevation from Shuttle Radar Topography Mission Digital Elevation Database v4.1.

Data ARC-Lake Surface Water Temperature

In this study, the ARC-Lake data set (version 3.0) retrieved from the ERS and Envisat satellites during 1996–2011 was used as a reference (MacCallum & Merchant, 2012) to evaluate and validate the FLake model. The LSWTs are estimated only in water pixels with clear-sky conditions. When the lake was close to freezing or frozen, the temperature values were replaced by 0°C. The data set contains daytime and nighttime observations on a 0.05° grid and on a lake-by-lake basis. Many studies have shown that ARC-Lake observations compare well with in situ observation records, with a small bias (<1°C) (MacCallum & Merchant, 2012). The ARC-Lake data set has been widely used to validate lake simulations (Thiery et al., 2015; Verseghy & MacKay, 2017). We used the all-day mean LSWT, calculated as the average of daytime and nighttime observations. The lake mean LSWT with bimonthly temporal resolutions was used to optimize the model parameters in our study. Daily reconstructed ARC-Lake observations were used to demonstrate the simulation performance of the ice-free season.

Geomorphometric Characteristics of Lakes

HydroLAKES data set version 1.0 estimates the volumes of all global lakes with surface areas greater than 10 ha using a geostatistical model based on surrounding topography information (Messager et al., 2016). The data set also provides other additional lake characteristics required by FLake, including lake boundary, lake mean depth, lake surface area and latitude/longitude of lake center coordinates. A basic assumption of the geostatistical model is that the topography around the lake is a reliable predictor variable of its bathymetry. The model has been trained and validated on 12,150 existing records of lake depth obtained from 14 national and international datasets. Overall, the lake depths from HydroLAKES match well with reference depths (Messager et al., 2016).

Moderate Resolution Imaging Spectroradiometer (MODIS) Albedo Products

We used Moderate Resolution Imaging Spectroradiometer (MODIS) BRDF/albedo products (MCD43B3, version 5) during the period of 2001–2005 as criteria to correct the albedo parameterization scheme in FLake. The products combine the albedo retrievals from the Aqua and Terra satellites and provide albedo values for different bands, which agree with the in situ observations (He et al., 2014; Schaaf et al., 2002). In this study, we used the white-sky and black-sky albedos for shortwave broadband with a spatial resolution of 1 km and an interval of 8 days. We carefully selected the pixels with quality flags between 0 and 2 to ensure relatively high data quality. The diffuse skylight ratio was calculated using the monthly diffuse and direct downward shortwave radiation products from the National Center for Environmental Forecasting and was introduced to calculate blue-sky albedo. The detailed method is presented in Li et al. (2016).

MODIS Water Clarity Product

Wang et al. (2020) developed a robust water clarity algorithm based on two water color parameters retrieved from the MODIS surface reflectance product MOD09. The light extinction coefficient (K_d) is closely correlated with water clarity. We obtained the K_d from the water clarity product published by Wang et al. (2020). Maps of Secchi desk depth (Z_SD) for 153 large lakes across China during the period of 2000–2017 were generated for the first time. The Z_SD product showed good consistency with in situ measurements, demonstrating reliability for use. The multiyear average Z_SDs of 84 lakes are available for our study. Following the empirical equations recommended by Layden et al. (2016), Z_SD was converted to K_d by using K_d = (0.757/Z_SD) + 0.07 m⁻¹ when Z_SD ranged from 2 to 10 m. K_d = 1.7/Z_SD was applied for lakes outside the former range. For the 10 lakes without water clarity values, K_d was set to a constant value of 3 m⁻¹.

China Meteorological Forcing Data Set (CMFD)

The FLake forcing variables include near-surface air temperature, specific humidity, wind speed, surface pressure, downward shortwave radiation and downward longwave radiation. We obtained 3-h forcing variables during the period of 1996–2011 from the China Meteorological Forcing Data set (CMFD). The CMFD is a high spatial resolution (0.1°) forcing data set developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (He et al., 2020; Yang et al., 2010) and is the amalgamation of a variety of data sources, including satellite products, reanalysis data sets and meteorological station measurements. This data set has good accuracy for surface meteorological variables and has been widely used in land surface process and hydrological modeling in China (Huang et al., 2017; Lazhu et al., 2016). The forcing data were extracted by the boundaries of lakes and then calculated through spatial averaging. However, considering that this data set are produced by merging various records of land stations, wind speed data are not suitable for representing the actual wind speed over lakes. Scaling of wind speed equation (Equation 1) recommended by Hsu (1988) was adopted to correct the bias. [Image Omitted. See PDF]Here, U_lake is the wind speed (m s⁻¹) over the lake, and U_land (m s⁻¹) is the wind speed over land.

Model Description

FLake is a one-dimensional freshwater lake model consisting of two vertical water layers (Mironov, 2008; Mironov et al., 2010). FLake divides the water column into an upper mixed layer and an underlying thermocline layer. The upper layer is well mixed and vertically uniform. The vertical temperature profile of the thermocline layer, extending down to the bottom, is parameterized based on the self-similarity theory. Moreover, FLake also describes the thermal structure of one ice layer, one snow layer on lake ice, and two thermally active layers for lake sediments. These layers are also parameterized by the same self-similarity theory. A prognostic entrainment equation for convective conditions and a relaxation-type equation for stable conditions are used to calculate the mixed-layer depth (Mironov, 2008). The temperature profile in the thermocline layer is computed through the self-similarity theory, that is, [Image Omitted. See PDF]Here, t is time (s). z is depth (m). ${\theta }_s\left(t\right)$ is the water temperature in mixed layer (K). h(t) is the depth of mixed layer (m). Δθ(t) = θ_s(t) − θ_b(t) is the temperature difference (K) across the thermocline with the depth of Δh(t), and θ_b(t) is the water temperature (K) at the bottom of the thermocline layer. ${\mathrm{\Phi }}_{\theta }$ is a dimensionless “universal” shape function of the relative (dimensionless) depth $\mathit{\zeta }=[z-h\left(t\right)]/\mathrm{\Delta }h\left(t\right)$, and ${\mathrm{\Phi }}_{\theta }$ satisfies two boundary conditions: ${\mathrm{\Phi }}_{\theta }\left(0\right)=0$ and ${\mathrm{\Phi }}_{\theta }\left(1\right)=1$.

The dimensionless shape function ${\mathbf{\Phi }}_{\mathit{\theta }}$ is obtained by [Image Omitted. See PDF]Here, C_θ is the shape factor and is assumed to vary with time, that is, [Image Omitted. See PDF]In Equation 4, t_rc is the relaxation time (s), which is parameterized as the time of the thermocline shape curve varying from one to the other. $C_{\theta }^{\max }$($C_{\theta }^{\min }$) is an empirical constant and represents the maximum (minimum) value of shape factor C_θ.

The water albedo is set to a constant value (0.07) in FLake, while ice albedo is assumed to be an empirical formulation of lake ice surface temperature: [Image Omitted. See PDF]In Equation 5, α_max is the albedo of white ice, set to 0.6. α_min is the albedo of blue ice, set to 0.1. T_f is the freezing temperature (K), and T_s is the ice surface temperature (K). The snow albedo is set equal to the ice albedo. The Beer-Lambert decay law is used to describe the vertical heat transfers in water and ice/snow. The fluxes of momentum and sensible and latent heat fluxes at the lake surface are estimated by a parameterization scheme, which comprehensively considers the relationship between roughness length, potential temperature, specific humidity and wind velocity.

Owing to the simple but realistic representation of major physics processes in lakes, FLake has a high computational efficiency compared with other lake models and has been extensively shown to have good accuracy in reproducing lake surface temperature changes in lakes of different morphometries and climate conditions (Balsamo et al., 2012; Huang et al., 2019; Lazhu et al., 2016; Lv et al., 2019). However, FLake has a limited ability to simulate the thermodynamic processes of deep lakes with depths greater than 60 m due to the absence of monimolimnion layer representation (Balsamo et al., 2012). However, this fact should have minor impacts in our study, as only two lakes in our study have depths greater than 60 m.

Lake-Specific Parameters

Because of the simplicity of the parameterization scheme and the few external lake-specific parameters used, it is relatively easy to tune FLake to specific applications. These external lake-specific parameters are albedo, lake mean depth, K_d, fetch, relaxation coefficient (c_relax_C), depth of the thermally active layer of bottom sediments (depth_bs), and temperature at the outer edge of the thermally active layer of bottom sediments (T_bs) (Layden et al., 2016; Mironov, 2008; Mironov et al., 2010). We evaluated and modified the former three parameters by referring to the model-tuning algorithm of Layden et al. (2016) to improve the model's capabilities in LSWT simulations during the period of 1996–2011. Fetch and c_relax_C were set to the recommended values from Layden et al. (2016). The bottom sediment module was switched off in the simulations. Namely, the heat flux from the bottom sediment was set to zero flux, and the values of depth_bs and T_bs had no effect on the simulations. Detailed modifications of the parameters are presented in Sections 2.4.1–2.4.3, and the sensitivity experiments are shown in Table 1. Additionally, to obtain physically reasonable initial conditions, we repeatedly ran the model with a cycled one-year atmosphere forcing of 1996 until an equilibrium state was obtained. Then, the initial conditions for 1997 were determined using the values of the prognostic variables at the end of the final year.

Table 1 Experimental Design of Freshwater Lake Model (FLake)

Experiments	Ice albedo	Water albedo	Depth (m)	Kd (m−1)
CTRL	Default	0.07	10	3
Exp_albedo	Modified	MCD43B3	10	3
Exp_depth	Modified	MCD43B3	HydroLAKES	3
Exp_Kd	Modified	MCD43B3	HydroLAKES	f(water clarity)a
Exp_tuning	Modified	MCD43B3	HydroLAKES*factor1 factor1:0.25–4.0	f(water clarity)a*factor2 factor2:0.25–4.0
Exp_optimization	Modified	MCD43B3	HydroLAKES* effective factor1	f(water clarity)a*effective factor2

Note. K_d, light extinction coefficient.

^aThe empirical formulations recommended by Layden et al. (2016) were used to calculate K_d based on a satellite-derived water clarity product (Wang et al., 2020).

Layden et al. (2016) used satellite-observed LSWT from ARC-Lake as the criterion to tune key lake properties (albedo, lake mean depth, and K_d) to make FLake reproduce simulations as closely as possible to the observed variations in LSWT. This model-tuning algorithm avoided the need for a detailed description of lake characteristic information, which is usually not available in many lakes. However, an increasing amount of lake characteristic information can now be obtained from satellite products, such as albedo, lake depth, and water clarity. We improved the tuning process by acquiring these parameters, which were more efficient and physically reasonable for representation of lake properties. Similarly, some necessary adjustments for lake mean depth and K_d were adopted in our experiments to achieve the best performance in simulating LSWT.

Albedo

As mentioned above, the ice albedo in FLake is calculated based on an empirical function of lake surface ice temperature. The scheme originates from a sea model (Mironov, 2008) and can describe seasonal variations in ice albedo to some extent. In China, most seasonally ice-covered lakes are in the TP region. However, some studies pointed out that FLake seriously overestimated the ice albedo in TP lakes (Lang et al., 2018; Li et al., 2018). For example, Lang et al. (2018) compared the simulated albedo against field measurements in Ngoring Lake in 2017 and found that the maximum albedo bias of FLake reached 0.48. The overestimated ice albedo indicates that more solar radiation is reflected off the lake, which affects the melting of lake ice and results in a delay in ice-off date. In sensitivity experiments, we set a series of combinations of white ice albedo (0.2–0.8) and blue ice albedo (0.1–0.6) to tune these two empirical constants until the modeled ice albedo closely approximated the MODIS albedo values. For water albedo, the default values of water albedos in the model were replaced by the averaged MODIS albedos during the open-water period.

Lake Mean Depth

Lake depth directly represents the thermal capacity of the lake and controls how long it takes to cool the lake in autumn and winter, thus controlling the ice-on date. FLake recommends using mean depth as an input parameter. In field observations, it is difficult to directly obtain the mean depth of lakes, especially for large lakes. We obtained lake mean depths from the HydroLAKES data set (Messager et al., 2016) and used them in sensitivity experiments. We also estimated the effective lake depth (multiplying the mean depth by a scale factor) to closely reproduce the observed LSWT dynamics by minimizing the root mean square error (RMSE) between the simulations and observations. The scale factor ranges from 0.25 to 4.0 in the tuning process. The effective lake depth is often not the same as the measured lake mean depth. Lake depth is not the only factor that affects the heat exchange between the lake surface and atmosphere. Other factors include lake bottom sediment, topography, altitude, and lake area. The effective lake depth will compensate for the influences of these factors on the surface heat exchange. Additionally, FLake does not consider the monimolimnion for deep lakes; thus, an artificial depth (60 m) is adopted for deep lake simulation (Le Moigne et al., 2016; Perroud et al., 2009).

Light Extinction Coefficient

The absorption and attenuation of solar radiation within water bodies are the main processes controlling water surface temperature. K_d was used to represent the amount of attenuation of solar radiation within the water column. Generally, a smaller K_d means that a greater amount of surface heat goes into deep water and results in a lower surface temperature of the lake. Conversely, a higher K_d leads to a higher surface water temperature. Existing field observations of K_d are scarce. Some constant values of K_d (for example, 3 m⁻¹) are usually used to drive the model. In the default configurations, K_d was set to 3 m⁻¹. The water transparency of lakes in China has large spatial heterogeneity, subject to the intensity of human activity and phytoplankton density (Feng et al., 2019). The recently published water clarity product (Wang et al., 2020) was used to estimate the K_ds of lakes in China. Similar to tuning the lake mean depth, K_d is multiplied by a series of scale factors (0.25–4.0) in the tuning process to obtain the effective K_d for each lake. The model was tuned using the optimal combination of depth and K_d factors, and 100 possible combinations for each lake were used in the tuning process.

Metrics Used for Parametrization

A model-tuning algorithm was used to estimate white/blue ice albedo, effective depth and K_d to make the simulations as close to the observations as possible, namely, reaching a minimum RMSE between the modeled and observed albedo or LSWT. The tuning process is widely used in lake model simulations (Layden et al., 2016; Li et al., 2019; Woolway & Merchant, 2019) and can further improve the LSWT simulation.

Three other statistics, including the mean bias error (BIAS), Pearson correlation (R), and Taylor score (TS), were also used to evaluate model skill. The formulas to calculate these statistics are shown below: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where N is the number of discrete points (in time or space). S_i and O_i are the simulation and observation at point i, respectively. $\overline{S}\left(\overline{O}\right)$ is the average simulation (observation) over N sample points. σ is the standard deviation of the simulation divided by that of the observation. R₀ is the achievable maximum correlation coefficient and is set to 1. Both lower RMSE and lower BIAS represent a better performance of the simulation. R shows the temporal/spatial variation similarity between the simulation and observation. TS can simultaneously assess models in amplitude and temporal/spatial patterns (Taylor, 2001). A higher TS implies better model performance to match the observation (Huang et al., 2016).

Considering that ice temperatures were excluded in the ARC-Lake product, model evaluations were performed during the ice-free season. Layden et al. (2015) showed that the first day with lake mean LSWT above/below 1°C was an effective indicator of ice-off/on date for large lakes. The lake ice phenology defined by Layden et al. (2015) has good consistency with in situ observations of ice-on and ice-off days for 21 Eurasian and North American lakes and has been widely adopted in research (Woolway & Merchant, 2019). We followed the method of Layden et al. (2015) and used 1°C as the threshold to define the ice-on/off date and ice-free season.

Results Evaluation of FLake Performance in the Control Experiment

In the control experiment (CTRL), ice albedo was calculated based on the default ice parameterization in FLake. Water albedo was set to a reference value of 0.07. Lake mean depth was set to 10 m. K_d was set to 3 m⁻¹. Since the purpose of the ARC-Lake data set is to monitor the surface water temperature of the lake, excluding the ice temperature, the evaluations were only carried out during the ice-free season.

Table 2 lists four statistical coefficients to demonstrate the average performances of FLake in simulating the LSWT in six subregions using CTRL settings. Overall, FLake roughly reproduced the seasonal cycle of LSWT, with R ranging from 0.89 to 0.98, and generally underestimated LSWT in the majority of lakes (Figure S1). However, the performances of CTRL showed large differences among the six subregions. CTRL clearly underestimated the average monthly LSWT of lakes in the TP region, XJ region, MP region, and NE region, slightly underestimated that of lakes in the YG region and slightly overestimated that of lakes in the EP region. CTRL has the worst representation for lakes in the TP region, but it has reasonable representation for lakes in the EP and YG regions. These results indicate that the LSWT simulation of FLake over seasonally ice-covered lakes is poorer than that over non-ice-covered lakes. The RMSEs of non-ice-covered lakes ranged from 1.44°C to 2.27°C, while the RMSEs of seasonally ice-covered lakes ranged from 1.72°C to 8.13°C (Figure S1). To further clarify the reasons for this difference, we examined the monthly RMSE and BIAS between the CTRL simulations and observations in six subregions to determine the months that were most responsible for the discrepancy (Figure 2). We found that the seasonally ice-covered lakes usually had the largest RMSEs and BIASs in winter and spring, reflecting the biases of the ice-off/on date and duration of ice cover (Figures S2 and S3). The seasonally ice-covered lakes in the TP and XJ regions showed the same bias pattern, with a positive bias in summer and negative biases in other seasons (Figure 2). For these lakes, CTRL underestimated the LSWT by a maximal BIAS in winter, and BIAS gradually decreased in spring and turned to positive values in June. The overestimations could last three to four months and then change to underestimations again in October. For lakes in the MP and NE regions, the overestimations continued even longer until November. The simulation performances of non-ice-covered lakes were quite different from those of seasonally ice-covered lakes. The RMSEs of lakes in the EP and YG regions did not show obvious seasonal variations, with a range of 0.97°C–2.53°C. The variations in BIASs for non-ice-covered lakes were largely random and lake specific, suggesting that parameters other than ice albedo may be more responsible for the model biases in lakes over the EP and YG regions.

Table 2 The Average Pearson Coefficient (R), Root Mean Square Error (RMSE), Mean Bias Error (BIAS), and Taylor Score (TS) of Lake Surface Water Temperature (LSWT, °C) Between Simulations From the Control Experiment (CTRL) and Observations From ARC-Lake of Lakes in Six Subregions During the Open-Water Period From 1996 to 2011

Subregion	R	RMSE	BIAS	TS
TP	0.90 ± 0.07	4.17 ± 1.39	−1.60 ± 1.12	0.74 ± 0.13
XJ	0.89 ± 0.08	4.51 ± 1.42	−1.44 ± 0.75	0.80 ± 0.14
MP	0.94 ± 0.01	4.03 ± 0.25	−1.25 ± 0.39	0.93 ± 0.03
NE	0.94	4.11	−1.16	0.94
EP	0.98 ± 0.00	1.91 ± 0.25	0.20 ± 0.50	0.99 ± 0.00
YG	0.93 ± 0.06	1.74 ± 0.45	−0.27 ± 0.66	0.93 ± 0.07

Note. ± represents 1 standard deviation.

Figure 2. The average monthly root mean square error (RMSE) and mean bias error (BIAS) between the simulated lake surface water temperature in the control experiment (CTRL) and observations from ARC-Lake of lakes in six subregions during the open-water period from 1996 to 2011. The bars represent the RMSE, and the red lines represent the BIAS. The black lines are used to distinguish the positive and negative BIASs. Evaluations over the ice-covered period are omitted.

Trial Lake Simulated With the Improved Lake Scheme Control Experiment for the Trial Lake

Qinghai Lake was used as a trial lake to exhibit the performance of CTRL in simulating LSWT. CTRL captured the seasonal variations in LSWT well (Figure 3). However, CTRL significantly overestimated the ice-covered period and underestimated the open-water period. The modeled ice-off date was ∼28 days later than the observed ice-off date, while the modeled ice-on date was ∼24 days earlier than the observed ice-on date. Additionally, the simulated LSWTs were higher than the observations in summer, while the simulated LSWTs were generally lower than the observations in other seasons (Table 3). This result was consistent with the simulations of Su et al. (2019). The simulated LSWT had a positive BIAS of 1.89°C in summer, a negative BIAS of −2.70°C in spring, and a negative BIAS of −4.01°C in autumn and early winter (Table S2). CTRL underestimated the annual average LSWT during the open-water period with a BIAS of −1.90°C and an RMSE of 1.97°C. Moreover, CTRL showed similar performances in simulating daytime and nighttime LSWT (Figures S4 and S5; Table S3).

Figure 3. The multiyear mean time series of lake surface water temperature (LSWT) from the ARC-Lake observations and the simulations in the control experiment (CTRL) for Qinghai Lake during the period of 1996–2011. The shaded area and error bar represent one standard deviation. ARC-Lake observations below 1°C represent the ice-covered period of Qinghai Lake. Inserted plots are the comparisons of ice-off (ice-on) dates between CTRL simulations and observations. The first day of LSWT > 1°C in spring/summer indicates the ice-off date. The first day of LSWT [less than] 1°C in winter indicates the ice-on date. Error bars represent one standard deviation of ice-off/on dates.

Modification of Albedo

The large negative biases between the modeled and observed LSWTs in spring and winter were caused by the biases in simulating ice phenology. Ice albedo is thought to be the main lake parameter that affects the ice-off date. During the period 2001–2005, the average ice albedo was 0.51 in FLake, much higher than the 0.21 in the MODIS albedo products (Figure S6). The overestimation of ice albedo can attenuate absorbed solar radiation and eventually result in postponement of the ice-off date. Hence, the ice albedo parameterization scheme was first modified in our study.

We used the MODIS albedo products to tune the empirical constants within the albedo parameterization scheme. Smaller albedo constants were determined, 0.3 for white ice albedo and 0.1 for blue ice albedo. The average water albedo observed (0.017) during the open-water period replaced the default water albedo (0.07) in the albedo experiment (Exp_albedo). Although large deviations in albedo still existed in some months, the modified albedo was generally close to the observed ice albedo (Figure S6). Exp_albedo clearly displayed improved skill in simulating LSWT in spring (Figure S7). The underestimation of spring LSWT and delay in ice-off date were both significantly reduced. The BIAS and RMSE between the simulated and observed LSWT in spring decreased from −2.70°C and 2.84°C to 0.25°C and 1.33°C, respectively (Table 3). Moreover, the overestimation in summer LSWT increased in Exp_albedo, with the BIAS increasing from 1.89°C to 2.50°C. Since the observed water albedo was lower than the default value, more solar radiation heated the surface water. Similarly, due to the decreasing water albedo, the negative BIAS in autumn and early winter was reduced (Table S4).

Modification of Lake Mean Depth

Reproducing the ice-on dates as realistically as possible is another difficulty in simulating lake thermodynamics. Lake depth plays a determinant role in controlling the ice-on dates. Lakes usually begin to freeze in shallows along the shore and then expand toward the center. Shallow waters have lower thermal inertia, requiring less time to lose heat to warm air and reach thermal equilibrium with the air, thus forming lake ice earlier. In the depth experiment (Exp_depth), the mean depth of Qinghai Lake was set to a HydroLAKES depth of 16.8 m. Figure S8 shows that the simulated dates of ice-on occurred later and were closer to the dates of ARC-Lake. The BIAS in autumn and early winter significantly decreased to −1.89°C. The differences between the modeled and observed LSWTs were also reduced in other seasons (Table S5). The RMSE (BIAS) of spring LSWT was 1.10°C (−0.34°C). The overestimation of summer LSWT was slightly mitigated (Table 3). The deeper the lake is, the more solar radiation is needed for the water column to fully mix in spring. Under the same radiation conditions, the surface water temperature of deep lakes will be lower than that of shallow lakes, resulting in a reduced spring LSWT in Qinghai Lake. The decrease in summer LSWT may be caused by deepening of the mixed layer in Exp_depth. The RMSE (BIAS) of the annual average LSWT covering the open-water period is 0.56°C (−0.32°C). Overall, the simulated LSWT was closer to the observed LSWT when the lake mean depth was adjusted from 10 to 16.8 m.

Table 3 The Monthly Root Mean Square Error (RMSE) and Mean Bias Error (BIAS) of Lake Surface Water Temperature (LSWT, °C) Between the Simulations From Five Experiments and Observations From ARC-Lake for Qinghai Lake During the Open-Water Period From 1996 to 2011

Month	RMSE					BIAS
	CTRL	Exp_albedo	Exp_depth	Exp_Kd	Exp_optimization	CTRL	Exp_albedo	Exp_depth	Exp_Kd	Exp_optimization
Apr	3.93	2.31	2.07	1.98	1.77	−3.87	−2.11	−1.93	−1.81	−1.55
May	2.09	2.56	1.31	0.81	1.22	−1.75	2.37	0.98	0.11	−0.79
Jun	2.70	3.63	2.91	1.34	1.01	2.54	3.55	2.78	1.17	0.18
Jul	2.12	2.57	2.20	1.19	0.78	2.05	2.49	2.09	1.11	0.38
Aug	1.25	1.52	1.73	1.22	0.88	0.96	1.36	1.57	1.08	0.68
Sep	1.20	0.85	0.73	0.79	1.00	−1.00	−0.6	0.10	0.42	0.67
Oct	2.59	2.08	1.58	0.77	0.52	−2.55	−2.01	−1.49	−0.40	0.21
Nov	3.08	2.49	1.57	0.59	0.78	−3.05	−2.43	−1.43	−0.17	0.60
Dec	9.95	7.56	5.12	4.12	0.77	−9.89	−7.51	−4.99	−3.82	−0.14

Note. Evaluations over the ice-covered period (January-March) are omitted.

Modification of the Light Extinction Coefficient

Solar radiation is the fundamental energy source of the whole lake system. K_d is a key parameter that determines the ability of solar radiation to penetrate water bodies, thus affecting the vertical thermal structure of lakes. In general, the more transparent the lake is, the easier it is for solar radiation to penetrate deeper into the water, resulting in a lower LSWT. Qinghai Lake is a typical dimictic lake with two turnovers in spring and autumn and one stratification in summer. The lake mixes completely from the surface to the bottom during the turnover period. Thus, the effect of K_d on LSWT is more notable in the stratification period. In CTRL, K_d was set to a constant of 3 m⁻¹. Qinghai Lake is in the TP region, with fewer aquatic plants and pollution discharges. The default value of K_d was too high for Qinghai Lake. K_d in the Exp_Kd experiment was tuned to a lower value of 0.31 m⁻¹, which was derived from a water clarity product (Wang et al., 2020). Exp_Kd produced a significantly decreased summer LSWT (Figure S9). BIAS (RMSE) in summer was decreased to 1.13°C (1.21°C). Moreover, by reducing K_d, RMSE decreased in both spring and autumn (Table S6). The higher lake transparency led to more rapid warming in spring, a deeper mixed layer and a warmer thermocline in summer, which made the stratification more stable and ultimately led to slower fall cooling. The RMSE and BIAS for annual LSWT during the open-water period were 0.60°C and −0.27°C, respectively.

Joint Optimization of the Parameters

After gradually adopting satellite-derived albedo, lake mean depth and K_d, the simulations became increasingly closer to the observations (Figure 4). Although Exp_Kd greatly improved the simulation performance, it still overestimated summer LSWT and estimated ice-on dates too early. Tuning lake mean depth and K_d can theoretically optimize the simulations. The underlying mechanism of depth optimization is to compensate for the influences of other factors on heat/water exchange at the lake-air interface. K_d was estimated from a Z_SD product based on two empirical formulas, with large uncertainties. Tuning K_d aims to reduce the errors and uncertainties in the empirical formulas. We conducted a joint tuning process to adjust the lake mean depth and K_d in the tuning experiment (Exp_tuning). The depth factor and K_d factor both ranged from 0.25 to 4.0. We designed 100 combinations of depth factors and K_d factors. The factors were determined based on the minimum RMSE between the simulated and observed LSWT during 1996–2011.

Figure 4. The multiyear mean time series of lake surface water temperature from the ARC-Lake observations and simulations in five experiments for Qinghai Lake during the period of 1996–2011. The shaded area and error bar represent one standard deviation. ARC-Lake observations below 1°C represent the ice-covered period of Qinghai Lake.

For Qinghai Lake, the effective depth factor is 1.5, and the K_d factor is 0.75. The experiment using optimal effective factors is called the optimization experiment (Exp_optimization, Figure S10). RMSE decreased by more than 57% compared with CTRL (Figure 5; Table S7). The declines in RMSE were most obvious from October to December, reaching 82%, while the declines were less obvious in August and September, reaching 23%. The most pronounced LSWT improvement in spring was in Exp_albedo and that in summer was in Exp_Kd. Exp_albedo, Exp_depth and Exp_Kd showed similar improvements in the LSWT simulations in autumn and early winter. Additionally, as the time series of 16 years was relatively short, R and TS in tuning experiments were prone to be affected by extreme values or null values, which produced large variability in monthly R and TS (Figure S11). The BIAS of ice phenology was reduced to 9.6 days for ice-off dates and 0.2 days for ice-on dates (Figure S12). Figure 6 shows the overall performances of standard deviation, R and RMSE compared with observations at annual scales using the Taylor diagram. The results showed that these experiments had obvious differences in simulating LSWT. Exp_optimization yielded the best simulation ability and outperformed other experiments for the whole open-water period. This result implied that the optimized FLake can realistically reproduce the LSWT in Qinghai Lake both seasonally and interannually.

Figure 5. Histograms of root mean square error (RMSE) and mean bias error (BIAS) of lake surface water temperature between the simulations from five experiments and the observations from ARC-Lake for Qinghai Lake during 1996–2011.

Figure 6. Taylor diagram of monthly lake surface water temperature of Qinghai Lake in the open-water period for five experiments.

Application of the Improved Lake Scheme in All Lakes

After adjusting these three parameters, FLake achieved good simulations in Qinghai Lake. Following the same strategy, we optimized FLake performances over all 94 large lakes across China. The ice albedo, water albedo, effective lake depth factor and K_d factor, which were finally applied in Exp_optimization for all lakes, are listed in Table S1. Table 4 shows the average RMSEs of subregions under five experiments.

Table 4 The Average Root Mean Square Error (RMSE) Between the Monthly Simulated Lake Surface Water Temperature (LSWT, °C) From Five Experiments and Monthly Observed LSWT From ARC-Lake of Lakes in Six Subregions During the Open-Water Period From 1996 to 2011

Subregion	RMSE
	CTRL	Exp_albedo	Exp_depth	Exp_Kd	Exp_optimization
TP	4.17 ± 1.39	3.18 ± 1.20	2.57 ± 0.85	2.46 ± 0.80	1.91 ± 0.62
XJ	4.51 ± 1.42	4.01 ± 1.66	3.71 ± 1.13	3.54 ± 0.81	3.02 ± 0.44
MP	4.03 ± 0.25	3.71 ± 0.51	3.33 ± 1.12	3.33 ± 1.12	2.04 ± 0.38
NE	4.11	3.62	3.11	3.11	3.07
EP	1.91 ± 0.25	1.91 ± 0.25	1.77 ± 0.25	1.78 ± 0.26	1.59 ± 0.23
YG	1.74 ± 0.45	1.73 ± 0.42	1.29 ± 0.12	1.29 ± 0.11	1.05 ± 0.17

Note. ± represents 1 standard deviation.

Figure 7 exhibits the spatial patterns of RMSE, BIAS, R and TS in Exp_optimization. The results showed that the RMSEs of 70 lakes were less than 2°C, accounting for 74% of the studied lakes. Ninety lakes with BIAS values within ±1°C accounted for 96% of the studied lakes. Moreover, almost all lakes had R and TS values above 0.85. Compared with the CTRL, the improvements were more significant in seasonally ice-covered lakes. For seasonally ice-covered lakes, the BIAS decreased from −1.56 ± 1.05°C to 0.03 ± 0.62°C, and the RMSE decreased from 4.21 ± 1.35°C to 2.07 ± 0.70°C. For non-ice-covered lakes, BIAS and RMSE decreased slightly from 0.14 ± 0.54°C and 1.89 ± 0.27°C to −0.08 ± 0.34°C and 1.52 ± 0.29°C, respectively. Overall, the optimization increased the R and TS values nationally, particularly in the TP region. Exp_optimization also greatly improved the simulation performance of ice phenology. The R-square of the simulated and observed ice-on dates reached 0.74, and the R-square of the simulated and observed ice-off dates was 0.61 (Figure S13). Figure 8 shows the interannual variation in the regional average LSWT during the open-water period from 1996 to 2011. The results showed that the averaged LSWTs of the subregions between Exp_optimization simulations and ARC-Lake observations were highly consistent. In summary, Exp_optimization can reconstruct the dynamic changes in surface water temperature well in large lakes across China.

Figure 7. Spatial patterns of root mean square error (RMSE), mean bias error (BIAS), Pearson coefficient (R) and Taylor score (TS) of lake surface water temperature between the simulations from the optimization experiment (Exp_optimization) and the observations from ARC-Lake during the open-water period from 1996 to 2011. The inset in map shows the frequency distribution of the statistic.

Figure 8. Time series of the yearly average lake surface water temperature during the open-water period from the optimization experiment (Exp_optimization) and from ARC-Lake in subregions from 1996 to 2011. Red and black lines represent simulations and observations, respectively. The red shade and gray shade represent one standard deviation of simulations and observations, respectively.

The largest improvements in Exp_optimization were in spring and winter for seasonally ice-covered lakes (Figure 9). The dominant parameters of LSWT vary with climate, season and lake type (Figure S14). Albedo dominated the LSWT simulation for seasonally ice-covered lakes in spring and winter. Non-ice-covered lakes were more sensitive to lake depth and K_d. In the albedo experiment, the RMSE in the TP region decreased by 76% in January, and the RMSE in the XJ region decreased by 52% in December. The LSWT in the MP region was very sensitive to lake depth in spring. In the depth experiment, the RMSE of LSWT in the MP region decreased by 43%–57% in spring. In Exp_depth and Exp_Kd, the RMSE in the YG region decreased by 18% and 16%, respectively. Notably, the RMSE of the NE region slightly increased in November after tuning the lake depth. This result occurred in Chagan Lake, the only large lake in ARC-Lake for northeastern China. The few effective observations of Chagan Lake in winter made the RMSE assessment less reliable. Taylor diagrams comprehensively show the simulation performances of five experiments in six subregions (Figure 10). Overall, Exp_optimization performed best and can realistically represent the LSWT temporal variations in all subregions. In addition, the improvements in FLake were significant in the TP region, XJ region and NE region.

Figure 9. Comparisons of seasonal lake surface water temperature between observations from ARC-Lake and simulations from the optimization experiment (Exp_optimization) for six subregions during 1996–2011. The black line is a linear fit line. The green line is the 1:1 line.

Figure 10. Taylor diagram of monthly lake surface water temperature of subregions in the open-water period for five experiments.

Discussion

Using the ARC-Lake satellite observations of lakes with surface areas larger than 100 km² from 1996 to 2011, the one-dimensional freshwater lake model FLake was evaluated and calibrated to simulate the LSWTs in 94 large lakes across China. Based on satellite observations, we tuned three lake-specific parameters (albedo, lake mean depth and K_d) in FLake to improve the ability to represent seasonal and interannual variations in LSWT. Compared with a previous study (Layden et al., 2016), our observation-driven tuning strategy improved the efficiency of the tuning process and reduced the risk of equifinality. The tuning process demonstrated the potential of FLake to reproduce realistic LSWT simulations across China, and FLake can be considered an appropriate tool for investigating and predicting LSWT changes in large lakes across China.

The correct lake-specific parameters greatly improved the simulation performance of FLake. We found that these parameters were consistent with observations or previous calibration studies. For example, the effective lake depths of Qinghai Lake and Nam Co Lake are 25 and 44 m in the study, respectively, which were very close to the previous results of 21 (Su et al., 2019) and 40 m (Lazhu et al., 2016). The K_d values of Poyang Lake and Liangzi Lake were 0.91 and 0.75 m⁻¹, respectively, which were similar to the results of 0.8 and 1.1 m⁻¹ by Li et al. (2019).

Regional differences were apparent in the ice albedo, lake mean depth and K_d of large lakes in China (Table S1). For seasonally ice-covered lakes, ice albedos were lower in the TP region and higher in other regions. In addition, lakes in eastern China were shallow and low transparent, while lakes in western China were deep and high transparent. White ice albedo represents the largest ice albedo in seasonally ice-covered lakes. In tuning the white/blue ice albedo, the white ice albedos were generally not higher than 0.3 in the TP region, while the white ice albedos in other regions were generally between 0.4 and 0.6. This finding implied that the default ice albedo in FLake seemed too high for the majority of large lakes in China. In the tuning process of effective lake depth, we found that lakes with tuning to a shallower effective depth were dominant, accounting for ∼62% of the total lakes. Additionally, 28% of lakes tended toward a deeper effective depth. Seventy-seven percent of the lakes with deeper tuning were less than 20 m deep (Figure 11). The bottom sediment module was switched off in our study, and the deeper effective depth can compensate for the influences of bottom sediment on surface temperature. A deeper effective depth has the same effect on the heat storage capacity of lakes as the heat retention in the sediment (Layden et al., 2016). In western China, the K_d values were mostly between 1.5 and 4.0 m⁻¹. Most K_ds in eastern China were lower than 1.5 m⁻¹. In general, the effective K_d factor decreased with increasing K_d in both western lakes and eastern lakes (Figure 11). Across all studied lakes, 72% were tuned to be more transparent. Twenty-five percent of lakes were tuned to be more turbid. Lakes with tuning that was more transparent were widely distributed in China. Lakes with more turbid tuning were distributed in some lakes with high transparency values (K_d < 0.25 m⁻¹) in the southern TP region, MP region and YG region (Table S1). The optimal lake-specific characteristics provide guidance for the improvement of LSWT simulations in China in the future.

Figure 11. (a) Comparison of lake mean depth and modeled effective lake mean depth. (b) Comparison of the light extinction coefficient (Kd) and modeled effective light extinction coefficient. Color circles represent the effective depth (Kd) factors. The insert is the box plot of modeled effective lake mean depth (Kd) factors at different lake mean depths (Kds).

Although the FLake model has been shown to be an effective model for simulating lake surface temperature, we recognize that some uncertainties and limitations still exist in our research. First, accurate simulation representation strongly depends on the reliability of the meteorological forcing data. The forcing data used in the study are an assimilation product from land meteorological observations. There are some biases between meteorological variables near lakes and those near land. For example, due to the small surface roughness of the water surface, the wind speeds over lakes are significantly higher than the wind speeds over surrounding land (Lazhu et al., 2016; Su et al., 2019). The underestimation of the 10 m wind speed can lead to a decreasing upward heat loss and therefore an overestimation of LSWT. The relationship between lake wind and land wind varies with lake and season (Phillips & Irbe, 1978; Quinn, 1979). No one wind speed adjustment is most suitable for all lakes. Layden et al. (2016) compared three kinds of wind speed adjustments in 35 globally distributed large lakes and found that Hsu's (1988) equation was most applicable to seasonally ice-covered lakes. Since most large lakes in China are seasonally ice-covered, we applied Hsu's (1988) equation to correct these wind speeds. The equation does introduce some biases and uncertainties in our study. More lake wind observations are needed in the future.

In addition, most lakes in the TP region are saline lakes (Song et al., 2013). The salinity of lake water strongly affects the maximum water density. Huang et al. (2019) found that adjusting the temperature of maximum density in Nam Co Lake can significantly enhance the convective mixing induced by the vertical density gradient and improve the simulation of the vertical thermal structure during cold seasons. Su et al. (2019) suggested that salinity affected the timing of overturn and duration of lake stratification in Qinghai Lake. These results implied that salinity indirectly affected the variation in surface temperature by affecting the vertical temperature structure of the lake. Moreover, other physical properties of lakes, such as freezing temperature, thermal conductivity, and saturated vapor pressure at the lake-air interface, are also related to the salinity of the lake (Jackett et al., 2006; Sun et al., 2008; Wen et al., 2015). These processes are not described in FLake and will introduce some errors to the simulations of saline lakes in the TP region.

Additionally, we used the ARC-Lake satellite observations to evaluate and calibrate FLake. Generally, satellite-derived LSWT is slightly lower than buoy observations. The satellite product retrieves the surface skin temperature, while the buoy is placed in the bulk layer deeper in the water, which represents the water temperature of this layer. The discrepancy between skin and bulk temperature is known as the cool skin effect (Donlon et al., 2002; Minnett, 2003; Robinson et al., 1984). The widely positive biases between simulated LSWT and ARC-Lake observations across China during the open-water period may be partly attributed to this phenomenon. Future studies should utilize more in situ observations to further confirm the robustness of lake temperature simulations.

Acknowledgments

This study is supported by the Second Tibetan Plateau Scientific Expedition and Research Program (Grant No. 2019QZKK0405) and the NSFC project Basic Science Center for the Tibetan Plateau Earth System (BCTES, Grant No. 41988101). The authors thank the two anonymous reviewers and the Editor C. Gentemann for their insightful comments.

Data Availability Statement

Satellite-derived LSWTs used in the study are available to download from http://www.laketemp.net/home_ARCLake/data_access.php (ARC-Lake, version 3.0). The geomorphological characteristics of the lakes, including latitude, longitude, surface area, and mean depth, used in this study were obtained from the HydroLAKES database (version 1.0, https://www.hydrosheds.org/pages/hydrolakes). The MODIS albedo product MCD34B3 (version 5) is available from the Land Processes Distributed Active Archive Center (LP DAAC; https://lpdaac.usgs.gov/data_access). Water clarity data in China are provided by Wang et al. (2020): Changes of water clarity in large lakes and reservoirs across China observed from long-term MODIS (https://doi.org/10.1016/j.rse.2020.111949). The China Meteorological Forcing Data set (CMFD) can be obtained from http://data.tpdc.ac.cn/en/data/8028b944-daaa-4511-8769-965612652c49/. The source code for FLake model is available at http://www.flake.igb-berlin.de/site/download. All the links were last accessed on June 27, 2021.