Observation and Analysis of Water Temperature in Ice-Covered Shallow Lake: Case Study in Qinghuahu Lake

: Water temperature serves as a key environmental factor of lakes and the most basic parameter for analyzing the thermal conditions of a water body. Based on the observation and analysis of the water temperature of Qinghuahu Lake in the Heilongjiang Province of China, this paper analyzed the variation trend of the heat ﬂux, effective thermal diffusivity of the icebound water, and revealed the temporal and spatial variation law of the water temperature and the transfer law beneath the ice on a shallow lake in a cold region. The results suggested a noticeable difference existing in the distribution of water temperature beneath the ice during different periods of ice coverage. During the third period, the water temperature vertically comprised three discrete layers, each of which remained unchanged in thickness despite the alternation of day and night. Sediment– water heat ﬂux and water–ice heat ﬂux both remained positive values throughout the freezing duration, averaging about 3.8–4.1 W/m 2 and 9.8–10.3 W/m 2 , respectively. The calculated thermal diffusivity in late winter was larger than molecular, and the time-averaged values increased ﬁrst and then decreased with water depth, reaching a maximum at a relative depth of 0.5. This research is expected to provide a reference for studies on the water environment of icebound shallow lakes or ponds in cold regions.


Introduction
The water temperature in lakes is fundamental to the understanding of their various physical and chemical processes and dynamic phenomena. It represents an important environmental factor of water ecosystems and also a vital factor influencing the metabolism and productivity of aquatic organisms [1,2]. Previous observational studies on the water temperature of lakes and reservoirs across the globe were principally carried out in temperate and tropical regions, which yielded fruitful results [3][4][5][6]. Freezing is a natural phenomenon that water in a cold region normally encounters in winter. Observational studies on water temperature beneath the ice are scarce due to the complex process of ice regime and difficulties in fieldwork [7]. The formation of ice sheets weakens the energy exchange between air and water and mitigates the disturbances caused by the wind to a lake. It also reduces the penetration depth of solar radiation passing vertically through a water body and slows the drop in water temperature beneath the ice [8]. Such changes in physical conditions result in differences in thermal conditions and temperature variations between lake water in cold regions and that in temperate and tropical zones [9]. The resulting unique water environment beneath the ice then exerts an influence on the biochemical process in water. Consequently, differences in the distribution and transformation of substances are visible between lakes in cold regions and those in temporal and tropical regions [10].
There have been many studies on the thermal conditions of deep or shallow lakes in cold regions before [11][12][13][14][15]. Zhou et al. [12] used a one-dimensional water temperature and negative temperature of Qinghuahu Lake in winter is between about −1458 and −2053 • C. The duration of ice coverage lasts 5 to 6 months, during which frost depth of the ground averages 1.63 m and stands at 1.86 m at a maximum, and the maximum ice thickness is between 1.0 and 1.3 m. Additionally, Qinghuahu Lake originates in a low-lying region, where drainage is impeded and rivers meander. It appears shallow, and the slopes of its lake basins are gentle. These hydraulic characteristics and ice regimes are typical features shared by the cluster of natural lakes in the Northeast China Plain. In consequence, research designed to observe and analyze the water temperature of Qinghuahu Lake is representative of studies on the temperature of shallow lakes in northeast China.
In situ observation of water temperature in this research was completed near the west bank of the pond that lies in the middle of Qinghuahu Lake, see Figure 1a. Water at the measurement point is around 1.87 m deep. Qinghuahu Lake starts to freeze in late October and early November every year and its ice sheets begin to melt in mid to late April in the following year. When the flood season arrives in July and August, the water level increases in this lake, which is one of the places in Daqing City that regulates excess water to control floods. Nevertheless, in winter, there is no river flowing into or out of the lake; therefore, the water level in this lake basically remains stable, and the water body creates a lentic environment that provides ideal conditions for in situ observation. duration averages 2840 h, and the average daily temperatures from October to the following April ranges from -18.9 to -12.3 °C. The cumulative negative temperature can directly reflect the cold situation of the local winter.
Selecting all days with a daily mean temperature below 0 °C and summing up these daily mean temperatures can result in a cumulative negative temperature. The cumulative negative temperature of Qinghuahu Lake in winter is between about −1458 and −2053 °C. The duration of ice coverage lasts 5 to 6 months, during which frost depth of the ground averages 1.63 m and stands at 1.86 m at a maximum, and the maximum ice thickness is between 1.0 and 1.3 m. Additionally, Qinghuahu Lake originates in a low-lying region, where drainage is impeded and rivers meander. It appears shallow, and the slopes of its lake basins are gentle. These hydraulic characteristics and ice regimes are typical features shared by the cluster of natural lakes in the Northeast China Plain. In consequence, research designed to observe and analyze the water temperature of Qinghuahu Lake is representative of studies on the temperature of shallow lakes in northeast China.
In situ observation of water temperature in this research was completed near the west bank of the pond that lies in the middle of Qinghuahu Lake, see Figure 1a. Water at the measurement point is around 1.87 m deep. Qinghuahu Lake starts to freeze in late October and early November every year and its ice sheets begin to melt in mid to late April in the following year. When the flood season arrives in July and August, the water level increases in this lake, which is one of the places in Daqing City that regulates excess water to control floods. Nevertheless, in winter, there is no river flowing into or out of the lake; therefore, the water level in this lake basically remains stable, and the water body creates a lentic environment that provides ideal conditions for in situ observation.

