Study of freeze-thaw cycle and key radiation transfer parameters in a Tibetan Plateau lake using LAKE2.0 model and field observations

Abstract

The Tibetan Plateau (TP) lakes are sensitive to climate change due to its seasonal ice cover, but few studies have paid attention to the freeze-thaw process of TP lakes and its key control parameters. By combining 216 simulation experiments using the LAKE2.0 model with the observations, we evaluated the effects of ice and snow albedo, ice (K_di) and water (K_dw) extinction coefficients on the lake ice phenology, water temperature, sensible and latent heat fluxes. The reference experiment performs well in simulating the lake temperature, with a small positive bias increasing with depth, but it underestimates the ice thickness. The increase of ice albedo, snow albedo and K_di induce a significant decrease in water temperature. Compared with the latent heat, the sensible heat flux is more sensitive to these three parameters. The ice thickness increases almost linearly with the increase of ice albedo but decreases with the increase of K_di. The ice thickness and frozen days vary little with K_dw, but increasing K_dw can decrease the water temperature. Compared with the ice albedo, the K_di and snow albedo have a large effect on the number of frozen days. This study brings to light the necessity to improve the parameterizations of the TP lakes freeze-thaw process.

1. Introduction

There are numerous lakes in the Tibetan Plateau (TP), with a total area of 5.0 × 10⁴ km² (Zhang and others, 2019). Since 1990, the area of TP lakes has expanded significantly due to the increase of glacial meltwater and precipitation (Zhang and others, 2017; Yang and others, 2018). A lot of observational and simulation studies have been carried out on the TP lakes (Biermann and others, 2014; Li and others, 2015, 2017; Wen and others, 2016; Lei and others, 2017; Wang and others, 2017; Huang and others, 2019) and preliminarily clarified the characteristics of water-heat exchange of lakes during the ice-free period. However, there are few reported studies focused on the lake frozen period. Most of the TP lakes show a warming trend (Zhang and others, 2014; Wan and others, 2018), which results in a shorter frozen period (Gou and others, 2017). Especially, the frozen lakes start to melt earlier (Cai and others, 2017). The earlier melted lakes promote the earlier establishment of the thermocline in summer and accelerate the warming of the upper waterbody (Austin and Colman, 2007; Hampton and others, 2008).

