Thermal Responses of the Largest Freshwater Lake in the Tibetan Plateau and Its Nearby Saline Lake to Climate Change

Abstract

There are thousands of lakes in the Tibetan Plateau (TP), and most are saline. However, little is known about the responses of TP lakes to climate change, especially saline ones. We investigated the thermal responses of the largest freshwater lake (Ngoring Lake) in the TP and its nearby small saline lake (Hajiang Salt Pond) to climate change using the improved lake scheme in the Community Land model (CLM4-LISSS), in which we primarily developed the salinity parameterizations previously evaluated in the Great Salt Lake in USA and further considered the effect of salinity on the temperature of the maximum density of saline water in the present study. The improved lake model with salinity parameterizations was first applied to a saline lake in the TP, where saline lakes make up the majority of water bodies. The CLM4-LISSS model could effectively simulate lake surface water temperature (LSWT), lake water temperature (LT) and ice thickness in Ngoring Lake. Additionally, the model including our salinity parameterizations significantly improved simulations of LSWT and LT in Hajiang Salt Pond, especially in winter. The LSWT of the two completely opposite lakes were warming in the simulations at a rate above 0.6 °C/decade. Meteorological forces were the main driving factor, with increasing downward longwave radiation, air temperature and air humidity, as well as weakening winds contributing to LSWT increase. Compared to a hypothetical shallow freshwater lake, the greater depth of Ngoring Lake made its surface warm faster, and salinity slightly accelerated the warming of Hajiang Salt Pond. Monthly mean LSWT differences between the two lakes were induced by salinity effects in cold periods and lake depth in the unfrozen period. In response to a warming climate, the LSWT in the ice-free Hajiang Salt Pond rapidly increased from January to April due to the warming climate, whereas the LSWT of Ngoring Lake increased faster in the first and last month of the ice-cover period due to later ice-on and earlier ice-off. This study will provide a useful tool for saline lakes in the TP and help deepen our knowledge about the responses of TP lakes, especially the saline lakes, to climate change, as well as response differences between freshwater and saline lakes and the reasons for these differences.

1. Introduction

The Tibetan Plateau (abbreviated as TP, see Abbreviation List at the end of the text) is known as “the Third Pole”, with an average altitude of 4000 m above sea level. The TP, with a total area of about 50,000 km², contains 1424 lakes (≥1 km² each) [1], most of which are saline, accounting for more than half of the total lake coverage in China. The specific TP climatic environment (low air density, pressure, and temperature, all-year intensive solar radiation) creates unique lake–atmosphere interactions [2,3,4,5,6,7,8]. TP lakes significantly influence the local and regional climate by heat and mass exchanges between lakes and the atmosphere, and resonate with the adjacent and remote regions [9,10].

Lakes are sentinels of large-scale climate variability which interact strongly with the atmosphere and respond fast and widely to climate change, especially in the TP. The TP is influenced by elevation-dependent intensive warming, at up to three times the global warming rate [8,11]. The lake surface water temperature (LSWT) has rapidly increased globally, with a mean increasing trend of 0.34 °C/decade in summer averaged over 235 lakes worldwide between 1985 and 2009 [12]. However, TP lakes have shown an overall warming trend of 0.37 °C/decade, based on data from 374 inland lakes [13]. This rate was slightly higher than the global mean, because the TP climate and the warming of TP lakes are highly heterogeneous [14]. The majority of the TP lakes are warming at a higher rate of 0.76 °C/decade primarily due to the increasing air temperature, downward longwave radiation, and decreasing wind speeds, while some lakes are cooling due to glacier meltwater inflow or reduced salinity [3,13,14,15,16,17]. Changes in thermal conditions profoundly influence a lake’s biological and chemical processes [18,19,20]. These processes may undergo substantial alterations, even with relatively small changes in lake temperature [21]. Moreover, the changing thermal characteristics of lakes further modulate local air–lake interactions, with significant impacts on the local climate. Therefore, a comprehensive investigation of the response of the thermal structure in TP lakes to climate change is needed to predict changes in lake ecosystems and the regional climate.