Observational Devices and Methods
This research analyzed previously collected data on local air temperatures and ice regimes and observed in situ the water temperature and meteorological factors of Qinghuahu Lake from 1 October 2017 to 1 May 2018, a period that covered the entire process from the ice forming to melting. The purpose was to explore the thermodynamic process and its corresponding mechanisms in this lake under the condition of ice coverage in winter. Factors observed included ice thickness, and the temperature of air, water, ice, and mud at the lake's bottom. Meteorological data of the observation location were accessed from the A753 WS, an automatic wireless weather station. Temperatures were measured through the PT100 platinum resistance temperature sensor (RH-8068 thermal

Observational Devices and Methods
This research analyzed previously collected data on local air temperatures and ice regimes and observed in situ the water temperature and meteorological factors of Qinghuahu Lake from 1 October 2017 to 1 May 2018, a period that covered the entire process from the ice forming to melting. The purpose was to explore the thermodynamic process and its corresponding mechanisms in this lake under the condition of ice coverage in winter. Factors observed included ice thickness, and the temperature of air, water, ice, and mud at the lake's bottom. Meteorological data of the observation location were accessed from the A753 WS, an automatic wireless weather station. Temperatures were measured through the PT100 platinum resistance temperature sensor (RH-8068 thermal resistance) developed by Roham Company, with a measurement accuracy of ±0.05 K. A thermistor chain, to which PT100 resistance temperature detectors were attached 5 cm apart from each other, was laid vertically before observation from lake surface to 10 cm below sediment, to automatically monitor variations in ice, water, and mud temperatures (see Figure 1b). The chain was held in place with a wooden frame of low thermal conductivity, and its top part was locked in ice as the water body gradually froze. Temperature data were transmitted to a CR1000X data logger, and the time taken by each data acquisition was 30 min. Ice thickness was measured with a sliding staff through a hole that an ice auger bored.

Calculation Method of Heat Flux and Thermal Diffusivity
The thermal conditions of the icebound water and thermal evolution of ice thickness were both related to water-ice heat flux Q wi at the upper boundary of water and sediment-water heat flux Q sw at the lower boundary. Therefore, it was necessary to analyze their seasonal variations according to the temperature data obtained during the observation process.

Sediment-Water Heat Flux
A considerable amount of heat that a lake absorbs from the atmosphere in spring and summer is stored in the thermally active layer at the top part of the bottom sediment and released back to water in autumn and winter. Heat transferred upward from sediment remains little in an icebound deep lake due to a profundal zone [28], whereas in a shallow lake, its bottom sediment is a vital source of heat. In 1927, Birge et al. [29] maintained that sediment at a shallow lake's bottom is crucial to its energy balance, and subsequent studies [30,31] drew similar conclusions. Results, however, varied substantially from research on sediment-water heat flux at the bottom of ice-covered shallow lakes under different geographical, meteorological, and hydrological conditions.
For the direct measurement of heat flux, the most commonly used measurement method is the use of thermal resistance heat flux sensor, the principle of which is based on Fourier's law, that is, to calculate the heat flux using the temperature gradient, and this is also the expression of heat flux. In this observation, since the temperature sensors had been installed, we chose to use the observed temperature data to obtain the heat flux indirectly through Fourier's law. According to Fourier's law, sediment-water heat flux can be calculated through the following formula: where λ w and λ s are the thermal conductivity of water and sediment measured, respectively (W/(mK)); dT w /dz and dT s /dz are the water temperature gradient and the sediment temperature gradient measured near the sediment-water interface, respectively (K/m); and T s1 and T s2 are the mud temperatures measured by utilizing the two nearest sensors below the sediment-water interface ( • C). Notably, the Fourier law is the law of heat conduction, when heat conduction is the only form of heat transfer, the molecular thermal conductivity of water can be plugged into the first formula on the right side of Equation (1) to calculate heat flux. Nevertheless, for water above the sediment-water interface, the heat transfer process includes not only heat conduction, but also probably heat convection, because solar radiation penetrating into water will intensify temperature and density heterogeneity among different layers, probably generating buoyancy-driven natural convection, which will transfer heat vertically among water layers. Therefore, it is necessary to firstly determine whether a heat convection process has taken place in the water.
The Grashof number Gr, a dimensionless number derived from the momentum equation of convective heat transfer, can serve as a determinant of the occurrence and transformation of thermal convection. Its expression is as follows: where g is the gravitational acceleration (m/s 2 ); α is the coefficient of volumetric expansion of water, which amounts to 2.563 × 10 −8 K −1 at 3.98 • C; ∆T is the temperature difference between the upper and lower boundaries of water, and here ∆T = 3.98°C, which is obtained with the corresponding water temperature at its maximum density (3.98 • C) minus the freezing point (0 • C); L is the characteristic length, which herein represents the water depth beneath the ice during the period of stable freezing, here the water depth corresponding to the maximum ice thickness was adopted as the L value, the non-freezing water depth at the measurement point is 1.87 m, and the maximum ice thickness is 1.02 m (see Section 3.1, paragraph 3); therefore, L = 1.87 m-1.02 m = 0.85 m; ν is the kinematic viscosity of water, amounting to 1.792 × 10 −6 m 2 /s at 0 • C. The identification of heat transfer in icebound water can be simplified as the same matter in a horizontal interlayer within a confined space. When Gr < 2430, heat in the interlayer is principally transferred by thermal conduction, and when Gr ≥ 2430, natural convection comes into existence in the interlayer and gradually becomes the major driving force of heat transfer [32]. By substituting the above values into Expression (2), we obtain Gr = 1.907 × 10 5 2430, indicating that there is apparently natural convection inside the layers of icebound water. If so, the thermal conductivity of water λ w in Equation (1) is its effective thermal conductivity [33] that takes natural convection heat transfer into account, rather than its molecular thermal conductivity λ wm [34]. The effective thermal conductivity of water λ w during natural convection proves difficult to determine due to its complex changes. Nevertheless, Equation (1) relies on the assumption that the heat flux at the interface is a continuous property across the interface; therefore, the heat-flux values of heat flux at the two sides are equal. Therefore, as a comparison, we applied the temperature gradients on both sides of the interface to calculate the heat fluxes.
When the temperature gradient on one sediment side was used for calculation, the thermal conductivity of the sediment λ s was calculated using the following formula [34]: where λ d is the thermal conductivity of soil matrix, set as 2.0 W/(m·K); λ wm is the molecular thermal conductivity of water, set as 0.59 W/(m·K); and θ sat is the moisture content of surficial sediment, which is normally set as 50% in a shallow lake [35].