The formation of ice is mainly determined by the heat storage of the lake, but the melting of ice mainly depends on the solar radiation (Efremova and Pal'shin, 2011), especially the radiation absorbed by the ice and water, which is closely related to the albedos and extinction coefficients. It has been reported that the melting date is delayed with the increase of the ice albedo (Kirillin and others, 2012). The snow cover also causes the ice cover to expand slightly and thus a later melting date (Martynov and others, 2010). Previous studies found that the ice and snow albedos range from 0.10 to 0.70 and from 0.50 to 0.90, respectively (Grenfell, 1979; Leppäranta and others, 2010), and the albedo parameterization for modeling applications is still not well solved. In our previous studies (Lang and others, 2018; Li and others, 2018), the ice albedo from observations and simulations of the Weather Research and Forecasting Model (WRF) (Wang and others, 2014), Community Land Model 4.5 (Oleson and others, 2013), FLake (Mironov and others, 2010) and WRF–FLake (Mallard and others, 2014) were compared, in which all models overestimated the ice albedo of the TP lakes.

Sunlight transmission through snow, ice and water, in terms of irradiance, is described using the extinction coefficient, which describes the rate at which downwelling irradiance decreases with depth due to scattering and absorption (Lei and others, 2011). Low values of the water extinction coefficient indicate a clear lake where light can penetrate deep into the water column, and high values indicate a turbid lake. The lake model is particularly sensitive to the water extinction coefficient below 0.5 m⁻¹ (Heiskanen and others, 2015). Compared with larger lakes, in small lakes, the water extinction coefficient shows a greater influence on the vertical thermal structure (Woolway and Merchant, 2019). From a purely mechanistic perspective, the water extinction coefficient affects lake hydrodynamics (Persson and Jones, 2008), even within 24 h (Woolway and others, 2015). Guerrero and others (2017) assessed the sensitivity of simulated fluxes to changes in the water extinction coefficient using a 1-D lake model and revealed that the uncertainty in fluxes’ estimations is strongly related to the accuracy with which the extinction coefficient is determined.

For low-transparent ice with a high extinction coefficient, sunlight is absorbed by the surface or the shallow layer. In contrast, for transparent ice, sunlight can reach the water to force convective heating with a fraction of heat returning up to the ice bottom (Leppäranta and others, 2019). However, few studies considered the extinction coefficient of ice in the simulation. Owing to a large variation of the radiation transfer parameters of the natural ice, the precise parameterization of the freeze-thaw process is difficult in the present lake models. Taking the CLM4-LISSS model as an example (Subin and others, 2012; Lawrence and others, 2018), for frozen lakes with snow <40 mm (no snow or snow depth <40 mm), the fraction of shortwave radiation absorbed at the lake surface (β) is set equal to the near-infrared fraction of the shortwave radiation reaching the surface and the remainder of nonreflected shortwave radiation (1-β) is absorbed in the top lake body layer. Usually, the thickness of the top lake body layer is set at 0.1 m in the model. This is unbefitting to lakes with transparent ice and easily leads to large bias under conditions of no or thin snow cover. In the CoLM lake model (Dai and others, 2018), the lake ice is opaque to solar radiation. However, some studies found that the extinction coefficient of the lake ice varies greatly (Lei and others, 2011; Leppäranta and others, 2019). Therefore, it is necessary to carefully consider the model extinction coefficient of lake ice for TP lakes because the snow cover does not always exist and the ice is relatively transparent. Snow cover is one of the most variable elements in the cryosphere and is difficult to simulate its effects during the frozen period (Yang and others, 2013). In the TP with strong solar radiation and thin snow cover, the effects of snow cover on the lake freeze-thaw process have not been fully understood. For improving the simulation, in addition to obtaining accurate parameters through observations, it is also critical to know the sensitivity of parameters in the lake simulation.

In this study, the 1-D model LAKE2.0, developed by Stepanenko and others (2016) is used to quantitatively evaluate the effect of key parameters (ice albedo, snow albedo, ice extinction coefficient, and water extinction coefficient) on the lake ice phenology, the water temperature and the turbulent flux simulation in the TP lakes. The data that drive and validate the model come from the field observations collected from the Ngoring Lake basin from 1 July 2011 to 31 December 2016. This study is structured as follows. Section 2 describes the study area, model, data collection and simulation experiments. The analysis and discussion of the simulation results are presented in Section 3. Finally, Section 4, presents the conclusions.

2. Study area, model, data and numerical experiment design

2.1. Study area

Ngoring Lake (4274 m AMSL) is located in the source region of Yellow River of the eastern TP (Fig. 1e), with a mean depth of 17 m and a surface area of 610 km², which is the highest large fresh water lake in China and regulates the runoff of the source of the Yellow River. As an important habitat of birds, Nogring Lake was listed in the list of wetlands of international importance in 2005. Moreover, the Ngoring Lake basin became one of the core areas of China's first national park (Three–River–Source National Park) in 2018. A cold and semi-arid continental climate prevails in this area, the air temperature varies from 7.8°C in July to −16.2°C in January, averaging −3.5°C (1953–2019) and the average annual precipitation is 326.8 mm (data from Madoi weather station of China Meteorological Administration). From December to April, the lake is usually covered with ice. Our observations show that the ice thickness can reach ~0.7 m which typically happens in late February (Fig. 2a). Around the lake, there are low hills covered with alpine meadow, and there is a small seasonal variation of vegetation height (~5–10 cm).

Fig. 1. Photos of the observation station and the map of the study area created from Landsat image. (a) Tower site at lakeside (TS), (b) the radiation and turbulent flux measurement system in TS, (c) Lake surface site (LS), (d) Grassland site (GS), (e) the map of the study area and the location of observation sites.

Fig. 2. Observed and simulated lake ice thickness (a) and snow depth (b) in a group of experiments with varying ice albedo (A _s = 0.70, K _di = 2.0 m⁻¹ and K _dw = 0.15). Lines of different colors represent different ice albedos. The reference experiment is with A _i = 0.2, which is the red line on the plot.

2.2. LAKE2.0 model and radiation transfer parameters

LAKE2.0 is a 1-D lake model for solving horizontally averaged equations for heat, gas and momentum transport (Stepanenko and others, 2016), which includes the processes of vertical heat transfer and considers the propagation of shortwave radiation in the layers of water, ice, snow and bottom sediments (Stepanenko and Lykossov, 2005; Stepanenko and others, 2011). In ice and snow, a coupled transport of heat and liquid water is reproduced (Volodina and others, 2000; Stepanenko and others, 2019). The model considers snow gravitational compaction and liquid water transport, among others. The mathematical description of these processes closely follows the formulation from Volodina and others (2000). The heat-balance equation at the top of snow includes heat and radiation fluxes, and the bottom boundary condition at the snow/ice interface is continuity of heat flux and temperature (Stepanenko and Lykossov, 2005).

The ice is considered fresh and its thermophysical properties are homogeneous vertically. The increment of ice thickness (dh _i) during time (dt) consists of the ice increment at the upper (dh _i0) and lower boundary (dh _il):

(1)$${\rm d}h_{\rm i} = {\rm d}h_{{\rm i0}} + {\rm d}h_{{\rm il}}$$

For the phase transitions inside the ice, only fresh water is involved in the phase change at both borders of the ice cover. The thickness increment at the lower boundary is the thickness of the formed or thawed layer of ice. The volume of fresh ice ($V_{{\rm ice}}^{{\rm fr}}$) formed during the freezing of volume of water ($V_{\rm w}^{{\rm fr}}$) is

(2)$$V_{{\rm ice}}^{{\rm fr}} = V^{\rm fr}_{\rm w} \rho _{{\rm w}0}/\rho _{{\rm i0}}$$

On top of the ice layer, the ice melting rate depends on the heat balance on its surface. Inside the ice, the heat exchange coefficient is assumed to be constant. The melting point temperature at the ice/water interface equals 0°C, and melted water is added to the water layer (Stepanenko and Lykossov, 2005). More detailed physical processes of ice are described by Stepanenko and others (2019). At the top boundary of sediment, the temperature and heat flux are continuous. All properties of sediments are assumed to be horizontally homogeneous. At the lower boundary, the heat and moisture fluxes are usually set to zero (Stepanenko and others, 2016).

The universal radiation concepts in the transfer through ice-covered lake and related parameters, based on a previous study by Leppäranta (2014), are shown in Figure 3. The albedo regulates the surface energy budget and the extinction coefficient controls the vertical distribution of radiation energy in snow/ice/water layers and has a big annual variability (Iakunin and others, 2020). The extinction coefficient of snow (K _ds) in LAKE2.0 model decreases slightly as the snow density (ρ _s) increases (Stepanenko and others, 2016; Iakunin and others, 2020). The snow cover is almost opaque to solar radiation and even a thin snow cover can completely prevent solar radiation from reaching the ice surface. In Ngoring Lake, the nonreflected radiation can penetrate through water and ice under snow-free conditions because of the high transparency. Moreover, the uncertainty in the water albedo is relatively small. Therefore, only the albedos of ice and snow, as well as the extinction coefficients of ice and water are selected to be analyzed in this study.

Fig. 3. Schematic picture of radiation concepts in the transfer through ice cover lake and related parameters. The prototype of this picture is derived from the study by Leppäranta (2014). The bold grey font represents the main process in solar radiation transfer through snow, ice and water. The blue and red fonts represent the related parameters in these processes.

Before the simulation, we modify the thresholds of a few input variables (surface air pressure and short-wave radiation) in the LAKE2.0 model code to make them suitable for the low air pressure and strong radiation of TP. The surface air pressure is ~600 hpa under the altitude of 4274 m condition (this study). In the TP, the surface downward short-wave radiation is even higher than the solar constant and the maximum value can exceed 1400 W m⁻² (Yang and others, 2010). Therefore, the lower limit of air pressure in the LAKE2.0 code is modified from 800 to 500 hpa and the upper limit of the short-wave radiation is modified from 1400 W m⁻² to 1450 W m⁻². In the model, several different turbulent mixing parameterizations can be selected. We use a standard k–ε (K-epsilon) scheme to calculate eddy diffusivity. The salinity needs to be set, but the effect of salinity on growth of ice is not considered (Stepanenko and others, 2019). In this study, the salinity in the mixed layer and at the bottom is set at 2.7 × 10⁻⁴ kg kg⁻¹ based on a previous study (Shen and others, 2012). The lake area is set at 610 km² and the lake morphometry curve is drawn from the depth maps of Ngoring Lake by Shen and others (2012). The lake depth is set at 25 m based on the water depth in 2015/16 and the lake water is divided into 40 layers in the vertical direction. There is no ice and snow at the initial moment, but the ice contains five layers in the model and the number of snow layers depends on the snow depth. At the lake bottom, there are five layers of soil and the soil type is silt loam. The inflow and outflow of the lake are disregarded because the lake depth and area do not change much in Ngoring Lake and thereby do not affect the thermal balance.

2.3. Data

In the Ngoring Lake basin, there are four observation sites. In the tower site (TS) (97°34′12″E, 34°54′24″N), a 20 m tower was built in October 2010, which was located 30 m west of the lakeshore and equipped with various meteorological measurement sensors (e.g. temperature, humidity, wind speed, precipitation) (Fig. 1a). A radiation and turbulent flux measurement system was established at the lakeside of TS later (Fig. 1b). Grassland site (GS) (97°33′16″E, 34°54′51″N), standing on a flat land surface, is located 1.5 km west of the lakeshore and 1.8 km from the TS site (Fig. 1d). GS was built in August 2011 and has been observing ever since. GS observes the same variables as TS. Lake surface site (LS) (97°38′59″E, 35°01′28″N) was built on a submerged rock and located 200 m to the northwest lakeshore (Fig. 1c). LS was located to the north-west from the center of the lake. Observations of LS were only carried out in the ice-free period of 2011–13 and the observed variables included conventional meteorological elements and turbulent fluxes. From 22 September 2015 to 21 September 2016, a full year lake temperature profile was observed in Ngoring Lake, and the water depth of this site (WS) was ~23–25 m.

The observations used in this study are listed in Table 1. In LAKE2.0, the forcing data, ranging from 1 July 2011 to 31 December 2016, consist of the following: air temperature, air specific humidity, air pressure, meridional and zonal wind speeds, downward short-wave and long-wave radiation fluxes, precipitation rate. Therein, except for the radiation variables, all the other forcing data come from TS. The downward short-wave and long-wave radiation fluxes from July 2011 to October 2012 are mainly obtained from GS and the missing data (~8.5% of the total) during the same period are replaced by the data from LS. Since November 2012, most of the radiation forcing data are derived from TS, and a few missing values are replaced by the data from GS and LS (~3% of total). The information of sampling frequency is listed in Table 1. The original sampling data with high frequency were uniformly processed into 30 min average data by the data collectors (e.g. CR3000, CR1000) (Campbell Sci., USA). The ice thickness was measured manually only in two frozen periods (2011/12 and 2015/16) near the TS, with a water depth of ~5 m. For the snow depth on the lake surface, there were no observations. The numbers of frozen days (the frozen area ⩾80% of lake area), estimated artificially from MODIS daily snow cover products (MOD10A1 and MYD10A1) combined with EOSDIS Worldview images, are ~164, 134, 141, 142 and 121 from 2011 to 2016, respectively.

Table 1. The information of observation data used in this study

2.4. Numerical experiment design

In this study, with a variation of four key parameters (ice albedo, snow albedo, ice extinction coefficient and water extinction coefficient), we designed 216 numerical experiments. Based on different snow albedo, ice and water extinction coefficients, these experiments were combined in three groups (Exice, Exwat and SnowA) (Table 2). In these experiments, the ice albedo (A _i) ranged from 0.10 to 0.80 with an increment step of 0.10, the snow albedo (A _s) ranges from 0.50 to 0.90 with an increment step of 0.05. These ranges were chosen according to the studies by Grenfell (1979) and Leppäranta (2014). There are no observations for the ice extinction coefficient (K _di) in the TP lakes, but one survey shows that this value varies from 1 to 5 m⁻¹ under snow-free conditions in the northern Europe lakes (Lei and others, 2011), which can be used for reference. In this study, the ice extinction coefficient varied from 1 to 5 m⁻¹ with an increment step of 0.5 m⁻¹. For the water extinction coefficient (K _dw), Zolfaghari and others (2017) found that the FLake model is particularly sensitive to water extinction coefficient values below 0.5 m⁻¹. A recent survey shows that it varies from 0.11 to 0.67 m⁻¹ in a few TP lakes (Shang and others, 2018). Therefore, in this study, the water extinction coefficient is only one-tenth of K _di and varies from 0.1 to 0.5 m⁻¹ with an increment step of 0.05 m⁻¹.

Table 2. Numerical experiment design

Explanations for parameter symbols are given in Section 2.4.

The initial lake temperature profile is derived from observations in Ngoring Lake in July 2011. To reduce the impact of the initial field, model simulations are repeated ten times using the data of the first year before the experiments. The time step of forcing and output data for simulation is set to be 1 hour in this study. Therefore, the observation data with 30 min average for every two groups are arithmetically averaged to generate the forcing data with an hourly average that meets the input of the model.

3. Results

3.1. Forcing data

The time series of hourly average variables in the forcing data are presented in Figure 4. From July 2011 to December 2016, the air temperature ranges from −31.9 to 19.9°C (Fig. 4a), the wind speed ranges from 0 to 20.4 m s⁻¹, in which the wind speed exceeding 6 m s⁻¹ accounts for 25%, and most of them occur in winter and spring (Fig. 4b). From 2011 to 2016, the wind speed in the frozen period, especially the events with high wind speed, shows a decreasing trend. The peak of downward short-wave radiation is 1215.5 W m⁻² (Fig. 4c) and the daily peak values can exceed 1000 W m⁻² from early March of each year. Peak values of downward long-wave radiation (>300 W m⁻²) in each year usually appear about a month and a half later than short-wave radiation. For the precipitation during the cold period (16 November–30 April), they are 39.6, 9.0, 33.5, 36.5 and 39.6 mm in each cold period from 2011 to 2016, respectively. Especially, precipitation from 2012 to 2013 is extremely little and is only a quarter of those in the other years. Although precipitation in the other years is similar in value, its distribution with time varies dramatically (Fig. 4d). Most of the precipitation in the 2013/14 cold period are concentrated in March and April, while 61.9% of precipitation in the 2015/16 happens in November and February. This difference may have a significant impact on the growth and melting of lake ice.

Fig. 4. Time series of hourly average (a–c) and daily accumulated (d) values of a few variables in the forcing data. (a) air temperature at 2 m height, (b) wind speed at 10 m height, (c) downward short-wave (SR) and long-wave (LR) radiation fluxes, (d) precipitation. For the precipitation, the warm period values from 1 May to 15 November are omitted.

3.2 Results in the reference experiment

3.2.1 Water temperature

Based on previous investigations (Lei and others, 2011; Li and others, 2018; Shang and others, 2018), the simulation with the parameters A _i = 0.20, A _s = 0.70, K _di = 2.0 m⁻¹ and K _dw = 0.15 m⁻¹ is chosen as the reference experiment. For the ice albedo, the observations show that 93.4% of the samples at noon are below 0.16 in Ngoring Lake (Li and others, 2018). The albedo during the day has a U-shaped diurnal cycle affected by the solar zenith angle. In Li and others (2018), the daily average albedo was calculated as the average of instantaneous albedos weighted by the instantaneous downward short-wave radiation flux at each hour interval during the daytime following the method by Wang and others (2015). As a result, a daily average ice albedo defined as 0.2 is appropriate. The value of K _dw is based on observations of Nam Co Lake (Nima and Zhuo, 2012) and Pangong Tso Lake in the TP (Shang and others, 2018) because there is no information available about Ngoring Lake. The ice extinction coefficient of TP lakes is difficult to estimate accurately due to the lack of observations. However, a study in Finland has found that the extinction coefficient of bare ice mostly varies from 1 to 2.5 m⁻¹ (Lei and others, 2011), and that, there is a positive correlation between the extinction coefficient and albedo for the lake ice. Our observations have shown that the ice albedo is small in Ngoring Lake (Li and others, 2018), thus, the ice extinction coefficient is likely also small.

As shown in the observations (Fig. 5a), there is a uniform distribution for the water temperature in the vertical direction from late September to the following March, and the lowest temperature appears in mid-December. From the middle of March, a weak stratification appears in the waterbody, with a warming of the upper layer. After the ice melting, the water density increases due to the cooling of the upper layer, which induces the lake water to turn over and mix again. In mid-June, the temperature of the middle and upper lake increases significantly, and the water temperature begins to stratify. It continues until September, and the maximum temperature over a year occurs from mid-August to early September. Compared with observations, except for the slightly higher minimum water temperature in December, the simulated water temperature shows a good agreement from late September to the following February. In March (Fig. 5b), the lake is stratified in the model run, but the stratification disappears in mid-April. After the ice melts, the warming rate of simulated temperature is faster than that observed. Besides, the simulated temperature is 1–2°C higher than the observed value. In later June, an obvious stratification appears and it is relatively deep. The simulated occurrence and disappearance dates of the stratification are consistent with the observed. In terms of water temperature, LAKE2.0 presents good performance. For the whole water column, the observed and simulated annual average temperatures are 6.3°C and 6.9°C, respectively, and the bias and the root mean square error (RMSE) are 0.6°C and 1.4°C. The simulated temperature is higher than the observed value and presents a larger error at the deep layer. For example, the bias and RMSE of the 15 m depth water temperature are 0.8°C and 1.8°C respectively, while that of the 3 m depth is 0.4°C and 1.2°C.

Fig. 5. Observed (a) and simulated in the reference experiment (b) water temperature from September 2015 to September 2016 with the ice albedo of 0.20, the snow albedo of 0.70, the ice extinction coefficient of 2.0 m⁻¹ and water extinction coefficient of 0.15 m⁻¹.

3.2.2 The ice thickness and frozen period

The time series of observed and simulated daily average ice thickness in the reference experiment is shown in Figure 2a. Therein, the maximum observed thickness is 0.73 m in 2 years of observations. The peaks of simulated ice thickness are less than the observations (2011/12 and 2015/16), with a 0.2–0.3 m underestimation. In the model, the lake water depth is set to 25 m, much deeper than that of the ice thickness measuring location (~5 m). In a single lake, the ice thickness is usually inversely correlated with the water depth. As a result, the ice thickness in the deep-water area should be less than the observed value, but the exact difference is unknown due to no observations. In terms of frozen period, the simulated results are in good agreement with the observations in the years 2012/13, 2014/15 and 2015/16, while the simulated periods are clearly shorter in another 2 years (2011/12 and 2013/14). The simulation error in the melting date is greater than that in the freezing date, and the largest error in the melting date appears in 2011/12, with 60 days earlier than the observation. This error may be a result of a combination of factors, and improvements in the ice melting simulation are still needed for this model in the future.

3.3. Influence of key parameters on water temperature

From Figures 6 and 7, it can be found that the temperature differences, derived from simulations with different ice extinction coefficients, ice and snow albedos, mainly appear in the frozen period and the following 2 months. The most significant differences occur before ice melting. The simulated water temperature is systematically higher than the observed value from ice melting to middle October, while the simulated value is lower from late October to early December.

Fig. 6. Difference between simulated and observed water temperature (simulation-observation) at 3 m depth from September 2015 to September 2016 in Exice experiments. Lines of different colors represent different ice extinction coefficients (K _di).

Fig. 7. Difference between simulated and observed water temperature (simulation-observation) at 3 m depth from September 2015 to September 2016 in SnowA experiments. Lines of different colors represent different snow albedos (A _s).

In Exice experiments (Fig. 6), the simulated water temperature at K _di = 1.0 m⁻¹ in the frozen period is significantly higher than others when the ice albedo is 0.1. With increasing K _di, the water temperature drops significantly. However, when K _di ⩾ 3.5 m⁻¹, the water temperature before and after ice melting is significantly lower than the observed value, and the pattern is different. With increasing A _i, the water temperature in the frozen period significantly decreases but does not show a linear trend. For example, when the ice albedo is 0.4, the water temperature is relatively high. The bare ice albedo mainly varies from 0.05 to 0.25 in Ngoring Lake during the whole frozen period (Li and others, 2018). Taking an ice albedo of 0.2 as an example, the pattern of simulated temperature is quite close to the observations at 1.5 ⩽ K _di ⩽ 3.0 m⁻¹, which indicates that it is reasonable to set K _di = 2.0 m⁻¹ in the reference experiment. When the ice albedo is 0.8 (obviously not reasonable) and the K _di is small, the simulated temperature is also close to the observed value. This means that the large ice albedo can compensate an effect of the high ice transparency on the simulated water temperature at 3 m depth.

In Exwat experiments (figure not shown), the simulated water temperatures are quite close to each other with low ice albedo. Similarly, at A _i = 0.20 and K _dw = 0.15 m⁻¹, the water temperature during the melting period can be well reproduced. With increasing ice albedo, simulated results of K _dw = 0.10 and 0.15 m⁻¹ are higher than those at other values. Different from Exice experiments, besides the frozen period, the simulated water temperatures with different K _dw also present obvious difference (0.5–1°C) in July and September.

Within the group of SnowA experiments (Fig. 7), the reference experiment with the ice albedo of 0.2 and the snow albedo of 0.7, the simulated value is in good agreement with the observed value in water temperature. In these experiments, if the snow albedo is not <0.7, the simulated water temperature is lower in the frozen period, especially with high ice albedo. With the same A _i, the difference between simulated and observed peaks of the temperature before melting appear later with the increase of A _s, which is mainly related to the longer frozen days.

Besides 3 m water temperature, the simulation bias in the water temperature at 15 m and the whole water column is also presented in Figure 8. In Exice experiments (Figs 8a–c), the bias is stable when K _di ranges from 1.5 to 3.5 m⁻¹. When K _di is <1.5 m⁻¹, the bias increases significantly but decreases clearly or even becomes a negative value when larger than 3.5 m⁻¹ (Figs 8a–c). However, for the water column temperature, the simulation bias of low A _i is larger than that of high A _i when K _di is <3.5 m⁻¹ (Fig. 8c). At K _di ⩾ 3.5 m⁻¹, the pattern is almost the opposite (Fig. 8c). However, at a certain range, the bias does not vary linearly with the A _i, whether in the shallow (3 m) or deep layer (15 m). For example, when 1.5 ⩽ K _di ⩽ 3.5 m⁻¹, the bias is the smallest with A _i of 0.7, followed by 0.5 (Fig. 8a), while the bias is the same at the ice albedo of 0.6 and 0.2. For Exwat experiments (Figs 8d–f), the difference among experiments is relatively small, but the bias still decreases with the increase of the water extinction coefficient. The relation between ice albedo and bias is similar to Exice. In SnowA experiment (Figs 8g–i), when the snow albedo varies from 0.60 to 0.75, the bias remains stable. Moreover, the variation of bias with ice albedo is similar to Exice.

Fig. 8. Mean bias in the water temperature (simulation-observation) calculated from September 2015 to September 2016. (a–c) Exice experiments, (d–f) Exwat experiments, (g–i) SnowA experiments, (a, d, g) at 3 m depth, (b, e, h) at 15 m depth, (c, f, i) for the whole water column. Lines of different colors represent different ice albedos (A _i).

3.4. Influence of key parameters on ice thickness and the number of frozen days

The ice thickness and the number of frozen days have a significant impact on the lake surface energy budget in winter and spring. For the simulation results, only when the ice thickness was larger than 0.05 m for 7 consecutive days or more, it was recorded into the frozen days, because it is difficult to form a stable ice sheet if the lake ice is too thin, and a few cold air events can also cause lakes to freeze intermittently outside the frozen period.

Figure 2 represents a group of experiments with varying ice albedo, which contains also the reference experiment (A _s = 0.70, K _di = 2.0 m⁻¹ and K _dw = 0.15). From this figure, in years with more snowfall, the difference in simulated maximum ice thickness of different albedo experiments is up to 0.2–0.3 m, while in 2012/13 with little snowfall, the difference in ice thickness caused by ice albedo increases significantly and reaches 0.44 m. In terms of the frozen period, the differences among different experiments are relatively small. In the cold period of 2013/14, the snowfall (33.5 mm) is only slightly less than in other years, but the difference in simulated ice thickness is the largest (0.6 m). As can be seen from Figure 2b, there is little snow on the lake in the early and middle frozen period (from December 2013 to February 2014), and the ice albedo plays a leading role in regulating the surface energy budget. Taking the frozen period of 2011/12, the simulated ice thickness is closest to the observed value with A _i of 0.2 in the early frozen period (December); after that, the simulated value is closest to the observed value when the ice albedo ranges from 0.6 to 0.7. This may be related to the huge difference in lake depth between the ice thickness measuring point and the model, as described in Section 3.2.2. Meanwhile, the simulated melting date is earlier than the observed. At the melting stage, all the lake ice simulated by different albedo experiments melt rapidly, with little date difference. In snowy years (e.g. 2011/12 and 2014/15), the maximum ice thickness generally appears later than in less snow years (e.g. 2012/13 and 2013/14).

In Exice experiments (Fig. 9), when the K _di value is reasonable (2.0 m⁻¹), simulated number of frozen days and ice thickness at small ice albedo are closer to the observations in 2012/13 and 2015/16. The number of frozen days tends to increase with decreasing K _di, especially at low K _di. However, the variation rate varies significantly from year to year, which may be related to the snow cover days on the ice surface. In the cold period of 2012/13, the observed snowfall is the lowest, but the simulated snow days are not (figure not shown). The winter period with the fewest snow days is from 2013 to 2014. In these two periods, the number of frozen days only increases slightly with decreasing K _di, which means that the effect of K _di on the number of frozen days is little.

Fig. 9. Simulated number of frozen days, mean and maximum ice thickness in Exice experiments every winter period. (a, f, k) 2011/12, (b, g, l) 2012/13, (c, h, m) 2013/14, (d, i, n) 2014/15, (c, j, o) 2015/16. Lines of different colors represent different ice albedos.

The maximum ice thickness shows a similar variation trend with the mean thickness. With a low K _di, solar radiation can penetrate through lake ice and reach the ice bottom or even the waterbody. As a result, the radiation is easily absorbed by the waterbody with higher specific heat capacity and the melting rate of ice will be slowed down. At the same time, the simulated ice thickness varies greatly in different years, but the observed interannual variation is small (Fig. 2a).

In SnowA experiments (Fig. 10), the number of frozen days and ice thickness increase with increasing snow albedo. Like Exice experiments, the maximum ice thickness under low ice albedo condition is closer to the observations in less snow year (2012/13). In snowy years, the result becomes more complex.

Fig. 10. Simulated number of frozen days, mean and maximum ice thickness in SnowA experiments every winter period. (a, f, k) 2011/12, (b, g, l) 2012/13, (c, h, m) 2013/14, (d, i, n) 2014/15, (c, j, o) 2015/16. Lines of different colors represent different ice albedos.

3.5. Influence of key parameters on sensible and latent heat fluxes

The 5-year period (from 1 December 2011 to 30 November 2016) average sensible and latent heat fluxes of Exice and SnowA experiments are shown in Figure 11. In Exice experiments (Figs 11a, b), the H is more sensitive to changes in K _di and A _i than the LE. Maximum H (35.1 W m⁻²) is 2.3 times the minimum value (15.2 W m⁻²), while LE during the same period is only 1.3 times. It can in part be explained by the different ways H and LE respond to varying lake temperature. Variations in lake temperature caused by K _di and A _i directly change the H by affecting the lake-air temperature difference. However, the influence of lake temperature on LE is indirect because the lake temperature affects the LE by changing the specific humidity on the lake surface. With the same K _di, the H decreases with increasing A _i. With the same ice albedo, H increases with the increase of K _di, mainly because the frozen days are shortened when K _di is large. H over a water surface is usually larger than over an ice surface. The observations show that H over the TP lakes is usually positive (from the lake surface to atmosphere) even at the stage when the ice has just melted (Li and others, 2015), which is different from the lakes at high latitudes (Verburg and Antenucci, 2010). When the A _i is not more than 0.6, H, A _i and K _di show a linear trend except for a few cases. Simulated LE during the frozen period is usually ~2–4 times larger than H. The difference between H and LE becomes larger with increasing frozen days, mainly because the simulated daily average H on the ice surface is mostly negative (from the atmosphere to lake surface).

Fig. 11. Simulated sensible (a, c) and latent (b, d) heat fluxes in experiments Exice (a, b) and SnowA (c, d) averaged for the period from 1 December 2011 to 30 November 2016. The unit for H and LE: W m⁻².

In the SnowA experiments, both H and LE decrease with increasing ice and snow albedos, because the increase in albedo causes the frozen days longer. Like Exice experiments, the maximum H is 3 times the minimum value, while the difference of LE is only 1.5 times. At the same time, H and LE are more susceptible to variation in snow albedo than in ice albedo. This may be related to the long-term snow cover on the ice during the frozen period, which weakens the influence of ice albedo. Meanwhile, with increasing snow albedo, both H and LE are more sensitive to ice albedo, especially for H. In this study, the snow cover is intermittent and thin. Taking the reference experiment as an example, the snow days with snow depth (⩾3.3 mm) only accounts for <20% of the frozen period in 5 years. Therefore, the simulation bias caused by inaccurate ice albedo cannot be ignored during most of the frozen period. The sum of H and LE can represent the intensity that the lake surface heats or cools the atmosphere. When the ice and snow albedos remain the minimum, the sum of H and LE reaches 113.8 W m⁻², while when both of albedos are the maximum, the sum is only 66.7 W m⁻². Based on the area of 610 km² of Ngoring Lake, the annual difference in the flux exchange at the air/lake interface is ~9.1 × 10¹¹ MJ.

3.6. Statistics of ice thickness, snow depth, the number of frozen days and water temperature depending on ice and snow radiation transfer parameters

The lake ice and snow isolate the lake water from the atmosphere, which makes the mass and energy exchange during the frozen period more complicated. The respective effects of four radiation transfer parameters (ice albedo, snow albedo, ice and water extinction coefficient) have been discussed in the previous sections. In this section, we compare the statistics of ice thickness and snow depth, the number of frozen days combined with 3 m water temperature depending on four radiation parameters in the cold period (16 November–15 May) from 2011 to 2016.

Variations of simulated ice thickness, snow depth, number of frozen days and water temperature along with ice albedo and extinction coefficient are shown in Figure 12. The ice thickness increases almost linearly with increasing ice albedo (Fig. 12(a)-2), with a maximum difference of 0.23 m. Conversely, the ice thickness approximately decreases linearly with the increase of ice extinction coefficient (Fig. 12(a)-3) and the effects of two parameters on ice thickness are almost a trade-off. The snow depth almost remains unchanged with A _i and only slightly reduced when K _di value is larger than 4.0 m⁻¹ (Fig. 12b). The effect of K _di on number of frozen days is significantly greater than that of A _i (Fig. 12c). The former can last up to 54 days, while the latter has only 14 days. Moreover, the number of frozen days is most sensitive to K _di <2.0 m⁻¹ but approximatively remains unchanged at A _i > 0.4. For the water temperature (Fig. 12d), the influences of K _di and A _i show good consistency and the water temperature also increases dramatically with the decrease of K _di when K _di ⩽ 2.0 m⁻¹. Taking 2.0 m⁻¹ of K _di for example, the solar radiation can penetrate the lake ice into the waterbody if the ice thickness is <0.5 m. When the K _di is <2.0 m⁻¹, more radiation energy is used to heat water rather than affect the lake ice as the K _di decreases. The low K _di reduces the energy of heating the ice, which is conducive to the ice growth and delays the melting of ice. However, it also increases the energy incident into the water and makes the water warmer, which is conducive to the melting of the bottom ice and shortens the frozen period. In this study, the number of frozen days also increases rapidly with the decrease of K _di when K _di ⩽ 2.0 m⁻¹, which means the melting of the bottom ice due to increased water temperature is not important.

Fig. 12. Variations in arithmetic mean values of simulated ice thickness (a), snow depth (b), number of frozen days (c) and water temperature at 3 m depth (d) along with ice albedo and extinction coefficient in Exice experiments in five cold periods (16 November–15 May) of years 2011–16. Taking figure (a) as an example (other graphs are similar), the X-axis and Y-axis represent ice albedo and ice extinction coefficient and the color lines represent the ice thickness in the figure (a)-1. The variations of the simulated mean ice thickness along with the ice albedo (figure (a)-2) and ice extinction coefficient (figure (a)-3) are shown with the line chart, respectively.

In Exwat experiments (Fig. 13), the snow depth almost remains unchanged with the water extinction coefficient (K _dw) and the ice thickness and number of frozen days increase very slightly when K _dw decreases. The increase in water temperature with decreasing K _dw is less than that of Exice experiments, but it is also larger than 1°C (Fig. 13(d)-3). This suggests that the increase in water temperature is mainly caused by the thermal change of lake water induced by the decrease of K _di, rather than the change in ice thickness and frozen days.

Fig. 13. Similar to Figure 12, but in Exwat experiments. The X-axis and Y-axis represent ice albedo and water extinction coefficient in figure (a)-1.

In SnowA experiments (Fig. 14), when grouped by ice albedo, the results are like those of Exice and Exwat experiments. While different, the ice albedo equal to 0.3 is a special turning point for the number of frozen days and water temperature in Exice experiments (Figs 12(c)-2 and 12(d)-2). This phenomenon is not found in SnowA experiments. When grouped by the snow albedo, the ice thickness, snow depth and the number of frozen days all significantly increase with increasing A _s. It is worth noting that the sensitive range for snow albedo varies from variable to variable. For example, the ice thickness varies little when 0.6 ⩽ A _s ⩽ 0.8, but this interval only ranges from 0.6 to 0.7 for the number of frozen days. When A _s > 0.7, the number of frozen days increases rapidly with increasing A _s, but for the snow depth, this turning point is A _s = 0.8 (Fig. 14(b)-3). Compared with the ice albedo, the snow albedo has a more significant effect on the number of frozen days. When the snow albedo ranges from 0.90 to 0.5, the differences in the number of frozen days is 50 days (Fig. 14(c)-3). Meanwhile, the maximum difference in the number of frozen days due to different ice albedo is only 20 days (Fig. 14(c)-2). Both the increasing ice and snow albedos reduce the water temperature and show a fair amount of influence (Fig. 14d).

Fig. 14. Similar to Figure 12, but in SnowA experiments. The X-axis and Y-axis represent ice albedo and snow albedo in figure (a)-1.

4. Conclusions and discussion

In this study, a 1-D lake model LAKE2.0 is used to quantitatively evaluate the effect of four radiation transfer parameters (ice albedo, snow albedo, ice and water extinction coefficients) on the lake ice phenology, water temperature and turbulent flux simulation in a TP lake. The main conclusions are as follows:

From the middle of March to ice melting, the water temperature increases in the upper waterbody and a weak stratification can be observed. In the reference experiment, the simulated water temperature is higher than the observed value (with an average bias of 0.6°C) and presents a larger bias at the deep-water layer. The peaks of simulated ice thickness are less than the observed, with a 0.2–0.3 m underestimation. Three groups of experiments (Exice, Exwat and SnowA) were performed. In Exice experiments, the influences of K _di and A _i on the water temperature show good consistency. When K _di ranges from 1.5 to 3.5 m⁻¹, the temperature bias is stable. When K _di is <1.5 m⁻¹, the bias increases significantly and decreases clearly when larger than 3.5 m⁻¹. Compared with the LE, the H is more sensitive to changes in K _di and A _i. The ice thickness increases almost linearly with an increase of A _i but decreases with an increase of K _di, and the effects of two parameters on ice thickness are almost a trade-off. The effect of K _di on number of frozen days is significantly greater than that of A _i. Moreover, the number of frozen days is most sensitive to K _di <2.0 m⁻¹ but approximatively remains unchanged at A _i > 0.4.

Different from Exice experiments, besides the frozen period, the simulated water temperatures with different K _dw also show obvious differences in July and September in Exwat experiments. The water temperatures are close to each other with low ice albedo. Three variables (ice thickness, snow depth, frozen days) vary little with K _dw and the decreasing rate of water temperature with the increase of K _dw is less than that of Exice experiment.

In SnowA experiments, when grouped by the snow albedo, the ice thickness, snow depth and the number of frozen days all significantly increase with the increase of snow albedo, but the sensitive ranges for snow albedo vary from variable to variable. For example, the ice thickness varies rapidly when A _s < 0.6 or A _s > 0.8, but this range is A _s > 0.7 for the number of frozen days. Compared with ice albedo, the snow albedo has a more significant effect on the number of frozen days. Both the increasing ice and snow albedos reduce the water temperature and show a fair amount of influence.

In this study, the influence of the ice extinction coefficient is quantitatively evaluated. This parameter is seldom paid attention to in studies of this kind, while it is valuable for improving lake models. However, limitations of 1-D offline model such as the absence of feedback processes and horizontal energy transfer may reduce the scientific nature of results to some extent. For example, if the feedback processes between lake and atmospheres are considered, the atmospheric variables should vary with snow and ice conditions in different experiments, especially the air temperature, which in turn will affect the surface water-heat exchange. However, the feedback processes will also increase the uncertainty of simulations, because the atmosphere and lake models are far from perfect. In Section 3.6, the average over five cold periods can present the general patterns, but to a certain extent, it also masks the differences at annual and monthly scales. Although our analysis shows that such differences have difficulties affecting the major patterns (figure not shown), it will still be interesting to understand the causes of these differences in future work. Therefore, both more adequate parameterization improvements in the offline model and an in-depth study using the coupled model are necessary in the future.