Most previous studies about the thermal responses of TP lakes were based on remote sensing data, which only reflected LSWT changes and the correlation between LSWT and possible driving factors. Additionally, results were mainly derived from statistical methods. However, this approach does not reveal the changes of internal phenomena in lakes, the quantitative contribution of driving factors, and the detailed mechanisms in lake processes. Numerical simulations appear to be the efficient method to reveal these key processes. A series of lake models, such as the Lake model, Flake (Freshwater lake) model, WRF (Weather Research and Forecast)-Lake, CLM (Community Land model)-Lake, CLM4-LISSS (the Lake, Ice, Snow, and Sediment Simulator), and the General Lake Model have been applied to studies of the TP lakes [2,4,7,15,16], with results showing that the vertically integrated mean lake water temperature (MLT) has been consistently changing corresponding to the increasing LSWT, while the bottom lake temperature (BLT) has varied in different ways depending on the lake depth [16,18]. However, with a scarcity of data, the development of TP lake models and numerical studies about the thermal responses of lakes to climate change have been focused solely on several large lakes, such as Nam Co Lake, Qinghai Lake, Ngoring Lake (NL) and Gyraing Lake [3,4,15], etc. In the absence of sufficient observational data and accurate forcing datasets, a previous long-term NL study employed the NCEP (National Centers for Environmental Prediction) and ERA-Interim (European Centre for Medium-Range Weather Forecasts Re-Analysis) data, in which the solar radiation was too great and decreasing quickly compared to the observations [22]. This resulted in predictions of insignificant NL warming in the simulations, whereas NL was actually warming, as shown by remote sensing data [14,23]. As such, the response of NL to climate warming should be restudied based on more accurate forcing data.

Most TP lakes are saline, but their responses to climate warming have been less studied than those of large lakes in the TP, mostly because of their small areas and the scarcity of observational data. Studies of saline lakes should be strengthened. The Hajiang Salt Pond (HSP, rich in soluble salts) is a paleo-saline lake which was formed by the joint action of traceability development of the Yellow River and climate change due to evaporation and condensation resulting from strong wind and sun. The pond is only about 11.2 km from NL and provides an ideal contrast as a saline lake to a freshwater lake NL with a similar climate.

The effects of salinity on responses of saline lakes to climate warming were poorly understood because of the lack of salinity parameterizations in commonly used models. Salinity could affect lake temperature, evaporation and ice appearance, etc. As such, a lake model considering salinity parameterizations is necessary for numerical studies of the majority of TP lakes.

Therefore, the CLM4-LISSS lake model, parameterized with salinity effects on several lake water characteristics and developed by ourselves and applied to the Great Salt Lake in USA, was introduced. Our previous saline lake model ignored the temperature of the maximum density (Tmaxd) of saline water that decreases with increased salinity and could affect the vertical thermal structure during the cold season in the TP [2,24]. To further improve our lake model, Tmaxd was further parameterized.

In the present study, we applied the CLM4-LISSS lake model developed with salinity parameterizations, in situ lake data, remote sensing data, and an assimilated meteorological dataset to study the thermal response of the largest freshwater lake NL in the TP and its nearby saline lake to climate change. The aim of the study is:

• to further improve our developed lake model with salinity parameterizations, to apply it to a TP saline lake, and to evaluate the model performance in the TP;
• to study the long-term warming trends in a deep freshwater lake and a saline lake in the TP; and
• to quantify the contribution of meteorological factors and salinity effects on the thermal regime changes of a deep freshwater lake and a saline lake.

2. Study Area, Data and Methods

2.1. Study Area

2.1.1. Freshwater Ngoring Lake

Ngoring Lake (NL), with a surface area of 610 km² and mean depth of 17 m, is the largest freshwater lake in the TP (Figure 1, 97.5~97.92°E, 34.75~35.08°N, 4274 m a.s.l.). Mineralization is low. A cold, semi-arid continental climate prevails in the NL basin. The monthly mean air temperature varies from 11.6 °C (August 2016) to −26.6 °C (January 1978), the annual average air temperature is −3.5 °C (1953–2016), and the annual precipitation is 322.4 mm (Data from China’s National Climate Center) at Maduo meteorological station (Figure 1, 34.91°N, 98.22°E, 4272 m a.s.l.). The lake is covered with ice from early December to early April.

2.1.2. Hajiang Salt Pond

Hajiang Salt Pond (HSP, 97.88–97.92°E, 35.02–35.05°N) is a small and shallow saline lake with less than 1 m depth and about 220 g L⁻¹ salinity [25]. The freezing point caused by the salinity is low enough to prevent the lake from freezing normally. It is located approximately 11.2 km east from NL at an altitude of 4240 m. The lake developed from the large Hajiang paleo-lake and currently covers an area of about 10 km².

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.

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,31,32,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]:

Tf = −0.0575 s

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

Rsvp = exp (−2/55.51 × (s/(1 − s/1000)/58.44 + 0.77))

(2)

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

c_psw = 4.188 − 4.4 s/1000

(3)

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

λ_sw = λ_fw [1.0 − 0.22 s/1000 + 0.1(s/1000)²]

(4)

where λ_fw is the freshwater thermal conductivity, equal to 0.6 W m⁻¹·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]:

Tmaxd = 3.98 − 0.216 s

(5)

2.3.2. Model Parameters for the Freshwater NL

In the model settings for NL (

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

	Experiment	S-NL	S-D1F	S-HSP
Parameter
Lake		NL	D1F	HSP
Salinity (‰)		0	0	220
Lake depth (m)		17	1	1
Albedo		Equation (7)	Equation (7)	0.15
Parameter of the light extinction coefficient η0		1.1925	1.1925	2 × 1.1925
Fraction of absorbed surface solar radiation β		0.5	0.5	0.6
Meteorological forcing		Bias-corrected ITPCAS	Bias-corrected ITPCAS	Bias-corrected ITPCAS

), 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].

The light extinction coefficient η was modelled as

η = η₀ d^−0.424

(6)

where d is the lake depth (m). The η₀ 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]:

α = α_max − α_max x + α_min x, x = exp (−95.6 (Tf − LSWT)/Tf)

(7)

where α_max and α_min are the max and min values of the lake ice albedo, respectively. In CLM4-LISSS, there are different (α_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 yet 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 η₀ were set to 0.6–0.8 and three times the freshwater lake value η₀ 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 η₀ 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. S-(Lake)-d(meteorological variable) sensitivity experiments.

Experiment	Lake	Forcing
S-NL-dTa	NL	Same as S-NL except that Ta was detrended
S-NL-dWS	NL	Same as S-NL except that WS was detrended
S-NL-dQ	NL	Same as S-NL except that Q was detrended
S-NL-dSWD	NL	Same as S-NL except that SWD was detrended
S-NL-dLWD	NL	Same as S-NL except that LWD was detrended
S-NL-dTa&LWD	NL	Same as S-NL except that Ta and LWD were detrended
S-D1F-dTa	D1F	Same as S-D1F except that Ta was detrended
S-D1F-dWS	D1F	Same as S-D1F except that WS was detrended
S-D1F-dQ	D1F	Same as S-D1F except that Q was detrended
S-D1F-dSWD	D1F	Same as S-D1F except that SWD was detrended
S-D1F-dLWD	D1F	Same as S-D1F except that LWD was detrended
S-D1F-dTa&LWD	D1F	Same as S-D1F except that Ta and LWD were detrended
S-HSP-dTa	HSP	Same as S-HSP except that Ta was detrended
S-HSP-dWS	HSP	Same as S-HSP except that WS was detrended
S-HSP-dQ	HSP	Same as S-HSP except that Q was detrended
S-HSP-dSWD	HSP	Same as S-HSP except that SWD was detrended
S-HSP-dLWD	HSP	Same as S-HSP except that LWD was detrended
S-HSP-dTa&LWD	HSP	Same as S-HSP except that Ta and LWD were detrended

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

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 (α, η₀, β, Rsvp, Tmaxd and Tf), i.e., S-HSP-(parameter) in

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

Experiment	Model Setting
S-HSP-α	Same as S-HSP except that albedo was set to that in S-NL
S-HSP-η0	Same as S-HSP except that the parameter of the light extinction coefficient η0 was set from 2 × 1.1925 to 1.1925 as in S-NL
S-HSP-β	Same as S-HSP except that the fraction of absorbed surface solar radiation β was set from 0.6 to 0.5 as in S-NL
S-HSP-Rsvp	Same as S-HSP except that the ratio of the saturated vapor pressure Rsvp over the saline water to that over the fresh water in Equation (2) was set to 1 as in S-NL
S-HSP-Tmaxd	Same as S-HSP except that the temperature of the maximum density Tmaxd was changed from minus to 3.98 °C as in S-NL
S-HSP-Tf	Same as S-HSP except that the freezing temperature Tf was changed from −12.65 °C to 0 °C as in S-NL

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

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

$$\mathrm{Bias}={\displaystyle \sum }_{i=1}^n(S_i-O_i)/n$$

(8)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \sum }_{i=1}^n{\left(S_i-O_i\right)}^2/n}$$

