Summary Study on Temperature Calculation Method of Water Accumulation in Permafrost Regions

1. Introduction

Water accumulation is a common engineering phenomenon in cold regions, and water accumulation gradually increases with climate changing and permafrost degenerating [1]. Research indicated water accumulation along subgrade sides had increased more than 10-fold in recent years [2]. The forming and development of water accumulation is usually caused by environment factors such as air temperature, geographical environment, soil and so on, and water accumulation can also change the local climate [3]. As a kind of exchange boundary between air and ground, the temperature change of water accumulation is different from permafrost, and water accumulation usually makes a thermal effect on permafrost and engineering in cold regions [4]. Moreover, water accumulation along embankment will also threaten the stability and bearing capacity of subgrade [5].

As a kind of water accumulation, thermokarst lake has been studied in permafrost [6,7,8]. Thermokarst lake always has a thermal impact on the permafrost and engineering around to increase and decline in the permafrost temperature and permafrost table, respectively [9,10]. Therefore, the thermokarst lake can make the permafrost degenerate and destroy the permafrost engineering to make it thaw settlement and bearing capacity decline [11,12,13,14]. The previous study show that the temperature of water is key factors to affect permafrost for thermokarst lake. Moreover, water depth is also an important factor to analyze the thermal effect of lake on permafrost, which can affect the lake bottom temperature. When the water depth is smaller than the frozen depth, temperature at the bottom of the water will be below 0℃in winter. While, when the water depth is greater than or equal to the freezing thickness, the temperature of the underlying soil will continue to be above 0℃ to make the permafrost thaw [15,16,17,18]. Therefore, the depth of the water and heat transfer in it cannot be neglected.

The recent study of water accumulation in cold regions is mainly focus on the thermokarst lake. It always makes a thermal analysis on the permafrost and embankment, but the heat transfer and temperature change in the lake being ignored. Therefore, this method is not applicant. The water should be as a boundary and consider the heat flux between the air and water and between the water and ground below in order to make a research about the thermal effect of pond in permafrost regions from the perspective of mechanism.

This paper reviews previous calculation methods of water temperature in water accumulation. Various kinds of analysis theory, calculation model and method of internal temperature in water accumulation are obtained, which provides the analysis scheme and basis for the calculation of water temperature in permafrost region, and then provides reference for the calculation and analysis of pond problems in projects.

2. Thermal surface boundary of water accumulation

The surface temperature of water accumulation is determined by the heat exchange of water surface layer, which affects the temperature distribution in the water accumulation. The surface of water in permafrost region is divided into unfrozen and frozen.

2.1. Unfrozen surface boundary

For unfrozen surface of water, surface temperature calculation has two kinds, one is using the fitting of measured temperature to get periodic function or empirical formula about water temperature at the surface over time [1]. The other is calculated by energy balance of the water at the surface according to meteorological elements including air temperature, wind speed and radiation, etc. [19].

Water surface temperature calculation used energy balance is always used in the calculation model of lake. The equation of heat balance proposed by Stepanenkoetal at the surface of pond based on meteorological data can be used to calculate the surface temperature of small lakes [19].

$$S\left(1-\alpha \right)+E-E-H-LE--\frac{\lambda }{h}\frac{\partial T}{\partial \xi }$$

(1)

Where,$S\left(1-\alpha \right)$ is solar radiation, $E_a$ is atmospheric long-wave radiation, $E_s$ is surface reflection, $H$ is sensible heat and $LE$ is latent heat, $-\frac{\lambda }{h}\frac{\partial T}{\partial \xi }$ is surface temperature considering the change of surface by evaporation.

When the change of surface by evaporation is ignored, the equation above can be expressed as follows.

$$S(1-\alpha )+E_a-E_s-H-LE=Q=-(k+K)\frac{\partial T}{\partial z}$$

(2)

Figure 1. Engergy exchange at the surface of water.

Figure 1. Engergy exchange at the surface of water.

Each heat value in the left side of the equation (2) can be obtained directly by observations and calculated by formula and different models of lake [20].

$$H_0=K(T_e-T)A_s$$

(3)

