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

**Jiahe Lang 1,2, Yaoming Ma 2,3,4,\*, Zhaoguo Li <sup>1</sup> and Dongsheng Su 1,2**


**Abstract:** Increasing air temperature is a significant 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 flux (LE), sensible heat flux (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 significantly different from December to July.

**Keywords:** Tibetan Plateau; climate warming; lake surface temperature; heat exchange; lake ice phenology

## **1. 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–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

**Citation:** Lang, J.; Ma, Y.; Li, Z.; Su, D. The Impact of Climate Warming on Lake Surface Heat Exchange and Ice Phenology of Different Types of Lakes on the Tibetan Plateau. *Water* **2021**, *13*, 634. https://doi.org/10.3390/ w13050634

Academic Editor: Thomas Meixner

Received: 26 December 2020 Accepted: 23 February 2021 Published: 27 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and 2009. In many cases, the observed rate of lake warming is more rapid than the air temperature increases seen at high latitudes, which may indicate a close relationship between these [6–9]. Lake surface temperature is an important indicator of the lake state and the evaporation process occurs only on the skin layer (~10 µm) of the water body [9,10]. The lake surface temperature is affected not only by the interactions between the atmosphere and the skin layer, but also between the skin layer and the water column [11]. The former interaction is governed by four primary meteorological variables: solar radiation, atmospheric humidity, air temperature, and wind speed; while the latter is controlled by the energy stored in the water body [12,13]. Trends towards earlier break-up and later freezeup dates for lake ice have been observed for most lakes in Canada between 1951 and 2012, with the latter showing a lower degree of temporal coherence than the former [14]. Studies have shown that a shortening of the frozen period can lead to an earlier establishment of the summer thermocline in lakes, which accelerates the warming of the upper lake water [9,15]. The freezing date of a lake is significantly affected by the lake's heat storage, air temperature and other climatic variables such as weed speed, while the melting of lake ice is mainly determined by air temperature and solar radiation [16,17]. Several additional lake responses to climate change have also been widely observed globally, such as more stable stratification and shallower thermocline depths [18–20].

The Tibetan Plateau (TP) is a large distinct geographic region, with a mean elevation that exceeds 4000 m. The TP is sometimes called the Asian Water Tower, and is more sensitive to climate change than global land surfaces are (0.46 ◦C per decade, 1984–2009) [21–23]. The TP has experienced rapid climate change, with surface air warming and moistening, solar dimming, and wind stilling. Precipitation is increasing over the central TP, and is decreasing in the Himalayan region [23,24]. The environmental changes experienced on the TP are mostly associated with rapid surface warming [23]. The annual mean air temperature between 1980 and 2018 obtained from the 95 China Meteorological Administration (CMA) weather stations on the TP is 4.1 ◦C [25]. For a global average temperature increase between 1.5 ◦C and 2.0 ◦C at the end of the century (2100), the increases in the maximum temperatures on the TP are projected to be between 2.34 ◦C and 3.20 ◦C under RCP4.5, and between 2.34 ◦C and 3.14 ◦C under RCP8.5, respectively [26]. These projections are helpful for the sensitivity experiments in our study.

The TP includes 57.2% of all lakes in China, and these account for ~1.9% of the total global lake surface area [27–29]. The lakes on the TP are completely covered by ice for 5–7 months each year, with a total ice area of approximately 5 <sup>×</sup> <sup>10</sup><sup>4</sup> km<sup>2</sup> in 2018 [29,30]. Mixing in ice-covered lakes can be caused by through-flow currents, by oscillations in the ice cover, by heat flow from the sediments, or by solar radiation penetrating the ice (influenced by the ice albedo). In contrast to high-latitude areas, solar radiation is very strong on the TP and there is little snow, and there is almost none on lake surfaces, where snow is blown by the wind, meaning that convection, driven by penetrating solar radiation, is more effective, particularly in the ice melt period [31]. The impact of the rising air temperature on lake ice phenology may be clearer for lakes on the TP. Studies have shown that the area, level and volume of lakes on the TP decreased slightly from 1976 to the mid-1990s, and then increased rapidly [25]. These lakes play an important role in the global water cycle and the Asian monsoon system. They are highly sensitive to global climate change and are therefore an effective indicator of the absence of direct anthropogenic influences [24,32–36]. Lakes can affect local land–atmosphere interactions and regional heat and water budgets, and can impact on local atmosphere boundary layer processes, as shown in numerical climate model simulations [25,37]. The characteristics of the surface energy balance and turbulent exchange differ from lake to lake, depending on the lake area and depth, as well as on meteorological conditions [38,39]. Long-term changes in lake evaporation, and in the latent heat fluxes (LE) and sensible heat fluxes (H), have been shown for Lake Nom Co and Lake Ngoring [40–44].

Almost all the previously mentioned research on lakes on the TP has focused either on lake–atmosphere interactions at a single lake in ice-free conditions by seasonal in

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.

#### **2. Study Area, Data and Methods**

#### *2.1. 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<sup>2</sup> , 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<sup>2</sup> 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].

#### *2.2. Data*

#### 2.2.1. 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–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. *2.2. Data* **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.

2.2.1. 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–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<sup>2</sup> , 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) 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<sup>2</sup> , 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 CO2/H2O 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.