Water-Ice Heat Flux
Heat flux at the water-ice interface, which occurs at the freezing front, can be calculated through the following Fourier law formula [36]: where λ i is the thermal conductivity of ice, standing at 2.14 W/(mK); dT w /dz and dT i /dz, respectively, represent the water temperature gradient and the ice temperature gradient near the freezing front measured in K/m; and T i1 and T i2 are temperatures in • C measured by the sensors closest to the freezing front. Similar to Equation (1), given that thermal conductivity of water λ w remains uncertain under the possible condition of natural convection heat transfer, the water-ice heat fluxes were calculated using the temperature gradients of both ice and water.

Effective Thermal Diffusivity
Thermal diffusivity is a measure of the rate at which heat diffuses from the hot side to the cool side in a substance. Essentially, it is an indicator used in the study of heat transfer to describe a material's ability to reach thermal equilibrium inside itself. It is also known as thermometric conductivity since the greater the thermal diffusivity is, the more rapidly temperature changes and spreads. It is of paramount importance to understand the properties of heat transfer in an unsteady state.
According to the definition formula of thermal diffusivity, thermal diffusivity is calculated by the following [36]: where λ is a substance's thermal conductivity (W/(mK)); ρ is density (kg/m 3 ); c is specific heat capacity (J/(kgK)); and ρc is the heat per unit of volume of a substance required to increase the temperature by 1 K. Under the circumstances where heat transfer among parts in a substance is not subject to relative displacements, the molecular thermal diffusivity under such specific conditions would be obtained when the λ substituted into the above formula is the molecular thermal conductivity of the substance. Conversely, the actual thermal diffusivity of icebound water, in which natural convection is involved, is called effective thermal diffusivity in this study. It can be obtained based on variations in the vertical profiles of water temperature within a fixed time frame (see Figure 2), even though direct calculation through Equation (5) would be in vain.
the properties of heat transfer in an unsteady state. According to the definition formula of thermal diffusivity, thermal diffusivity is calculated by the following [36]: where λ is a substance's thermal conductivity (W/(mK)); ρ is density (kg/m 3 ); c is specific heat capacity (J/(kgK)); and ρc is the heat per unit of volume of a substance required to increase the temperature by 1 K. Under the circumstances where heat transfer among parts in a substance is not subject to relative displacements, the molecular thermal diffusivity under such specific conditions would be obtained when the λ substituted into the above formula is the molecular thermal conductivity of the substance. Conversely, the actual thermal diffusivity of icebound water, in which natural convection is involved, is called effective thermal diffusivity in this study. It can be obtained based on variations in the vertical profiles of water temperature within a fixed time frame (see Figure 2), even though direct calculation through Equation (5) would be in vain. Next, we tried to establish an expression for the efficient thermal diffusivity. The starting point of the consideration is the heat flux balance of a horizontally homogeneous layer, given by the following first law: where the quantity HS denotes the sensible heat flux in units of W m −2 , the vertical coordinate z and the heat flux HS count positive in nadir direction. Averaging Equation (6) over the time interval 2 1 t t t Δ = − yields the following: Next, we tried to establish an expression for the efficient thermal diffusivity. The starting point of the consideration is the heat flux balance of a horizontally homogeneous layer, given by the following first law: where the quantity H S denotes the sensible heat flux in units of W m −2 , the vertical coordinate z and the heat flux H S count positive in nadir direction. Averaging Equation (6) over the time interval ∆t = t 2 − t 1 yields the following: The integration of Equation (7) over a layer of thickness ∆z = z B − z > 0 yields the following: Equation (8) relies on the assumption that the height dependence of the mass density and specific heat of water can be neglected. This assumption is safely satisfied for a shallow lake. To specify the boundary conditions, the sensible heat flux at the bottom of the lake is taken to be the heat flux across the sediment-water interface, i.e., H S t (z B ) = −Q sw . The negative sign considers the fact that the positive sediment-water heat flux Q sw = 3.8W m −2 > 0 is directed from the sediment to the water body aloft, i.e., in negative z direction. The sensible heat flux at the reference level z is parameterized employing the small-eddy approximation underlying the classical down gradient ansatz as follows: By inserting H S t (z B ) and H S t (z ) into Equation (8), one arrives at the following equation: By approximating the integral in the numerator of Equation (10) by the sum and inserting Equation (11) into Equation (10), we finally arrive at an equation for the determination of effective thermal diffusivity A eff :

