The Impact of Climate Warming on Lake Surface Heat Exchange and Ice Phenology of Different Types of Lakes on the Tibetan Plateau

: Increasing air temperature is a signiﬁcant feature of climate warming, and is cause for some concern, particularly on the Tibetan Plateau (TP). A lack of observations means that the impact of rising air temperatures on TP lakes has received little attention. Lake surfaces play a unique role in determining local and regional climate. This study analyzed the effect of increasing air temperature on lake surface temperature (LST), latent heat ﬂux (LE), sensible heat ﬂux (H), and ice phenology at Lake Nam Co and Lake Ngoring, which have mean depths of approximately 40 m and 25 m, respectively, and are in the central and eastern TP, respectively. The variables were simulated using an adjusted Fresh-water Lake (FLake) model (FLake_ α _ice = 0.15). The simulated results were evaluated against in situ observations of LST, LE and H, and against LST data derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) for 2015 to 2016. The simulations show that when the air temperature increases, LST increases, and the rate of increase is greater in winter than in summer; annual LE increases; H and ice thickness decrease; ice freeze-up date is delayed; and the break-up date advances. The changes in the variables in response to the temperature increases are similar at the two lakes from August to December, but are signiﬁcantly different from December to July.


Introduction
Climate change has received much attention in recent decades. The Intergovernmental Panel on Climate Change (IPCC) Working Group I's contribution to the IPCC Fifth Assessment Report (WGI AR5) [1], specifically the chapter "Observations: atmosphere and surface", shows that the global average surface temperature warmed by 0.85 • C (0.65 to 1.06 • C) between 1880 and 2012. Global surface temperature is likely to change by more than 1.5 • C by the end of the 21st century under almost all representative concentration pathway (RCP) scenarios, and the change is likely to exceed 2 • C under RCP8.5 [1]. Lakes have a near-global distribution. Differences between land and lake surfaces mean that lakes play a unique role in determining local and regional climate, for example through their low albedo, small roughness length, and high heat capacity [2][3][4]. They are considered to be sentinels of climate change [5]. Recent studies have found significant warming for lakes throughout the world, with a mean increasing trend of 0.34 • C per decade between 1985 situ observation, or on the response of lake area, level and volume to climate change. The responses of lake surface heat exchange and ice phenology to increasing air temperature on the TP have rarely been addressed [25].
In this study, our goal is to investigate the response of lake ice phenology, lake surface temperature (LST), sensible heat flux (H) and latent heat flux (LE) to rising air temperature for two different lakes (Lake Nam Co and Lake Ngoring) on the TP. Our results were simulated using the Fresh-water Lake (FLake) model, developed by Mironov [45], and the model was driven by a dataset of long-term in situ observations. We adjusted the model parameterization scheme for lake ice albedo and improved the accuracy of the winter simulations. The simulated results were evaluated against in situ observed LST, H and LE data, and against LST data derived from Moderate Resolution Imaging Spectroradiometer (MODIS) observations. Then, lake ice phenology, LST, H and LE were investigated under different air temperature scenarios. Finally, we analyzed the maximum possible impact of future air temperature increases on these two lakes on the TP.

Study Area
The Lake Nam Co and Lake Ngoring basins are characterized by a cold and semi-arid continental climate, and the thermal structure of the lakes means that they belong to the dimictic lake type. However, the two lakes differ in surface area, depth, latitude, and altitude. Lake Ngoring is the highest large freshwater lake in China, with a mean depth of 17 m and a surface area of 610 km 2 , and is located in the Yellow River source region on the eastern TP (34.46-35.4 • N, 97.3-97.55 • E; 4274 m a.m.s.l.; Figure 1). The average precipitation between early December and early April is only 28.16 mm (1954-2014) [46]. The minimum and maximum air temperatures occurred in January (−14.2 • C) and August (9.0 • C), respectively, and the annual mean air temperature was approximately −1.9 • C (2011-2016). Lake Ngoring is usually completely ice-covered from early December to early April [15]. The thickest ice appears in late February, and was~0.7 m in 2013 and 2016 [46]. Lake Nam Co is the highest large lake on Earth, with an area of 2021.3 km 2 as of 2010 [47], and its mean depth is approximately 40 m. It is located on the southern part of the TP (30.5-30.95 • N, 90.2-91.05 • E; 4710 m a.m.s.l.; Figure 1). The average annual precipitation observed at Nam Co Station amounts to more than 400 mm, and mainly occurs from May to October (Figure 1) [48]. The minimum and maximum air temperatures occurred in January (4 • C) and July (9.3 • C), respectively, and the annual mean air temperature was approximately 0.5 • C (2011-2016). Lake Nam Co is usually completely ice-covered from early January to late March. Since 1978, the persistence of full ice cover for Nam Co Lake has decreased by 19 to 20 days [34].

In Situ Measurements
There are four sets of observations (two weather stations and two monitoring stations) available for the two lakes. The station data (2011)(2012)(2013)(2014)(2015)(2016) are used as long-term forcing to drive the FLake models, which were derived from the Lake Nam Co station and Lake Ngoring station on the lakeshore. The monitoring stations' data (2015)(2016) are used to evaluate the model results. The observation site at Nam Co provides LST, H and LE data (2015)(2016), and the observation site at Ngoring provides H and LE data (2015-2016), as well as lake ice albedo data (2014).
Water 2021, 13, x FOR PEER REVIEW 4 of 24 Figure 1. The topography around Lake Nam Co and Lake Ngoring, and the location of Nam Co station (upper right), and of Ngoring station (lower right). The location of the monitoring station at Lake Ngoring for ice albedo, and at Lake Nam Co for lake surface temperature (LST), sensible heat flux (H) and latent heat flux (LE) is also marked.

