Mechanisms and effects of under-ice warming water in Ngoring Lake of Qinghai–Tibet Plateau

The seasonal ice cover in lakes of the Qinghai–Tibet Plateau is a transient and vulnerable part of the cryosphere, whose characteristics depend on the regional climate: strong solar radiation in the context of the dry and cold environment because of the high altitude and relatively low latitude. We use the first under-ice temperature observations from the largest Tibetan freshwater lake, Ngoring Lake, and a one-dimensional lake model to quantify the mechanism of solar thermal accumulation under ice, which relies on the ice optical properties and weather conditions, as well as the effect of the accumulated heat on the land–atmosphere heat exchange after the ice breakup. The model was able to realistically simulate the feature of the Ngoring Lake thermal regime: the “summer-like” temperature stratification with temperatures exceeding the maximum density point of 3.98 ^∘C across the bulk of the freshwater column. A series of sensitivity experiments revealed solar radiation was the major source of under-ice warming and demonstrated that the warming phenomenon was highly sensitive to the optical properties of ice. The heat accumulated under ice contributed to the heat release from the lake to the atmosphere for 1–2 months after ice-off, increasing the upward sensible and latent surface heat fluxes on average by ∼ 50 and ∼ 80 W m⁻², respectively. Therefore, the delayed effect of heat release on the land–atmosphere interaction requires an adequate representation in regional climate modeling of the Qinghai–Tibet Plateau and other lake-rich alpine areas.

1 Introduction

Seasonal lake ice is a part of the cryosphere, gaining recent attention from researchers due to its sensitivity to climate change (Kirillin et al., 2012; Sharma et al., 2020). The duration of ice cover affects the stability and vertical mixing of lakes, as well as the lake–atmosphere matter and energy exchange (Rösner et al., 2012; Efremova et al., 2013; Ramp et al., 2015). Ice cover regulates lake biochemical indicators, such as the concentration of dissolved oxygen, nitrogen, and phosphorus, changing the biochemical reaction rate and affecting the water quality and distribution of aquatic organisms (Weitere et al., 2010; Dokulil, 2013; G. C. Li et al., 2015; Hardenbicker et al., 2016). Shortening of the ice season has been observed worldwide (Sharma et al., 2019; Dauginis and Brown, 2021) and attributed to anthropogenic warming (Grant et al., 2021). Future climate predictions indicated the accelerated reduction in seasonal lake ice, especially pronounced in the lake-rich Arctic regions (Brown and Duguay, 2011). Global assessment of seasonal lake ice changes requires quantification of the major heat sources and sinks on seasonal to climatic timescales. While the major prerequisite for the ice cover development is a sufficiently long season with air temperature below the freezing point of water, the heat budget of ice-covered lakes varies with latitude and altitude, depending strongly on the available solar radiation, the latter being the major source of heat for under-ice lake water (Kirillin et al., 2012). During the polar night in the Arctic and temperate lakes covered by snow, the solar heating is minor, and the bottom sediment is the main heat source (Winter I according to Kirillin et al., 2012); at later stages of the ice season (Winter II), as the snow melts, solar radiation becomes the main heat source governing thermal stratification and mixing under ice and the melting process at the ice base (Kirillin et al., 2018, 2020). Further, lakes with seasonal ice cover can be divided into cryomictic and cryostratified according to their maximum depth, surface area, and wind speed (Yang et al., 2021). In dry and cold areas with little snow, Winter II can occupy the entire ice-covered period (Kirillin et al., 2012), making solar radiation the major factor affecting the lake ice regime. The situation is relevant to the alpine lakes.

In particular, the largest alpine lake system of the Qinghai–Tibet Plateau (TP) is located not only at the highest plateau on Earth, with an average altitude of 4000–5000 m, but also at the relatively low latitude of 26–39^∘ N in the TP, ensuring a high amount of solar radiation. The TP is covered by more than 1400 lakes with areas larger than 1 km², and the total lake area is more than 5 × 10⁴ km², accounting for 57.2 % of that in China (Wan et al., 2016; Zhang et al., 2019). Recent studies reported the first observations from ice-covered Tibetan lakes, indicating the major role of solar radiation in their thermal regime (Wang et al., 2021). Water temperatures in lakes Bangong Co and Nam Co constantly increased during the ice-covered period, with a stronger increase in shallower Bangong Co (Lazhu et al., 2021). Observations in meromictic Dagze Co demonstrated that stable temperatures in the early ice-covered period start warming only in the late ice-covered period, conditioned by the high water salinity (Wang et al., 2014; Lazhu et al., 2021). Salinity has a strong influence on the temperature and mixing regime of all three abovementioned lakes, by altering their density stratification and vertical heat transport. Among freshwater lakes in the TP, Ngoring Lake is the largest one (Kirillin et al., 2017; Wen et al., 2022). Kirillin et al. (2021a) found strong solar radiation under ice cover heating the entire lake water column to the maximum freshwater density temperature (T_md; ∼ 3.98 ^∘C) more than a month before the ice breakup; this situation is never found in lowland freshwater lakes. As a result, strong heat release from water to the ice base turned into the major factor governing the ice melt, with the water temperature under ice reaching 6 ^∘C. This radiation-dominated regime differs dramatically from the typical heat budget known from earlier studies on ice-covered lakes. Although it belongs to the classification that Winter II occupies the whole ice-covered period, the under-ice lake temperature exceeds 4 ^∘C in the late stage due to the strong solar radiation (SR) on the TP. Quantification of the resulting heat balance and thermal stratification characteristics of alpine conditions is the subject of the present study.