(9)

$$\mathrm{R}={\displaystyle \sum }_{i=1}^n((S_i-\overline{S})\left(O_i-\overline{O}\right))/\left(\sqrt{{\displaystyle \sum }_{i=1}^n{(S_i-\overline{S})}^2}\sqrt{{\displaystyle \sum }_{i=1}^n{(O_i-\overline{O})}^2} \right)$$

(10)

where $O$_i represents the observations, n is the total number of observations, and S_i 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.

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.

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.

The long-term LSWT trends of the three lakes were positive for all months. In non-freezing 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

Annual MLT

The annual vertically averaged lake water temperature (MLT) was 5.11, 3.55, and 3.35 °C in NL, D1F, and HSP respectively (Figure 4b). The annual MLT of HSP was the lowest because of the stronger heat loss from the open surface of HSP in winter compared to the ice-covered NL and D1F.

The annual NL MLT increased by 0.19 °C/decade, i.e., much slower than the annual LSWT in NL (0.68 °C/decade) and slower than the annual MLT trends in D1F (0.54 °C/decade) and HSP (0.56 °C/decade) (Figure 4). Apparently, the MLT of the deep lake increased more slowly due to its large heat capacity and thermal inertia. Salinity differences between shallow lakes D1F and HSP did not have a strong influence on the observed MLT trend.

Monthly MLT

The monthly means of HSP MLT in winter were the lowest, reaching around −10 °C in January (Figure 5b). In turn, because of the shallow depth, the MLT of HSP and D1F were the highest in summer, reaching about 14 °C in July. Hence, lake depth was the main factor controlling the magnitude of the annual and summer monthly MLT in the three simulations. In contrast, salinity determined the winter MLT minima by preventing ice cover formation in saline lakes.

In NL, the long-term trends of monthly MLT were less than 0.6 °C/decade, and variations therein were much smaller than those of monthly mean LSWT (0.2–1.6 °C/decade) (Figure 5). In NL, the monthly MLT showed a strong positive trend from April to November and stayed almost stable during the ice period, while the monthly LSWT became higher throughout the whole year, especially in freezing and breakup months.

In shallow lakes (HSP and D1F), the trends of monthly MLT in summer coincided with monthly LSWT trends (0.4–0.8 °C/decade); the difference (<0.3 °C/decade) between MLT and LSWT in shallow lakes was only seen in winter and at times of ice-formation and melting. Still, their differences were much smaller than in NL.

Akin to the surface temperatures, the MLT of HSP in winter (February to April) increased at the highest rate among the three simulations, i.e., at about 0.6 °C/decade, compared to around 0.2 °C/decade in D1F and <0.1 °C/decade in NL.

3.2.3. Bottom Temperature and the Variation Trends

Annual BLT

An annual BLT of 3.63 °C in HSP was the lowest among the three experiments, with 4.36 °C and 5.33 °C in NL and D1F, respectively (Figure 4). The simulated BLT of the two shallow lake experiments, i.e., in S-HSP and D1F, were warming by 0.52 and 0.44 °C/decade, while that of NL tended to become slightly cooler (−0.03 °C/decade) without passing the significance test. This is consistent with results from other deep dimictic lakes [16,50,51].

Monthly BLT

In NL, the monthly BLT became obviously lower (about −0.4 °C/decade) from July to September with a 0.01 significance level (Figure 5c). Its warming surface intensified the stability in the stratification period of the deep lake and resulted in less heat being transferred to the bottom. Changes in NL BLT in other months were not insignificant.

Owing to the shallow depth of D1F and HSP, there were no big differences between monthly LSWT, MLT, and BLT in the two runs during ice-free periods (Figure 5). The shallow lake depth and high salinity made the winter BLT increase faster in HSP than in other runs.

3.3. Effects of Local Climate Drivers on the Lake Warming

3.3.1. Correlation between LSWT and Meteorological Forcing

The study region was experiencing the same rapid climate change as most of the TP during the study period [52]. All trends of ITPCAS meteorological variables (Figure 6) over the study region passed the 0.01 significance level except for solar radiation, with a 0.1 significance level. Ta increased by 0.49 °C/decade, downward longwave radiation LWD increased at a rate of 3.40 W m⁻²/decade, and the specific humidity Q grew at 0.16 g kg⁻¹/decade.