In situ Measurements
There are four sets of observations (two weather stations and two monitoring stations) available for the two lakes. The station data (2011)(2012)(2013)(2014)(2015)(2016) are used as long-term forcing to drive the FLake models, which were derived from the Lake Nam Co station and Lake Ngoring station on the lakeshore. The monitoring stations' data (2015)(2016) are used to evaluate the model results. The observation site at Nam Co provides LST, H and LE data (2015-2016), and the observation site at Ngoring provides H and LE data (2015-2016), as well as lake ice albedo data (2014).
The weather station at Lake Nam Co was established on the eastern shore of the lake, about 1.5 km from the shoreline, in 2005 (green five-pointed star in Figure 1). The automatic weather station (AWS) tower records wind speed, wind direction (WD), air temperature, humidity, pressure, precipitation, and downward short-and long-wave radiation. The in situ observation site for Lake Nam Co is situated on an island (an area of approximately 0.18 km 2 , shown as a black triangle in Figure 1). An eddy covariance (EC) system (4.5m above the island surface) and AWS tower (1.52m and 9.52m above the land surface) were established on the island, which is about 10 m from the shore, on July 28th, 2015. The EC system has been applied to all kinds of lakes including several lakes of TP to measure turbulent fluxes (momentum, sensible heat and latent heat flux) [40,49]. The EC system consisted of an open-path CO2/H2O infrared gas analyzer (LI-7500A, LI-COR Biosciences) Figure 1. The topography around Lake Nam Co and Lake Ngoring, and the location of Nam Co station (upper right), and of Ngoring station (lower right). The location of the monitoring station at Lake Ngoring for ice albedo, and at Lake Nam Co for lake surface temperature (LST), sensible heat flux (H) and latent heat flux (LE) is also marked.
The weather station at Lake Nam Co was established on the eastern shore of the lake, about 1.5 km from the shoreline, in 2005 (green five-pointed star in Figure 1). The automatic weather station (AWS) tower records wind speed, wind direction (WD), air temperature, humidity, pressure, precipitation, and downward short-and long-wave radiation. The in situ observation site for Lake Nam Co is situated on an island (an area of approximately 0.18 km 2 , shown as a black triangle in Figure 1). An eddy covariance (EC) system (4.5 m above the island surface) and AWS tower (1.52m and 9.52m above the land surface) were established on the island, which is about 10 m from the shore, on July 28th, 2015. The EC system has been applied to all kinds of lakes including several lakes of TP to measure turbulent fluxes (momentum, sensible heat and latent heat flux) [40,49]. The EC system consisted of an open-path CO 2 /H 2 O infrared gas analyzer (LI-7500A, LI-COR Biosciences) and a three-dimensional sonic anemometer (CSAT3, Campbell Scientific, Inc.). Temperature, humidity and three-dimension wind speeds are measured at a frequency of 10 Hz by the gas analyzer and ultrasonic anemometer, respectively. A water temperature profiler, measuring to a depth of 0.5 m, was installed (90.7979 • E, 30.8107 • N) from July 28th to November 19th in 2015, and from July 7th to November 18th in 2016. The temperature at the 0.5 m depth is used as the LST in our study (Ts, • C).
The lake station at Ngoring was installed in June 2011. Initially, an AWS tower and an EC observation system were situated in the northwestern part of the lake, standing on a submerged rock about 200 m from the shore (35.038 • N, 97.658 • E) [37,41]. These were damaged by ice in the winter of 2011-2012. Another two systems were located on the southwestern lake shore (yellow five-pointed star in Figure 1, 34.918 • N, 97.558 • E) to ensure continuity with records from 2012, and these then developed into Ngoring Station. Air temperature and humidity were measured with a temperature and humidity probe (HMP45C, Vaisala) at a 3 m height. The incoming and outgoing shortwave and longwave radiations were measured with a net radiometer (CNR-1, Kipp and Zonen) 1.5 m above the lake surface. The sensor signals were recorded by a data logger (CR5000, Campbell Scientific, Inc.) at a frequency of 10 Hz. The AWS tower measures air temperature, humidity, precipitation, pressure, wind speed and direction, and downward short-and long-wave radiation. In this study, the AWS data are used as long-term forcing data to drive the lake model for Lake Ngoring. The EC observation system measures H and LE and the data are used to evaluate the model results.
Half-hourly data for H and LE at Nam Co and Ngoring lakes were calculated by processing the high-frequency EC system observations using the "Turbulence Knight 3" software (Windows TK311) (https://zenodo.org/record/20349#, accessed on 27 February 2021). All the relevant corrections (spike removal, time lag compensation, spectral correction, planar fit rotation, and Webb-Pearman-Leuning density correction) were included [40,50]. Sporadic missing data were replaced through interpolation combined with other meteorological variables such as radiation, wind speed, and temperature (specific humidity) differences between the lake surface and air [49]. Sensible heat flux (H) and latent heat flux (LE) are given by Equations (1) and (2), respectively.
Here, ρ a is the air density, c p is the specific heat of air at a constant pressure, L v is the latent heat of vaporization, w is fluctuation of the vertical wind component, and T and q are the temperature and specific humidity fluctuations. To ensure the data quality of the EC observations, we considered many criteria [43,49]. Biermann et al. [13] showed the in situ-observed turbulent heat flux for conditions with wind direction (WD) from the lake surface. Footprint analysis [51] was used to identify observations collected when Lake Ngoring and Lake Nam Co Lake were the dominant source areas. So, we discarded the turbulent heat flux data when wind direction (WD) criteria were not met (Nam Co WD < 135 • and WD > 270 • , Ngoring 35 • < WD < 215 • ), since these fluxes are contaminated with land. In this study, EC system observation data from 2015-2016 are used to validate our lake model results.
In situ observations of ice albedo were recorded over west Ngoring Lake from 3 to 6 January 2014 (red triangle in Figure 1). The instrument used for this observation is a Kipp & Zonen 4-Component Net Radiometer (CNR4) (1.20 m above the ice surface), in which the pyrgeometer and pyranometers measure longwave (4500-42,000 nm) and shortwave (300-2800 nm) infrared radiation, respectively [15]. The albedo (α) measured by the pyranometers over the lake ice surface is obtained from the ratio: α = R S _dw/R s _up, where R S _dw and R S _up are the measured upward and downward shortwave radiances, respectively. In this study, the data were used with the MODIS-observed data to modify the parameterization scheme used to represent lake ice albedo.

MODIS Lake Surface Temperature, Ice Albedo and Snow/Ice Cover Ratio
The MODIS albedo dataset, MCD43A3, from October 2014 to June 2015, was used to adjust the parameterization scheme for lake ice albedo, using the in situ observation data in this study. The MODIS data include albedo measurements at a 500 m spatial resolution and are updated every 8 days. The two lakes are large enough (the area of Lake Nam Co is 2021.3 km 2 while Lake Ngoring is 610 km 2 ) that multiple MODIS footprints can fit. The pixels that cover the two lakes are carefully selected (two pixels within the lake boundary are removed) to ensure that land contamination is not an issue. The information in this product is detailed in Ref. [52]. The MCD43A product used in this study is the White Sky Albedo (WSA), from the short waveband product. Studies have shown that the snow and ice albedo obtained from MCD43A3 have good accuracy for individual regions, with mean biases ranging from 0.06 to 0.07, when compared with in situ observations [15,53]. MCD43A2 is the quality assurance product for MCD43A3, and it records whether each pixel has snow-cover or not (ice-cover). The MCD43A2 data for winters in 2012-2013, 2013-2014, and 2014-2015 are used to analyze the snow/ice cover ratio for Lake Ngoring and Lake Nam Co during the frozen period.
In this study, the MODIS LST product (MOD11A2) is used to evaluate the results simulated by the lake model for 2015-2016. The MOD11A2 product includes two instantaneous observations each day (approximately 11:00 and 21:00 local time), with a spatial resolution of 1 km. Many studies have shown that the MOD11A2 product has good accuracy for individual regions, with mean biases ranging from 0.2 • C to 1.05 • C [12,54,55]. The water surface temperature from the MODIS observations of the skin layer is generally lower than the in situ observation of the mixed layer [56][57][58], because of the cool skin effect [12,59].

Lake Model and Numerical Experiment Design
One-dimensional (1-D) lake models have demonstrated their ability of simulating the lake thermodynamics [60,61]. The driving assumption of the 1-D lake model is horizontal homogeneity [62]. The one-dimensional lake surface energy balance is given by where Q * is net radiation, Q LE is the latent heat flux (evaporative heat flux), Q H is the sensible heat flux, and Q GL is the heat flux across the lake bottom. The heat storage in the lake is determined as Q S = C W ∆T W /∆t, where C W is the heat capacity of water. The temperature change with time, ∆T W /∆t, is integrated over the total depth of the lake. Certain 1-D lake models with different degrees of sophistication in physical processes have been widely developed, including the 2-layer FLake model based on similarity theory [62,63]; turbulence lake models [64,65]; and radiation-diffusion lake models [11]. Many studies have conducted a series of intercomparisons of the available 1-D lake models' performances in simulating the thermal features of different lakes during the past several years [66][67][68][69][70]. Based on these studies, much considerable progress has been made to improve the 1-D lake models [61,[70][71][72][73].
Compared to the relatively comprehensive studies on the evaluation and development of the 1-D lake models over lowlands and wet regions, few annual modeling studies being have previously been done for TP lakes because of the harsh environment conditions in winter. Recently, the lakes on TP have received increasing attention [37,43,70,74,75]. The competences of the FLake, WRF-Lake, and CoLM-Lake models in simulating the thermal features of Lake Nam Co have been evaluated, and FLake performs the best in simulating the temporal evolution and intensity of temperature in the shallow layers [70]. Huang et al. [70] also adjusted three key parameters (temperature of maximum water density, light extinction coefficient, surface roughness lengths) within FLake and improved the model results. Li et al. [75] evaluated the effects of ice and snow albedo, ice and water extinction coefficients on the lake ice phenology, water temperature, and sensible and latent heat fluxes using the LAKE2.0 model. The computation of the FLake is efficient, due to its relatively simple construction, and it performs reasonably across lake categories in predicting both surface temperatures and ice characteristics [61]. In terms of accuracy, previous studies have found small positive biases in the H and LE simulated in the FLake model, as well as good correlations between the LST from the FLake model and in situ observations, and between FLake LST and MODIS observations, for the ice-free period [76,77]. Results from simulations for the freezing period show that the lake temperature in winter has a negative bias [15,42,67,70].
FLake has been applied to typical temperature freezing lakes of North America and Alaska (Bear Lake) for several sensitivity analyses of lake depth, water transparency, explicit snow and snow/ice albedo, snow density and heat conductivity [68,78]. FLake has also been driven by regional climate scenarios, applied to small lakes in Berlin, Germany, for modeling the impact of global warming on water temperature and seasonal mixing regimes [79,80]. However, FLake has never previously been applied to Lake Nam Co and Lake Ngoring for assessing air temperature sensitivity to compare the effects of different warming degrees on the lake. The lake model used in this study is the 1-D Mironov FLake model, which is a freshwater model, and is suggested to be appropriate for lakes with depths of less than 50 m because of its relatively simple stratification scheme. The model requires a small number of lake parameters to be specified: the lake depth and the optical characteristics of the lake. The FLake model is based on the concept of self-similarity (assumed shape) for the vertical temperature profile in the thermocline, and simulate lakes' temperature profiles and surface heat fluxes [62,81,82]. The essence of the concept of self-similarity of the temperature profile θ(z, t) in the thermocline is that the dimensionless temperature profile in the thermocline can be parameterized through a "universal function of dimensionless depth" with reasonable accuracy, that is [62]: , that satisfies the boundary conditions Φ θ (0) = 0 and Φ θ (1) = 1. A snow layer, a lake ice layer, an upper mixed layer, a thermocline layer and a lake sediment layer are considered in the FLake model, and the concept of self-similarity is applied to all layers except the upper mixed layer. The water temperature of the upper mixed layer is nearly vertically uniform. The depth of the mixed layer is described with an entrainment equation for convective conditions and a relaxation-type equation for stable conditions. The water surface albedo, with respect to solar radiation, is set to 0.07 in the default configuration, and the albedo of the ice surface is given by [15]: where α max is white ice albedo, α max = 0.6. The α min is blue ice albedo, α min = 0.1. T f 0 is the freezing temperature (273.15 K), T i is the ice surface temperature. As the ice surface temperature (T i ) approaches the freezing temperature (T f 0 ), the albedo approaches 0.10. The solar radiation transfer between the water and the snow or ice is calculated using a one-band exponential approximation of the Beer-Lambert decay law, with an extinction coefficient of 3 m −1 for water, and 1.0 × 10 7 m −1 for both ice and snow, in the default configuration. The parameterized scheme for the turbulent fluxes of momentum, and sensible and latent heat at the lake surface is adopted in the FLake model [62].

Modification of the Lake ice Albedo Parameterization Scheme in the FLake Model
Previous studies have found small positive biases in H and LE simulated in the FLake model, as well as good correlations between the LST from the FLake model and in situ observations, and between FLake LST and MODIS observations, for the ice-free period [76,77]. Results from simulations for the freezing period show that the lake temperature in winter has a negative bias [15,42,67,70]. Since our study focuses on changes to lake ice phenology in different air temperatures, the accuracy of the model LST and ice phenology in winter is very important. Since there is almost no snow on lake surfaces on the TP in winter, and albedo is a key parameter for calculating lake ice phenology, we modified the parameterization scheme for lake ice albedo to improve the negative model bias for LST in winter. The MODIS-observed average ice surface albedo is 0.15 for lakes on the TP [15]. Based on results from our previous study, we calculated the lake ice albedo from a new albedo parameterization scheme (Equation (4)), and tested the value for α max by replacing the constant value set for the maximum ice albedo (α max = 0.6, Equation (4)) to 0.4 and 0.15 for the two lakes. We conducted a further two tests of Equation (4) by replacing the maximum ice albedo (α max = 0.6, Equation (4)) with 0.15, and setting α ice to 0.15. The in situ observation data for the ice albedo were used to evaluate the results from each test for daily changes at Lake Ngoring (Figure 2a), and the MODIS observation data were used to evaluate the results of the two tests for the frozen period at Lake Ngoring ( Figure 2b) and at Lake Nam Co (Figure 2c). It can be seen from Figure 2 that the albedo parameterization scheme in the FLake model greatly overestimated the ice albedo, and did not describe daily changes well. It is very difficult to accurately describe daily changes in ice albedo. Since the lake ice albedo values from the in situ observations are concentrated around 0.15 for most of the daytime, the ice albedo had a large positive bias when α max = 0.4 (open down triangle in Figure 2), or α max = 0.6 (solid up triangle in Figure 2). The results of the parameterized scheme show a small negative bias when α max = 0.15 (asterisk in Figure 2), compared with the MODIS observation data. The results when α ice = 0.15 (red line in Figure 2a, red dot in Figure 2b,c) are the most accurate, and we therefore selected this scheme as the final albedo parameterization scheme and used this in all following simulations, after evaluating the model results for LST, H and LE for this scheme against the in situ and MODIS observations. parameterization scheme for lake ice albedo to improve the negative model bias for LST in winter. The MODIS-observed average ice surface albedo is 0.15 for lakes on the TP [15]. Based on results from our previous study, we calculated the lake ice albedo from a new albedo parameterization scheme (Equation (4)), and tested the value for α max by replacing the constant value set for the maximum ice albedo (α max = 0.6, Equation (4)) to 0.4 and 0.15 for the two lakes. We conducted a further two tests of Equation (4) by replacing the maximum ice albedo (α max = 0.6, Equation (4)) with 0.15, and setting α ice to 0.15. The in situ observation data for the ice albedo were used to evaluate the results from each test for daily changes at Lake Ngoring (Figure 2a), and the MODIS observation data were used to evaluate the results of the two tests for the frozen period at Lake Ngoring ( Figure 2b) and at Lake Nam Co (Figure 2c). It can be seen from Figure 2 that the albedo parameterization scheme in the FLake model greatly overestimated the ice albedo, and did not describe daily changes well. It is very difficult to accurately describe daily changes in ice albedo. Since the lake ice albedo values from the in situ observations are concentrated around 0.15 for most of the daytime, the ice albedo had a large positive bias when α max = 0.4 (open down triangle in Figure 2), or α max = 0.6 (solid up triangle in Figure 2). The results of the parameterized scheme show a small negative bias when α max = 0.15 (asterisk in Figure 2), compared with the MODIS observation data. The results when α ice = 0.15 (red line in Figure 2a, red dot in Figure 2b and c) are the most accurate, and we therefore selected this scheme as the final albedo parameterization scheme and used this in all following simulations, after evaluating the model results for LST, H and LE for this scheme against the in situ and MODIS observations.

Numerical Experiment Design
In this study, we applied the FLake model to evaluate the influence of air temperature on simulations for two different lakes: Lake Ngoring and Lake Nam Co. To evaluate the influence of different degrees of air temperature rise (1.5 • C to 3.5 • C) on LST, H, LE, and ice phenology at the two lakes, we first carried out a control experiments (CTRL) by running the FLake model with the default model configuration. We then calculated an additional simulation to test the modified ice albedo parameterization scheme (same in both of the two lakes). Based on the results from the previous tests of the ice albedo parameterization scheme, we conducted a series of sensitivity experiments by increasing the air temperature of the input data by different amounts (1.50 • C, 1.75 • C, 2.00 • C, 2.25 • C, 2.50 • C, 2.75 • C, 3.00 • C, 3.25 • C, and 3.50 • C) to investigate the influence of increasing air temperature on the two lakes. The model lake depth is set to equal the mean depths of the lakes (25 m for Lake Ngoring, 40 m for Lake Nam Co). The initial water temperature for the bottom of Lake Nam Co was set to 276.65 K. The value is close to the temperature that has been observed at a 44 m depth in other years. Other parameters that it is not essential to specify will be set as default with no corresponding observations. The forcing data are derived Water 2021, 13, 634 9 of 22 from the observations made at the stations at Lake Ngoring and Lake Nam Co. The main input variables include the surface air temperature, relative humidity, wind speed, air pressure, downward shortwave and longwave radiation, and precipitation from 1 July 2011 to 1 January 2017. The simulations began in July 2011 and ended in December 2016. The integration time step was 30 min. The simulations began in July 2011, instead of January 2012, to allow for model spin-up.

Evaluation of the Simulated Results
Observations of LST, H and LE from 2015 to 2016 were used to validate the results from both the original FLake model and for the FLake_α ice = 0.15 model, for Nam Co and Ngoring, as shown in Figure 3 (Nam Co) and Figure 4 (Ngoring). The results show that the FLake model's performance is sensitive to the ice albedo parameterization scheme, especially in winter. Setting the ice albedo to 0.15, instead of 0.6, in the FLake model consistently led to the improved simulation of H and LE during winter and the ice melt period for both lakes, and improved the simulated LST throughout the year at Lake Nam Co. Reducing the ice albedo in the FLake model results in fewer frozen days and a smaller maximum ice thickness, leading to a better agreement with the in situ observations reported in other studies [33,37,83]. Further intercomparison indicates that FLake_α ice = 0.15 performs better than the standard FLake model in simulating lake ice phenology for both of the two lakes, and in simulating LST for Lake Nam Co; it also simulates LST, H and LE better than the FLake model during the ice melt period for Lake Ngoring.    The simulated LSTs from FLake_α ice = 0.15 are higher in winter, and lower in summer, at both lakes, relative to the original FLake model. The FLake_α ice = 0.15 model simulates LST better than the FLake model at Lake Nam Co, with a mean bias of 2.22 • C, relative to the bias in the FLake model bias of 2.69 • C. The root-mean-square error (RMSE) values for the LSTs simulated at Lake Nam Co in the FLake_α ice = 0.15 and FLake models are 2.70 • C and 3.25 • C, respectively. At Lake Ngoring, the mean LST bias is 3.04 • C for the original FLake model, and 1.94 • C for FLake_α ice = 0.15, and the RMSEs are 3.74 • C and 2.42 • C, respectively. One possible reason for the warm biases simulated for Lake Ngoring in winter is that the ice at Lake Ngoring is thicker than at Lake Nam Co, and the frozen period is longer. This effectively inhibits energy exchange between the lake water and the atmosphere. In the model, all the solar radiation entering the lake via the lake ice remains in the ice layer, resulting in the model's overestimation of the lake surface temperature in winter (ice layer surface temperature). Another possible reason is that Lake Nam Co is deeper than Lake Ngoring. Mixing in ice-covered lakes is caused by many factors, and convection driven by penetrating solar radiation is one effective driver (Bengtsson, 1996). The shallower the lake, the greater the impact of the through-ice solar radiative flux on mixing. Lake Ngoring is much shallower than Lake Nam Co, so it is more sensitive to changes in ice albedo.

Latent Heat Flux and Sensible Heat Flux
The model can simulate seasonal variations in H and LE at the lake surface. The simulated H and LE in FLake_α ice = 0.15 are higher in winter and spring for both lakes (Lake Nam Co: January to May; Lake Ngoring: December to June) than in the original FLake model (Figures 3 and 4). The FLake_α ice = 0.15 model simulates LE and H better than the FLake model during the melt period (March to April) at Lake Ngoring, with RMSEs of 31.76 W m −2 and 23.34 W m −2 , respectively. The RMSEs for the FLake model simulations of LE and H are 35.14 W m −2 and 28.25 W m −2 , respectively. There are no in situ EC observation data for the corresponding ice melt period at Lake Nam Co to facilitate a comparison. Ref. [40] indicated that sublimation during winter is non-zero over Lake Nam Co, as LE was observed to be very high during the ice-covered period in the EC data. In simulations from the original FLake model, LE is zero throughout frozen period, and simulations from FLake_α ice = 0.15 significantly improve on this bias.

Lake Ice Phenology
The simulated number of frozen lake ice days is lower for both lakes in the FLake_α ice = 0.15 simulations than in the original FLake simulations. At Lake Nam Co, the number of frozen days in simulations from FLake_α ice = 0.15 is 82, and there are 101 frozen days in simulations from FLake. At Lake Ngoring, the number of frozen days is 139 in FLake_α ice = 0.15, and 154 in FLake. The simulated lake ice freeze-up date is several days earlier in FLake than in FLake_ α ice = 0.15, and the onset of ice melt occurs half a month earlier in FLake_ α ice = 0.15 than in the original FLake model for both lakes. In situ observations show that freeze-up generally occurs in mid-January at Lake Nam Co, and that ice break-up occurs in early April [33,83]. The average duration of complete freezing (i.e., complete ice cover) at Nam Co is 58.5 days [84]. The lake ice thickness at Lake Nam Co peaks in February, and the maximum thickness in the in situ observations is about 0.4 m, which occurs relatively near to the Nam Co station [83]. The in situ observations of ice thickness at Lake Ngoring show that it remained above 0.6 m from mid-January to early March, with a maximum (0.73 m) measured in late February [37]. Lake ice thickness at Ngoring was measured from December 2012 to March 2013 near the lakeshore, very close to the AWS. In Figure 5, the simulated lake ice phenology is shown to agree well with the in situ observations reported in other studies. to the AWS. In Figure 5, the simulated lake ice phenology is shown to agree well with the in situ observations reported in other studies.
In summary, the comparison of the simulations against in situ and MODIS observations shows that the FLake_α ice = 0.15 model can approximately reproduce the observed lake ice phenology, and improves upon the original FLake model for the freezing and melting periods, but there is no significant improvement in the simulation results for the ice-free period. The simulated LST, H and LE values at Ngoring are higher than the observed values, with mean biases of 3.04 °C , 17.61 W m −2 , 45.17 W m −2 , respectively.

The Influence of Rising Air Temperature at the Two Different Lakes
All the above-mentioned positive biases in LST and LE for Lake Ngoring are present in both the FLake and the FLake_α ice = 0.15 simulations, and have relatively little influ- In summary, the comparison of the simulations against in situ and MODIS observations shows that the FLake_α ice = 0.15 model can approximately reproduce the observed lake ice phenology, and improves upon the original FLake model for the freezing and melting periods, but there is no significant improvement in the simulation results for the ice-free period. The simulated LST, H and LE values at Ngoring are higher than the observed values, with mean biases of 3.04 • C, 17.61 W m −2 , 45.17 W m −2 , respectively.

The Influence of Rising Air Temperature at the Two Different Lakes
All the above-mentioned positive biases in LST and LE for Lake Ngoring are present in both the FLake and the FLake_α ice = 0.15 simulations, and have relatively little influence on seasonal variations for the lake. FLake_α ice = 0.15 simulates the lake ice freezing and melting periods at Lake Ngoring, and the LST and frozen period at Lake Nam Co, better than the original FLake model. We therefore assume that the FLake_α ice = 0.15-simulated LST, H, LE and ice phenology data can be used to analyze changes in seasonal variations at the two lakes.
Although global climate models (such as CMIP5 and CMIP6) can provide climatic variables which can be used as a future climate scenario, it is well known that the simulation errors of global climate system models for the TP are large [23,25]. The forcing data of the control experiment in our study are in situ-observed climatic variables data; these data can provide a more realistic driving field of the model. Because the environmental changes experienced on the TP are mostly associated with rapid surface warming [23], we use offline simulation and only change the air temperature, mainly to compare the effects of different warming degrees on the lake. Here, we use simulations from FLake_α ice = 0.15 to analyze the effects of rising air temperatures on LST, H, LE and ice phenology.  H (Figure 6a,b) is the most sensitive parameter to changes in air temperature, with a decreasing trend for every month as air temperature rises, except for February at Lake Nam Co and December in Lake Ngoring. The difference in H simulated under different air temperatures reaches a maximum in November (before freeze-up), and suddenly decreases when the two lakes begin to freeze. The minimum decreases in H at Lake Nam Co reach −16.45 Wm −2 and −14.94 Wm −2 in November and December, respectively. After Lake Nam Co has frozen, the minimum decrease is in January, and is −4.46 Wm −2 . The minimum decreases at Lake Ngoring reach −11.97 Wm −2 and −14.34 Wm −2 in October and November, respectively. After Lake Ngoring has frozen, the minimum decrease is in December, and is −2.49 Wm −2 . LST rises rapidly in winter with the increasing air temperature, and the change in H at Lake Nam Co in February (Lake Ngoring in December) is opposite to that for other months.

Seasonal Variations in the Effect of Rising Air Temperature on the Lake Surface Temperature and on Heat Fluxes
With increasing air temperature, LE (Figure 6c,d) generally decreases from the end of the frozen period to a period after ice break-up (March to June at Lake Nam Co, April to June at Lake Ngoring), and increases in the other months, especially during the frozen period. The maximum decreases at Lake Nam Co reach 16.00 Wm −2 and 13.45 Wm −2 in May and June, respectively, and the maximum increase reaches 25.89 Wm −2 in February. The maximum decrease at Lake Ngoring reaches 4.18 Wm −2 in April, while the maximum increase reaches 14.88 Wm −2 in December.
LST (Figure 6e,f) increases every month as the air temperature rises, and two relatively large peaks appear in winter and summer, with the largest increase in winter. The maximum increases at Lake Ngoring reach 3.02 • C and 3.40 • C in December and January, respectively, while the maximum increase at Lake Nam Co reaches 2.64 • C in January. Figures 7 and 8 show the impact of different temperature increases on the seasonal changes in LE, H, LST, and lake mixed layer temperature (MLT). Here, we use the differences between the different sensitivity tests (line in Figures 7 and 8) to analyze the influence of different increase intervals on seasonal variations in the four variables. The differences before and after the frozen period fluctuate significantly. The variability in the sensitivities for the four variables in the ice-free period is roughly the same at Lake Ngoring and at Nam Co, but the magnitude of the sensitivities is small. Since the frozen period is longer at Lake Ngoring, we can see differences between the variations in the sensitivities at the two lakes during the frozen period. In the middle of the frozen period, there is a relatively regular one-month (February) change in LST, LE and H with increasing air temperature at Lake Ngoring, and we do not see this at Nam Co. The thicker ice at Lake Ngoring means that the MLT is unaffected by the increase in air temperature for the three months after freeze-up. However, the MLT at Lake Nam Co changes significantly once the air temperature rises beyond 2.75 • C. At both lakes, there is a clear difference between the simulations with unchanged air temperature, and the simulations when the air temperature has been increased by 1.5 • C, shown with a black line in Figures 7 and 8. The next obvious difference is between the simulations calculated with temperature increases of 2.75 • C and 3.00 • C at Lake Nam Co, shown with blue line in Figures 7 and 8 (2.25 • C and 2.50 • C at Lake Ngoring, shown with a red line in Figures 7 and 8).
Differences between the sensitivity tests can be used to analyze the rate of change for the lake surface and MLT for different intervals of increasing temperature. The larger the positive (negative) value, the faster the LST rises (drops) over the temperature interval. The variations in the differences for LST and MLT are the same in the ice-free period, with relatively regular changes except for the months before and after the ice has frozen. After freeze-up, the thick ice layer isolates the lake body from the air, and the responses of the mixed layer to different air temperature increases become increasingly similar, until they become the same, i.e., the sensitivity to the magnitude of the change in air temperature disappears. This happens at Lake Ngoring, but differences appear again between the mixed layer responses to different air temperature increases at Lake Nam Co once the air temperature rises beyond 2.75 • C. When the air temperature rises beyond 2.75 • C at Lake Nam Co, the lake ice is thin enough, and solar radiation will be absorbed by the lake water through the ice layer, and then affect the mixed layer temperature.  , d) at Lake Ngoring. The plotted differences are the differences between the sensitivities found for different air temperature changes, as detailed in the legend, i.e., (+1.5 °C )-(+0 °C ) refers to the difference between the sensitivity found for a temperature increase of 1.5 °C , and for a temperature increase of 0. The y coordinates for each line show how much the variable has changed within this air temperature increase interval. The gray-shaded area and the left-slash-filled area show the simulated values for no change in air temperature, and for an air temperature increase of 3.5 °C , respectively. Differences between the sensitivity tests can be used to analyze the rate of change for the lake surface and MLT for different intervals of increasing temperature. The larger the positive (negative) value, the faster the LST rises (drops) over the temperature interval. The variations in the differences for LST and MLT are the same in the ice-free period, with relatively regular changes except for the months before and after the ice has frozen. After freeze-up, the thick ice layer isolates the lake body from the air, and the responses of the mixed layer to different air temperature increases become increasingly similar, until they become the same, i.e., the sensitivity to the magnitude of the change in air temperature disappears. This happens at Lake Ngoring, but differences appear again between the mixed layer responses to different air temperature increases at Lake Nam Co once the air temperature rises beyond 2.75 °C. When the air temperature rises beyond 2.75 °C at Lake Nam Co, the lake ice is thin enough, and solar radiation will be absorbed by the lake water through the ice layer, and then affect the mixed layer temperature. Seasonal variations in lake ice thickness (shading), and differences between the sensitivity of H and LE to different air temperature increases (lines), averaged from January 2012 to December 2016 (a,c) at Lake Nam Co and (b,d) at Lake Ngoring. The plotted differences are the differences between the sensitivities found for different air temperature changes, as detailed in the legend, i.e., (+1.5 • C)-(+0 • C) refers to the difference between the sensitivity found for a temperature increase of 1.5 • C, and for a temperature increase of 0. The y coordinates for each line show how much the variable has changed within this air temperature increase interval. The gray-shaded area and the left-slash-filled area show the simulated values for no change in air temperature, and for an air temperature increase of 3.5 • C, respectively.  , d) at Lake Ngoring. The plotted differences are the differences between the sensitivities found for different air temperature changes, as detailed in the legend, i.e., (+1.5 °C )-(+0 °C ) refers to the difference between the sensitivity found for a temperature increase of 1.5 °C , and for a temperature increase of 0. The y coordinates for each line show how much the variable has changed within this air temperature increase interval. The gray-shaded area and the left-slash-filled area show the simulated values for no change in air temperature, and for an air temperature increase of 3.5 °C , respectively.  Figure 7.
Differences between the sensitivity tests can be used to analyze the rate of change for the lake surface and MLT for different intervals of increasing temperature. The larger the positive (negative) value, the faster the LST rises (drops) over the temperature interval. The variations in the differences for LST and MLT are the same in the ice-free period, with relatively regular changes except for the months before and after the ice has frozen. After freeze-up, the thick ice layer isolates the lake body from the air, and the responses of the mixed layer to different air temperature increases become increasingly similar, until they become the same, i.e., the sensitivity to the magnitude of the change in air temperature disappears. This happens at Lake Ngoring, but differences appear again between the mixed layer responses to different air temperature increases at Lake Nam Co once the air temperature rises beyond 2.75 °C. When the air temperature rises beyond 2.75 °C at Lake Nam Co, the lake ice is thin enough, and solar radiation will be absorbed by the lake water through the ice layer, and then affect the mixed layer temperature. Changes in H are closely related to changes in LST. The LST for Lake Nam Co from the end of July to mid-January follows a generally increasing trend, when the air temperature increases from 2.00 • C to 2.75 • C, and differences in LST are more obvious with temperature increases of between 0 • C and 1.5 • C. When the air temperature rises from 1.50 • C to 2.00 • C, and above 2.75 • C, the LST decreases with the increasing air temperature for three months after the ice break-up, and the most significant decrease occurs in the last month. When the lake begins to freeze-up in January, the differences appear to be irregular, and fluctuate until the end of February. The LST at Lake Ngoring from the end of June to the end of November follows a steadily increasing trend for air temperature increases from 0 • C to 3.5 • C, and the differences in LST are more obvious for a temperature increase of 0 • C to 1.5 • C. When the air temperature rises from 1.5 • C to 2.25 • C, and from 2.75 to 3 • C, LST decreases with increasing air temperature during two months after ice break-up. When the lake begins to freeze-up in November, the differences appear to fluctuate until the end of March. By February, the rising LST trend has stabilized, and then there is a rapid change in March, when the ice is about to begin break-up.
3.2.2. The Effect of Rising Air Temperature on Lake Surface Temperature and Ice Phenology Figure 9 shows the impact of the change in air temperature on the simulated LST in summer and winter, and Table 1 shows the impact on the simulated lake ice phenology. As expected, the LST rises with increasing air temperature. When the air temperature rises, the number of frozen days decreases and the maximum ice thickness decreases, meaning the freeze-up date of lake ice will be delayed, and the break-up date will be advanced.   The LST in summer and winter increases with the increases in air temperature. However, the LST rises much faster in winter than in summer, and the increase is greater at Lake Ngoring than at Lake Nam Co in both seasons. When the air temperature rises by 3.50 • C, the LST at Lake Nam Co rises by 1.90 • C in winter, and by 0.57 • C in summer, while the LST at Lake Ngoring rises by 2.88 • C in winter, and by 0.92 • C in summer. The pattern with which LST rises with increasing air temperature is different in winter and summer. Once the winter air temperature rises beyond 2.75 • C at Lake Nam Co, LST rises faster with increasing air temperature. The same thing occurs at Lake Ngoring when the winter air temperature rises beyond 2.25 • C. However, when the air temperature increases from 2.75 • C to 3.00 • C at Lake Nam Co, the LST decreases from 0.76 • C to 0.50 • C, and at the same time, the air temperature begins to affect the temperature of the mixed layer. Table 1. Changes in lake ice thickness, the number of frozen days, and dates for the onset of freeze-up and melt when the air temperature increases by between 1.5 • C and 3.5 • C at Lake Nam Co and Lake Ngoring (averaged from January 2012 to December 2016).