Where, $H_0$ is the heat received at surface of water in unit time, $K$ is the surface heat exchange coefficient, $T_e$is the equilibrium temperature; $T$ is surface water temperature; $A_s$ is the surface area of reservoir. $K$ and $T_e$ is expressed as:

$$K=15.7+\left(\beta +0.26\right)W_f$$

(4)

$$T_e=T_D+{\beta }_1{H}_s/K$$

(5)

Where, $\beta$ is the curve slope of water vapor pressure-temperature,$\beta =\left(E_s-E\right)/(T_A-T_D)$, $E_s$ and $E$ are saturation and actual vapor pressure respectively,

$T_A$and $T_D$ are air temperature and dew point temperature respectively; $W_f$ is the wind function related to evaporation; ${\beta }_1$ is the absorption coefficient of water surface to solar radiation; $H_s$ is the sun radiation reaching water surface.

The thermal balance at surface is also affected by flow conditions of water. Due to the density of water being a function of temperature as the Figure 2, according to the difference of water temperature, different density of water flows into different position. New energy is brought to the surface layer with the water flowing to the surface.

2.2. Frozen surface boundary

Ignoring the snow, there are two kinds of method to calculate thermal boundary at the frozen surface of water including formula and energy balance. Ashton [21] thought ice surface temperature is same as air temperature. Sydor [22] simulated ice growth through taking the measured water temperature under the ice as the thermal boundary.

In addition, more methods using energy balance method to calculate the frozen surface of pond are widely used.

Water is gradually frozen from the surface to lower as the surface reaches freezing temperature. And ice begins to melt, when the accumulated heat of ice reaches above 0℃. When water accumulation is completely frozen, heat is transferred from air and soil under water to ice. When water accumulation is not completely frozen, heat is transferred from high temperature and water to ice. So there is usually freezing in one direction from up to down and melting in two direction from up and down to middle position in permafrost regions [23].

Heat exchanges between ice-gas and ice-water is important during the process of freezing and melting. These heat exchanges above and the change of water temperature are affected by meteorological factors such as accumulated negative temperature, sunshine and wind speed, and controlled by convection and phase transition of water [24]. So icing time is related with accumulated air temperature, freezing time, and the area and depth of water. At the early stage of ice growth, the heat exchange between ice and air temperature is relatively large and the growth rate is fast. As ice layer becomes thicker, the heat exchange between air and ice are decreased, so the growth rate is stable. The thickness of ice is related to the average daily temperature and freezing time in cold season. Burn obtained the relationship between ice thickness and mean temperature in winter is linear by the measured data of Todd Lake in Canada [25].

3. Water temperature

3.1. Empirical law

The empirical methods of water temperature calculation can be divided into two types, which is the calculation method based on the boundary temperature of water and the vertical average temperature of water. The vertical temperature distribution of water is shown as Figure 3. Zhang calculated the vertical water temperature each month, according to the water temperature distribution of reservoir surface and bottom [26].

$$T_y=(T_0-T_b){e^{(-y/x)}}^n+T_b+T_y+T_0+T_b$$

(6)

Where, $T_y$ is the vertical water temperature related to depth y, $T_0$ is the reservoir surface temperature, $T_b$is the reservoir bottom temperature, x, n is a function related to month.

Chen [27] observed that the water temperature of reservoir mainly changed in the vertical direction, especially in summer, the temperature difference in vertical can be more than 10℃ based on the temperature measurement data of Miyun Reservoir in Beijing for two years.

But the empirical formula of water temperature obtained from measured data has particularity, ignoring the factors such as water flow, local meteorological conditions, wind mixing and so on. However, the referring factors can be taken into mathematical model by theoretical way, since its calculation results is more accurately, the mathematical model method may be more practical and universal.

3.2. Mathematical model method

The research on mathematical model of water accumulation temperature originated in 1930s. At that time, a large number of measured data of reservoir water temperature were analyzed in USA, and the stratification of water temperature was mainly reflected in the vertical direction. The diffusion model based on convection-diffusion is based on this assumption, the longitudinal and transverse temperature changes are ignored, and vertical water temperature stratification is considered, and all the elements of inflow in the same water layer are basically the same. Based on the convection-diffusion theory, WRE and MIT one-dimensional mathematical models were put forward successively in the 1960s. In the model, it is assumed that the reservoir is composed of mixed homogeneous horizontal isothermal thin layers, and the heat transfer only occurs between the adjacent horizontal thin layers, and the external flow flows into the water layer unit with the same pre support density [29,30].