$$\mathbf{H} = \rho\_a c\_p \overline{w \prime T \prime} \tag{1}$$

$$\text{LE} = L\_{\upsilon} \rho\_{a} \overline{w!q!} \tag{2}$$

Here, *ρ<sup>a</sup>* is the air density, *c<sup>p</sup>* is the specific heat of air at a constant pressure, *L<sup>v</sup>* is the latent heat of vaporization, *w*0 is fluctuation of the vertical wind component, and *T*0 and *q*0 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: *α* = *RS*\_dw/*R<sup>s</sup>* \_up, where *R<sup>S</sup>* \_dw and *R<sup>S</sup>* \_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.

#### 2.2.2. 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<sup>2</sup> while Lake Ngoring is 610 km<sup>2</sup> ) 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–58], because of the cool skin effect [12,59].

#### *2.3. 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 *Q*<sup>∗</sup> = *Q*LE + *Q*<sup>H</sup> + *Q<sup>S</sup>* + *Q*GL, where *Q*<sup>∗</sup> is net radiation, *Q*LE is the latent heat flux (evaporative heat flux), *Q*<sup>H</sup> 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<sup>S</sup>* = *CW*∆*TW*/∆*t*, where *C<sup>W</sup>* is the heat capacity of water. The temperature change with time, ∆*TW*/∆*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–70]. Based on these studies, much considerable progress has been made to improve the 1-D lake models [61,70–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.

#### 2.3.1. FLake Model

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]:

$$\frac{\theta\_{\mathbf{s}}(t) - \theta(z, t)}{\Delta \theta(t)} = \Phi\_{\theta}(\zeta) \text{ at } h(t) \le z \le h(t) + \Delta h(t), \tag{3}$$

where *t* is time, *z* is depth, *θs*(*t*) is the temperature of the upper mixed layer of depth *h*(*t*), ∆*θ*(*t*) = *θs*(*t*) − *θb*(*t*) is the temperature difference across the thermocline of depth ∆*h*(*t*), *θb*(*t*) is the temperature at the bottom of the thermocline, and Φ*<sup>θ</sup>* ≡ [*θs*(*t*) − *θ*(*z*, *t*)]/∆*θ*(*t*) is a dimensionless "universal function of dimensionless depth", *ζ* ≡ [*z* − *h*(*t*)]/∆*h*(*t*), that satisfies the boundary conditions Φ*<sup>θ</sup>* (0) = 0 and Φ*<sup>θ</sup>* (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]:

$$a\_{ice} = \left[a\_{max} - (a\_{max} - a\_{min}) \exp\left(-95.6 \left(T\_{f0} - T\_i\right) / T\_{f0}\right)\right] \tag{4}$$

where *αmax* is white ice albedo, *αmax* = 0.6. The *αmin* is blue ice albedo, *αmin* = 0.1. *T<sup>f</sup>* <sup>0</sup> is the freezing temperature (273.15 K), *T<sup>i</sup>* is the ice surface temperature. As the ice surface temperature (*T<sup>i</sup>* ) approaches the freezing temperature (*T<sup>f</sup>* <sup>0</sup> ), 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−<sup>1</sup> for water, and 1.0 <sup>×</sup> <sup>10</sup><sup>7</sup> <sup>m</sup>−<sup>1</sup> 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].

#### 2.3.2. 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. 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.

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

*Water* **2021**, *13*, x FOR PEER REVIEW 8 of 24

**Figure 2.** Lake ice albedo from the in situ and MODIS observations, and calculated using the different albedo parameterization schemes, at Lake Ngoring and Lake Nam Co for (**a**) 4 to 5 January 2014, and (**b**, **c**) 1 February to 31 March 2015. Time is shown in local time. **Figure 2.** Lake ice albedo from the in situ and MODIS observations, and calculated using the different albedo parameterization schemes, at Lake Ngoring and Lake Nam Co for (**a**) 4 to 5 January 2014, and (**b**,**c**) 1 February to 31 March 2015. Time is shown in local time.

#### 2.3.3. Numerical Experiment Design 2.3.3. 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

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.

#### **3. Results**

#### *3.1. 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. *Water* **2021**, *13*, x FOR PEER REVIEW 10 of 24