Lake
Increasing When the air temperature rises from 0 • C to 3.50 • C, the number of frozen days at Lake Nam Co decreases from 83 days to 62 days, while the number of frozen days at Lake Ngoring decreases from 133 days to 104 days; the maximum ice thickness at Lake Nam Co reduces from 0.3 m to 0.1 m while the maximum ice thickness at Lake Ngoring reduces from 0.76 m to 0.56 m; ice freeze-up is delayed by 9 days, and the onset of ice melt advances by 12 days at Lake Nam Co, while freeze-up is delayed by 15 days and the onset of ice melt advances by 14 days at Lake Ngoring. When the air temperature rises from 2.25 • C to 2.50 • C, the frozen period at Lake Ngoring changes more obviously than it does between cooler temperatures. Specifically, when the air temperature at Lake Ngoring rises from 0 • C to 2.25 • C, and from 2.50 • C to 3.50 • C, the number of frozen days decreases by 7 days and 4 days, respectively, and when the air temperature rises from 2.25 • C to 2.50 • C, the number of frozen days decreases by 18 days. There is no obvious rule to relate the change in the duration of the frozen period at Lake Nam Co with rising air temperature. The processes that determine the duration of the frozen period are relatively complicated, particularly around the time of freeze-up and break-up. Lake Ngoring has a relatively long and stable freezing period, and we speculate that temperature increases must reach a threshold of around 2.25 • C before they influence the duration of the frozen period at Lake Ngoring.