The energy equation is

$$\frac{\partial T}{\partial t}+\frac{\partial }{\partial \mathrm{z}}\left(\frac{T{Q}_V}{A}\right)=\frac{1}{A}\frac{\partial }{\partial z}\left(A{D}_Z\frac{\partial T}{\partial z}\right)+\frac{B}{A}\left(u_i{T}_i-u_{0}T\right)+\frac{1}{\rho AC}\frac{\partial \left(A{\phi }_Z\right)}{\partial z}$$

(7)

Where, T is the water temperature of unit layer $\left(\mathrm{℃}\right)$;$T_i$ is the inflow temperature $\left(\mathrm{℃}\right)$; A is the unit level area $\left(m^2\right)$; B is the average width of the unit layer $\left(m\right)$; $D_Z$ is the vertical temperature diffusivity $\left(m^2/s\right)$; $\mathsf{\rho }$ is the density of water $\left(kg/m^3\right)$, is a function of temperature; C is the Specific Heat Capacity of water $\left[J/\left(kg\mathrm{℃}\right)\right]$; ${\phi }_Z$ is the solar radiation flux $\left(W/m^2\right)$; $u_i$ is the inflow velocity $\left(m/s\right)$;$u_0$ is the outflow velocity $\left(m/s\right)$;$Q_V$ is the vertical flow through the upper boundary of the element $\left(m^3/s\right)$; N is the surface unit; $V_N$ is the volume of surface unit $\left(m^3\right)$; $Q_{V,N-1}$ is the vertical flow through the interface between layer N and layer N-1 $\left(m^3/s\right)$;${\phi }_N$ is the heat absorbed by the surface water through the water air interface $\left(W/m^2\right)$.

In the middle and late of 1970s, the distribution of vertical water temperature was proposed considering the inflow and outflow of reservoir, and the influence of wind based on the total energy. In 1975, Stefan and Ford proposed MLTM model based on this theory [31]. In 1978, Imberger proposed the DYRESM model [32] suitable for simulating the temperature and salinity changes of reservoirs, and introduced the concept of mixed layer for the first time, which can simulate the distribution of reservoir inflow and outflow, water temperature and salinity at the same time. Since 1980, the mixed layer model has been widely used and improved because of its good characteristics of simulating water temperature and water quality. Some models have been applied to the study of water temperature in lakes and other ponding types. Although the mixed layer model increases the migration of turbulent kinetic energy compared with the previous diffusion model, the study of diffusion in the lower layer of ponding is fuzzy and has poor universality. It can be seen from the equation that the key to the water temperature distribution is the selection of vertical diffusion coefficient. Therefore, the improvement of water temperature numerical model is embodied in the treatment of vertical diffusion coefficient. Some models set the vertical diffusion coefficient as constant according to the focus of research problems, and more empirical formulas of vertical diffusion coefficient have been proposed, and various improved water temperature models have been applied.

With the gradual development of one-dimensional water temperature model, researchers have successively proposed two-dimensional or even three-dimensional water temperature model. The two-dimensional model integrates the input in term of N-S equation transversely and obtains the longitudinal and vertical water temperature distribution. In 1975, Edinger developed the earliest two-dimensional model LARM [33]. After that, two-dimensional water temperature models are proposed. Based on the energy balance and the vertical one-dimensional hybrid model. Farrell and Heinz [34] introduced k-ᵋ model into the simulation of reservoir density flow (k-ᵋ model is a mathematical treatment of water momentum equation), which can well simulate the characteristics of reservoir density flow, such as submergence, vertical vortex and temperature stratification. Young [35] proposed a two-dimensional model of turbulent diffusion of momentum and heat under different vertical empirical formulas. By using numerical software, a universal multi-dimensional prediction model of reservoir water temperature was established, which considered the factors of water surface heat exchange, reservoir inflow, reservoir discharge, reservoir water level and storage capacity change greatly. The three-dimensional model is based on the turbulence model, which considers the temperature changes in the transverse, longitudinal and vertical directions at the same time, and the heat transport equation and turbulence equation are coupled to solve. Due to the relatively large workload of 3D model calculation, it is not widely used in practice. By comparing the simulated and measured results of water temperature in the water accumulation, it is concluded that the variation of water accumulation temperature is related to its internal diffusion coefficient, and indirectly affected by the density of water under the action of temperature. In winter, the bottom temperature of deep water is high, and the heat is transferred to the upper part; in summer, the heat is transferred from the surface to the bottom, and the temperature from the surface to the bottom shows a linear or nonlinear trend. In the process of melting and freezing, the internal temperature stratification is obvious.