Due to the harsh environment of the TP and difficulties in collecting field observations, numerical models are often used to investigate phenomena and mechanisms of TP lakes. At present, the widely used lake models are the FLake model and the lake scheme coupled in the CLM (Community Land Model), CoLM (Common Land Model), and WRF (Weather Research and Forecasting Model) (Lazhu et al., 2016; Wen et al., 2016; Fang et al., 2017; Dai et al., 2018; Huang et al., 2019; Song et al., 2020; Wu et al., 2021). However, for computational efficiency, winter dynamics in these highly parameterized lake models are represented in a rather simplified way, lacking the detailed mechanisms of heating by radiation and resulting vertical heat transports across the water column (Lazhu et al., 2016; Wen et al., 2016; Huang et al., 2019). As an alternative, we adopt for this study a “classical” two-equation turbulence modeling approach proving its reliability in decades of studies on the environmental turbulent fluid dynamics. The one-dimensional model LAKE implements the approach in application to lake dynamics and was applied previously to different lakes (Stepanenko et al., 2011, 2016; Guseva et al., 2016). We combine modeling with in situ observations from Ngoring Lake, data on weather forcing, and remote sensing to (i) test the ability of a one-dimensional lake model LAKE to simulate temperature and stratification driven by intense solar heating in ice-covered Lake Ngoring, (ii) conduct series of sensitivity experiments aimed at revealing the role of meteorological forcing and ice optical properties in lake temperature and mixing regime, and (iii) reveal the effects of temperature distribution before ice breakup on lake heat storage and lake–atmosphere heat transfer.

2 Study area and data

2.1 Study area

Ngoring Lake (34.76–35.08^∘ N, 97.53–97.90^∘ E; Fig. 1) is located in the western valley of Maduo County on the eastern TP, with an average lake surface elevation of 4274 m a.s.l. It is the largest freshwater lake in the Yellow River source region, with a salinity of about 0.27 g kg⁻¹ (Shen et al., 2012). It has a surface area of 610 km², a maximum depth of 32 m, and an average depth of 17 m. The pH is 8.49, and there are very few fish in the lake. Aquatic plants grow only in the riparian area. The lake thermally stratifies in summer and is covered by ice from late November or early December to late April (Wen et al., 2016). According to observational data from 1953 to 2016 at Maduo station (34.9^∘ N, 98.2^∘ E) of the China Meteorological Administration, the average annual precipitation was 322.4 mm, mostly concentrated from May to September. The average annual air temperature was −3.53 ^∘C. The maximum air temperature was 24.3 ^∘C and occurred on 20 July 2006, and the minimum air temperature was −48.1 ^∘C and occurred on 2 January 1978.

Figure 1 (a) Location and (b) bathymetry of Ngoring Lake. Panel (b) is adapted from Kirillin et al. (2021a). The pentagrams denote the lake border station (LBS) and water temperature measurement site (WS).

2.2 Data

2.2.1 Observations: LBS station and WS site