Ice Formation and Melting
The temperature remains a principal factor contributing to the formation and growth of ice cover in the water of cold regions. Autumn and winter reveal a gradual decrease in water temperature as the water body constantly transfers heat to the low temperature atmosphere. A thin ice layer starts to form when the lake's surface cools to 0 • C. Ice needles and flowers then aggregate on the surface, and ultrathin ice cover emerges. As the air temperature continues to decline, the water body beneath the ice continuously loses heat to maintain a thermal equilibrium and eventually changes into ice, thickening the existing ice cover.
Take Qinghuahu Lake in the winter of 2017-2018, for example. The thin ice that formed for the first time on the lake's surface on the evening of 22 October 2017 disappeared when the temperature increased soon after the sun rose the next morning, a phenomenon that reoccurred for the next three consecutive days. Early on the morning of 25 October, the stable ice belt took shape along the bank and persisted in the daytime, despite the rising temperature. The ice coverage gradually extended to the surface of the area where the water went deeper. On 26 October, the ice cover at the surface closed down except for a few strips of water. The next day, a stable ice cover formed that coated the entire surface of the lake. Figure 3 demonstrates the processes of ice formation and melting and variations in the average daily temperatures in Qinghuahu Lake during the winter (2017-2018). The period of rapid growth in ice cover lasted from 27 October to early February of the next year, during which ice cover grew quickly thanks to the fast heat loss from water that was attributed to the constant decrease in air temperature and thin ice cover. Meanwhile, the freezing front moved downward rapidly and continuously. The increase in ice thickness retarded heat loss from the icebound water. Between early February and early March, the ice cover grew slowly. The ice thickness reached its maximum value on 9 March, amounting to 102 cm. Starting from mid-March 2018, the ice cover began to melt as the average daily temperatures rose above 0 • C, and the ice melting accelerated, particularly after 2 April. By 21 April, the entire ice cover disappeared. The growth and thaw period spanned 132 d and 44 d, respectively, resulting in an ice-formation period of 176 d. The ice growth and melting rates averaged 0.79 and 2.36 cm/d, respectively. One caveat was that the snow cover has an important influence on the ice process. During the winter of observation, there was little snowfall, and the light snow was often accompanied by strong winds; therefore, it was difficult for snow to stay and accumulate on the ice cover. During the observation period, no snow cover of a sufficient area and thickness was formed.

Temporal Variations of Water Temperature
Temporal variations of water temperature are subject to the equilibrium of all the factors contributing to the heat budget. The variations in the water temperature beneath the ice in the first three days of ice formation are shown in Figure 4. Diurnal variations in water temperature, as it displays, were evident with daily changes in solar intensity and air temperature. The water temperature peaked at 5-6 p.m. every day and hit bottom at 6-8 a.m. The magnitude of temperature variations proved greater in surface water than any other part and decreased at increasing depths. Meanwhile, phase transitions lagged behind changes in air temperature. The effects of air temperature on the water temperature beneath the ice declined during the first few days of ice formation. Despite that, the hypolimnion temperature, which was subject to the influence of original thermal storage and ground temperature, apparently rose by 1-2 K to above 1.5 °C. The water temperature from the surface to the bottom ranged from 0.5 to 1.8 °C shortly after thin ice appeared on the surface of the lake. Three days later, the near-bottom temperature rose to nearly 4 °C. The gradual increase in ice thickness weakened the water-ice heat exchange and solar radiation that penetrated into the water body. Consequently, a gradual decline in the magnitude of diurnal temperature variations occurred. When the ice thickness exceeded 30 cm, apparent diurnal variations in the water temperature beneath the ice were no longer visible. Instead, temperature inversion in a steady state occurred in icebound water, a phenomenon of water temperature that increased with depth.

Temporal Variations of Water Temperature
Temporal variations of water temperature are subject to the equilibrium of all the factors contributing to the heat budget. The variations in the water temperature beneath the ice in the first three days of ice formation are shown in Figure 4. Diurnal variations in water temperature, as it displays, were evident with daily changes in solar intensity and air temperature. The water temperature peaked at 5-6 p.m. every day and hit bottom at 6-8 a.m. The magnitude of temperature variations proved greater in surface water than any other part and decreased at increasing depths. Meanwhile, phase transitions lagged behind changes in air temperature. The effects of air temperature on the water temperature beneath the ice declined during the first few days of ice formation. Despite that, the hypolimnion temperature, which was subject to the influence of original thermal storage and ground temperature, apparently rose by 1-2 K to above 1.5 • C. The water temperature from the surface to the bottom ranged from 0.5 to 1.8 • C shortly after thin ice appeared on the surface of the lake. Three days later, the near-bottom temperature rose to nearly 4 • C. The gradual increase in ice thickness weakened the water-ice heat exchange and solar radiation that penetrated into the water body. Consequently, a gradual decline in the magnitude of diurnal temperature variations occurred. When the ice thickness exceeded 30 cm, apparent diurnal variations in the water temperature beneath the ice were no longer visible. Instead, temperature inversion in a steady state occurred in icebound water, a phenomenon of water temperature that increased with depth. peared on the surface of the lake. Three days later, the near-bottom temperature rose to nearly 4 °C. The gradual increase in ice thickness weakened the water-ice heat exchange and solar radiation that penetrated into the water body. Consequently, a gradual decline in the magnitude of diurnal temperature variations occurred. When the ice thickness exceeded 30 cm, apparent diurnal variations in the water temperature beneath the ice were no longer visible. Instead, temperature inversion in a steady state occurred in icebound water, a phenomenon of water temperature that increased with depth.  In spring, due to the rising air temperature and solar radiation, the ice thickness decreased rapidly, leading to more solar radiation penetrating through ice cover into the In spring, due to the rising air temperature and solar radiation, the ice thickness decreased rapidly, leading to more solar radiation penetrating through ice cover into the water body. Eventually, the resulting increase in water temperature boosted the waterice heat flux, which accelerated the melting of the ice cover. The variations in the water temperature during the late period of ice melting presented in Figure 5 suggests that temperature inversion persisted during the late winter despite an evident increase in the water temperature near the freezing front. A rapid increase in water temperature occurred, and evident diurnal variations resumed when the complete disappearance of ice cover on 21 April made a direct energy exchange between the water and the air possible. The temperatures of water layers at the depths of 150 cm and above were vertically distributed in a normal mode, decreasing with an increased depth, yet the temperature at a depth of 170 cm near the lake's bottom was clearly higher than that at a depth of 160 cm. The reasons may be twofold. One is that the two parameters of water transparency and salt content in the near-bottom layers were higher than that in the water above, which may weaken the development of convection. Then, the water temperature at a depth of 160-170 cm was a remnant of the winter quiescent layer that had not been disturbed by radiatively driven convection. The other is that the warmer weather led to higher ground temperatures and increased heat flux between the bottom and the water body, which contribute to the rise in bottom water temperatures. water body. Eventually, the resulting increase in water temperature boosted the water-ice heat flux, which accelerated the melting of the ice cover. The variations in the water temperature during the late period of ice melting presented in Figure 5 suggests that temperature inversion persisted during the late winter despite an evident increase in the water temperature near the freezing front. A rapid increase in water temperature occurred, and evident diurnal variations resumed when the complete disappearance of ice cover on 21 April made a direct energy exchange between the water and the air possible. The temperatures of water layers at the depths of 150 cm and above were vertically distributed in a normal mode, decreasing with an increased depth, yet the temperature at a depth of 170 cm near the lake's bottom was clearly higher than that at a depth of 160 cm. The reasons may be twofold. One is that the two parameters of water transparency and salt content in the near-bottom layers were higher than that in the water above, which may weaken the development of convection. Then, the water temperature at a depth of 160-170 cm was a remnant of the winter quiescent layer that had not been disturbed by radiatively driven convection. The other is that the warmer weather led to higher ground temperatures and increased heat flux between the bottom and the water body, which contribute to the rise in bottom water temperatures.