The increasing Ta, LWD, and Q had positive effects on lake warming, with around 0.8 correlation with LSWT in the three experiments (Figure 7a). Wind speed decreased by −0.23 m s⁻¹/decade and was negatively correlated with the LSWT (r = −0.6). The downward shortwave radiation SWD decreased at −1.51 W m⁻²/decade, acting in the opposite direction, with a correlation coefficient of −0.15.

3.3.2. Quantified Contribution of Individual Meteorological Forcing to Lake Warming

The above correlation analysis showed the effects of meteorological variables on lake warming, but their quantified contributions were still unknown. Therefore, detrended sensitivity experiments, referred to as S-(Lake)-d (variable) in Table 2, were performed to attempt to answer this question, although the sensitivity experiments were a little artificial. The increase of atmospheric longwave radiation contributed the most to the warming of TP lakes, causing 30–40% of the annual LSWT change (Figure 7b). Increasing Ta induced a 30% increase of LSWT. The combined increase of atmospheric longwave radiation and Ta could explain almost the all of the observed LSWT warming. The decreasing WS caused a 20–30% increase in LSWT, especially in HSP. Air humidity increase accelerated the lake surface warming by 20% and, consequently, should not be ignored, as it was in a study of the Nam Co Lake in TP [15]. The decrease of SWD decelerated warming by about 10%.

3.4. Effects of Salinity Parameters on the Lake Warming

The cumulative effect of salinity on the lake water properties caused the simulated lake surface to be 2.6 °C warmer and experience 0.02 °C/decade faster warming in S-HSP than in its freshwater counterpart (Figure 4).

The salinity effect on Tf and Rsvp were simulated, making the saline lake surface 2.34 °C and 0.9 °C warmer (Figure 7c) in S-HSP than in the two sensitivity experiments of salinity effects (S-HSP-Tf and S-HSP-Rsvp), inducing about 90% and 31% annual differences between the saline lake and the freshwater one. The lower freezing point also significantly accelerated the increase of simulated LSWT in saline lake warming to 0.07 °C/decade (3.5 times the warming rate difference between the saline lake and the hypothetical freshwater lake with the same 1 m depth), while the salinity effect on the saturation water vapor pressure had no impact on the long-term temperature trend.

The differences between the S-HSP and S-HSP-η₀ simulations showed that the lower transparency of salt water with higher η₀ increased the annual LSWT by only 0.02 °C and slowed the warming by 0.01 °C/decade.

The difference between the S-HSP and S-HSP-α simulations showed that the higher α of salt water cooled the lake surface by 1.24 °C annually (about −47% of the annual LSWT difference between the saline lake and the freshwater one) and slowed the warming rate by 0.02 °C/decade (the same magnitude as the warming rate caused by all the salinity effects).

Changes in the temperature of the maximum density (Tmaxd) of saline water and the absorption of solar radiation (β) by the lake surface caused by salinity did not affect the annual mean temperature or the changing rate.

4. Discussion

4.1. Salinity Effects and Parameterizations

Most TP lakes are saline, but existing numerical studies have focused exclusively on several large lakes with small amounts of salt because of the lack of observations and salinity parameterizations in lake models. Based on previous salinity parameterizations coupled in CLM applied to the Great Salt Lake in USA and the significant impacts of Tmaxd on density convection and thermal stratification, the effects of Tmaxd were further parameterized in our lake model. The improved model was first applied to the TP saline lake and significantly reduced errors in the simulation of LSWT and LT in the saline lake, especially in winter. The salinity-extended lake model will be an efficient tool for studying saline lakes in the TP.

Salinity parameterization of Tmaxd had no obvious impacts on the warming of a small saline lake, mainly because of the shallow lake depth, absence of ice cover and the strong effects of wind in terms of turbulence mixing in HSP. However, Tmaxd could alter the thermal structure of a lake with a certain depth [2,24]. Therefore, considering the salinity effect on Tmaxd makes the developed saline lake model more accurate.

Salinity will play a major role in terms of the impact of climate change on TP lakes in future. Salinity decreases due to increased precipitation and inflow of glacial meltwater, but it could also increase due to increased evaporation. Since many lakes are fairly small, even small changes in the water balance can be important. Also, salinity evolution influences the vertical stratification of lakes and, consequently, the water temperature structure. Variations in salinity will be addressed in future according to the mass balance.