**Figure 3.** (**a**–**c**) Comparison and (**d**–**i**) evaluation of simulated LST, H and LE with in situ and MODIS observation data for Nam Co: (**a, d, e**) LE, (**b, f, g**) H, (**c, h, i**) LST. A fitted curve is shown for each variable, where the x represents the in situ observations and the y represents the simulated results. The comparison period is from January 2015 to December 2016. **Figure 3.** (**a**–**c**) Comparison and (**d**–**i**) evaluation of simulated LST, H and LE with in situ and MODIS observation data for Nam Co: (**a,d,e**) LE, (**b,f,g**) H, (**c,h,i**) LST. A fitted curve is shown for each variable, where the x represents the in situ observations and the y represents the simulated results. The comparison period is from January 2015 to December 2016.

Nam Co: (**a, d, e**) LE, (**b, f, g**) H, (**c, h, i**) LST. A fitted curve is shown for each variable, where the x represents the in situ observations and the y represents the simulated results. The comparison period is from January 2015 to December 2016.

**Figure 4.** (**a**–**c**) Comparison and (**d**–**i**) evaluation of simulated LST, H and LE with in situ and MODIS observation data for Ngoring: (**a,d,e**) LE, (**b,f,g**) H, (**c**,**h**, **i**) LST. A fitted curve is shown for each variable, where the x represents the in situ observations and the y represents the simulated results. The comparison period is from January 2015 to December 2016.

#### 3.1.1. Lake Surface Temperature

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.

#### 3.1.2. 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−<sup>2</sup> and 23.34 W m−<sup>2</sup> , respectively. The RMSEs for the FLake model simulations of LE and H are 35.14 W m−<sup>2</sup> and 28.25 W m−<sup>2</sup> , 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.

## 3.1.3. 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. *Water* **2021**, *13*, x FOR PEER REVIEW 12 of 24 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.

**Figure 5.** Seasonal variations in lake ice thickness from this study, from Qu et al. 2012 [83] (averaged over 2006–2011), and from Wen et al. 2016 [37] (averaged over December 2012 to March 2013) (**a**) for Lake Ngoring and (**b**) for Lake Nam Co. **Figure 5.** Seasonal variations in lake ice thickness from this study, from Qu et al. 2012 [83] (averaged over 2006–2011), and from Wen et al. 2016 [37] (averaged over December 2012 to March 2013) (**a**) for Lake Ngoring and (**b**) for Lake Nam Co.

variations at the two lakes.

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

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.

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−<sup>2</sup> , 45.17 W m−<sup>2</sup> , respectively.

#### *3.2. 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.

3.2.1. Seasonal Variations in the Effect of Rising Air Temperature on the Lake Surface Temperature and on Heat Fluxes

Figure 6 shows seasonal variations in LST, H and LE after air temperatures have risen by 1.50 ◦C, 2.25 ◦C, 3.00 ◦C, and 3.50 ◦C at Lake Nam Co and Lake Ngoring, from 2012 to 2016. In ten groups of sensitivity tests, we selected four modes (1.50 ◦C, 2.25 ◦C, 3.00 ◦C, 3.50 ◦C), for which there were clear changes, to analyze the response of the seasonal variability at the two lakes to the increase in air temperature. *Water* **2021**, *13*, x FOR PEER REVIEW 14 of 24

**Figure 6.** Seasonal variations in LST, H and LE after the air temperature has risen by 1.5, 2.25, 3 and 3.5 °C (averaged from January 2012 to December 2016) (**a, c, e**) at Lake Nam Co and (**b, d, f**) at Lake Ngoring. **Figure 6.** Seasonal variations in LST, H and LE after the air temperature has risen by 1.5, 2.25, 3 and 3.5 ◦C (averaged from January 2012 to December 2016) (**a,c,e**) at Lake Nam Co and (**b,d,f**) at Lake Ngoring.

2.50 °C at Lake Ngoring, shown with a red line in Figures 7 and 8).

Figures 7 and 8 show the impact of different temperature increases on the seasonal

ence 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

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 <sup>−</sup>16.45 Wm−<sup>2</sup> and <sup>−</sup>14.94 Wm−<sup>2</sup> in November and December, respectively. After Lake Nam Co has frozen, the minimum decrease is in January, and is <sup>−</sup>4.46 Wm−<sup>2</sup> . The minimum decreases at Lake Ngoring reach <sup>−</sup>11.97 Wm−<sup>2</sup> and <sup>−</sup>14.34 Wm−<sup>2</sup> in October and November, respectively. After Lake Ngoring has frozen, the minimum decrease is in December, and is <sup>−</sup>2.49 Wm−<sup>2</sup> . 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.

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−<sup>2</sup> and 13.45 Wm−<sup>2</sup> in May and June, respectively, and the maximum increase reaches 25.89 Wm−<sup>2</sup> in February. The maximum decrease at Lake Ngoring reaches 4.18 Wm−<sup>2</sup> in April, while the maximum increase reaches 14.88 Wm−<sup>2</sup> 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. *Water* **2021**, *13*, x FOR PEER REVIEW 15 of 24 *Water* **2021**, *13*, x FOR PEER REVIEW 15 of 24