The long-term automatic lake border station (LBS; 34.91^∘ N, 97.55^∘ E; Fig. 1) was installed in October 2012, at an altitude of 4282 m a.s.l., providing meteorological forcing data: wind speed at 10 m and air temperature, specific humidity, and air pressure at 2 m as well as downward shortwave (SR) and longwave radiation (LR) at 1.5 m from September 2015 to September 2016 (Li et al., 2020). The detailed information about site configuration and measured quantities are referred to in Z. G. Li et al. (2015) and Wen et al. (2016). The precipitation was obtained from the daily value data set (V3.0) at Maduo station of Chinese surface climate data (http://data.cma.cn, last access: 27 August 2022).

The water temperature measurement site (WS; 35.03^∘ N, 97.70^∘ E; Fig. 1) was located in the north of Ngoring Lake, where the total water depth was about 26.5 m. The multi-layer water temperature observation system consisted of 16 self-recording RBR SOLO water temperature probes with a precision of 0.01 ^∘C. The sampling distance and time intervals were 1 m and 10 min, respectively.

2.2.2 MODIS lake surface temperature

The 8 d L3 global lake surface temperature product (MYD11C2) was derived from the data of Moderate Resolution Imaging Spectroradiometer (MODIS) and was used to evaluate the simulated results. MODIS offers long-term daily global coverage data with high spatial resolution. This product provides an 8 d combined radiative surface temperature at approximately 10:30 and 22:30 LT (local time), which is the satellite transit time. The resolution is 0.05^∘ latitude–longitude (5600 m at the Equator) for the Climate Modeling Grid (CMG) (https://ladsweb.nascom.nasa.gov/search, last access: 27 August 2022) (Wan et al., 2004).

2.2.3 ERA5-Land data

ERA5-Land is produced as an enhanced global data set for the land component of the fifth generation of European Reanalysis (ERA5) by the European Centre for Medium-Range Weather Forecasts (ECMWF), framed within the Copernicus Climate Change Service (C3S) of the European Commission. It is available for the hourly record for about 40 years, from 1981 to the present. Expediently, ERA5-Land has an enhanced horizontal resolution of 9 km (∼ 0.08^∘) compared to ERA5 (31 km) and ERA-Interim (80 km) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=form, last access: 27 August 2022) (Hersbach et al., 2020; Muñoz-Sabater et al., 2021).

ERA5-Land data are applied for a comparative analysis of warming mechanisms and thermal conditions in Tibetan ice-covered lakes against those in the Arctic. The reanalysis forcing data for the geographical position 69.05^∘ N, 20.83^∘ E, were adopted as “typical” arctic weather conditions; that is the SR is 0 W m⁻² during the winter polar night period. Northern Fennoscandia is covered by several lakes characterized by the longest ice-covered period in western Europe. The largest of these lakes, Kilpisjärvi, has a similar morphometrical feature to Ngoring (average depth of 19.5 m, maximum depth of 57 m, surface area of 37 km²). However, they receive different SR because of different latitudes. The lake has been intensively studied in the last decades (Kirillin et al., 2015, 2018; Leppäranta et al., 2017, 2019). Its under-ice water temperature remained stable during winter from 1992 to 1993 (Tolonen, 1998). In the following, model experiments forced by the ERA5 weather data (1992–1993) for the Arctic refer to “Kilpisjärvi” runs.

3 Methods

3.1 LAKE model

The one-dimensional model LAKE, simulating thermodynamic, hydrodynamic, and biogeochemical processes, is used to solve the horizontally averaged transfer of gases, heat, salts, and momentum in an enclosed water body (Stepanenko et al., 2011, 2016). The vertical heat diffusion is simulated, and the penetration of solar radiation into the water ice, snow, and bottom sediment layers (Heiskanen et al., 2015; Cao et al., 2020) is taken into account. The exchange between the water and the inclined bottom is modeled explicitly because the model equations have been averaged over horizontal sections of the water body. The second-order κ–ε parametrization of turbulence is applied (Stepanenko et al., 2016).

3.1.1 Heat transfer in water body

The water temperature is calculated according to the one-dimensional thermal diffusion equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{c}_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\partial T_{\mathrm{w}}}{\partial t}}& =-c_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\mathrm{1}}{A}}{\int }_{{\mathrm{\Gamma }}_A}{T}_{\mathrm{w}}\left({\mathit{u}}_{\mathrm{h}}\cdot \mathit{n}\right)\mathrm{d}l\\ & +{\displaystyle \frac{\mathrm{1}}{A{h}^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial \mathit{\xi }}}\left(A_{\mathrm{w}}{K}_{\mathrm{T}}{\displaystyle \frac{\partial T_{\mathrm{w}}}{\partial \mathit{\xi }}}\right)-{\displaystyle \frac{\mathrm{1}}{Ah}}{\displaystyle \frac{\partial AS}{\partial \mathit{\xi }}}\\ & +{\displaystyle \frac{\mathrm{1}}{Ah}}{\displaystyle \frac{\partial A}{\partial \mathit{\xi }}}\left[S_b\left(\mathit{\xi }\right)+F_{iz,\mathrm{b}}\left(\mathit{\xi }\right)\right]+{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}{\displaystyle \frac{\mathit{\xi }}{h}}{\displaystyle \frac{\partial T_{\mathrm{w}}}{\partial \mathit{\xi }}},\end{array}\end{array}$$

where c_w is water specific heat, ρ_w is water density, T_w is water temperature, h(t) is lake depth, t is time, $\mathit{\xi }=z/h$ is a normalized vertical coordinate $\left(z\in \left[\mathrm{0},h\right]\right),z=\mathrm{0}$ is located at the free water surface of the lake, S is downward shortwave radiation, A_w is the z-dependent cross-sectional area of water, K_T is the thermal diffusivity coefficient equal to the sum of molecular and turbulent diffusivities, S_b(ξ) is shortwave radiation flux, F_iz,b is soil heat flux at the level z, n is an outer vector normal to the boundary Γ_A of the horizontal cross-section A, and u_h is a horizontal vector in water (Stepanenko et al., 2016; Guseva et al., 2016).

3.1.2 Heat transfer in ice cover

When the air temperature decreases below 0 ^∘C, and the surface water temperature drops to the freezing point, the initial ice cover forms. When the net radiation of the lake is positive, the ice melts continuously until the ice thickness declines to zero. The general heat conduction equation in ice cover follows the equation

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{c}_{\mathrm{i}}{\mathit{\rho }}_{\mathrm{i}}{\displaystyle \frac{\partial T_{\mathrm{i}}}{\partial t}}& =c_{\mathrm{i}}{\mathit{\rho }}_{\mathrm{i}}{\displaystyle \frac{\mathit{\xi }}{h_{\mathrm{i}}}}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{i}}}{\mathrm{d}t}}{\displaystyle \frac{\partial T_{\mathrm{i}}}{\partial \mathit{\xi }}}-c_{\mathrm{i}}{\mathit{\rho }}_{\mathrm{i}}{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{i}}}}{\displaystyle \frac{\mathrm{d}{h}_{i\mathrm{0}}}{\mathrm{d}t}}{\displaystyle \frac{\partial T_{\mathrm{i}}}{\partial \mathit{\xi }}}-{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{i}}}}{\displaystyle \frac{\partial S}{\partial \mathit{\xi }}}\\ & +{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{i}}{h_{\mathrm{i}}}^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial \mathit{\xi }}}\left(A_{\mathrm{i}}{\mathit{\lambda }}_{\mathrm{i}}{\displaystyle \frac{\partial T_{\mathrm{i}}}{\partial \mathit{\xi }}}\right)+{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{i}}{h}_{\mathrm{i}}}}{\displaystyle \frac{\partial A_{\mathrm{i}}}{\partial \mathit{\xi }}}{F}_{\mathrm{T},\mathrm{b}}\\ & -L{\mathit{\rho }}_{\mathrm{i}}{\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}},\end{array}\end{array}$$