Vertical Profiles of Water Temperature
Epilimnion differed from hypolimnion in phases induced by changes in water temperature, generating vertical spatial differences in temperature. Meanwhile, temporal differences also existed in the vertical distribution of water temperature during different periods of ice coverage due to variations in meteorological factors and ice thickness. The duration of ice coverage in this research comprised three periods. The first was the initial stage of ice formation, which lasted from the formation of ice cover on the lake's surface to the moment of ice cover reaching 30 cm in thickness. Stable freezing, the second period, commenced when the ice thickness exceeded 30 cm and ended once constant melting of the ice cover begun. Lastly, the late winter spanned from the end of the second period to the complete disappearance of the ice cover. Figure 6 shows vertical profiles of the water temperature during the three periods of the ice-cover duration, in which observa-

Vertical Profiles of Water Temperature
Epilimnion differed from hypolimnion in phases induced by changes in water temperature, generating vertical spatial differences in temperature. Meanwhile, temporal differences also existed in the vertical distribution of water temperature during different periods of ice coverage due to variations in meteorological factors and ice thickness. The duration of ice coverage in this research comprised three periods. The first was the initial stage of ice formation, which lasted from the formation of ice cover on the lake's surface to the moment of ice cover reaching 30 cm in thickness. Stable freezing, the second period, commenced when the ice thickness exceeded 30 cm and ended once constant melting of the ice cover begun. Lastly, the late winter spanned from the end of the second period to the complete disappearance of the ice cover. Figure 6 shows vertical profiles of the water temperature during the three periods of the ice-cover duration, in which observational data on water temperature were acquired at 10 a.m. every day. It should be noted that since the ice thickness changed over time, and the icebound water depth was also time-varying, the relative water depth (ratio of local depth and maximum water depth) was adopted for convenience and for consistency in the analysis.  During the first period, temperature inversion occurred vertically in the water on 28 October when the ice cover took shape. Nevertheless, the temperature gradient appeared small, standing at around 0.8 K/m, and the water temperature was uniformly distributed. As the ice cover gradually thickened, the bulk temperature of the water went up and the temperature gradient widened. Meanwhile, the gradient in the hypolimnion, up to 7.6 During the first period, temperature inversion occurred vertically in the water on 28 October when the ice cover took shape. Nevertheless, the temperature gradient appeared small, standing at around 0.8 K/m, and the water temperature was uniformly distributed. As the ice cover gradually thickened, the bulk temperature of the water went up and the temperature gradient widened. Meanwhile, the gradient in the hypolimnion, up to 7.6 K/m at a maximum, was significantly larger than that in the epilimnion.
In a stable freezing period, the water temperature was linearly distributed in the vertical direction, while the temperature gradient averaged around 5.3 K/m. Ice cover was relatively thicker, making it more difficult for the sensible atmospheric heat and solar radiation to be transferred to the water. The heat budget of the water simply comprised the heat flux of the lake bottom Q gw (incoming heat) and the water-ice heat flux Q wi (outgoing heat). Previous research suggested that the difference between Q sw and Q wi in this period was low [37]. Consequently, the water's heat budget maintained an equilibrium and the distribution of water temperature remained stable.
For the late winter, the intensified solar radiation caused by rising temperatures penetrated through the ice cover into the water body. The variations were complex in the vertical distribution of the water temperature, which specifically featured three layers [38]. The range from the ice-water freezing front to 5-10 cm downward is the conduction layer (CL), where the water temperature rises rapidly with increasing depth and the temperature gradient stands at 20-40 K/m. Beneath the CL is the mixed layer (ML), where water temperature is uniformly distributed vertically. The reason is that the intensified solar radiation during the late winter increases the temperature and density heterogeneity among the different layers of water. The resulting natural convection driven by buoyancy facilitates heat exchange and mass transfer between the upper and lower layers, leading to the uniform distribution of a vertical water temperature [39]. The hypolimnion beneath the ML constitutes the quiescent layer (QL), in which the water temperature is linearly distributed in a stable state since influences imposed by solar radiation on the water body are minimized by the great depth. Figure 6c displays that water temperature, particularly in the CL and ML, displaying an overall increase overtime during the late winter. Variations in the thickness of the ML and QL appeared significantly, while the temperature gradient of the CL widened constantly. For instance, the water temperature of the ML on the fringe of the CL rose by about 1.9 K between 30 March and 10 April. The thickness of the CL remained unchanged, whereas that of the ML increased from 14 to 61 cm, which accounted for over two thirds of the icebound water. The temperature gradient of the CL rose from 40 to 64 K/m over time. In contrast, the gradient of the QL remained at around 4 K/m. Figure 7 demonstrates diurnal variations in the vertical distribution of the average water temperature on 6 and 7 April during the late winter, and the parts between the dotted and solid lines in this figure represent the standard deviations of the average water temperature. The diurnal variations in the thickness of water layers were considerably small. The average water temperature in the daytime was evidently higher than that in the nighttime, while the maximum variation in the diurnal temperature, which stood at 1.4 • C, occurred in the transitional zone between the CL and the ML. The temperature gradient of the CL in the daytime was around twice that at nighttime. The standard deviation of the daytime water temperature was obviously higher than that of the nighttime at a given depth, indicating that fluctuations in daytime temperatures were more drastic than in nighttime temperatures. erably small. The average water temperature in the daytime was evidently higher than that in the nighttime, while the maximum variation in the diurnal temperature, which stood at 1.4 °C, occurred in the transitional zone between the CL and the ML. The temperature gradient of the CL in the daytime was around twice that at nighttime. The standard deviation of the daytime water temperature was obviously higher than that of the nighttime at a given depth, indicating that fluctuations in daytime temperatures were more drastic than in nighttime temperatures.    (1) and (3) based on observational temperature data. As a comparison, Figure 8 shows the heat fluxes calculated by using both sediment temperature gradients and water temperature gradients. It can be seen from the figure that the two kinds of results both show that the sediment-water heat flux throughout the whole freezing duration remained as a positive value, implying that the sediment was releasing heat into the water. The sediment-water heat fluxes calculated using water temperature gradients are distributed in the range 3.4~5.3 W/m 2 , and the average value is 4.1 W/m 2 during the whole freezing period, while the fluxes calculated by sediment temperature gradients are distributed in the range 3.5~4.6 W/m 2 , with an average of 3.8 W/m 2 . The average value calculated using the water temperature gradient is slightly larger than that by sediment temperature gradient, and the difference is 0.3 W/m 2 , which is equivalent to 7.3 and 7.9% of them, respectively.   (1) and (3) based on observational temperature data. As a comparison, Figure 8 shows the heat fluxes calculated by using both sediment temperature gradients and water temperature gradients. It can be seen from the figure that the two kinds of results both show that the sediment-water heat flux throughout the whole freezing duration remained as a positive value, implying that the sediment was releasing heat into the water. The sediment-water heat fluxes calculated using water temperature gradients are distributed in the range 3.4~5.3 W/m 2 , and the average value is 4.1 W/m 2 during the whole freezing period, while the fluxes calculated by sediment temperature gradients are distributed in the range 3.5~4.6 W/m 2 , with an average of 3.8 W/m 2 . The average value calculated using the water temperature gradient is slightly larger than that by sediment temperature gradient, and the difference is 0.3 W/m 2 , which is equivalent to 7.3 and 7.9% of them, respectively.  Figure 9 demonstrates the temporal variations in water-ice heat flux at the water-ice interface calculated through Equation (4) based on the observational data temperature. As a comparison, Figure 9 also shows the heat fluxes calculated by using both ice temperature gradients and water temperature gradients. It can be seen from the figure that the two kinds of results both show that the water-ice heat flux throughout the whole freezing duration remained as a positive value, indicating that the water body was emitting heat into the ice cover. In addition, the two results both show the rising trend as time passed. The water-ice heat fluxes calculated using the water temperature gradients are distributed in the range 7.3-13.1 W/m 2 , and the average value is 10.3 W/m 2 during the whole freezing period, while the fluxes calculated by ice temperature gradients are distributed in the range 6.9-12.1 W/m 2 , with an average of 9.8 W/m 2 . The average heat flux calculated using the water temperature gradient is also larger than that using the ice temperature gradient, and the difference is 0.5 W/m 2 , which is equivalent to 4.7 and 5.1% of them, respectively. The average value of water-ice heat flux, 9.8 or 10.3 W/m 2 , is greater than some results obtained in studies on deep lakes [40,41] (5-7 W/m 2 ). Therefore, the calculation of ice thickness in a shallow lake would be far from reasonable if the wa-  Figure 9 demonstrates the temporal variations in water-ice heat flux at the water-ice interface calculated through Equation (4) based on the observational data temperature. As a comparison, Figure 9 also shows the heat fluxes calculated by using both ice temperature gradients and water temperature gradients. It can be seen from the figure that the two kinds of results both show that the water-ice heat flux throughout the whole freezing duration remained as a positive value, indicating that the water body was emitting heat into the ice cover. In addition, the two results both show the rising trend as time passed. The water-ice heat fluxes calculated using the water temperature gradients are distributed in the range 7.3-13.1 W/m 2 , and the average value is 10.3 W/m 2 during the whole freezing period, while the fluxes calculated by ice temperature gradients are distributed in the range 6.9-12.1 W/m 2 , with an average of 9.8 W/m 2 . The average heat flux calculated using the water temperature gradient is also larger than that using the ice temperature gradient, and the difference is 0.5 W/m 2 , which is equivalent to 4.7 and 5.1% of them, respectively. The average value of water-ice heat flux, 9.8 or 10.3 W/m 2 , is greater than some results obtained in studies on deep lakes [40,41] (5-7 W/m 2 ). Therefore, the calculation of ice thickness in a shallow lake would be far from reasonable if the water-ice heat transfer was neglected and the lower boundary of the ice cover was regarded as thermally isolated or observed values of the water-ice heat flux from a deep lake were adopted. Overall, the average relative difference of heat flux calculated by using the temperature gradients on two sides of either interface is less than 7.9%. As it happens, the temperature observation operation was not very strictly controllable due to the limitation of natural conditions, and the average relative difference within 7.9% can be regarded as an insignificant statistical difference for the prototype observation. Therefore, the natural convection heat transfer process in water cannot be arbitrarily judged by the difference between the results of the two calculation methods. For the sediment-water heat flux, the average relative difference between the values calculated using the water side and either side is 7.9%, larger than that of water-ice heat flux, which might be due to the uncertainty of some of the parameters in the calculation of sediment thermal conductivity.