**Figure 7.** 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. **Figure 7.** 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. **Figure 7.** 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.

different temperature increases (lines), averaged from January 2012 to December 2016 (**a,c**) at Lake Nam Co and (**b,d**) at Lake Ngoring. The details of the legend are the same as in Figure 7. **Figure 8.** Seasonal variations in the lake ice thickness (shading), and the difference between LST and MLT sensitivities to different temperature increases (lines), averaged from January 2012 to December 2016 (**a,c**) at Lake Nam Co and (**b,d**) at Lake Ngoring. The details of the legend are the same as in Figure 7. **Figure 8.** Seasonal variations in the lake ice thickness (shading), and the difference between LST and MLT sensitivities to different temperature increases (lines), averaged from January 2012 to December 2016 (**a,c**) at Lake Nam Co and (**b,d**) at Lake Ngoring. The details of the legend are the same as in Figure 7.

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

Nam Co, the lake ice is thin enough, and solar radiation will be absorbed by the lake water

temperature rises beyond 2.75 °C. When the air temperature rises beyond 2.75 °C at Lake

Differences between the sensitivity tests can be used to analyze the rate of change for

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. *Water* **2021**, *13*, x FOR PEER REVIEW 17 of 24

**Figure 9.** Plot of changes in LST in summer and winter for air temperature increases of between 1.5 °C and 3.5 °C at Lake Nam Co and Lake Ngoring (summer data are averaged from June, July, **Figure 9.** Plot of changes in LST in summer and winter for air temperature increases of between 1.5 ◦C and 3.5 ◦C at Lake Nam Co and Lake Ngoring (summer data are averaged from June, July, and August 2012–2016, and winter data are averaged from December, January, and February 2012–2016).

and August 2012–2016, and winter data are averaged from December, January, and February 2012–2016). **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 Air Temperature (**℃**) Freezing Date Melting Date Frozen Days Maximum Ice Thickness (m)** 0 09 January 1 April 83 0.300 0.150 11 January 23 March 72 0.226 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.

> 0.175 17 January 22 March 65 0.216 0.200 18 January 21 March 63 0.197 0.225 11 January 21 March 70 0.195 0.250 11 January 21 March 70 0.192 0.275 11 January 21 March 70 0.196 0.300 11 January 20 March 69 0.159

> 0 06 December 17 April 133 0.761 0.150 11 December 15 April 126 0.678 0.175 11 December 15 April 126 0.656 0.200 13 December 15 April 124 0.645 0.225 11 December 15 April 126 0.638 0.250 20 December 06 April 108 0.607 0.275 20 December 05 April 107 0.602 0.300 21 December 03 April 104 0.580 0.325 21 December 04 April 105 0.573

Lake Nam Co

Lake Ngoring


**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).

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

3.2.3. 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−<sup>2</sup> , which occurs at the beginning of February, and the maximum for Lake Ngoring is 47.37 W m−<sup>2</sup> , which occurs at the end of December. The maximum decreases occur at the end of May and at the beginning of April, reaching <sup>−</sup>38.22 W m−<sup>2</sup> and <sup>−</sup>23.95 W m−<sup>2</sup> , 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−<sup>2</sup> and 24.73 W m−<sup>2</sup> , respectively. In contrast to the LE responses, the maximum decreases in H occur at the beginning of March and in mid-November, reaching <sup>−</sup>22.98 W m−<sup>2</sup> and <sup>−</sup>21.39 W m−<sup>2</sup> 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. *Water* **2021**, *13*, x FOR PEER REVIEW 19 of 24

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

**4. Conclusions**

2016): (a, b) LE, (c, d,) H, and (e) lake ice thickness.

the former), for H, LE and ice thickness for Lake Nam Co and Lake Ngoring (averaged from January 2012 to December

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

#### **4. Conclusions 4. 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.

**Author Contributions:** Y.M. conceived and designed the experiments; J.L., Z.L. and D.S. performed the experiments; Y.M. and J.L. analyzed the data; J.L. wrote the main manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (2019QZKK0103), the Strategic Priority Research Program of Chinese Academy of Sciences (XDA20060101), the National Natural Science Foundation of China (91837208, 42075089), the Youth innovation Promotion Association CAS (QCH2019004).

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

**Acknowledgments:** We acknowledge the hard work of Binbin Wang, Shaobo Zhang and Xiaohui Peng in the field observations.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**