where c_i is ice specific heat, ρ_i is ice density, T_i is ice temperature, λ_i is ice thermal conductivity, h_i is ice thickness, $\frac{\mathrm{d}{h}_{i\mathrm{0}}}{\mathrm{d}t}$ is the increment of ice thickness on its surface, F_T,b is the heat flux at the ice–sediment boundary, A_i is the z-dependent cross-sectional area of the ice cover determined by the basin morphometry, L is the latent heat of water, and p is ice porosity (Stepanenko et al., 2019). The last term to the right-hand side presents heat of phase transition of salty water in ice pores.

The penetration of solar radiation into the medium is calculated using the Beer–Lambert law (Stepanenko and Lykossov, 2005; Stepanenko et al., 2019):

$$\begin{array}{}\text{(3)}& S\left(\mathit{\xi }\right)=S\left(\mathrm{0}\right)\exp \left(-a_{\mathrm{e}}h\mathit{\xi }\right),\end{array}$$

where a_e is the medium extinction coefficient. To solve the temperature in Eqs. (1) and (2), it is necessary to specify the top and bottom boundary conditions and provide the method to calculate the heat flux at each depth z. The atmospheric turbulent heat flux schemes are based on the Monin–Obukhov similarity theory (Stepanenko et al., 2016).

When the lake is covered by ice, the temperatures of the bottom layer of ice and the top layer of water are equal and fixed to the melting/freezing point temperature (Stepanenko et al., 2019), which is calculated by the following formula:

$$\begin{array}{}\text{(4)}& T_{\mathrm{mp}}=-C\cdot \left|\partial T_{\mathrm{mp}}/\partial C\right|,\end{array}$$

where T_mp is the melting/freezing point temperature (^∘C), C is salinity at the water–ice interface, and $\left|\partial T_{\mathrm{mp}}/\partial C\right|=\mathrm{66.7}$ and is assumed constant.

Based on the study by Leppäranta (2015), the albedo regulates the surface energy budget, and the extinction coefficient controls the vertical distribution of radiation energy in the medium. In the LAKE model, the albedo of water (A_w) is 0.06, and the snow extinction coefficient (E_s) decreases with increased snow density. Snow accumulation in the Ngoring Lake area is almost zero. Therefore, only A_i, E_i, and E_w are analyzed in this study. Version 2.3, called LAKE2.3, is used in this article.

3.2 Methods to evaluate the model accuracy

The indexes to evaluate accuracy of the model are the root mean square error (RMSE), bias, and correlation coefficient (CC):

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{j=\mathrm{0}}^n{\left(m_j-o_j\right)}^{\mathrm{2}}},\text{(6)}& {\displaystyle }\mathrm{Bias}=\stackrel{\mathrm{‾}}{m}-\stackrel{\mathrm{‾}}{o},\text{(7)}& {\displaystyle }\mathrm{CC}={\displaystyle \frac{\mathrm{Cov}\left(M,O\right)}{\sqrt{\mathrm{Var}\left(M\right)\mathrm{Var}\left(O\right)}}},\end{array}$$

where m_j and o_j represent the simulations and observations; $\stackrel{\mathrm{‾}}{m}$ and $\stackrel{\mathrm{‾}}{o}$ are the corresponding average values; Var(M) and Var(O) are the variances of observed and simulated values, respectively; and Cov(M,O) is the covariance.

3.3 Calculation method of heat storage

The heat storage evolution in water is calculated by the following formulation:

$$\begin{array}{}\text{(8)}& Q=c_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{w}}{\sum }_{k=\mathrm{1}}^n{T}_k\mathrm{\Delta }{z}_k,\end{array}$$

where c_w=4192 J kg⁻¹ K⁻¹ and ρ_w=10³ kg m⁻³, n is the layer number, Δz_k is depth interval between two successive layers, and T_k (K) is the average temperature in layer k (Nordbo et al., 2011; Gan and Liu, 2020).

4 Characteristic analysis of water temperature and local climate in two lakes

4.1 Characteristics of water temperature