Heat Flux at Interfaces
In addition, we found that the average value of sediment-water heat flux calculated by the sediment temperature gradient is 3.8 W/m 2 , and the average value of water-ice heat flux calculated using the ice temperature gradient is 9.8 W/m 2 , with a difference of 6 W/m 2 (if the water temperature gradient is adopted, the difference is 6.2 W/m 2 ). From the relationship between the two average values, the heat flux released by the water body from the water-ice interface to the ice is greater than the heat flux obtained by the water from the sediment-water interface. If only these two parts of the heat budget existed in water, the water temperature should continue to decrease with time, but as it can be seen from Figure 6, the water temperature increases gradually with the passage of time during the whole freezing duration. Therefore, there must be other heat sources to make up for the difference between Qsw and Qwi and supply additional heat to the water, so that the water temperature continues to rise. Obviously, the heat source is the solar radiation penetrating into the water. Next, the solar radiation penetrating into water during the freezing duration will be evaluated to verify the above hypothesis.
The solar radiation flux is related to latitude and date. According to the relevant meteorological research results [42], we obtained the variation process of the solar radiation flux in the prototype observation area, whose latitude is 45 °N. Through a fitting regression, we developed an empirical expression of the solar radiation flux with the ordinal number of dates. The solar radiation will be reflected on the ice surface; therefore, the incident solar radiation flux from the ice surface is as follows: Overall, the average relative difference of heat flux calculated by using the temperature gradients on two sides of either interface is less than 7.9%. As it happens, the temperature observation operation was not very strictly controllable due to the limitation of natural conditions, and the average relative difference within 7.9% can be regarded as an insignificant statistical difference for the prototype observation. Therefore, the natural convection heat transfer process in water cannot be arbitrarily judged by the difference between the results of the two calculation methods. For the sediment-water heat flux, the average relative difference between the values calculated using the water side and either side is 7.9%, larger than that of water-ice heat flux, which might be due to the uncertainty of some of the parameters in the calculation of sediment thermal conductivity.
In addition, we found that the average value of sediment-water heat flux calculated by the sediment temperature gradient is 3.8 W/m 2 , and the average value of water-ice heat flux calculated using the ice temperature gradient is 9.8 W/m 2 , with a difference of 6 W/m 2 (if the water temperature gradient is adopted, the difference is 6.2 W/m 2 ). From the relationship between the two average values, the heat flux released by the water body from the water-ice interface to the ice is greater than the heat flux obtained by the water from the sediment-water interface. If only these two parts of the heat budget existed in water, the water temperature should continue to decrease with time, but as it can be seen from Figure 6, the water temperature increases gradually with the passage of time during the whole freezing duration. Therefore, there must be other heat sources to make up for the difference between Q sw and Q wi and supply additional heat to the water, so that the water temperature continues to rise. Obviously, the heat source is the solar radiation penetrating into the water. Next, the solar radiation penetrating into water during the freezing duration will be evaluated to verify the above hypothesis.
The solar radiation flux is related to latitude and date. According to the relevant meteorological research results [42], we obtained the variation process of the solar radiation flux in the prototype observation area, whose latitude is 45 • N. Through a fitting regression, we developed an empirical expression of the solar radiation flux with the ordinal number of dates.
Q rad = −0.0001d 3 + 0.0233d 2 + 0.6326d + 112.53 (13) where Q rad is the solar radiation flux, W/m 2 ; d is the ordinal number of the date, d on 1 January is 1, and d on 31 December of the previous year is −1; the application scope of Equation (13)  The solar radiation will be reflected on the ice surface; therefore, the incident solar radiation flux from the ice surface is as follows: (14) where Q rad-ice is the incident solar radiation flux from the ice surface, W/m 2 ; A is the albedo of ice, the value is about 0.5. As part of Q rad-ice is absorbed by the ice, the solar radiation flux penetrating into the water is as follows: where Q rad-water is the solar radiation flux penetrating into the water, W/m 2 ; A is the albedo of the ice, the value is about 0.5; γ is the extinction coefficient of ice; and z is the ice thickness, m. In combination with Equations (13)- (15), the Q rad-water values during the freezing duration were calculated and are listed in Table 1. It can be seen from Table 1 that in the whole freezing duration, the solar radiation flux penetrating into water is distributed in the range of 7.34-162.42 W/m 2 . Except for the large flux values in the initial period of ice formation and late ablation period, the average value of solar radiation flux in the stable sealing period (when the ice thickness is greater than 30 cm) is 12.48 W/m 2 , which is greater than the average value of the difference between Q wi and Q sw of 6-6.2 W/m 2 ; therefore, the water temperature will continue to rise in the whole freezing duration. Figure 10 presents the calculated values of effective thermal diffusivity A eff varied with water depths and time. It shows that the A eff in water during the late winter was in the range of 0.32-1.37 × 10 −6 m 2 /s, which was 2.3-9.8 times the molecular thermal diffusivity, amounting to 0.14 × 10 −6 m 2 /s. A eff varied with the position of relative water depth ( Figure 10). Figure 11 is drawn based on the time-averaged values of A eff at different water depths in late winter. It displays that the average value of A eff in this period first increased and then decreased with relative water depth. The A eff value in the middle relative depth was significantly larger, and the maximum value, 0.97 × 10 −6 m 2 /s, occurred approximately at the relative depth of 0.5. It indicated that the water temperature in the mixed layer area tended to be more uniformly distributed since the natural convection there proved to be more drastic and water temperature spread more quickly between the different layers. In contrast, the water near the freezing front (conduction layer) and near sediment (quiescent layer) had larger vertical temperature gradients and the temperatures were more inhomogeneously distributed in comparison to the mixed layer, due to the small effective thermal diffusivities and slow spread of water temperature. average value of solar radiation flux in the stable sealing period (when the ice thickness is greater than 30 cm) is 12.48 W/m 2 , which is greater than the average value of the difference between Qwi and Qsw of 6-6.2 W/m 2 ; therefore, the water temperature will continue to rise in the whole freezing duration. Figure 10 presents the calculated values of effective thermal diffusivity Aeff varied with water depths and time. It shows that the Aeff in water during the late winter was in the range of 0.32-1.37 × 10 −6 m 2 /s, which was 2.3-9.8 times the molecular thermal diffusivity, amounting to 0.14 × 10 −6 m 2 /s. Aeff varied with the position of relative water depth ( Figure 10). Figure 11 is drawn based on the time-averaged values of Aeff at different water depths in late winter. It displays that the average value of Aeff in this period first increased and then decreased with relative water depth. The Aeff value in the middle relative depth was significantly larger, and the maximum value, 0.97 × 10 −6 m 2 /s, occurred approximately at the relative depth of 0.5. It indicated that the water temperature in the mixed layer area tended to be more uniformly distributed since the natural convection there proved to be more drastic and water temperature spread more quickly between the different layers. In contrast, the water near the freezing front (conduction layer) and near sediment (quiescent layer) had larger vertical temperature gradients and the temperatures were more inhomogeneously distributed in comparison to the mixed layer, due to the small effective thermal diffusivities and slow spread of water temperature.   average value of solar radiation flux in the stable sealing period (when the ice thickness is greater than 30 cm) is 12.48 W/m 2 , which is greater than the average value of the difference between Qwi and Qsw of 6-6.2 W/m 2 ; therefore, the water temperature will continue to rise in the whole freezing duration. Figure 10 presents the calculated values of effective thermal diffusivity Aeff varied with water depths and time. It shows that the Aeff in water during the late winter was in the range of 0.32-1.37 × 10 −6 m 2 /s, which was 2.3-9.8 times the molecular thermal diffusivity, amounting to 0.14 × 10 −6 m 2 /s. Aeff varied with the position of relative water depth ( Figure 10). Figure 11 is drawn based on the time-averaged values of Aeff at different water depths in late winter. It displays that the average value of Aeff in this period first increased and then decreased with relative water depth. The Aeff value in the middle relative depth was significantly larger, and the maximum value, 0.97 × 10 −6 m 2 /s, occurred approximately at the relative depth of 0.5. It indicated that the water temperature in the mixed layer area tended to be more uniformly distributed since the natural convection there proved to be more drastic and water temperature spread more quickly between the different layers. In contrast, the water near the freezing front (conduction layer) and near sediment (quiescent layer) had larger vertical temperature gradients and the temperatures were more inhomogeneously distributed in comparison to the mixed layer, due to the small effective thermal diffusivities and slow spread of water temperature.