4. Phase change in pond

4.1. Empirical law

The empirical method is based on the measured data, which is from field observation and indoor test data. The method was improved continuously and according to the theory of water flow, thermal conduction, and the ice forming [36]. From the measured reservoir temperature data, it showed that the ice temperature near the ice surface is the same as the total radiation on the ice surface, and the closer the water temperature is to the ice bottom, the faster the water temperature decreases [37]. Temperature of static water condition is the main factor affecting ice condition, ice thickness is linearly related to accumulated hourly negative temperature. And water depth affects the start time of initial ice formation, but has little effect on the thermal change of ice. In the early stage of ice formation, the surface water body loses heat quickly and falls to freezing point to form ice needle ice. During freezing period, the bottom layer of ice layer crystallizes continuously under the heat conduction in the ice and the heat balance in the water to realize the of ice growth [38]. In the growth period, the temperature in the ice layer is basically linear distribution, and the slope of ice temperature increases with the decrease of temperature. The distribution of water temperature under ice during melting ice is basically unchanged, and the vertical temperature difference is mainly concentrated in the surface water body. It is shown that the water body is in a static state in the indoor experiment, and the heat exchange between the ice surface and the atmosphere is mainly affected by the temperature and relative humidity, and the icing only occurs at the interface of ice water. In view of the limitations of field measurement and model test, it can not be used directly and widely, and the water temperature calculation is nonlinear during phase transition. More mathematical models of water phase transition icing have been proposed one after another.

4.2. Mathematical model method

In the ice mathematical model, the development of river ice is better than that of reservoir ice [39]. The first one-dimensional vertical Chen-Orlob water temperature model considering icing was proposed in Carlson 1977, However, due to the lack of stability, it is difficult to simulate the of icing and melting process [40]. CE-QUAL-R1 model was proposed in Ashton 1982 to simulate the water temperature distribution under the ice layer [41]. In 1984, one dimensional unsteady river ice model was proposed, this model takes into account the variation of ice and density along the process [42,43]. The existing 1 and 2 dimensional ice mathematical model has developed relatively mature, Classic RICE models. Lal and Shen [44] improved and validated and applied in practice, Shen also proposed to replace the critical velocity in the RICE model with ice transport capacity, RECEN model for simulating ice transport and accumulation in ice cap. The temperature distribution models of the water layer and the lower soil layer after freezing under the heat exchange of each meteorological factor and the bottom of the water accumulation were established [45]. The empirical formulas and parameters in the above model were refined in Fang and Stefan in 1994 [46]. Fang [47] assumed that each water layer produces a heat exchange with the bottom of the water in 1996, the water temperature and ice growth were simulated for a long time. With the development of one-dimensional model, the two-dimensional icing model has been developed accordingly. Shen [48] proposed two-dimensional river ice mathematical models considering three modules: hydrodynamic force, ice transport and ice mass conservation. Although the research of river ice model has made great progress, but because of the complexity of the mechanics of the river itself, in the process of calculating temperature, various physical mechanism processes, such as hydrodynamic, thermal and phase transition, Therefore, there will still be more problems in the development of river ice model. The key to simulate the icing of rivers and lakes lies in the consideration and treatment of ice water and ice gas interface.