It is pointed out that the under-ice water temperature from 2015 to 2016 in Ngoring Lake rose continuously during the entire ice-covered period according to observations (Wang et al., 2021; Kirillin et al., 2021a). In November, because of the strong wind and wide surface, the lake mixed evenly with slight oscillation (<1 ^∘C between 2 and 22 m), and water temperature decreased gradually until the lowest point of 0.47 ^∘C at 2 m on 12 December, after which the lake froze up completely (Fig. 2a). Meanwhile, the 2 m air temperature fell to −7.79 ^∘C. Ngoring Lake is mostly covered only by bare ice in winter due to drought and less precipitation and snow. In the early ice-covered phase (from 12 December to 7 March), the whole lake mixed completely because solar radiation penetrated ice and heated the upper water, which was warm (<T_md), heavy, and sinking (Fig. 2b) (Kirillin et al., 2012). In parallel, water temperature continued to warm until it reached T_md on 7 March (Fig. 2a).

In the late ice-covered stage (from 7 March to 18 April), the lake stratified, on the one hand owing to strengthened solar radiation and on the other hand since radiation absorption of water decayed with depth based on the Beer–Lambert law. Water temperature increased at the rate of 0.052 ^∘C d⁻¹ in the layers from 2 to 6 m, which was more rapid than the early stage of 0.035 ^∘C d⁻¹. On 18 April, the ice melted entirely, and water temperature rose to 5.83 ^∘C at 2 m while remaining at T_md below 9 m. After that, full mixing took place rapidly because the lake warmed, gaining heat from the sun and atmosphere as a result of ice breaking up (Fig. 2b).

Figure 2 (a) The daily average water temperature observations of Ngoring Lake at the surface (Ts​​​​​​​) and at 2, 9, and 22 m from 1 November 2015 to 1 June 2016. Ts is MODIS lake surface temperature. The gray reference lines denote 3.98 and 0 ^∘C, respectively. The shaded pink area denotes the ice-covered period. The water temperature profile (b) observed and (c) simulated in CTL. The ice-covered period is represented between the two dashed red lines.

4.2 Local climate: Tibet vs. Kilpisjärvi

The daily averages of meteorological variables between Tibet and Kilpisjärvi are shown in Fig. 3, and the ranges and averages of those meteorological variables during the ice-covered period (from 12 December to 18 April based on Ngoring Lake) are compared (Table 1). The average differences in air temperature, specific humidity, and downward LR were −0.42 ^∘C, −0.38 g kg⁻¹, and 41.9 W m⁻², respectively. The wind speed of Tibet was 1.7 times that of Kilpisjärvi, and the downward SR in Tibet of 199.41 W m⁻² was stronger than that in Kilpisjärvi of 40.46 W m⁻². In Tibet, the precipitation was a multiple of 0.037 less than that in Kilpisjärvi.

On the whole, there were few differences in air temperature, specific humidity, and downward LR in the two regions. Nevertheless, there was much lower precipitation and much higher downward SR and wind speed in the high-altitude and relatively low-latitude Tibet. Surface pressure was not considered in this study due to the minor effect on water temperature.

Figure 3 Comparison of daily average values of the meteorological variables for Tibet from 2015 to 2016 and for Kilpisjärvi from 1992 to 1993. (a) Precipitation, (b) downward SR, (c) wind speed, (d) downward LR, (e) specific humidity, and (f) air temperature. “Ngoring” represents the Tibet region, and “Kilpis” represents Kilpisjärvi region.

Table 1 Ranges and averages of the meteorological variables of Tibet (2015–2016) compared with Kilpisjärvi (1992–1993) during the ice-covered period (12 December–18 April).

5 Simulation setup

To reveal the mechanism of water temperature rising during the ice-covered period in Ngoring Lake and its further influences, one control simulation (CTL) and 27 sensitivity simulations (SIMs) depending on CTL were set in this study (Table 2).

5.1 Setup in CTL

The depth was set as 26.5 m, measured at the WS point, and divided vertically into 35 layers. The simulation period was from September 2015 to September 2016. The initial vertical profile of water temperature, mixed layer, and the bottom temperature was set following the observations (Fig. 2b). The albedos of snow and ice and the extinction coefficients of ice and water were set as A_s=0.7, A_i=0.25, E_i=2.5 m⁻¹, and E_w=0.15 m⁻¹ based on previous investigations (Lei et al., 2011; Li et al., 2018, 2020; Shang et al., 2018). The inputs driving meteorological variables were air pressure, wind speed, specific humidity, air temperature, precipitation, and downward SR and LR. The forcing data and model run intervals were 30 min and 15 s respectively.

5.2 Setup in SIMs

To explore the influence of a single meteorological variable, SIM_^* simulations were set up. The symbol ^* is SR, Precip, LR, U, T_air, or q in Kilpisjärvi. These scenarios were quite artificial because these variables are closely correlated. Despite that, these sensitivity simulations can shed light on the influence of local climate on lake temperature evolution during the ice-covered period.

To discuss the effect of the main physical parameters, SIM_# simulations were set up. The sign # represents the values of A_i, E_i, or E_w. SIM_E^* (^* equal to 1, 2, or 3) is set for exploring the effects of three different initial water temperature profiles before ice breakup on the lake heat storage and heat fluxes.

Table 2 Names, explanations, and numbers of all experiments.

6 Simulation results

6.1 Model validation

The CTL simulation was relatively consistent with the observations (Fig. 2b, c), even though the whole ice season of the CTL started about 15 d earlier than the observed. Water temperature was a little higher in CTL than that in the observations from mid-March to late May. The deviation was greater in the deep water. The simulated temperature warmed faster and was higher by 1 ^∘C than the observed value after ice melted.