Conclusions
In this study, a field observation of Qinghuahu Lake, a typical ice-covered shallow lake of Heilongjiang Province, was carried out to clarify the thermal conditions of icebound water. The temporal variations and vertical profile patterns of water temperature throughout the ice-cover duration were analyzed, and heat fluxes at interfaces and effective thermal diffusivity inside the water were, thus, evaluated based on temperature data.
The temporal variation and vertical distribution of the water temperature in the frozen shallow lake are obviously affected by the process of the ice's growth and decay. The whole ice duration can be divided into three different periods according to ice thickness to analyze the water temperature, i.e., ice formation period, stable freezing period, and late winter. During the initial stage of the ice formation period, the vertical gradient of the water temperature is small. With the gradual thickening of the ice cover, the overall water temperature and gradient increase. In the stable freezing period, vertical water profile is approximately linear. During late winter, the vertical profile of water temperature presents a typical three-layer structure, namely, a conduction layer, a mixed layer, and a quiescent layer. The thickness of the three layers varies with time. The maximum diurnal difference of water temperature, amounting to 1.4 K, appears in the transition zone between the conduction layer and the mixed layer.
Based on the observed temperature data, the heat fluxes at the water-ice interface and the sediment-water interface were evaluated. The sediment-water heat flux and water-ice heat flux both remain in positive values throughout the freezing duration, and the average relative difference of heat flux calculated by using the temperature gradients on both sides of either interface is less than 7.9%. The average difference between the water-ice heat flux and the sediment-water heat flux is about 6~6.2 W/m 2 , while the averaged flux of solar radiation penetrating into the water is at least 12.48 W/m 2 , which is larger than 6~6.2 W/m 2 ; therefore, the water temperature continues to rise in the whole freezing duration. The results also demonstrate that the thermal structure of extremely shallow lakes follows the same seasonal pattern of water temperature dynamics as in lakes with a relatively larger depth.
A model solving effective thermal diffusivity based on temperature data was established, and the calculated effective thermal diffusivity is about 2.3-9.8 times the molecular thermal diffusivity, and the time-averaged values first increase and then decrease with the relative water depth, reaching a maximum at the relative depth of 0.5, which indicates that the water temperature of a mixed layer tends to be more uniformly distributed than that of a conduction layer and a quiescent layer.