Evaluating and Improving the Performance of Three 1-D Lake Models in a Large Deep Lake of the Central Tibetan Plateau.

The ability of FLake, WRF-Lake, and CoLM-Lake models in simulating the thermal features of Lake Nam Co in Central Tibetan Plateau has been evaluated in this study. All the three models with default settings exhibited distinct errors in the simulated vertical temperature profile. Then model calibration was conducted by adjusting three (four) key parameters within FLake and CoLM-Lake (WRF-Lake) in a series of sensitive experiments. Results showed that each model's performance is sensitive to the key parameters and becomes much better when adjusting all the key parameters relative to tuning single parameter. Overall, setting the temperature of maximum water density to 1.1 °C instead of 4 °C in the three models consistently leads to improved vertical thermal structure simulation during cold seasons; reducing the light extinction coefficient in FLake results in much deeper mixed layer and warmer thermocline during warm seasons in better agreement with the observation. The vertical thermal structure can be clearly improved by decreasing the light extinction coefficient and increasing the turbulent mixing in WRF-Lake and CoLM-Lake during warm seasons. Meanwhile, the modeled water temperature profile in warm seasons can be significantly improved by further replacing the constant surface roughness lengths by a parameterized scheme in WRF-Lake. Further intercomparison indicates that among the three calibrated models, FLake (WRF-Lake) performs the best to simulate the temporal evolution and intensity of temperature in the layers shallower (deeper) than 10 m, while WRF-Lake is the best at simulating the amplitude and pattern of the temperature variability at all depths.

Introduction

As one of the typical land surface types, there are more than 1,200 inland lakes with an individual surface area more than 1 km2 in the Tibetan Plateau, which is the region with the highest lake density in China (Ma et al., 2011; Song et al., 2013; Zhang et al., 2014). The total lake surface area on the Tibetan Plateau is over 47,000 km2 and contributes more than 52% of the total lake coverage in China (Zhang et al., 2014). As situated in the alpine and semiarid to arid climate zone, the water levels, surface areas, and thermal conditions of these lakes are very sensitive to climate change (Lei et al., 2014; Liao et al., 2013; Lu et al., 2017; Phan et al., 2012; Song et al., 2014; Wu & Zhu, 2008; Zhu et al., 2010). The lakes in the Tibetan Plateau have become sensitive indicators of climate and environment change (Liu et al., 2010; Neckel et al., 2014; Zhang et al., 2013).

Compared to the other land surface types, those lakes have much lower albedo and roughness, larger heat capacity, and thermal conductance and are larger sources of moisture for the lower atmosphere (Biermann et al., 2014; Blanken et al., 2011; Huang et al., 2010; Nordbo et al., 2011; Xiao et al., 2016; Xu et al., 2009, 2011). Lakes can exert important impacts on the surface energy exchange (Gu et al., 2016; Li, 2015; Rouse et al., 2005; Wang et al., 2015) and therefore the atmospheric circulations (Eerola et al., 2014; Gerken et al., 2013, 2014), and further affect the weather and climate at local to regional scales (Gu et al., 2016; Leah & Steenburgh, 2017; Martynov et al., 2012; Steenburgh & Campbell, 2017; Veals & Steenburgh, 2015; Welsh et al., 2016; Wen, Lv, et al., 2015; Xue et al., 2017; Yeager et al., 2013).

Climate models are the primary and most efficient tools for the study of the lake‐air interactions and lake climatic effect (Mackay et al., 2009; Notaro, Holman, et al., 2013; Notaro, Zarrin, et al., 2013). Many earlier studies have shown that inclusion of lake effect leads to improved performance of climate models (Mallard et al., 2014; Martynov et al., 2013; Xue et al., 2017; Zhao et al., 2012). With the consideration of the air‐lake interaction in climate models, such as idealized lakes within global climate models (Balsamo et al., 2012; Bonan, 1995; Le Moigne et al., 2016; Lofgren, 1997; Subin et al., 2012) and one‐dimensional (1‐D) or three‐dimensional (3‐D) lake models in regional climate models (Bennington et al., 2014; Gula & Peltier, 2012; Long et al., 2007; Martynov et al., 2010; Mironov et al., 2010; Samuelsson et al., 2010; Xue et al., 2017), the errors in the simulations of the surface air temperature, evaporation, thermal conditions, convection, and precipitation over the lake‐rich regions can be remarkably reduced (Balsamo et al., 2012; Le Moigne et al., 2016; Mallard et al., 2014; Martynov et al., 2013; Salgado & Le Moigne, 2010; Xiao et al., 2016).

Despite the fact that a few attempts have recently been made to fully couple 3‐D lake models with regional climate models (Long et al., 2007; Song et al., 2004; Xue et al., 2017), the 3‐D air‐lake coupling would add much more complexity to the climate model and require considerable computational power (Martynov et al., 2010). In addition, the 3‐D lake models generally have much finer horizontal resolution (around 2–2.5 km) than most current regional climate models and 1‐D lake models are commonly used to be coupled with regional climate models (Bennington et al., 2014; Gu et al., 2015; Le Moigne et al., 2016; Martynov et al., 2012; Notaro, Holman, et al., 2013; Xiao et al., 2016; Xu et al., 2017).

During the recent several decades, 1‐D lake models with different degrees of sophistication in physical processes have been widely developed to simulate the lake thermodynamics (Mironov et al., 2010; Stepanenko et al., 2013). These 1‐D lake models include FLake model based on similarity theory (Kirillin et al., 2011; Mironov et al., 2010; Stepanenko et al., 2013; Stepanenko et al., 2014), eddy‐diffusive lake models (Hostetler et al., 1993; Oleson et al., 2010; Subin et al., 2012), and k‐ɛ turbulence closure lake models (Jöhnk et al., 2008; Perroud et al., 2009; Stepanenko et al., 2013).

Some of these lake models have been coupled with weather or climate models. For example, the FLake model has been adopted by UK Met Office (Rooney & Jones, 2010), Meteo‐France (Salgado & Le Moigne, 2010), Centre National de Recherches Météorologiques (CNRM; Le Moigne et al., 2016), Swedish Meteorological and Hydrological Institute (Samuelsson et al., 2010), Finnish Meteorological Institute (Eerola et al., 2010), and European Centre for Medium‐Range Weather Forecasts (ECMWF; Dutra et al., 2010). The FLake model has also been coupled with the Weather Research and Forecasting (WRF) model (Mallard et al., 2014) and the fifth generation Canadian Regional Climate Model (CRCM5) model (Martynov et al., 2012). The eddy‐diffusive lake models have been recently coupled with the WRF model (Gu et al., 2015, 2016; Xiao et al., 2016; Xu et al., 2017), the fourth version of the International Centre for Theoretical Physics (ICTP) Regional Climate Model (RegCM4; Bennington et al., 2014; Notaro, Holman, et al., 2013; Xue et al., 2017), and CRCM5 model (Huziy & Sushama, 2017; Martynov et al., 2012).

