*2.2. Data*

### 2.2.1. Observations Data in NL and HSP

The NL lake temperature below the surface was observed (at a distance of ~2 km from the shore) using a Campbell 109 L logger sensor from June to October 2012, and from May to September 2013, and using RBR SOLO sensor from September 2015 to September 2016. The lake ice thickness was manually measured at irregular intervals from December 2012 to March 2013 and from December 2015 to March 2016 near the shore. Water temperature in HSP was measured (~1 km from the shore) under the lake surface from September 2015 to April 2016 using HOBO Water Temperature Pro v2 Data Loggers U22-001. These data were used to evaluate the performance of the lake model in the TP.

### 2.2.2. ITPCAS Data and Its Correction

Data for the lake model forcing were obtained from the China meteorological forcing dataset (1979–2018) (ITPCAS), developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences [26]. The data include air temperature and specific humidity at 2 m above the surface, wind speed, surface air pressure, precipitation as well as downward shortwave and longwave radiation. The spatial resolution is 0.1◦ and the temporal resolution is 3 h [27]. The dataset was produced by merging a variety of data sources, including Princeton reanalysis data, Global Land Data Assimilation System data, the Global Energy and Water Cycle Experiment-surface radiation budget shortwave radiation dataset, Tropical Rainfall Measuring Mission satellite precipitation analysis data and China Meteorological Administration (CMA) station data.

The trend of annual mean ITPCAS air temperature (Ta) in the period of 1979–2016 was 0.77 ◦C/decade (*p* < 0.01), i.e., higher than that measured at the Maduo meteorological station (0.55 ◦C/decade (*p* < 0.01)). A closer inspection of the data revealed that the warming rate, i.e., 0.38 ◦C/decade, was identical in the two datasets for 1979–1997, with a systematic cold bias in the ITPCAS data compared to observations (Figure 2). The difference reversed in 1998–2006 but returned in 2007. Consequently, the forcing data from ITPCAS was bias-corrected based on the monthly observational data from the Maduo Station from 1979 to 2016 before being used to drive the lake model.

**Figure 2.** Monthly air temperature difference between ITPCAS and observations at the Maduo station (ITPCAS-Observation).

### 2.2.3. MODIS Data

The global 8-day composite daytime and nighttime LSWT from MODIS AQUA product data MYD11C2 V006 (0.05◦ resolution) and MYD11A2 V006 (1 km resolution) from

2003–2016 were used to evaluate the results of the long-term simulations for NL and HSP, respectively. Albedo from MODIS MCD43B3 V006 with 1-km resolution was used to set up the lake model.

### *2.3. Lake Model and Setup*

In this study, we used an enhanced version of the Hostetler 1-D lake thermal model code, known as CLM4-LISSS. It was originally the lake scheme coupled in CLM, which is the land component of the Community Earth System Model (Details may be found in Subin et al., 2012). The lake scheme includes 0–5 snow layers, 10 lake liquid water and ice layers, and 10 sediment layers, which were also applied in WRF and the RegCM (Regional Climate Model) [28,29], etc.

CLM4-LISSS has been shown to provide reasonable performance in simulations of lake temperature, surface energy fluxes, and ice and snow thicknesses for several different-sized lakes around the world [30–33]. The model has also been improved in terms of ice albedo and mixing process, etc., for the TP [34,35]. It has been shown to effectively simulate the amplitude and pattern of temperature variability at all depths [2,7,36]. However, there is little progress for the saline lake model in the TP. Therefore, we first applied the developed CLM4-LISSS model with salinity parameterizations (evaluated in the Great Salt Lake in USA) to a simulation study of a saline lake (HSP) in the TP [17,30,37]. Further, the change of Tmaxd caused by dissolved salt in the water was parameterized in the current salinity scheme because it could significantly change the thermal processes of the lake water [2,24].

### 2.3.1. Lake Model with Salinity Parameterizations

The following equations, except for Tmaxd, were incorporated into the developed CLM4-LISSS lake model to evaluate for the salinity [17,37].

The dependence on salinity *s* (‰) of the freezing point (Tf, ◦C) was approximated by the seawater formula [38]:

$$\text{Tf} = -0.0575 \text{ s} \tag{1}$$

The effect of dissolved salts on evaporation was expressed by the ratio of the saturated vapor pressure Rsvp over saline water to that over freshwater as [39]:

$$\text{Rsvp} = \exp\left(-2/55.51 \times \text{(s/(1 - s/1000)/58.44 + 0.77)\right) \tag{2}$$

The specific heat capacity of saline water cpsw (kJ/kg/K) was determined as [17,37]:

$$x\_{\text{psw}} = 4.188 - 4.4 \,\text{s/1000} \tag{3}$$

The thermal conductivity of saline water λsw was determined as [17,37]:

$$
\lambda\_{\rm sw} = \lambda\_{\rm fw} \left[ 1.0 - 0.22 \,\mathrm{s/1000} + 0.1 \mathrm{(s/1000)^2} \right] \tag{4}
$$

where λfw is the freshwater thermal conductivity, equal to 0.6 W <sup>m</sup><sup>−</sup>1·K.

In addition to the previously used salinity parametrizations [37], the effect of salinity on the temperature of the maximum density (Tmaxd) of saline water was introduced in the study as Tmaxd decrease with increased salinity [40]:

$$\text{Transxd} = 3.98 - 0.216 \,\text{s} \tag{5}$$

2.3.2. Model Parameters for the Freshwater NL

In the model settings for NL (Table 1), the fraction of the net solar radiation absorbed near the lake surface was set to β = 0.5 in the absence of snow, and derived from the snow optics sub-model when snow was present [30].


**Table 1.** Numerical experiments S-(Lake) and model parameters.

The light extinction coefficient η was modelled as

$$
\eta = \eta\_{\rm{O}} \text{ d}^{-0.424} \tag{6}
$$

where d is the lake depth (m). The η0 parameter was set as 1.1925 in the lake scheme of the CLM model [30].

The lake albedo was fixed as 0.06 for open water conditions according to observations at the NL station without considering the diurnal change, as in CLM4-LISSS, and calculated for an ice-covered lake with the following function [30,41]:

$$\alpha = \alpha\_{\text{max}} - \alpha\_{\text{max}} \ge + \alpha\_{\text{min}} \ge \text{\text{\textdegree}} \left( -\theta 5.6 \left( \text{Tf} - \text{LSWT} \right) / \text{Tf} \right) \tag{7}$$

where αmax and αmin are the max and min values of the lake ice albedo, respectively. In CLM4-LISSS, there are different (<sup>α</sup>max, αmin) values for near infrared and visible radiation. Without making a distinction between the two radiation types in the study, αmax and αmin were set to 0.6 and 0.1, respectively [30,41].

In the lake model, NL depth was set to 17 m, i.e., the same as the mean lake depth. Variations of lake depth were not considered, for two reasons: (1) Some lake models, such as the General Lake Model [15], are capable of considering changes in lake depth, but they have not ye<sup>t</sup> been coupled with atmospheric models. In recent lake–air coupled simulation studies, the lake model CLM4-LISSS and Flake model are the two most commonly applied lake models in the TP [42,43]. Both models use a fixed lake depth. Returning to this study, long-term lake–air interactions will be further studied in the coupled atmospheric model, so the CLM4-LISSS lake model with a fixed lake depth was employed for the sake of consistency. (2) The NL lake level varies by less than 1 m per year, and varied by less than 3 m from 1985 to 2014 [44]. Such variations only induced small effects on the simulated LSWT [7].

### 2.3.3. Model Parameters for the Saline HSP

The salinity in HSP was set to 220‰ (Table 1) according to observations [25], and the lake depth was set to 1 m. Salinity variations were not considered because there was always insoluble salt at the bottom of HSP and the lake was shallow and well mixed.

Due to its shallowness, high salinity (corresponding Tf around −20 ◦C), and turbidity, HSP rarely freezes and has a higher albedo than the ice-free NL. The albedo of HSP was set to 0.15 (Table 1) as its annual mean albedo, as shown by MODIS.

The parameters β and η0 were set to 0.6–0.8 and three times the freshwater lake value η0 for the shallow turbid hypereutrophic Taihu Lake, respectively [45]. HSP has less phytoplankton and more transparent water than Taihu; as such, β was set to 0.6 and η0 was set to twice the freshwater value (Table 1).

### 2.3.4. Numerical Experiments Design

The bias-corrected ITPCAS data and the two lake configurations from Table 1 were used to simulate the temperature conditions in NL and HSP in two model runs henceforth referred to as S-NL and S-HSP, respectively (Table 1).

Additionally, to segregate the effects of lake depth and salinity on the lake heat budget, a sensitivity model run called S-D1F (Table 1) was performed for a hypothetical freshwater lake with a depth of 1 m under the same ITPCAS atmospheric forcing. In this way, the only difference between the S-D1F and S-NL configurations was the lake depth, and the differences between S-HSP and S-D1F runs were caused solely by the salinity.

To understand the effects of climate change on lake warming, we performed sensitivity experiments in which the monotonic trend in each meteorological forcing variable (air temperature Ta, wind speed WS, specific humidity Q, downward shortwave radiation SWD and downward longwave radiation LWD) was removed individually based on a linear regression analysis and control runs using S-NL, S-D1F and S-HSP for NL, D1F and HSP, respectively. Experiments were called S-(Lake)-d(meteorological variable), as shown in Table 2. Owing to the consistency and significant impacts on the simulated lake temperature of Ta and LWD, more experiments with detrended Ta and LWD together were run. Some of the above forcing meteorological variables could be interconnected, and the above sensitivity experiments were quite artificial in nature. However, these experiments shed light on the controlling factors of lake warming [15,19] and quantified their individual effects on lake warming rate.


**Table 2.** S-(Lake)-d(meteorological variable) sensitivity experiments.

To estimate the individual effect of salinity on lake heat budget, we ran an additional series of sensitivity experiments without considering the effects of salinity on each parameter ( α, η0, β, Rsvp, Tmaxd and Tf), i.e., S-HSP-(parameter) in Table 3. The salinity effects on simulated lake temperature caused by the specific heat capacity, and the thermal conductivity in saline lakes were not studied as they were negligible [17,37].


**Table 3.** S-HSP-(parameter of salinity effect) sensitivity experiments.

### 2.3.5. Model Performance Criteria

The performance of the model was tested against the observed temperature and heat fluxes using three model efficiency scores: bias, root mean square error (RMSE) and the correlation coefficient (R) [46]:

$$\text{Bias} = \sum\_{i=1}^{n} (\mathbf{S}\_i - \mathbf{O}\_i) / n \tag{8}$$

$$\text{RMSE} = \sqrt{\sum\_{i=1}^{n} (\mathcal{S}\_i - \mathcal{O}\_i)^2 / n} \tag{9}$$

$$\mathcal{R} = \sum\_{i=1}^{n} \left( \left( \mathcal{S}\_{i} - \overline{\mathcal{S}} \right) \left( O\_{i} - \overline{\mathcal{D}} \right) \right) / \left( \sqrt{\sum\_{i=1}^{n} \left( \mathcal{S}\_{i} - \overline{\mathcal{S}} \right)^{2}} \sqrt{\sum\_{i=1}^{n} \left( O\_{i} - \overline{\mathcal{D}} \right)^{2}} \right) \tag{10}$$

where *Oi* represents the observations, *n* is the total number of observations, and *Si* represents the simulated results.

### **3. Results and Analysis**

### *3.1. Performance of the Lake Model*

3.1.1. Performance on the Freshwater Lake (NL)

The model showed long-term seasonal variations with maximum and minimum values consistent with the MODIS LSWT (Figure 3a). The simulated LSWT forced by ITPCAS was slightly overestimated compared to the MODIS data with bias = 0.6 ◦C, RMSE = 3.2 ◦C, and R = 0.94. The simulated LSWT precision was close to that of Qinghai Lake and Nam Co Lake [4,15,16] using long-term ITPCAS. Except for the model errors, the simulated errors were from two sources: (1) compared to in situ observations, MODIS had the bias averaged from −1.5 ◦C to 0.2 ◦C and RMSE of around 2.0 ◦C owing to the cool skin effect [14,47,48]; and (2) ITPCAS was closer to the observations than other reanalysis data (e.g., NCEP and ERA) but still not very accurate, especially in lake basins with strong underlying heterogeneity. When the model was driven by observations from the NL lakeshore, the bias and RMSE between the simulated and observed LSWT were only −0.21 ◦C and 1.44 ◦C, respectively, in 2012 [49]. The model produced a better simulation with in situ observed forcing, but still could not reproduce long-term lake thermal conditions which were consistent with the assimilated meteorological dataset.

**Figure 3.** Observed and simulated LSWT (**a**), LT (**b**) and ice thickness (**c**) in S-NL, and LSWT (**d**) and LT (**e**) in S-HSP and S-D1F.

The model was able to accurately reproduce the lake temperature (Figure 3b), i.e., compared to the shallow-layer observations in the ice-free period with bias = 0.5 ◦C, RMSE = 2.8 ◦C, and R = 0.91. Similarly, the bias and RMSE between the simulated and observed lake temperature could be reduced to −0.25 ◦C and 0.41 ◦C, respectively, in 2016, when the model was driven by observations taken on the shore of NL [49].

The observed ice thickness remained above 0.6 m from the middle of January to early March, with the maximum value measured in late February (Figure 3c). With bias < 0.1 m and RMSE < 0.2 m, the model had a good ability to simulate the ice thickness.

### 3.1.2. Performance over a Saline Lake

The S-HSP numerical experiment with representative fixed salinity yielded more accurate LSWT simulations (Figure 3d) than the reference S-D1F experiment without the salinity effects, although both simulations were able to mimic the variations of LSWT. Compared to MODIS LSWT, the bias and RMSE of LSWT in S-D1F were 8.3 ◦C and 10.3 ◦C, respectively, while in S-HSP, they were reduced 3.4 ◦C and 7.0 ◦C, mainly because the application of a lower freezing point led to improved simulations in winter. The bias and RMSE from S-D1F to S-HSP in winter were reduced from 7.4 ◦C to 3.0 ◦C and from 9.3 ◦C to 3.1 ◦C, respectively. S-HSP could effectively reproduce the observed drop and increase in lake water temperature in winter, while S-D1F resulted in an ice cover and a nearly fixed water temperature below the ice. Factoring in salinity effects, the lake model was able to reflect the unfrozen state of the saline lake and the real variability of the lake temperature in winter, which is essential for studying the physics and chemistry of cold saline lakes.

### *3.2. Lake Temperature Variations and the Influence of Forcing Data* 3.2.1. LSWT and Variation Trends Annual LSWT

The annual mean LSWT of NL varied from −2 to 3 ◦C during 1979–2016, with a long-term average of 0.37 ◦C; the long-term mean LSWT of HSP was 2.95 ◦C, varying within 1–6 ◦C from one year to another (Figure 4a,c). The 2.6 ◦C lower mean LSWT in freshwater NL than in non-freezing saline HSP was conditioned by the ice cover reflecting solar radiation in winter. In the reference simulation, i.e., D1F with 1 m freshwater lake depth, the annual mean LSWT of 0.38 ◦C was close to that of NL and much smaller than that of HSP (Figure 4). Compared with salinity, lake depth did not have significant effects on the annual mean LSWT in the three simulations because the ice formation played the dominant role.

**Figure 4.** Annual LSWT (red), MLT (black) and BLT (blue) and their trends in S-NL (**a**), S-HSP (**b**) and S-D1F (**c**) Monthly LSWT.

The trends of annual LSWT, MLT and BLT in each of the three simulations except NL BLT surpassed 99% significance. The annual LSWT in NL increased by 0.68 ◦C/decade during 1979–2016, which was slightly faster than the 0.64 ◦C/decade simulated for HSP. The increasing rate of annual LSWT in S-D1F was 0.62 ◦C/decade, the lowest among the three simulations. Still, the difference between the three experiments was not significant as

long as the meteorological forcing was the same and remained the main driver responsible for the annual LSWT changing rate.

In monthly means, LSWT differences between the three experiments were mainly controlled by salinity during the cold periods, while lake depth was the primary factor during open-water periods. From November to May, the monthly LSWT of NL and D1F were similar to each other and about 4 ◦C lower than those in HSP (Figure 5a), because the high salinity of the latter prevented the development of ice cover. In turn, during the warm period from June to October, the LSWT of D1F and HSP closely followed air temperature variations; both lake surfaces were 2 ◦C warmer than that of the deeper NL in mid-summer (June and July), and both cooled faster than NL in autumn, with LSWT differences increasing gradually from 1 to 3 ◦C from September to October.

**Figure 5.** The monthly mean temperatures (◦C, lines) and their long-term trends (◦C/year, bars) in the model runs S-NL (red), S-HSP (blue) and S-D1F (black). Panels (**<sup>a</sup>**–**<sup>c</sup>**) correspond to LSWT, MLT and BLT, respectively. Solid and hollow points at the end of bars mean passing the significance test of *p* < 0.01 and *p* < 0.05, respectively.

The long-term LSWT trends of the three lakes were positive for all months. In nonfreezing HSP, the temperatures between January and April increased at 0.6 ◦C/decade (Figure 5a), with the largest increase taking place in February (1.0 ◦C/decade) in response to climate change. The LSWT of NL and D1F increased significantly in May and November because of the shortened ice period and ice-albedo feedback. In May, the LSWT in D1F warmed a little faster than in NL (1.2 ◦C/decade vs. 1.1 ◦C/decade), and in November, the

warming rate of NL LSWT (1.6 ◦C/decade) was the most acute, demonstrating that LSWT increase in months of delayed ice cover is faster in deeper lakes.

### 3.2.2. MLT and the Variation Trends