The Maximum
Possible Impact of Rising Air Temperature on the TP On Lake Nam Co and Lake Ngoring Figures 10 and 11 show the impact of air temperature rises of 3.5 • C at the two lakes. Here, we use the differences between responses to air temperature increases of 0 • C and 3.5 • C to analyze the influence of the most significant warming on seasonal variations in five variables (H, LE, LST, MLT, and ice thickness) for the two lakes. In general, when the air temperature rises by 3.5 • C, annual LE and LST increase and H decreases at both lakes (Figure 10b,c, Figure 11b,c). The increase in LE and LST at Lake Ngoring is higher than at Lake Nam Co, and the decrease in H at Lake Nam Co is greater than at Lake Ngoring. The ice thickness at Lake Ngoring decreases more than at Lake Nam Co. The changes in the variables in response to the temperature increases are similar at the two lakes from the end of July to the beginning of December. From mid-December to the end of July, the two lakes respond differently to the temperature increase, due to their different frozen periods and the continuous impact on lake surface energy exchange after the lake ice melts. The difference between the responses of H to the temperature change at the two different lakes is greater than the difference between the LE responses at the two lakes. The maximum value for the increase in LE at Lake Nam Co is 39.52 W m −2 , which occurs at the beginning of February, and the maximum for Lake Ngoring is 47.37 W m −2 , which occurs at the end of December. The maximum decreases occur at the end of May and at the beginning of April, reaching −38.22 W m −2 and −23.95 W m −2 , respectively. The maximum values for the increases in H at Lake Nam Co and Lake Ngoring are also at the beginning of February and at the end of December, when they reach 18.08W m −2 and 24.73 W m −2 , respectively. In contrast to the LE responses, the maximum decreases in H occur at the beginning of March and in mid-November, reaching −22.98 W m −2 and −21.39 W m −2 for Lake Nam Co and Lake Ngoring, respectively. Similar to H and LE, the maximum increases in surface temperature at Lake Nam Co and Lake Ngoring are at the beginning of February and the end of December, when they reach 5.06 • C and 6.58 • C, respectively. The maximum decrease in ice thickness at Lake Nam Co and Lake Ngoring occurs at the beginning of February and the end of March, reaching −0.24 m and −0.36 m, respectively. The possible reasons can be summarized as the following points. Firstly, the model lake depth is set to 25 m for Lake Ngoring and 40 m for Lake Nam Co. Secondly, the meteorological conditions of the two lakes are different, such as air temperature, precipitation and wind speed.   Figure 10. Comparison of changes in the differences, when the air temperature rises by 0 °C and by 3.5 °C (the latter minus the former), for H, LE and ice thickness for Lake Nam Co and Lake Ngoring (averaged from January 2012 to December 2016): (a, b) LE, (c, d,) H, and (e) lake ice thickness. Figure 11. Comparison of changes in the differences, when the air temperature rises by 0 °C and by 3.5 °C (the latter minus the former), for LST, MLT and ice thickness, for Lake Nam Co and Lake Ngoring (averaged from January 2012 to December 2016): (a, b) LE, (c, d,) H, and (e) for lake ice thickness. Figure 11. Comparison of changes in the differences, when the air temperature rises by 0 • C and by 3.5 • C (the latter minus the former), for LST, MLT and ice thickness, for Lake Nam Co and Lake Ngoring (averaged from January 2012 to December 2016): (a,b) LE, (c,d) H, and (e) for lake ice thickness.