The current 1‐D lake models have been developed for given environmental applications with different degrees of simplification, but none covers all processes to describe the lake‐air interactions in reality. This leads to different model results due to the different physical concepts among these lake models (Perroud et al., 2009; Stepanenko et al., 2010). Before coupling a lake model with a weather or climate model, it is necessary to validate and intercompare the performance of available lake models in reproducing adequately the behavior of surface conditions of different lakes within the simulation domain (Stepanenko et al., 2013; Stepanenko et al., 2014). With the special emphasis on the applications of lake models in climate simulation and weather forecast to simulate the lake‐air interaction, the Lake Model Inter‐comparison Project (LakeMIP, http://ww.unige.ch/climate/lakemip; Stepanenko et al., 2010) was launched in late 2008. This project was set to intercompare the performance of different 1‐D lake models and identify the key processes, such as lake physics, chemistry, hydrology, and biology, and further develop and improve their parameterization schemes.

During the past several years, many studies have conducted a series of validation and intercomparison of the available 1‐D lake models' performance in simulating the thermal features of different lakes around the word. These lakes include tropical (Thiery, Martynov, Darchambeau, Descy, Plisnier, Sushama, & van Lipzig, 2014; Thiery, Stepanenko, Fang, Jöhnk, Li, Martynov, Perroud, Subin, Darchambeau, Mironov & van Lipzig, 2014), subtropical (Lazhu et al., 2016), and temperate (Kheyrollah et al., 2012) large deep lakes; subtropical (Deng et al., 2013; Gu et al., 2013; Perroud et al., 2009; Wen et al., 2016) and temperate (Vörös et al., 2010) large shallow lakes; and temperate small shallow lakes (Martynov et al., 2010; Stepanenko et al., 2014; MacKay et al., 2017; Verseghy & MacKay, 2017; Yao et al., 2014). Based on the findings and knowledge obtained from these studies, considerable progresses have been made to improve the 1‐D lake models (Bennington et al., 2014; Gu et al., 2015; Subin et al., 2012; Xiao et al., 2016; Xu et al., 2017).

Compared to the relatively comprehensive studies on the evaluation and development of the 1‐D lake models over lowlands and wet regions, few such studies over the Tibetan Plateau have been conducted due to few systematic direct field observations limited by the remoteness and inaccessibility of the lakes on the Tibetan Plateau (Wang et al., 2015; Wen et al., 2016). Recently, with the gradual development of the observational networks for the lake internal thermal features and the air‐lake heat and water exchanges over the Tibetan Plateau (Biermann et al., 2014; Guo et al., 2016; Li, 2015; Wang et al., 2010, Wang et al., 2015, Wang et al., 2017; Wen, Lv, et al., 2015), the alpine lakes have received increasing attention, including the differences in the eddy co‐variance fluxes between lake surface and nearby grassland surface (Biermann et al., 2014), and the turbulent flux transport of heat and water between alpine lakes and overlying air (Lazhu et al., 2016; Li, 2015; Wang et al., 2015). The processes of air‐lake heat and water transfer in the lake models can strongly affect the simulation of the lake thermal structure (Subin et al., 2012; Xiao et al., 2016).

However, whether the current lake models are suitable for the simulation of the alpine lake thermal structure over the Tibetan Plateau or not is still unclear so far. These issues need to be deeply revealed to deepen our understanding of the key processes related to the simulations of the lake thermal features over the Tibetan Plateau. Based on 2 years' in situ observation in Lake Nam Co on the Central Tibetan Plateau, the performance of three 1‐D lake models (FLake, WRF‐Lake, and CoLM‐Lake) with default settings of some key model parameters in simulating the alpine lake thermal regime is systematically evaluated first. Then the ability of each lake model to simulate the thermal features of alpine lakes on the Tibetan Plateau is further calibrated and improved by conducting a series of sensitivity experiments. Findings of this study may provide valuable information for further improvements and applications of the current 1‐D lake models over the alpine and semiarid to arid areas.

Study Area

Lake Nam Co, which is located in the Central Tibetan Plateau (N30°30´–30°55´, E90°16´–91°03´, 4,718 m above sea level; Figure 1a), is the second largest lake in the Tibetan Autonomous Region (Lu et al., 2005; Lazhu et al., 2016). Its surface area is around 2,000 km2 (Huang et al., 2017; Zhang et al., 2014), and the maximum depth is over 100 m (Figure 1b). The mean depth of Lake Nam Co is about 40 m (Lazhu et al., 2016). Precipitation and glacier melt water are the main water supply for Lake Nam Co, in which the water loss is entirely contributed by the intensive evaporation and potential groundwater discharge due to the absent outflow (Liu et al., 2010). Precipitation can be stored in the snowpack in winter, thus delaying its arrival at the lake. In addition, the Lake Nam Co is a saline lake with the water salinity of 1.7 g/L (Wang et al., 2009), which can slightly affect the lake water freezing point (Wang et al., 2019).

Figure 1

(a) Distribution of the lakes (blue) with the area over 1 km2 in Tibetan Plateau and the location of Lake Nam Co; the solid red line indicates the terrain height of 2,500 m. (b) Bathymetry of Lake Nam Co and the locations of the water temperature profile observation site and weather station.

Observations suggest that Lake Nam Co is a dimictic lake characterized by thermal stratification during late June to early November and thorough mixing with a uniform temperature throughout the whole water body during winter and spring, and it is generally covered with ice from January to subsequent May (Huang et al., 2015; Lazhu et al., 2016).

Models, Data, Experiments, and Methods

Lake Models

Three 1‐D lake models including the FLake (Mironov, 2008; Mironov et al., 2010), WRF‐Lake (Gu et al., 2015, 2016; Subin et al., 2012; Xu et al., 2017), and CoLM‐Lake (Dai et al., 2003; Dai et al., 2018) are evaluated in current study. The detailed descriptions of each lake model are shown as follows:

FLake Model

The FLake model is developed based on the concept of self‐similarity (assumed shape) of the vertical temperature profile in the thermocline (Kirillin et al., 2011; Mironov, 2008; Mironov et al., 2010; Stepanenko et al., 2013, Stepanenko et al., 2014) and is a two‐layer integral or bulk fresh water model. FLake separates the lake water into two layers vertically, namely, upper mixed layer where the water temperature is vertically uniform and thermocline layer, which is described using the concept of self‐similarity of the evolving vertical water temperature profile. The water temperature profile in the thermocline can be fairly accurately parameterized by where is time (s), is the depth (m), is the upper mixed layer temperature (K), is the mixed‐layer depth (m), which is computed based on the convective entrainment or relaxation‐type equation in terms of wind mixing (Mironov, 2008). More details for the calculation of are given by Mironov (2008). Δ is the temperature differences across the thermocline with the depth ofΔ, is the lake depth (m), and is the temperature (K) at the bottom of the thermocline. Φ is a dimensionless universal function of the dimensionless depth , which satisfies the boundary conditionsΦ and Φ. Thus, the two‐layer parameterization of the vertical temperature profile at time is given by

The shape function Φ is calculated by where is the shape factor and can be computed by where is the relaxation time (s), which is the time of the evolution of the temperature profile in the thermocline from one limiting curve to the other, following the change of sign in . and are the minimum and maximum values of the shape factor.

In addition, one lake‐ice layer, one snow layer on the lake ice, and two layers in the thermally active layer of the lake sediments are also included in the FLake model. The same basic concept of self‐similarity is also applied to these layers. The water surface albedo with respect to the solar radiation is set to 0.07 in the default configuration, while the albedo of the ice or snow surface is given by , is the albedo of white ice and is the albedo of blue ice. In current study, and. . and are lake surface ice or snow temperature (K) and freezing temperature (K). The solar radiation transfers in water and snow or ice is calculated by a one‐band exponential approximation of the Beer‐Lambert decay law with an extinction coefficient of 3 m−1 for water and 1.0 × 10m−1 for both ice and snow in the default configuration. The parameterized scheme of the turbulent fluxes of momentum, sensible, and latent heat at the lake surface is adopted in the Flake model (Mironov, 2008). As a key parameter for the turbulent fluxes, the roughness lengths over the water or ice/snow surface in the FLake model are computed with respect to the wind velocity, potential temperature and specific humidity (Zilitinkevich et al., 2001; Andreas, 2002). Meanwhile, the parameterized scheme of lake ice and two‐layer parametric representation of the evolving temperature profile in the thermally active layer of bottom sediments proposed by Golosov et al. (1998) are adopted in the FLake model (Mironov, 2008).

The prominent merit of FLake is that it has a small number of specified parameters, namely, the lake depth, fetch, and optical characteristics of water (Lazhu et al., 2016). However, the FLake model is only suitable to simulate the vertical temperature structure and mixing conditions of the lakes with lake depths less than 50 m due to its simple stratification without considering the hypolimnion, which is usually present between the thermocline and the lake bottom in deep lakes (Mironov, 2008; Perroud et al., 2009; Stepanenko et al., 2013). The source codes of FLake and related external data are available at http://www.flake.igb‐berlin.de/sourcecodes.shtml. More detailed descriptions of the FLake model can be referred by Mironov (2008).

WRF‐Lake Model

The lake model called WRF‐Lake, which is derived from the lake component of the WRF model version 3.7 (Xu et al., 2017) and based on the original concept of Hostetler and Bartlein (1990) and Hostetler et al. (1993) (available at http://www2.mmm.ucar.edu/wrf/users/download/get_source.html), is also adopted in this study. It solves the 1‐D thermal diffusion equation by dividing the lake vertical profile into several discrete layers (Subin et al., 2012): snow layers on the lake ice (up to five layers based on the snow depth), water and ice layers (10 layers for global simulations and 25 layers for site simulations), and 10 soil layers in the lake bottom sediment (Gu et al., 2013, Gu et al., 2015, Gu et al., 2016; Wen et al., 2016; Xiao et al., 2016; Xu et al., 2017). The energy balance at the lake surface is given by where is the lake surface absorption fraction of the net surface solar radiation (+downwards). Based on Oleson et al. (2013), the visible wave band (<0.7 μm) penetrates and the near‐infrared wave band (≥0.7 μm) is absorbed by the surface layer. Here the WRF‐Lake assumes that 40% () of the net solar radiation is absorbed in the surface layer and 60% (1 − ) of the net solar radiation penetrates. is the net long wave radiation absorbed by the lake surface (+downwards). and are the sensible heat flux and water vapor flux from lake surface to the atmosphere (+upwards), and G is the ground heat flux from the surface into the lake (+downwards). All the fluxes in equation (5) are determined by the lake surface temperature. converts to an energy according to  K is the freezing temperature.

The lake temperature of each layer is calculated based on a Crank‐Nicholson thermal diffusion solution. The temperature governing equation for the WRF‐Lake model is given by Subin et al. (2012): where T is the temperature (K) at depth z (m), t is the time (s), is the volumetric heat capacity (J · m−3 · k−1) at depth z (the volume‐weighted sum of the respective heat capacities of the water, ice, and mineral constituents). (m2/s) is the molecular diffusion coefficient, (m2/s) is the eddy diffusion coefficient due to wind‐driven eddies (Hostetler & Bartlein, 1990), and (m2/s) is a modest enhanced diffusivity intended to represent the 3‐D mixing processes, which are not explicitly parameterized (Bennington et al., 2014; Fang & Stefan, 1996; Subin et al., 2012). is a factor that depends on the lake depth d (m) and increases the total diffusivity for large lakes to represent 3‐D mixing processes, such as caused by horizontal temperature gradients (Subin et al., 2012). When the water surface temperature is below the freezing point, the eddy diffusion coefficient is turned off in the model (Gu et al., 2015). is the solar radiation (W/m2) penetrating to the depth z. is set to 0.5 (0.6) m at the top of unfrozen lakes in which no additional radiation is absorbed for the lakes with the depths more (less) than 4 m. is the light extinction coefficient (m−1), which is a function of the lake depth d (m) and depends strongly on water quality (Gu et al., 2013; Xu et al., 2017).

The net solar radiation in equation (5) is given by where and are the incident direct beam and diffuse solar fluxes andΛindicate the visible and near‐infrared wave bands, which can be derived from the surface incoming solar radiation by the decomposition methods (Jacovides et al., 1996; Reindl et al., 1990; Zhang et al., 2018). and are the albedos for direct beam and diffuse solar radiation. The albedo for the direct solar radiation at unfrozen lake surface is given by where is the solar zenith angle, while the albedo for the diffuse radiation is set to 0.1. For frozen lake surface with snow depth less than 40 mm, the albedos of both direct and diffuse solar radiation are set to 0.6 for visible radiation and 0.4 for near‐infrared radiation (Subin et al., 2012). For frozen lake surface without resolved snow layers, the albedos for direct and diffuse radiation of ice surface are calculated by where for visible radiation and 0.4 for near infrared radiation. and are the lake surface ice or snow temperature (K) and freezing temperature (K). The albedos of ice surface are restricted to be no less than that defined by equation (9) (Subin et al., 2012; Xiao et al., 2016). For frozen lakes with resolved snow layers, the reflectance of the ice surface is fixed at , and the snow reflectance is computed as over nonvegetated surfaces. Both albedos are combined to obtain the snow‐fraction‐weighted albedo of ice‐snow mixed surface.

Energy exchanges including sensible, latent and radiation fluxes, and momentum exchange between the surface (water, ice, or snow) and overlying atmosphere are calculated during the solution of the surface temperature simultaneously. As key parameters for the calculation of the surface turbulent fluxes, the surface roughness lengths for momentum, heat, and water vapor are set in the default configuration (Gu et al., 2013; Xu et al., 2017) as follows: For unfrozen lake, the lake surface roughness length for the momentum (z 0m), heat (z 0h), and water vapor (z 0q) are set to 0.001 m. For frozen lake with resolved snow, z 0m = z 0h = z 0q = 0.0024 m. For frozen lake without resolved snow, z 0m = z 0h = z 0q = 0.004 m.

CoLM‐Lake Model

We also used the lake scheme called CoLM‐Lake, which is developed based on Zeng et al. (2002) and Subin et al. (2012), and has been adopted in the version 2014 of Common land surface model (CoLM; Dai et al., 2003; Dai et al., 2018; available at http://globalchange.bnu.edu.cn/research/model). The temperature governing equation, main physical processes, and vertical discretization of the snow‐lake‐sediment system in the CoLM‐Lake are similar to the WRF‐Lake model.

Different from the WRF‐Lake model, which uses a prescribed constant lake surface roughness, the CoLM‐Lake model adopts parameterized lake surface roughness lengths of momentum, heat, and water vapor given by Subin et al. (2012) as follows. For unfrozen lake, the lake surface roughness length of momentum (m), the roughness length of heat (m), and the roughness length of water vapor (m; Zilitinkevich et al., 2001). For frozen lake with resolved snow, m and (m). For frozen lake without resolved snow, m and . The variables used for the calculation of the lake surface roughness lengths are described as follows: is the surface friction velocity (m/s) and g = 9.80616 m/s2 is the acceleration of gravity. is the effective Charnock coefficient and and are the maximum and minimum Charnock coefficient, respectively. and define the fetch and depth limitation, respectively. and u is the atmospheric forcing wind. is the lake fetch (m) depending on the lake depth d (m). is set to 1 and can be adjusted as data are available. is the near‐surface atmospheric roughness Reynolds number. is the kinematic viscosity, , , , and is the air pressure at the atmospheric reference height.

In the CoLM‐Lake model, the fraction of nonreflected shortwave radiation absorbed at the surface, β, is set equal to the NIR fraction (≥0.7 μm) predicted by the atmospheric model or forcing data, which is typically ∼0.5 (Subin et al., 2012). Moreover, any shortwave radiation penetrating the snow is absorbed in the top layer of lake ice, 𝛽 = 1 (Subin et al., 2012). More details can be found in the paper of Dai et al. (2018).

Data

The 3‐hourly long‐term (1979–2014) data set with the horizontal resolution of 0.1° developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITPCAS; Lazhu et al., 2016), is adopted in this study to produce the forcing data for the three lake models. This data set includes 2‐m air temperature and specific humidity, 10‐m wind speed, surface pressure, precipitation, and downward solar and long wave radiations. The data can be freely downloaded from http://en.tpedatabase.cn/portal/. In current study, the three lake models are driven by the ITPCAS data corrected by Lazhu et al. (2016) based on the measurements at the weather station with a distance of about 1.5 km from the shoreline in the northeastern Lake Nam Co (Figure 1b).

To validate the model results, we used the daily lake water temperature profile data at the depths of 3, 6, 16, 21, 26, 31, and 36 m, which is observed at the station on the southeast part of Lake Nam Co (30°45.74′N, 90°46.83′E) with the water depth of 92 m (Figure 1b) during 1 January 2012 to 31 December 2013 (Lazhu et al., 2016). In addition, the observed lake surface temperature (LST) at approximately11:00 and 21:00 local time every day during 2012–2013 with a horizontal resolution of 1 km from the Moderate Resolution Imaging Spectroradiometer (MODIS) product (MOD11A1; Savtchenko et al., 2004; Wan et al., 2004) is also used for the model evaluation. We note that much rarer MODIS LST data are observed in cold seasons than in warm seasons due to wrong discrimination of clouds from lake ice/snow. Meanwhile, the bathymetry of Lake Nam Co (Wang et al., 2009) is adopted in current study and shown in Figure 1b.

Numerical Experiment Design

To evaluate the performance of each lake model in simulating the thermal structures of Lake Nam Co, we first carried out a series of control experiments (CTRL) at the temperature profile observation site (Figure 1b) by running each lake model with the default model configurations. Following Subin et al. (2012), we set 25 vertical layers of the lake body in both WRF‐Lake and CoLM‐Lake models, in which the lake depth was set to the real depth of 92 m at the temperature profile observation site (Figure 1b). FLake model is only suitable for the lakes with lake depths less than 50 m due to its simple stratification without considering the hypolimnion (Stepanenko et al., 2013). Thus, the lake depth for Flake is set to 40 m that equals to the mean lake depth in current study.

Based on the results of the previous studies (Wang et al., 2009; Lazhu et al., 2016; Xu et al., 2017; Dai et al., 2018) and the additional runs for testing the key model parameters in each lake model, we further conducted a series of sensitive experiments (Table 1) by tuning some key parameters or replacing some constant parameters (surface albedo and roughness) with parameterized schemes to calibrate and improve the model performance. The key parameters tuned in current study include the light extinction coefficient (0.07–0.17 m−1), temperature of maximum water density (Tdmax; 1.1~4 °C), and mixing factor (0.5–50).

Table 1

Experimental Design of the FLake, CoLM‐Lake, and WRF‐Lake Models

Models	Model parameters	CTRL	SenExp1	SenExp2	SenExp3	SenExp4
FLake	Extinction coefficient (m−1)	3	0.1	0.1	0.1
	Lake surface albedo	0.07	0.07	Subin et al. (2012)	Subin et al. (2012)
	Tdmax (°C)	4	4	4	1.1
CoLM‐Lake	Extinction coefficient(m−1)	1.1925d −0.424	0.1	0.1	0.1
	mixing Factor md	10	10	40	40
	Tdmax (°C)	4	4	4	1.1
WRF‐Lake	Extinction coefficient (m−1)	1.1925d −0.424	0.1	0.1	0.1	0.1
	Mixing factor md	10	10	40	40	40
	Tdmax (°C)	4	4	4	1.1	1.1
	surface roughness (m)	Constant	Constant	Constant	Constant	Subin et al. (2012)

For the FLake model, in addition to the CTRL experiment with the default settings of the three key model parameters including the light extinction coefficient, the lake surface albedo and Tdmax, we also conducted three sensitive experiments by gradually adjusting each parameter mentioned above (Table 1). In the CTRL experiment, the light extinction coefficient is set to 3 m−1, the lake surface albedo is a constant of 0.07, and the Tdmax is set to 4 °C. Because its clear water allows solar radiation to penetrate to deeper depths, the light extinction coefficient of 3 m−1 seems too high for Lake Nam Co. Wang et al. (2009) showed that the light extinction coefficient in Lake Nam Co ranges from 0.07 and 0.17 m−1 with the average of 0.12 m−1. According to Dai et al. (2018), the light extinction coefficient is set to 0.1 m−1 in the sensitive experiments in current study. Study of Wang et al. (2009) showed that the Tdmax is around 3.6 °C in Lake Nam Co with the surface pressure of 570 hPa and water salinity of 1.7 g/L. We empirically set the Tdmax to 1.1 °C after a number of tests to get optimal results in this study. From Figures 7b, 7c, 7e, 7h, and 7i, the simulation of the vertical water temperature profile at the thermocline in summer is slightly affected by the Tdmax. However, the Tdmax of 1.1 °C may not reflect the reality but be efficient to reproduce the vertical water temperature profile in winter. The current 1‐D lake models do not consider the processes such as the horizontal temperature advection, underflows, drag force of the sediment layer, and glacial melt water intrusion (Deng et al., 2013; Martynov et al., 2010; Perroud et al., 2009). Therefore, we set the Tdmax to 1.1 °C as an equivalent Tdmax to reflect the mixing processes not explicitly parameterized in this study. Meanwhile, we experiment with the formulation of Subin et al. (2012) as an alternative to a constant albedo of 0.07.

Figure 7

Time‐depth distribution of the daily mean modeled water temperature in each sensitive experiment of the FLake, CoLM‐Lake, and WRF‐Lake models at the observation site over 2012–2013.

As shown in the Table 1, both the CoLM‐Lake and WRF‐Lake models adopt the parameterized light extinction coefficient (m−1) relying on the lake depth d (m). In current study, the lake depth in the CoLM‐Lake and WRF‐Lake models is set to the real depth of 92 m, which is corresponding to the light extinction coefficient of around 0.175 m−1 default value. Similar to the FLake model, the light extinction coefficient is set to 0.1 m−1 and the Tdmax is set to 1.1 °C in the sensitive experiments of the CoLM‐Lake and WRF‐Lake models. Meanwhile, the default value of the mixing factor m in equation (7) is set to 10 for both CoLM‐Lake and WRF‐Lake models. As Lake Nam Co is a deep lake, the mixing factor of 10 is not large enough to well composite the vertical mixing, which can be enhanced by increasing the mixing factor (Gu et al., 2015). After a number of tests, the mixing factor is set to 40, which is comparable to Dai et al. (2018). Compared to the CoLM‐Lake model, the WRF‐Lake model adopts constant surface roughness lengths for the momentum, heat, and water vapor, which may not reflect the real lake surface conditions but significantly affect the surface turbulent fluxes (Charusombat et al., 2018). Following Subin et al. (2012), the constant lake surface roughness lengths are replaced by the parameterized scheme in the experiment SenExp4 for the WRF‐Lake model (Table 1).

In each experiment, the model simulation starts at 00:00 Beijing time 1 June 2008 and ends at 24:00 31 December 2013 with the integration time step of 10 min. Each lake model is driven by the forcing data with the time interval of 10 min derived from the 3‐hourly corrected ITPCAS data (Lazhu et al., 2016) using a linear interpolation method. Lake Nam Co is ice‐free, and the water is well mixed with almost uniform temperature vertically in early June (Lazhu et al., 2016), so the initial lake water temperature and ice fraction at each model layer can be simply set to 276.65 K and 0 according to the observation, respectively. Considering the effects of the model spin‐up time on the model results, we took the model results during 1 January 2012 to 31 December 2013 with the spin up time of three and half years for analysis.

Methodology

In this study, the temporal correlation (TC) is used to reveal the temporal variation similarity between the model simulation and observation. The root‐mean‐square error (RMSE) is adopted to evaluate the model errors in quantity. The formulas for these statistics (Huang et al., 2016) are given as follows: where is the simulation (observation) at the point in time. is the number of sample. ( ) is the mean value of the simulation (observation) with a sample size of N. High TC indicates large similarity in the temporal variation between the simulation and observation, while low RMSE suggests that the simulation has close agreement with the observation in quantity.

In addition, to validate the model simulation in the amplitude and pattern of variability simultaneously, a measure of skill score from Taylor (2001) is given by where is the TC given by the equation (11) between the observation and simulation. is the temporal standard deviation of the simulation normalized by that of the observation and reveals the similarity of the temporal variability between the simulation and observation. is an achievable maximum correlation (here set as 1). It is clear that the Taylor score (TS) ranges from 0 to 1 and better performance is indicated by higher TS (Huang et al., 2016; Kan et al., 2015).

Results

Validation of Each Lake Model With Default Model Configuration

To reveal the performance of each lake model with the default settings in simulating the thermal structure of Lake Nam Co, comparisons between the simulations of the CTRL experiment for each lake model and the observations are systematically conducted in this section. Figure 2 first gives the daily time series of the LST from the MODIS observations and the CTRL experiment simulations for each lake model in both daytime and nighttime. Compared to the observations, each lake model with the default settings (Table 1) produced distinct errors with relatively larger ones in nighttime than in daytime (Table 2). It is also noted that the LST simulated by the three lake models without any model tuning shows large differences especially during cold seasons. Although the FLake model clearly overestimated the LST in both daytime and nighttime from July to September, it can reasonably reproduce the temporal variation of the LST with a TC of 0.96 (0.85) in daytime (nighttime; Table 2). However, the CoLM‐Lake model apparently underestimated (overestimated) the LST in both daytime and nighttime by 2–8 °C during May to July (August to October) and produced comparable LST relatively to the observation during December (Figures 2a and 2b), implying later start of the thermal stratification and thereafter shorter duration of the thermal stratification compared to the observation (Figure 3b). The WRF‐Lake model shows relatively good ability to simulate the LST in quantity during April to early August. While pronounced underestimation can be noted during middle August to December, especially during the nighttime of November and December the LST is significantly underestimated by more than 10 °C, indicating earlier end of the thermal stratification and thereafter shorter duration of the thermal stratification than the observation (Figure 3c). Overall, among the three lake models with the default model settings, the FLake model performs the best to simulate the temporal variation and intensity of the LST in terms of TC and RMSE. However, the CoLM‐Lake model tends to show the best performance in simulating the amplitude and pattern of the LST variability in terms of TS (Table 2). Meanwhile, each lake model generally shows relatively better performance in daytime than in nighttime in terms of the evaluation statistics of TC, RMSE and TS (Table 2).

Figure 2

The daily time series of lake surface temperature (LST) from the Moderate Resolution Imaging Spectroradiometer (MODIS) observations and the simulations from the CTRL experiment for each lake model at the observation site during daytime and nighttime over 2012–2013.

Table 2

The Verification Statistics Including TC, RMSE, and TS Between the Modeled LST in the CTRL Experiment for Each Model Against the Observation Shown in Figure 2 During Daytime and Nighttime Over 2012–2013

	TC		RMSE (°C)		TS
Lake models	Daytime	Nighttime	Daytime	Nighttime	Daytime	Nighttime
FLake	0.96	0.85	1.66	3.18	0.76	0.62
CoLM‐Lake	0.87	0.84	3.21	3.34	0.89	0.86
WRF‐Lake	0.81	0.59	4.09	9.68	0.87	0.43

Note. The calculation is conducted using the data when the Moderate Resolution Imaging Spectroradiometer (MODIS) estimates are available. CTRL, control experiments; LST, lake surface temperature; RMSE, root‐mean‐square error; TC, temporal correlation; TS, Taylor score.

Figure 3

Time‐depth distribution of the daily mean water temperature from the observations and the simulations of the CTRL experiment for each lake model at the observation site over 2012–2013.

Figure 3 shows the modeled and observed thermal structure at the observation site in Lake Nam Co. From the observation (Figure 3d), Lake Nam Co is a typical dimictic and holomictic lake (Huang et al., 2015; Wang et al., 2019), which is well mixed and characterized by very small vertical temperature gradient during December to the subsequent early June. The annual minimum temperature of the entire water column occurs in February with the values less than 1 °C and then the water temperature starts to gradually increase in the entire water mass until the onset of the thermal stratification in late June in response to the enhanced heating from the surface (Huang et al., 2015). During July to late October, the lake was stratified with a warm upper mixed layer of 10 to 30 m and a cool lower layer. It can also be clearly seen that the thickness of the mixed layer is gradually increased during the thermal stratification period. Meanwhile, the deepening of the thermocline from late October until the water column is mixed again from top to bottom with a homogeneous thermal status in the entire water body due to the surface heat losses and increased vertical mixing associated with the enhanced winds (Huang et al., 2015).

All the three lake models with default model settings can reproduce the overall features of the observed time‐depth distribution of the daily water temperature (Figures 3a–3c). The mixed layer depth produced by CoLM‐Lake and WRF‐Lake increases until the entire column reaches Tdmax, even after the surface temperature begins to decrease, while FLake produces nearly the same time of maximum temperature at all depths and the mixed layer becomes shallower toward the end of the positively stratified season. However, compared to the observation (Figure 3d), apparent model errors can be noted, that is, all the three lake models produced an inverse thermal stratification with relatively larger vertical temperature gradients during December to the following May than the observation. Meanwhile, all the three lake models tend to produce a relatively shallower mixed layer during late September to October (Figures 3a–3c). The WRF‐Lake modeled water temperature in the mixed layer and thermocline is underestimated by 1–3 °C (Figure 3c).

To quantify the performance of each lake model with the default model settings in simulating the water temperature at different depths, Figure 4 further gives the vertical distribution of the TC, RMSE, and TS for the simulated water temperature in the CTRL experiment of each lake model against the observation. All of the three lake models generally show that the model performance in simulating the water temperature temporal variation (Figure 4a) and the amplitude and pattern of variability (Figure 4c) decreases at increased water depth. Rather small TC and TS can be noted at the depth deeper than 30 m, where the TC (TS) is less than 0.8 (0.6) in most cases. Meanwhile, the model's ability to simulate the water temperature magnitude with all of the RMSEs less than 2.5 °C, and the water temperature in the deep waters is relatively better simulated in most cases (Figure 4b). Overall, among the three lake models without any model tuning, the FLake model performs the best to simulate the water temperature at most depths in terms of the verification statistics including TC, RMSE, and TS (Figure 4).

Figure 4

The vertical distribution of the temporal correlation (TC), root‐mean‐square error (RMSE), and Taylor score (TS) for the daily mean water temperature simulated by the CTRL experiment of each lake model against the site observation over 2012–2013.

It should be pointed out that all of the three lake models with the default configuration show large model errors (Figures 2, 3, 4 and Table 2), especially in the thermocline and lower water layers where the verification statistics such as TC and TS are still rather low. As some key parameters in the lake models were estimated based on the observations in plains or wet regions, they may not be suitable for the lakes over the alpine and semiarid to arid areas. To expand the lake models' applicability in the alpine and semiarid to arid regions such as the Tibetan Plateau, there is still much room to calibrate and improve the performance of current lake models to simulate the lake thermal structures based on the field observations. In the following sections, we will systematically reveal the impacts of some key model parameters in each lake model on the model performance based on the simulations from a series of sensitive experiments, and then try to explore a way to improve the model performance in simulating the lake thermal structures over the Tibetan Plateau.

Evaluation of Each Lake Model With the Key Parameters Calibrated

Figure 5 gives the daily time series of the LST from the MODIS observations and the simulations from each experiment of each lake model at the observation site in Lake Nam Co. It is clear that the LST simulations are very sensitive to the key model parameters in each lake model. As shown in Figures 5a and 5b, the FLake model with the default configuration tends to overestimate the LST in both daytime and nighttime during July to October. This overestimation can be significantly reduced by adjusting the light extinction coefficient from 3 to 0.1 m−1 (Table 1). Simulations from the SenExp2 experiment by replacing the constant lake surface albedo with the parameterized scheme suggested by Subin et al. (2012) and from the SenExp3 experiment by setting the Tdmax to 1.1 °C instead of 4 °C are very comparable to those from the SenExp1 experiment. This can be confirmed by the verification statistics such as TC, RMSE, and TS shown in Figures 6a–6c. Compared to the lake surface albedo and Tdmax in the FLake model, the light extinction coefficient affects the vertical water temperature profile through redistributing the solar radiation absorption at different depths (Xu et al., 2017; Dai et al., 2018).

Figure 5

The daily time series of lake surface temperature (LST) from the Moderate Resolution Imaging Spectroradiometer (MODIS) observations and the simulations of each experiment of each lake model at the observation site during daytime and nighttime over 2012–2013.

Figure 6

The temporal correlation (TC), root‐mean‐square error (RMSE), and Taylor score (TS) for the modeled lake surface temperature in each experiment of each lake model against the Moderate Resolution Imaging Spectroradiometer (MODIS) observations in daytime and nighttime. The calculation is conducted when the MODIS data are available over 2012–2013.

From Figures 5c and 5d, the LST simulated by the CoLM‐Lake model shows large differences among the four experiments, the CTRL experiment with the default configuration tends to significantly underestimate (overestimate) the LST during middle June to late July (late August to early November) in both daytime and nighttime. With the light extinction coefficient, mixing factor and Tdmax tuned one by one in the experiments from SenExp1 to SenExp3, the CoLM‐Lake simulated LST gradually becomes more and more consistent with the observations (Figures 6d–6f).

The LST magnitude of underestimation by the CTRL experiment of WRF‐Lake model in both daytime and nighttime during September to December (Figures 5e and 5f) can be slightly reduced by decreasing the light extinction coefficient, increasing the mixing factor and reducing the Tdmax one by one in the experiments from SenExp1 to SenExp3 (Figures 6g–6i). However, significant improvements in the LST simulation during its decline period (late August to December) can be obtained by further adopting the parameterized lake surface roughness lengths instead of constant values in the SenExp4 experiment (Figures 5e and 5f and 6g–6i).

In addition to indicating the impacts of the key model parameters on the LST simulation, revealing how the key model parameters affect the simulations of the lake vertical thermal structure can further deepen our understanding the associated processes such as air‐lake energy exchanges, turbulent mixing, and convective overturning. Figure 7 shows the time‐depth distribution of the daily modeled water temperature in each sensitive experiment of the FLake, CoLM‐Lake, and WRF‐Lake models. Compared to the CTRL experiment (Figure 3a), setting the light extinction coefficient from 3 to 0.1 m−1 can significantly deepen the mixed layer and heating the thermocline during warm seasons by penetrating more solar radiation into the deep waters (Figure 7a). Further adopting the parameterized lake surface albedo instead of constant values and reducing the Tdmax lead to remarkable improvements in simulating the time‐depth distribution of water temperature in cold seasons (Figures 7b and 7c). In the other three experiments with the Tdmax of 4 °C (Figures 3a, 7a, and 7b), the relatively stronger inverse thermal stratification during cold seasons is attributed to the much weaker heat exchanges induced by the suppressed convective mixing from the bottom to the top due to the relatively stronger stability (Subin et al., 2012; Dai et al., 2018). Setting the Tdmax to 1.1 °C in the experiment SenExp3 can enhance the convective mixing induced by the vertical density gradient and further lead to a well‐mixed water column and a uniform temperature throughout the whole water body in cold seasons (Figure 7c), which is much closer with the observations (Figure 3d).

To quantify the performance of each experiment for the FLake model in simulating the water vertical thermal structure, Figure 8 gives the vertical distribution of the TC, RMSE, and TS for the modeled water temperature against the site observation and their percentage changes of each sensitive experiment relative to the CTRL experiment. From Figure 8a, the TC values produced by the four experiments ranging from 0.6 to 0.97 decrease with the water depth increased. And all the four experiments produced larger consistency of TC at the depth shallower than 18 m (Figure 8a) with the percentage changes less than 3% (Figure 8d), indicating that the light extinction coefficient, lake surface albedo, and Tdmax slightly affect the performance of the FLake model in simulating the temporal variation of water temperature at the depth less than 18 m. However, the TC percentage changes range from 8 to 45% at the depth deeper than 30 m with the largest values over 35% at the depth of 36 m (Figure 8d), suggesting that the temporal variation of water temperatures at deep waters is strongly affected by the three parameters. Meanwhile, the TC percentage changes for the SenExp1 and SenExp2 experiments are very comparable to those for the CTRL experiment, indicating that the light extinction coefficient remarkably affect the FLake model performance in simulating the temporal variation of the water temperature at most depths especially in the deep waters.

Figure 8

The vertical distribution of the temporal correlation (TC), root‐mean‐square error (RMSE), and Taylor score (TS) for (a–c) the daily mean water temperature simulated by each experiment of the FLake model against the site observation over 2012–2013 and (d–f) their percentage changes produced by each sensitive experiment relative to the CTRL experiment.

As shown in Figure 8b, the RMSE values of the water temperature modeled by the four experiments of Flake model at different depths range from 0.8 to 2.2 °C. Most experiments produce much larger RMSE values in the thermocline than at other depths, the SenExp1 and SenExp2 experiments produced very comparable RMSE at different depths, suggesting that the FLake model's performance can be slightly affected by adopting the parameterized scheme of the lake surface albedo (Subin et al., 2012) instead of the constant value of 0.07. Compared to the CTRL experiment, the RMSE of the modeled water temperature in the SenExp1 experiment with the adjustment of the light extinction coefficient from 3 to 0.1 m−1 can be significantly reduced (increased) by 15–45% at the layers shallower than 15 m (in the thermocline between 20 and 30 m; Figure 8e). The RMSE of the water temperature produced by the SenExp1 and SenExp2 experiments at the depth less (more) than 15 m can be slightly (significantly) reduced by further setting the Tdmax to 1.1 °C instead of 4 °C in the SenExp3 experiment with the percentage changes of about 2% (5%~45%; Figure 8e), suggesting that the water temperature simulation in the thermocline is much more strongly affected by the Tdmax compared to the other two parameters in Flake model.

From Figure 8c, the TS values at the depth deeper than 27 m in the CTRL experiment of FLake range from 0.4 to 0.8, which can be increased to over 0.8 by the SenExp1 experiment, suggesting that the amplitude and pattern of the temperature variability in the deep layers can be significantly improved by reducing the light extinction coefficient with the percentage changes ranging from 10% to 90% (Figure 8f). With the further adjustments of the lake surface albedo and Tdmax in addition to the light extinction coefficient, the TS for the water temperature modeled by Flake can also be slightly increased.

Overall, setting the light extinction coefficient to 0.1 m−1 instead of 3 m−1 tends to greatly improve the ability of FLake in simulating the water temperature temporal variation and the amplitude and pattern of the water temperature variability at the depth deeper than 25 m and the water temperature in quantity at the depth shallower than 15 m (Figure 8). However, the performance of the FLake model in simulating the water temperature in quantity at the layers deeper than 15 m can be further significantly improved by setting the Tdmax to 1.1 °C instead of 4 °C (Figures 8b and 8e). Meanwhile, among the four experiments, the SenExp3 experiment with all of the three parameters tuned displays the best skill in simulating the water temperature at almost all depths in terms of the verification statistics including TC, RMSE, and TS.

Compared to the observation (Figure 3d), the CoLM‐Lake with default configuration tends to produce much stronger inverse stratification in the cold seasons and much shallower mixed layer and colder thermocline during the warm seasons (Figures 3b and 3d). Reducing the light extinction coefficient and increasing the mixing factor in the CoLM‐Lake model lead to much deeper mixed layer and higher temperature at the deep waters in the warm seasons (Figures 7d and 7e), which is much closer to the observation than the CTRL simulation. In addition, the vertical thermal structure in the cold seasons can be further significantly improved by reducing Tdmax (Figure 7f). Possible reasons are similar to the situations of the FLake simulations.

As shown in Figure 9, with the three parameters gradually tuned one by one in the three sensitivity experiments of the CoLM‐Lake model, the TC and TS (RMSE) values gradually increase (decrease) at almost all depths compared to the CTRL experiment. From Figures 9d and 9f, reducing the light extinction coefficient produces much larger TC (TS) for the modeled temperature with the percentage changes ranging from 5 to 15% (5% to 160%) with relatively larger values at the layers between 15 and 30 m (the depth deeper than 15 m). Further increasing the mixing factor leads to the TS increased by 10% (10 to 80%) at the depth shallower (deeper) than 15 m. Compared to the SenExp1 and SenExp2 experiments, the SenExp3 experiment with reduced Tdmax produces relatively larger TC (TS) values at the depth shallower than 24 m (deeper than 15 m) with the percentage changes of ~5–7% (10–20%). As shown in Figure 9e, the RMSE produced by the CoLM‐Lake model with default settings can be clearly reduced by decreasing the light extinction coefficient with the percentage changes ranging from 6 to 20% at all depths. Meanwhile, increasing the mixing factor and reducing the Tdmax in addition to the deceased light extinction coefficient can further reduce the RMSE by~20–30%.

Figure 9

The vertical distribution of the temporal correlation (TC), root‐mean‐square error (RMSE), and Taylor score (TS) for (a–c) the daily mean water temperature simulated by each experiment of the CoLM‐Lake model against the site observation over 2011–2013 and (d–f) their percentage changes produced by each sensitive experiment relative to the CTRL experiment.

As a whole, among the four experiments of the CoLM‐Lake model, the SenExp3 experiment with all of the three parameters tuned produces the highest TC and TS and lowest RMSE at all depths, indicating the best skill in simulating the temporal variation and quantity of water temperatures and the amplitude and pattern of temperature variability at all depths. The temporal variation of water temperatures and the amplitude and pattern of the temperature variability at the depth deeper than 15 m are strongly affected by both the light extinction coefficient and turbulent mixing. In addition to reducing the light extinction coefficient and increasing the turbulent mixing, further reducing the Tdmax can lead to better simulated water temperature temporal variation at the depth less than 24 m (the amplitude and pattern of the temperature variability at the depth deeper than 15 m; Figures 9d and 9f). Meanwhile, the light extinction coefficient and Tdmax in the CoLM‐Lake model show much larger effect on the model performance in simulating the water temperature in quantity at almost all depths. However, the turbulent mixing mainly affects the CoLM‐Lake model's ability in simulating the water temperature in quantity at the up layers less than 20 m.

Similar to the CoLM‐Lake model, the WRF‐Lake model with default settings (Figure 3c) tends to produce much stronger inverse stratification during the cold seasons and much shallower mixed layer and cooler thermocline in the warm seasons compared to the observation (Figure 3d). The modeled thermal structure in warm seasons can be apparently improved by reducing the light extinction coefficient and increasing the vertical turbulent mixing (Figures 7g and 7h). Meanwhile, the modeled thermal structure in cold seasons can be further significantly improved by reducing the Tdmax, but the water temperature in warm seasons is slightly underestimated (Figure 7i). These errors can be significantly reduced by further adopting the parameterized lake surface roughness lengths (Subin et al., 2012) instead of constant values (Figure 7j).

As shown in Figures 10a and 10d, compared to the CTRL experiment, the SenExp1 experiment with much smaller light extinction coefficient produces larger (smaller) TC at the depth shallower (deeper) than 15 m. This indicates that the temporal variation of the modeled water temperature in the upper (lower) layers becomes relatively better (worse) by decreasing the light extinction coefficient in the WRF‐Lake model. In the experiments from SenExp2 to SenExp4 with the mixing factor, Tdmax, and the lake surface roughness lengths gradually tuned one by one (Table 1), the TC becomes larger and larger at all depths compared to the CTRL and SenExp1, suggesting better and better skill in simulating the water temperature temporal variation.

Figure 10

The vertical distribution of the temporal correlation (TC), root‐mean‐square error (RMSE), and Taylor score (TS)for (a–c) the daily mean water temperature simulated by each experiment of the WRF‐Lake model against the site observation over 2012–2013 and (d–f) their percentage changes produced by each sensitive experiment relative to the CTRL experiment.

Meanwhile, the TC percentage changes produced by the experiments from the SenExp2 to SenExp4 range from 5 to 40% and much larger values are located in the layers deeper than 20 m (Figure 10d). Compared to the other three parameters tuned in the WRF‐Lake model, adopting the parameterized scheme instead of constant surface roughness lengths tends to produce much larger percentage changes ranging from 20 to 40% at all depths (Figure 10d). This suggests that adopting the parameterized surface roughness lengths instead of constant values exhibits remarkable effects on the improvements of the model performance in simulating the water temperature temporal variation at all depths especially in the deep layers among the four parameters tuned in the WRF‐Lake model.

From Figure 10b, the SenExp1, SenExp2, and CTRL experiments of the WRF‐Lake model produce very comparable RMSE at all depths, indicating that the model performance in simulating the water temperature in quantity cannot be apparently improved by tuning the light extinction coefficient and mixing factor in the WRF‐Lake model (Table 1). However, the RMSE values at all depths can be remarkably reduced by further reducing the Tdmax in the SenExp3 experiment with the percentage changes ranging from 10 to 25% relative to the CTRL experiment (Figure 10e). Further replacing the constant lake surface roughness lengths with the parameterized scheme in the SenExp4 experiment (Table 1) can produce much smaller RMSE (less than 1 °C at all depths) than the other four experiments (Figure 10b).

Moreover, the RMSE values are significantly reduced by the SenExp4 experiment with the percentage changes more than 50% at all depths relative to the CTRL experiment (Figure 10e), suggesting that adopting the parameterized lake surface roughness lengths in the WRF‐Lake model can significantly improve the model performance in simulating the water temperature in quantity. Meanwhile, reducing the light extinction coefficient can significantly improve the model performance in simulating the amplitude and pattern of the temperature variability in the depth deeper than 20 m (Figure 10c) with the percentage changes relative to the CTRL experiment ranging from 20 to 90% (Figure 10f), and the improvements increase with the depths. With the other three parameters gradually tuned one by one, the skill becomes better and better with relatively larger improvements in the hypolimnion than in the epilimnion (Figures 10c and 10f). Overall, the SenExp4 with all of the four parameters tuned shows the best skill in simulating the temporal variation and quantity of water temperatures and the amplitude and pattern of temperature variability at all depths among the five experiments of the WRF‐Lake model.

Intercomparisons of the Three Improved Lake Models

As mentioned above, each lake model was calibrated and improved by tuning some key model parameters to simulate the thermal features of Lake Nam Co over the central Tibetan Plateau. In this section, we will further intercompare the performance of the three calibrated and improved lake models in the simulating the lake thermal structures and reveal the suitability of each lake model for the alpine lakes over Tibet. As shown in Figure 11, all the three improved lake models can reasonably reproduce the temporal variation of the LST with the TC values ranging from 0.90 to 0.98 and well simulate the amplitude and pattern of LST variability with the TS values ranging from 0.91 to 0.96 in both daytime and nighttime (Table 3). It is also noted that all the improved lake models produce much smaller RMSE values (Table 3), which are significantly reduced compared to the CTRL experiment of each model (Table 2). Overall, all of the three improved lake models show comparable skill in simulating the LST. Meanwhile, among the three lake models, the FLake model produces the largest TC and TS and smallest RMSE, indicating that the FLake model performs the best to simulate the LST in terms of temporal variation, quantity, and amplitude and pattern of variability, this is similar to the findings of Stepanenko et al. (2010).

Figure 11

The daily time series of lake surface temperature (LST) from Moderate Resolution Imaging Spectroradiometer (MODIS) observations and the simulations of each improved lake model at the observation site during daytime and nighttime over 2012–2013.

Table 3

The Verification Statistics Including TC, RMSE, and TS Between the LSTs Modeled by Each Improved Lake Model Against the Observation Shown in Figure 11 During Daytime and Nighttime Over 2012–2013

	TC		RMSE (°C)		TS
Lake models	Daytime	Nighttime	Daytime	Nighttime	Daytime	Nighttime
FLake	0.98	0.94	1.19	2.10	0.96	0.93
CoLM‐Lake	0.97	0.93	1.53	2.19	0.95	0.91
WRF‐Lake	0.95	0.90	1.69	1.92	0.93	0.91

Note. The calculation is conducted using the data when the Moderate Resolution Imaging Spectroradiometer (MODIS) estimates are available. LST, lake surface temperature; RMSE, root‐mean‐square error; TC, temporal correlation; TS, Taylor score.

From the time‐depth distribution of the daily water temperature from the observations (Figure 3d) and the simulations (Figures 7c, 7f, and 7j) of each improved lake model at the observation site over 2012 to 2013, it is clear that all of the three lake models well reproduce the time‐depth distribution of the water temperature despite of some errors in intensity. Among the three improved lake models, the FLake (WRF‐Lake) model tends to show the best skill in simulating the temporal evolution and quantity of the water temperatures at the depth shallower (deeper) than 10 m (Figures 12a and 12b); meanwhile, the WRF‐Lake performs the best to simulate the amplitude and pattern of the temperature variability simultaneously at all depths in terms of TS (Figure 12c).

Figure 12

The vertical distribution of the temporal correlation (TC), root‐mean‐square error (RMSE), and Taylor score (TS) for the daily mean water temperature simulated by each improved lake model against the site observation over 2012–2013.

Concluding Remarks and Discussions

Based on the observed data, three 1‐D lake models (FLake, CoLM‐Lake, and WRF‐Lake) have been systematically evaluated and calibrated to simulate the thermal features of Lake Nam Co in the Central Tibetan Plateau in this study. In addition, the performance of the three 1‐D lake models is intercompared and the suitability of each lake model to reproduce the evolution of the thermal properties in the Lake Nam Co is further revealed. Main findings are summarized as follows:

The three 1‐D lake models with default model configurations exhibit large errors in the simulated LST and evolution of the water temperature profile. Rather low TC (<0.87 in most cases) and TS (<0.89) and large RMSE (>3 °C in most cases) for the LST simulated by the three lake models indicate poor model performance in simulating the temperature temporal variation, intensity, and amplitude and pattern of temperature variability. Meanwhile, all the three lake models produced an inverse thermal stratification with relatively larger vertical temperature gradients during December to the subsequent May and a relatively shallower mixed layer and much colder thermocline during late September to October compared to the observation.

Studies on the sensitivity of the model performance to some key parameters have shown that the light extinction coefficient strongly affects the simulation of LST and vertical thermal structure during summer. However, the Tdmax (lake surface roughness lengths) in the CoLM‐Lake (WRF‐Lake) model exhibits significant effects on the LST simulation during May to August (August to December). Adjusting the Tdmax from 4 to 1.1 °C in the three lake models consistently leads to much higher skill in simulating the vertical thermal structure of Lake Nam Co during cold seasons. Reducing the light extinction coefficient from the default value of 3 to 0.1 m−1 in the FLake model leads to much deeper mixed layer and warmer thermocline in warm seasons with better agreement with the observations. Meanwhile, the vertical thermal structure modeled by the WRF‐Lake and CoLM‐Lake models can be distinctly improved by decreasing the extinction coefficient and increasing the vertical turbulent mixing during the warm seasons. In addition to reducing the light extinction coefficient, increasing the turbulent mixing, and reducing the Tdmax, further replacing the constant lake surface roughness lengths with a parameterized scheme (Subin et el., 2012) in the WRF‐Lake model can significantly improve the time‐depth distribution of the water temperature simulation in warm seasons.

Further intercomparison of the three calibrated and improved 1‐D lake models indicates that the FLake model performs the best to simulate the temporal evolution and quantity of temperature at the lake surface and epilimnion with the depth shallower than 10 m, while the WRF‐Lake model shows the best performance in simulating the temporal evolution and quantity of the water temperature at the layers deeper than 10 m and the amplitude and pattern of the temperature variability at all depths. In all, the FLake model is computationally the most efficient and shows the best skill in simulating the temporal variation and quantity of the lake near‐surface temperatures relative to the other models, suggesting that it is a good candidate for the applications when a quick and reliable estimation of lake near‐surface temperatures is important, that is, be coupled with weather forecast models or climate models. While, the eddy‐diffusive lake models such as the WRF‐Lake and CoLM‐Lake are suitable for revealing the lake processes or limnological applications due to their better ability to reproduce the amplitude and pattern of the temperature variability.

Overall, the three 1‐D lake models calibrated by adjusting some key model parameters can well reproduce the evolution of the thermal features in Lake Nam Co. Sensitivity analysis has deepened our understanding of some important physical processes related to the thermal vertical profiles, such as the convective mixing induced by the vertical density gradients in cold seasons, decay of the light penetrating into the lake water, the turbulent mixing induced by winds, and the heat and water transfer between lake surface and overlying atmosphere strongly affected by the surface roughness lengths etc. In current study, the thermal structure in cold seasons can be clearly improved due to the enhanced convective mixing via simply reducing the Tdmax in the three lake models.

However, the Tdmax may vary with local air pressure and water salinity (Boehrer & Schultze, 2008; Chen & Millero, 1986; Perroud et al., 2009), whether setting Tdmax to 1.1 °C is suitable for the other alpine lakes over the Tibetan Plateau or not is still unclear. To answer such questions, simulations of the thermal features in more lakes with different sizes, depths, and altitudes over the Tibetan Plateau should be carried out when more observed data are available in the future. Meanwhile, the water temperature profile in warm seasons can be clearly improved by reducing the light extinction coefficient, while we just simply set a constant light extinct coefficient of 0.1 m−1 in each lake model. However, the light extinction coefficient strongly depends on the water transparency (Gu et al., 2013; Xu et al., 2017), so a parameterized scheme of light extinct coefficient related to the water transparency should be developed based on the available intensive observations in the future.

In addition, the mixing factor in the eddy‐diffusive lake models such as the CoLM‐Lake and WRF‐Lake is another important parameter apparently affecting the simulation of the lake vertical thermal structure in the warm seasons (Dai et al., 2018; Gu et al., 2015). Many recent studies mainly relate it to the lake depth (Dai et al., 2018; Gu et al., 2015; Subin et al., 2012; Xu et al., 2017). However, the mixing factor may also be related to other aspects of lakes such as lake size, underflows, drag force of the sediment layer, and horizontal temperature gradients (Deng et al., 2013; Martynov et al., 2010; Perroud et al., 2009). Moreover, all the three 1‐D lake models are fresh water models, but most lakes over the Tibetan Plateau are saline lakes (Song et al., 2013). To expand the applicability of these lake models over the Tibetan Plateau, the impacts of salinity on the physical properties of lake water, such as specific heat capacity, thermal conductivity, freezing point, density, and saturated vapor pressure at the lake surface (Jackett et al., 2006; Sun et al., 2008; Wen, Nidhi, et al., 2015), should be well considered and parameterized within the lake models in the future.