The results were evaluated by calculating RMSE, bias, and CC between simulated and observed water temperature at the lake surface (Ts) and at 2, 9, and 22 m (Fig. 4). The CC in each layer was equal to or even greater than 0.95, and the CC at 2 and 9 m was as high as 0.98, even though RMSE and bias of lake surface were 3.25 and 1.42 ^∘C, respectively. The Ts RMSE was largely due to the uncertainty in MODIS lake surface temperature (Donlon et al., 2002; Tavares et al., 2019). The absolute bias values in the internal lake were less than 0.01 ^∘C, and RMSE was less than 0.95 ^∘C. More importantly, the under-ice water temperature warming phenomenon was reproduced reasonably.

Figure 4 The daily average water temperature observed and simulated in CTL of (a) the surface (Ts), (b) 2 m, (c) 9 m, and (d) 22 m in Ngoring Lake from November 2015 to June 2016. The dotted lines represent T_md 3.98 ^∘C.

6.2 Influences of local climate on water temperature

To explore the influences of local climate on water temperature, six SIM_^* (^* represents six meteorological variables; Table 2) simulations were designed. The 3 m water temperature was typically selected to analyze since water temperature at different depths varied consistently over time.

SIM_SR was the simulation when the Kilpisjärvi downward SR was substituted for that in Tibet. During the ice-covered period, the downward SR difference between CTL (199.41 W m⁻²) and SIM_SR (40.46 W m⁻²) was 158.95 W m⁻². In the sensitivity simulation SIM_SR, the 3 m water temperature was stable, keeping in the range of 0–0.1 ^∘C (Fig. 5a). The ice formation date was earlier, and the ice-breaking date was delayed, which led to the growth of the whole ice season. The mixed layer depth increased (Fig. 5d). Consequently, the strong downward SR on the TP generated the under-ice water warming in Ngoring Lake.

In the simulation SIM_Precip, the Tibet precipitation was replaced by that in Kilpisjärvi. In the sensitivity experiment SIM_Precip, the 3 m water temperature was fixed, then increased but did not exceed T_md in the early ice cover stage (Fig. 5a). The stratification and temperature maximum center disappeared in late March, and the lake was fully mixed (Figs. 2c, 5g). Because the average precipitation in SIM_Precip (0.12 mm h⁻¹) was approximately 30 times larger than that in CTL (0.0044 mm h⁻¹) during the ice-covered period, more solar radiation was reflected and absorbed by snow due to more snowfall accumulation. Thus, the high precipitation damped the water temperature rise.

In the SIM_LR simulation, the downward LR in Kilpisjärvi superseded that in Tibet. The average downward LR was 233.62 W m⁻² in SIM_LR, which was stronger than that in CTL (191.73 W m⁻²) during the ice-covered period. The 3 m water temperature still warmed, and the complete ice melting time of late February was ahead. The heat was transferred from lake to atmosphere because of lower air temperature after ice breakup. The water temperature underwent a cooling process (2 ^∘C) until reaching a new equilibrium with the atmosphere (Fig. 5b). Compared with the CTL, water mixing in the ice-covered period was more uniform, the stratification in late March was weakened, and the temperature maximum center advanced by about 15 d (Fig. 5e).

In the SIM_U simulation, the wind speed of Kilpisjärvi was substituted for that of Tibet. The wind speed in SIM_U (2.83 m s⁻¹) was weaker than that in CTL (4.93 m s⁻¹) for the ice-covered period. In the sensitivity experiment SIM_U, the 3 m water temperature kept rising, but it was about 3 ^∘C higher than that in CTL during the whole simulation period (Fig. 5b). Due to the decrease in wind speed, the mixed layer depth was reduced, and the lake stratification was more stable (Fig. 5h).

In the SIM_T_air simulation, the air temperature of Kilpisjärvi was replaced by that of Tibet. The average air temperature difference between SIM_T_air (−9.83 ^∘C) and CTL (−10.25 ^∘C) was negligible (0.42 ^∘C). In the sensitivity experiment SIM_T_air, the water temperature decreased more quickly, and in late October, the lake froze, no longer releasing heat into the atmosphere. The lake stratification was enhanced, and the water temperature maximum center was ahead by about 10 d (Fig. 5c, f).

In the SIM_q simulation, the specific humidity of Kilpisjärvi was substituted for Tibet. The difference in specific humidity between SIM_q and CTL was 0.38 g kg⁻¹ during the ice-covered period. In the sensitivity experiment SIM_q, the simulations were coincidental to that in CTL, and thus the specific humidity had little effect on the water temperature (Fig. 5c, i).

In conclusion, the stronger downward SR and lower precipitation in the high-altitude and relatively low-latitude TP played positive roles in the water temperature warming during the ice-covered period in Ngoring Lake. Less downward LR, lower air temperature, and larger wind speed did not change the warming trend but affected the warming amplitude and rate. Specific humidity had no significant influence.

Figure 5 The simulated 3 m daily average water temperature in (a), (d) SIM_SR, (a), (g) SIM_Precip, (b), (e) SIM_LR, (b), (h) SIM_U, (c), (f) SIM_T_air, and (c), (i) SIM_q sensitivity experiments from November 2015 to June 2016 is compared with the CTL and the observation, and the change in vertical stratification is shown. The dotted line represents 3.98 ^∘C.​​​​​​​