4.2. Simulated Warming Rates of LSWT in Different Studies

LSWTs globally have increased rapidly, with a mean trend of 0.34 °C/decade in summer between 1985 and 2009 [12]. Although the TP is warming at twice or even three times the global warming rate [8,11], TP lakes have overall been warming with a trend of 0.37 °C/decade, based on data from 374 inland lakes. They are warming slightly more rapidly than the global mean because the warming of TP lakes is highly heterogeneous [14]. Most TP lakes are warming with the higher 0.76 °C/decade rate, while some lakes are cooling due to glacier meltwater inflow or reduced salinity [3,13,14,15,16,17].

The simulated LSWT warming rate of NL was 0.68 °C/decade, as shown by remote sensing data that the lake was warming [14,23], and its warming rate was between the simulated rates of 0.52 ± 0.25 °C/decade in Nam Co by the GLM and 0.74 °C/decade in Qinghai Lake by Flake. The trend was not insignificant as in the previous NL simulation, in which the model was forced by the NCEP and ERA data where solar radiation was larger and decreased quickly compared to the observations [22]. Thus, insignificant warming in NL was concluded in the previous study. Although field work is hard, enough observations should be performed and accurate forcing datasets should be built for more accurate simulation studies.

Our results make it clear that the largest freshwater lake in the TP and a nearby small saline lake have indeed warmed over the last several decades and are warming faster owing to the amplification effect of their high altitude. The warming rate of LSWT in NL (0.68 °C/decade) significantly exceeded that of the regional air temperature (0.49 °C/decade) and was similar to Qinghai Lake and Lake Superior due to reduced ice cover [16,53]. While lakes in some temperate climate regions are warming in line with increased air temperatures [50,54], even the warming rates of tropical lakes are smaller than those of air [55,56]. Moreover, the bottom layers of NL were simulated to isolate from direct atmospheric influence, and tended to show long-term cooling at a rate −0.03 °C/decade on account of strengthening stratification. This result is similar to those reported from other stratified dimictic lakes, such as Qinghai Lake, Heiligensee Lake, and so on [16,50].

5. Conclusions

Compared to our observations, the salinity-extended lake model demonstrated a good ability to represent lake–air interactions and the thermal regime in both a freshwater lake and a saline lake. The newly introduced salinity parameterization significantly improved the model performance for a saline lake in winter.

The simulated long-term increasing rates of the annual LSWT in NL, HSP, and a shallower hypothetical freshwater lake amounted to more than 0.6 °C/decade, mainly due to meteorological forcing. Increasing LWD and Ta, weakening wind, and increased air humidity had positive effects on the warming trend of TP lakes in decreasing order, while solar radiation dimming counteracted the warming. The LWD and Ta contributed the most to lake warming in the sensitivity experiment; although it was overlooked in previous studies [3,15], increasing atmospheric humidity over TP should be considered a significant climatic factor.

The shallow lakes in experiments S-D1F and S-HSP lacked seasonal thermal stratification and were well-mixed vertically, revealing similar long-term warming trends across their depths. Comparing to a fresh water lake with 1 m depth, the 17-m-deep NL experienced a 0.06 °C/decade faster surface warming and a slower MLT rise of 0.35 °C/decade.

High salinity prevented ice cover formation in HSP and induced more heat release in winter and lower MLT and BLT than in freshwater lakes NL and D1F. However, the high salinity made the annual mean LSWT 2.6 °C higher and resulted in a 0.02 °C/decade stronger warming trend than in the freshwater lake with the same depth. The salinity effect on the freezing point contributed most to this difference, inducing a 90% higher LSWT compared to the freshwater D1F. The salinity effect on evaporation caused a 31% higher LSWT in HSP. The opposite salinity effect on the lake surface albedo cooled the lake surface and decelerated the warming trend.

The monthly mean LSWT differences between Ngoring Lake and the Hajiang Salt Pond were induced by salinity effects in cold periods and lake depth in the unfrozen period. The LSWT in ice-free Hajiang Salt Pond increased rapidly from January to April due to climate change, whereas the LSWT of Ngoring Lake increased faster in the first and last months of the ice-cover period due to later ice-on and earlier ice-off.