Conclusions
The model simulations show that LST rises with increasing air temperature, the rate is much faster in winter than in summer, and the increase is greater at Lake Ngoring than at Lake Nam Co. When the air temperature rises, the number of frozen days is reduced, the maximum ice thickness decreases, the freeze-up date for the lake ice is delayed, and the break-up date advances. From the end of July to the beginning of December, the variables changed similarly for the two lakes, and from mid-December to the end of July, due to their different frozen periods and the continuous impact on the lake surface energy exchange after ice melt, the two lakes respond significantly differently to changes in temperature. The difference in the response of H is greater than the difference in the response of LE between Lake Ngoring and Lake Nam Co.
It is interesting that LST increases in summer and winter with rising air temperature, but it increases much faster in winter than in summer, especially as the lake begins to freeze, and it increases more at Lake Ngoring than at Lake Nam Co. This is completely consistent with the fact that the increase in surface temperature on the TP is higher in winter than in summer. The mechanism for surface warming on the TP has not yet been resolved, so in future work, we will explore a plateau warming mechanism that includes multiple freezing and thawing processes, such as glaciers, frozen soil, and snow. In combination, these phenomena may explain the freezing period response of lake ice.
Another interesting phenomenon is that when the air temperature rises from 2.25 • C to 2.50 • C, the freezing period at Lake Ngoring changes more obviously than following similarly sized increases between cooler air temperatures. There is no obvious rule to relate the duration of the frozen period at Lake Nam Co with rising air temperature. The physical processes that determine the duration of the frozen period are relatively complicated, especially during the times for ice freeze-up and break-up. Lake Ngoring has a relatively long and stable freezing period, and we speculate that temperature increases must reach a threshold of 2.25 • C before they influence the duration of the frozen period at Lake Ngoring. Our future work will try to investigate this by carrying out model simulations for other lakes on the TP. We will also focus on a wider range of climate change indicators, including wind speed, precipitation, and humidity, and try to analyze the impact of future climate change on lakes on the TP.
One aspect of this work that could be improved is that the diurnal pattern in the simulated albedo does not match the typical U-shaped curve seen in the observations. This is because the albedo parameterization scheme in the FLake model was developed from the sea ice albedo parameterization scheme, and its application to lakes requires further improvement. Our future work will establish long-term observation sites at multiple lakes to obtain sufficient observation data to improve the parameterization scheme for ice albedo, and thereby make it more suitable for lakes on the TP. Another aspect of this work that could be improved is that we use offline simulation, and only change the air temperature, mainly to compare the effects of different warming degrees on the lake. This method is commonly used, even though it has some obvious flaws, but the uncertainty is relatively small. In our future study, we will further use reginal climate models to comprehensively explore the effects of the offline and coupled climate models on the simulation results of the method.

Data Availability Statement:
The input data used for FLake model simulations and the in situ observed data used for evaluate the simulated results presented in this study are available upon reasonable request from the corresponding author.