6.3 Influences of main physical parameters on water temperature

The radiation transfer, which depended on the albedo and extinction coefficient, played a decisive role in the water temperature. Only influences of A_i, E_i, and E_w on water temperature simulation were discussed with sensitivity experiments due to less snow in Tibet. According to previous observations, A_i as observed on the TP was mostly less than 0.12, and the albedo of clear blue ice was only 0.075 (Li et al., 2018). The range of A_i without snow cover was set as 0.1–0.8 with an interval of 0.1 in SIM_A_i experiments.

E_i has not been observed on the TP, but surveys in Finnish lakes show that the value of bare ice varies between 1–4 m⁻¹, while the value of snow-covered ice can reach 5 m⁻¹ (Lei et al., 2011). In SIM_E_i simulations, E_i was equal to 1–5 m⁻¹ with an interval of 1 m⁻¹.

For the E_w, Zolfaghari et al. (2017) found that the FLake model is particularly sensitive at E_w≤0.5 m⁻¹. Shang et al. (2018) observed that E_i varies from 0.11 to 0.67 m⁻¹ in a few TP lakes. Therefore, the sensitivity simulations SIM_E_w were designed so that the E_w varied from 0.1 to 0.5 m⁻¹ with an increment step of 0.1 m⁻¹. The experimental settings are shown in Table 3.

Table 3 Numerical sensitivity simulations of parameters affecting the radiative transfer.

In the SIM_A_i sensitivity experiment, the 3 m water temperature decreased by approximately 1 ^∘C when ice albedo increased by 0.1. When the albedo grew to 0.80, and the water temperature warming decreased from 4 to 2 ^∘C. The increase in ice albedo did not affect the ice formation date but remarkably delayed the ice melting time, accordingly prolonging the ice-covered period. When the albedo increased from 0.1 to 0.8, the ice-covered period was extended for 15–30 d for every increase of 0.1 (Fig. 6a).

In the sensitivity experiment SIM_E_i, the ice extinction coefficient changes did not all make a continuous rise in water temperature, but the 3 m water temperature decreased by 1–2 ^∘C when the ice extinction coefficient increased by 1 m⁻¹ (Fig. 6b). The ice absorbed more heat, and less heat entered the lake water under ice due to the larger ice extinction coefficient.

To further discuss influences of A_i and E_i on lake temperature during the ice-covered period, the period was divided into two stages: Period A and Period B in CTL and SIM_A_i. Period A ranged from freezing point to T_md, and Period B ranged from T_md to maximum temperature (T_m). The duration of Period A was longer than that of Period B, and the temperature heating rate in Period B (∼ 0.1 ^∘C d⁻¹) was 2.5 times greater than that of Period A (∼ 0.04 ^∘C d⁻¹). The reason was that the lake was completely covered by ice, and the inner lake was evenly mixed in Period A, while the ice thickness decreased, and the radiation absorbed by the ice decreased in Period B. The upper layer absorbed more heat than the deeper layer, and the upper water temperature increased rapidly. When A_i and E_i increased, the heating rate decreased, and the duration increased in Period A, while the maximum temperature decreased, and the heating rate and duration fluctuated in Period B. When A_i≥0.6, the heating rate during the ice-covered period decreased and did not rise to T_md.

In the SIM_Ew sensitivity experiment, the water extinction coefficient had just a little influence on the winter water temperature; 3 m water temperature decreased with the increase in Ew (Fig. 6c). When only the extinction coefficient of water changed, the solar radiation entering the water through the ice is unchanged, and so the heat storage of the lake was unaffected. It is just that the heat distribution in the vertical direction is changed. The higher the extinction coefficient of water was, the more heat was absorbed by the surface layer and the less heat reached the deep layer. The phenomenon that the 3 m water temperature decreases with increasing extinction coefficient becomes more and more obvious in the later stage of ice melting.

Figure 6 Comparison of the 3 m simulated daily average water temperature with the observed value under different (a) A_i, (b) E_i, and (c) E_w.

6.4 Influences of water temperature on lake–atmosphere exchange

The thermal conditions in an ice-covered lake just before ice melting have a significant influence on the air–lake energy exchange. To analyze the effects of lake temperature characteristics on the atmosphere at ice melting, three experiments – SIM_E1, SIM_E2, and SIM_E3 (Table 1) – were set up based on the CTL and the observed lake temperature profile on 25 March 2016, 5 d before ice completely melted (Fig. 7a). The characteristics of the initial water temperature profile were as follows.

• SIM_E1. The stratification was weak, the first layer temperature was at the melting point, and from the second layer down the water temperature was set as 2 ^∘C, corresponding to Bangong Co (Wang et al., 2014).

• SIM_E2. The lake was strongly stratified. The first layer was at the melting point, and the temperature increased linearly, reaching T_md at the bottom, corresponding to Valkea-Kotinen lake (Bai et al., 2016).

• SIM_E3. The temperature of the first layer was at the melting point, the temperature gradually increased with the depth from the second layer to the middle layer, and the temperature in the middle layer increased to T_md, corresponding to Thrush Lake (Fang and Stefan, 1996).