For the study of reservoir ice condition, it is still in the initial stage. At present, in the reservoir ice model, the Stefan classical assumption is used to determine the ice surface temperature. It is considered that the ice surface temperature is equal to the air temperature and the water temperature under the ice is linearly assumed [49]. In addition to meteorological conditions, the reservoir capacity and the flow into the reservoir are also the main factors affecting the water temperature of the reservoir. Hence, the ice condition in the reservoir is affected by many factors, such as meteorological conditions, reservoir morphology, operation mode and hydrodynamic conditions, and the icing mechanism is relatively complex [49]. Reservoir surface ice is mostly formed in winter night temperature difference, the surface is less than 0℃, ice formation begins on the shore, as the temperature decreases, more ice is gradually formed at the deep level of reservoir. Xing Fang [50] considering the influence of wind speed, a lake ice model was established to calculate the ice thickness by using the heat lost from water cooling.

5. Heat exchange between the bottom of water accumulation and the underlying soil

The research and development of water temperature calculation model shows that the heat situation in water accumulation depends on the energy exchange between the water and the environment. Therefore, the water, the upper atmosphere and a part of the adjacent soil layer should be taken to calculate the complete calculation model. Especially for the shallow water, the influence of the heat exchange of the bottom soil layer on the heat distribution in water accumulation can not be ignored. In the previous calculation methods of water temperature distribution in lakes and reservoirs, the thermal effect of soil layer was also considered. In 1972, the calculation model of temperature distribution in water and underwater sediment by finite difference method was established [51]. The thermal effects between the water accumulation and the soil layer are also mutual, especially when studying the influence of water accumulation on nearby projects, the soil layer at the water bottom should be considered in calculation and analysis. Stepanenko [52] used the plant soil model [53] to analyze the heat and water exchange between the bottom of the ponding and the soil. Xing Fang [54] proposed a modified ponding model considering the heat exchange between the ponding and the bottom sediments, indicating that the change of bottom water temperature can affect the depth of 10 meters.

Water depth and water contact of soil affect the degree of water temperature affected by external changes and the heat exchange between the bottom and soil. The measured bottom water temperature of lake shows that the lake edge temperature is lower than that of the inner deep water in the cold season, while the temperature at the lake edge is higher than that of the inner deep water in the warm season. The heat exchange between the lake bottom and the soil layer is greater than that of the internal deep water. The bottom temperature varies in the horizontal direction due to the influence of water depth and heat exchange between water and soil [55]. In the alpine environment, taking permafrost region as an example, the relationship between water depth and freezing depth will also affect the lower soil layer affected by the heat of ponding, and the temperature of the bottom layer of ponding determines the development of its lower melting zone [56,57]. When the depth of water is less than the freezing depth, the water bottom temperature is less than 0 ℃ [58]. However, when the water depth is greater than or equal to the freezing depth, the water accumulation in winter is not completely frozen or the bottom of the cold season is frozen, and the temperature is still higher than 0℃. The long-term effects of the above two scenarios will drive the bottom soil layer of water accumulation to gradually form a melting layer [59]. If the accumulated water does not exist for a long time and the temperature at the bottom does not change recently, the lower soil layer and the bottom of the ponding reach the thermal balance, and the lower melting layer will not develop any more [60,61,62]. However, if the accumulated water exists for a long time or the bottom temperature changes obviously due to the external temperature, such water will have a great impact on the development of the surrounding permafrost. By measuring and simulating the heat budget at different positions of the water accumulation, it is found that the lower part of the water has formed a through hot melt zone, and the surrounding frozen soil has degenerated [59]. The results show that the long-term existence of the thermokrast lake in cold region can be regarded as the heat source of the permafrost around it [63]. Therefore, the heat exchange between water and soil should be considered in the water temperature calculation model, and the heat exchange between water and soil can not be ignored when studying the heat effect of water accumulation, and the heat impact caused by water is also affected by the external environment, geological conditions, water accumulation type and other factors.

6. Discussion

There is vertical temperature stratification in the water accumulation, which can be divided into mixed type, stable type and transition type between them because of the degree of stratification. For the mixed type, water temperature gradient in vertical direction is small and the bottom temperature varies little with surface temperature. For the stable type, there are layers of variable temperature, oblique temperature and average temperature from the upper to the lower, and the surface water temperature changes faster than the middle and lower layers. The transition type has the temperature characteristic of the two types above. When the surface is frozen, the ice temperature is distributed with temperature gradient, and the water layer under ice still accords with the distribution characteristics of above water body.