In CTL, the first layer temperature was equal to the freezing/melting point, and the second layer reached the maximum temperature on 25 March. The temperature became lower with the deeper layer, until the temperature reached T_md.

With the different initial temperature profiles, the heat storage was different after ice breakup, and the difference persisted for about 2 months (Fig. 7b). In CTL, from 30 to 31 March, when ice melted completely, the lake heat storage ranged from 30893.02 to 30874.51 MJ m⁻², and the heat released was 18.51 MJ m⁻². In the three experiments, from 1 to 2 April, when ice melted completely, the lake heat storage changed from 30657.51 to 30651.67 MJ m⁻² in SIM_E1, from 30781.07 to 30769.91 MJ m⁻² in SIM_E2, and from 30833.28 to 30822.42 MJ m⁻² in SIM_E3, and the heat release was 5.84, 11.16, and 10.86 MJ m⁻², respectively (Fig. 7b).

The heat released was in the form of sensible heat and latent heat, accounting for 0.060 % (CTL), 0.019 % (SIM_E1), 0.036 % (SIM_E2), and 0.035 % (SIM_E3) of the ice-covered heat storage, respectively. As the initial lake temperature profiles were different before the ice completely melted, the ice melted earlier and faster with the higher lake temperature. The lake heat storage increased from 25 March to 24 May, and the heat release rate was different under different circumstances. After late May, the heat balance between the lake and the atmosphere was the same, and so the heat storage basically stayed equal after that.

Figure 7 (a) The initial water temperature profiles in the model are set on 25 March 2016, and the corresponding daily average (b) lake heat storage is simulated. SIM_E1, SIM_E2, and SIM_E3 are three different sensitivity simulations.

Figure 8 The difference in the daily average sensible and latent heat between the CTL and the three sensitivity simulations under the water temperature profiles in SIM_E1, SIM_E2, and SIM_E3.

The lake surface temperature also affected the sensible and latent heat release, whose differences were calculated between CTL and the three experimental simulations (Fig. 8). The influence of different initial water temperature profiles started on 31 March, which is when the ice had melted completely in CTL, and when the sensible and latent differences between CTL and three experimental simulations were less than 0.1 W m⁻² for 3 consecutive days, we judged that the influence had ended. The maximum differences in the sensible heat (51.0 W m⁻²) and latent heat (76.7 W m⁻²) between SIM_E1 and CTL appeared on 31 March and ended on 12 and 30 June, respectively (Fig. 8a). In SIM_E2 the corresponding numbers were 51.4 W m⁻² (31 March to 5 June) for sensible heat and 81.7 W m⁻² (1 April to 17 June) for latent heat (Fig. 8b), and in SIM_E3 they were 51.5 W m⁻² (31 March to 23 May) for sensible heat and 86.0 W m⁻² (1 April to 5 June) for latent heat (Fig. 8c). Compared with the three lake temperature experiments, the heating characteristics of Ngoring Lake made the heat release higher and faster during ice breakup. The duration of heat release difference was from 59 (to 23 May) to 97 (to 30 June) days, and for the latent heat release, the situation lasted about 12–18 d longer than for the sensible heat release.

7 Conclusions

The analysis demonstrates a significant increase in lake temperature during the ice-covered period in Ngoring, the largest freshwater lake on the Tibetan Plateau (TP), with water temperatures exceeding the freshwater maximum density value T_md. The heating is governed by strong solar radiation, the factor making alpine lakes on the high-altitude and relatively low-latitude TP differ from the high-latitude and low-altitude northern lakes with similar winter air temperature patterns. The one-dimensional lake model LAKE2.3 successfully captured the major mechanisms of warming and vertical thermal stratification during the ice-covered period. Compared with MODIS surface temperature data, the bias, RMSE, and CC were 1.42, 3.25 ^∘C, and 0.96, respectively. The absolute values of bias and RMSE were less than 0.1 and 1 ^∘C at 2, 9, 14, and 22 m. The CC of simulated and observed water temperature at 2, 9, and 14 m was as high as 0.98, and the CC of simulated and observed water temperature at 22 m was 0.95.

Sensitivity simulations with perturbed local climate data confirmed the decisive role of subsurface solar radiation in the water temperature rise and demonstrated strong negative feedback with winter precipitation amount. The downward longwave radiation, air temperature, and wind speed had only a minor influence on the water temperature.

The warming rates and, as a result, the duration of the ice-covered period were sensitive to the physical ice properties: ice albedo and light extinction coefficient both reduced the amount of the subsurface solar radiation. An increase in the albedo of ice reduced the rising trend of water temperature and prolonged the ice season. At the critical albedo of 0.6, the lake water warming decreased obviously, and temperature remained stable at no more than 3.98 ^∘C. The extinction coefficient of water had just a minor effect on water temperature under the ice.

An important consequence of the under-ice solar heat accumulation consisted of increased sensible and latent heat releases in the subsequent open-water phase. According to the model results, the effects on the surface fluxes of Ngoring Lake lasted for 59–97 d after the ice melt and increased the upward latent and sensible surface heat fluxes up to ∼ 80 and ∼ 50 W m⁻², respectively. Herewith, the phenomenon of under-ice solar heating may have a significant effect on the land–atmosphere interaction on regional scales and has to be accounted for in coupled climate models.