The main factors of water temperature stratification are surface absorption heat, wind speed, solar short wave radiation, heat conduction and buoyancy convection in the water accumulation, which are determined by the environment, climate, depth, scale and thermal parameters of the pond. The surface absorption heat is the source power of the water accumulation. The wind speed can affect the kinetic energy of surface layer diffusing heat to lower layer. The solar short wave radiation absorbed in every layer affects the internal energy exchange of water. Besides that, the density of water is a nonlinear function of temperature, and the density of about 4℃ is the largest. When there is an unstable temperature stratification in the water accumulation, the difference of density can form natural buoyancy flow, which causes the corresponding water flow and promotes the temperature change in the pond. The depth of water not only changes the influence of solar radiation, but also changes the effect degree by the surrounding soil. Shallow water is more affected by solar radiation and surrounding soil. The area of water accumulation can affect the flow state of water and the heat transfer ability of water and the degree of natural buoyancy convection. Moreover, the water temperature distribution is also affected by the inlet and outlet water temperature and velocity. The external forced inflow will cause additional energy and momentum exchange in the pond. Then the redistribution of internal energy of water accumulation is affected to change the state of water temperature distribution.

When the water is frozen, the latent heat of phase change between ice and water and the movement of interface will affect the distribution of water temperature. The ice layer is a simple solid heat transfer and the attenuation coefficient of solar short wave radiation is different from that of water layer. Therefore, the influence factors considered in the temperature distribution of water accumulation with icing process include ice melting period, ice thickness and thermodynamic parameters of ice. The energy transfer process and law after freezing should also consider the influence of the depth and area of the water accumulation. The shallow water is greatly affected by the heat exchange between the soil and water, and the water accumulation area will affect the icing sequence of the surface layer.

The water accumulation in permafrost region also involves phase transformation. So these factors above which can affect the water temperature change are also be considered in permafrost regions, in order to make the water temperature calculation more accurate. The key to the calculation of water accumulation in cold regions lies in Figure 4. They are including that the energy exchange at the interface of water and air, including the energy exchange at the interface of water and air, the heat transfer in water accumulation and the energy exchange between soil and water bottom. For the interface of water and air, it is important to calculate the surface energy balance of water considering the meteorological elements above the surface, and there is difference between the surface of water and ice. For the heat transfer in water accumulation, it is important to make an analysis of water temperature stratification taking the actual buoyancy flow into account and the solution of melting and freezing of ice. The heat transfer in ice layer is solid and that in water is fluid, and ice always floats on the surface of water. In addition, the solar short wave radiation will exist in each water layer as a heat source in the heat exchange inside the water accumulation. The heat exchange between water accumulation bottom and permafrost below pond is key to calculate water temperature.

7. Conclusion

As a kind of boundary of the earth, water accumulation changes the energy exchange between the earth and the atmosphere, affects the changes of the surrounding environment, and indirectly affects the stability of adjacent projects. In order to deeply analyze its impact on local environment and engineering, it is necessary to define the response of water accumulation to climatic conditions, that is, to calculate and analyze the internal temperature of water accumulation with time and changes of external conditions. Taking lakes and reservoirs as the research objects, this paper summarizes the methods, mechanisms and application conditions of surface and internal temperature calculation and analysis of water accumulation. However, for the special soil environment in permafrost regions, the calculation and analysis of water heat and its impact on frozen soil engineering are still in the initial stage, and there are many problems to be studied in deep.

1. It is more applicable to calculate the water surface temperature by energy balance method, which includes solar short wave radiation, net long wave radiation, sensible heat and latent heat. The mechanism and influence degree of various factors on the surface temperature of water accumulation can be discussed in depth, which is conducive to accurate prediction of water temperature.

2. The depth and area of water accumulation will affect the water temperature distribution by affecting the internal flow. Therefore, it is necessary to analyze the variation law of water temperature and its influence on the thermal state of the lower soil layer from the perspective of water accumulation scale.

3. The existence of water accumulation in special areas not only affects the environment, but also affects the stability of surrounding engineering structures. These impacts are also related to the condition and existence time of water accumulation. Therefore, it is necessary to analyze the impact of water accumulation on the underlying soil based on the response of water accumulation to climate under different existing time.