This paper presents a general approach to solving geophysical problems related to solar heating and the melting of ice on the water surface. The spring melting of lake ice is a relatively simple problem and this will be considered in more detail. After that, a general view of the more complex problem of melting of sea ice during the polar summer will be discussed. Ice melting on the Arctic Sea surface is really important for global climate change. What is less well known is that the opening from the ice of high-mountain lakes in Tibet, sometimes referred to as the Earth’s third pole, significantly affects the climate of not only central Asia. Recall that the Tibetan Plateau is home to about a thousand large lakes with a total area of about 15,000 km² (Su et al., 2020; Zhang and Duan, 2021). Interestingly, the early opening time of these lakes is an indicator of global warming (Su et al., 2019; Zhang and Duan, 2021).

It should be noted that the physical picture of heating and melting of the ice on the lake surface turns out to be quite different for the case when the ice is not covered with snow and in the case when a snow cover is present. However, the approach to solving problems of radiative transfer in snow and ice layers has a common methodological basis, and this can be used in theoretical studies of both lake and sea ice melting. In the spectral range of semitransparency of snow and ice, one should take into account the scattering of radiation either by ice grains in snow or by gas bubbles, which are usually contained in the ice. This means that one should focus on choosing a simple and sufficiently accurate method for solving the radiative transfer equation (RTE) (Dombrovsky and Baillis, 2010; Howell et al., 2021; Modest and Mazumder, 2021).

The main difficulties in radiative transfer modeling are caused by the RTE’s integral term, which takes into account the anisotropic scattering of radiation by particles or bubbles in the medium and contains a scattering phase function. Fortunately, we deal with multiple scattering of radiation in an optically thick medium when the transport approximation is quite sufficient. According to this approach, the scattering function is replaced by a sum of the isotropic and forward components. The resulting transport RTE looks like that for the hypothetic isotropic scattering but with the transport scattering coefficient (Dombrovsky, 2012; Dombrovsky, 2019). It should be recalled that the transport approximation has been successfully employed in diverse problems of radiative transfer in thermal engineering (Dombrovsky et al., 2017; Dombrovsky et al., 2020), biomedicine (Tuchin, 2007; Sandell and Zhu, 2011; Dombrovsky et al., 2012; Jacques, 2013; Eisel et al., 2018; Dombrovsky, 2022), and geophysics (Dombrovsky and Kokhanovsky, 2022; Dombrovsky and Kokhanovsky, 2023; Dombrovsky et al., 2019; Dombrovsky et al., 2022).

The linearity of the RTE makes possible another significant simplification of the problem. In the case of direct solar irradiation of a scattering medium, the radiation intensity can be represented in the form of two additive components: direct radiation and a diffuse component formed due to the scattering. In highly scattering media, the diffuse component of radiation intensity can be calculated using one of the simplest differential approximations: either the known P ₁ approximation of the spherical harmonics method or the two-flux approximation. These approximations make it sufficient to solve a boundary-value problem for a second-order ordinary differential equation instead of the RTE. The choice between these approximations is determined by the problem statement. Comparison with exact numerical solutions (Dombrovsky and Baillis, 2010) has shown that the two-flux method is preferable when solving one-dimensional problems typical of solar heating. In this case, the error of P ₁ is larger because this method does not take into account the discontinuity in the angular dependence of the radiation intensity on the illuminated surface of the medium.

Generally speaking, there are various thermal processes in a snowpack or scattering ice sheet and these processes should be involved in a complete heat transfer model. As an example, one can recall the ice sublimation and diffusion of water vapor through a snow layer. This may be important for the snow microstructure but the related effects are insignificant in our problem.

The 1D transient energy equation for temperature T t , z in a layer of snow or ice and the accompanying initial and boundary conditions can be written as follows: ρ c ∂ T ∂ t = ∂ ∂ z k ∂ T ∂ z + P , t > 0 , 0 < z < d th , (1a) T 0 , z = T init z , (1b) z = 0 , k ∂ T ∂ z = h T air − T − q rc , z = d th , ∂ T ∂ z = 0 , (1c) where z is the is the coordinate measured from the illuminated surface, ρ , c , and k are the density, the specific heat capacity, and the thermal conductivity of the medium, T air t and h t are the temperature of ambient air and convective heat transfer coefficient which depends on the wind velocity. The adiabatic condition at z = d th should be replaced by the condition of the first kind when the temperature of the lower interface is known (as in the case of the ice-water interface). The heat flux q rc is the mid-infrared cooling due to thermal radiation of the snow or ice surface in the atmospheric transparency window of λ w 1 < λ < λ w 2 ( λ w 1 = 8  µm, λ w 2 = 13  µm) (Raman et al., 2014; Chen et al., 2016; Hossain and Gu, 2016). Of course, the latent heat of ice melting, L = 0.34  MJ/kg, should be taken into account in the calculations. This can be done using an equivalent additional heat capacity in a narrow temperature interval near the melting temperature as it was done by Dombrovsky et al. (2019). This simple technique is not new and was successfully used by the author in the numerical solution of other problems to account for the latent heat of phase change of the first kind, including the general case when the so-called mushy zone is formed between solidus and liquidus (Dombrovsky et al., 2015; Roy et al., 2023).

It is usually difficult to choose a realistic initial profile of temperature, T init z , for the heat transfer calculations. Fortunately, the effect of this temperature profile decreases with time. As an example, it was shown by Dombrovsky et al. (2019) that the choice of this temperature profile makes no difference for the snow layer after about 4 h from the initial time moment.

In the case of ice not covered by snow, sunlight penetrates through the ice layer and leads to water heating. In limnology, the transfer of solar radiation in ice and water is described using the exponential Bouguer law, which is applicable only in the case of direct incident radiation and single scattering of light in the medium. This methodological drawback remained even in recent papers (Kirillin et al., 2012; Leppäranta, 2015; Kirillin et al., 2021). The error of this simplified model for the radiative transfer was partially compensated for by the selection of an extinction coefficient, which provided a satisfactory agreement between the calculations and the field measurements. This semiempirical approach to radiative transfer is incorrect. The multiple scattering of short-wave solar radiation by microcracks and gas bubbles leads to stronger solar heating of ice. In addition, the radiation transmitted through the ice layer contains both direct and diffuse components.

According to (Hale and Querry, 1973), the absorption index of water increases with the wavelength from κ w = 10 − 9 − 2 × 10 − 9 at λ = 0.4 − 0.5  µm to κ w = 10 − 5 at λ = 1.2  µm. The corresponding water absorption coefficient, defined as α λ w = 4 π κ w / λ , changes from α λ w = 0.025 − 0.063  m⁻¹ to α λ w = 100  m⁻¹. The latter means that visible radiation penetrates water to a depth of tens of meters, while infrared radiation with wavelength λ = 1.2  µm penetrates to a depth of less than 1 cm. In the wavelength range 0.6 < λ < 1.2  µm, the absorption index of ice differs slightly from that of water (Warren and Brandt, 2008), and a similar situation occurs for ice. The radiation with wavelength λ > λ * = 1.2  µm is absorbed in a relatively thin surface layer of ice. That is why the threshold value of λ * is considered as a conventional boundary of the range of semitransparency.

The influence of uncertainty in the experimental values of the ice absorption index on the computational data for solar heating of snow and ice has been analyzed by Dombrovsky and Kokhanovsky (2022), Dombrovsky and Kokhanovsky (2023). It was shown that the discrepancy between the data of (Warren and Brandt, 2008; Picard et al., 2016) in the λ < 0.6  μm range does not significantly affect the results of calculations. Therefore, the more complete spectral data of Warren and Brandt (2008) are used in the present work.

The gas bubbles in ice are assumed to be spherical and the radii of these bubbles, a , are much greater than the wavelength. The corresponding diffraction parameter x = 2 π a / λ is very large, and therefore the geometrical optics approximation can be used instead of the rigorous Mie theory (Bohren and Huffman, 1998). As a result, the absorption coefficient of ice containing bubbles can be calculated as follows (Dombrovsky, 2004): α λ = 1 − f v α λ 0 , (2) where f v is the volume fraction of bubbles, α λ 0 = 4 π κ ice / λ is the absorption coefficient of ice without bubbles, and κ ice is the index of absorption of ice. In most real cases, f v ≪ 1 and α λ ≈ α λ 0 . Not only the gas bubbles but also the average distance between them is very large compared to the wavelength of the radiation. Therefore, the hypothesis of independent scattering by single bubbles is true (Mishchenko, 2014; 2018).

According to (Kirillin et al., 2012), the gas bubbles are usually not uniformly distributed in the ice layer: there are more of them in the lower part of the layer. This can be taken into account in the calculations. However, the simplest model including only two bubble parameters is chosen in the present study: the average values of volume fraction and radius of bubbles. According to (Dombrovsky and Kokhanovsky, 2020b), the scattering parameter S = f v / a 32 , where a 32 is the Sauter’s mean radius of bubbles, is sufficient to describe light scattering by polydisperse bubbles. The transport scattering coefficient of ice with gas bubbles can be calculated as follows (Dombrovsky, 2004): σ λ tr = 0.675   n ice λ − 1   S , (3) where n ice is the refraction index of ice. The value σ λ tr is usually greater than the absorption coefficient in the entire spectral range under consideration. As a result, the extinction of sunlight in an ice layer is determined by the light scattering.

Following the recent study by Dombrovsky and Kokhanovsky (2023), the present paper is focused primarily on modeling the thermal regime of mountain lakes like those of the Qinghai–Tibet Plateau, located at a height of more than 4 km above sea level, for which the problem is somewhat easier: it can be assumed that the sky is clear and only a small part of the solar light is scattered in the atmosphere. As a result, it is sufficient to consider the transfer of direct solar radiation, the intensity and spectral composition of which change with the zenith angle of the Sun. Minor atmospheric precipitation and strong winds (Wang et al., 2022) result in the absence of snow cover on the ice surface. This leads to a significant solar heating of the ice and water in the lake.

The problem can be simplified because the refractive indices of water and ice differ only slightly in the wavelength range of λ uv ‐ vis < λ < λ * ( λ uv ‐ vis ≈ 0.4  µm). As a result, the reflection and refraction of light at the ice-water interface are small and can be neglected. The effect of solar light reflection at the ice-water interface is even smaller. It can be estimated from the Fresnel formulas for the reflection coefficient at the interface of transparent media with different refractive index values. At wavelength λ = 0.6  µm the refractive indices of ice and water are 1.31 and 1.33, respectively, and the reflection coefficient of light incident at an angle of 45^o is equal to 1.2%, and when illuminated along the normal–0.06%. The refraction of light at the ice-water interface is also negligibly small.

It is assumed that optical properties of ice do not change along the horizontal ice surface and the radiative transfer along the surface of an ice-covered lake may not be considered. The 1D model for the propagation of obliquely incident radiation in an ice layer containing gas bubbles has been described recently by Dombrovsky and Kokhanovsky (2022). It is assumed that there are no any bubbles in water, and scattering by a small amount of plankton is insignificant. In this case, the problem of radiative transfer in the ice layer can be solved independently of the propagation of light in water. It is convenient to write the transport RTE and the boundary conditions at an oblique illumination of the ice layer with thickness d in dimensionless variables (Dombrovsky and Kokhanovsky, 2023): μ ∂ I ¯ λ ∂ τ λ tr + I ¯ λ = ω λ tr 2 G ¯ λ , G ¯ λ τ λ tr = ∫ − 1 1 I ¯ λ τ λ tr , μ   d μ , (4a) I ¯ λ 0 , μ = r λ I ¯ λ 0 , − μ + 1 − r λ   δ μ j − μ , I ¯ λ τ λ , 0 tr , − μ = 0 , μ , μ j > 0 , (4b) where I ¯ λ = I λ / I λ inc , G ¯ λ = G λ / I λ inc , I λ inc is the spectral intensity of incident radiation in direction μ i = ⁡ cos ⁡ θ i ( θ i is measured from the external normal), r λ = r λ n ice , μ i is the Fresnel reflection coefficient, μ j = 1 − 1 − μ i 2 / n ice 2 is the cosine of the refraction angle, τ λ tr z = ∫ 0 z β λ tr z   d z is the optical depth and τ λ , 0 tr = τ λ tr d .

The intensity of radiation and the spectral irradiation are presented as follows: I ¯ λ = J ¯ λ + 1 − r λ E λ j δ μ j − μ , E λ j = ⁡ exp − τ λ tr / μ j , (5a) G ¯ λ = G ¯ λ dif + 1 − r λ E λ j , G ¯ λ dif = ∫ − 1 1 J ¯ λ   d μ . (5b)

The resulting equations for the diffuse component J ¯ λ are: μ ∂ J ¯ λ ∂ τ λ tr + J ¯ λ = ω λ tr 2 G ¯ λ , (6a) J ¯ λ 0 , μ = r λ J ¯ λ 0 , − μ , J ¯ λ d , − μ = 0 ,   μ > 0 . (6b)

The two-flux method gives the following boundary-value problem for the irradiation G ¯ λ , i dif τ λ tr : − ( G ¯ λ dif ) ″ + ξ λ 2 G ¯ λ dif = 4 ω λ tr 1 − r λ , i E λ j ,   ξ λ = 1 − ω λ tr , (7a) τ λ tr = 0 , ( G ¯ λ dif ) ′ = 2 γ G ¯ λ dif , τ λ tr = τ λ , 0 tr , ( G ¯ λ dif ) ′ = − 2 G ¯ λ dif , (7b) where r ¯ λ is the angle-averaged reflection coefficient of the ice layer surface and γ = 1 − r ¯ λ / 1 + r ¯ λ . This problem statement is correct for any variation of the optical properties across the ice layer and can be solved numerically (Dombrovsky and Baillis, 2010). In the case of uniform optical properties of ice, one can obtain an analytical solution to the problem.

The calculations made after publication of paper by Dombrovsky and Kokhanovsky (2023) showed that the reflection of sunlight from the illuminated ice surface is negligible and can be disregarded without any significant loss of calculation accuracy. At the same time, the refraction of light at an oblique illumination of ice should be taken into account and μ j > μ i . Therefore, below is the solution of the problem at r ¯ λ = 0 ( γ = 1 ). This simplified analytical solution at ξ λ ≠ ν j = 1 / μ j is as follows: G ¯ λ dif = 4 ω λ tr ξ λ 2 − ν j 2 E λ j + 2 − ν j 2 − ξ λ A E λ dif + B / E λ dif E λ , 0 dif 2 − ζ λ 2 ,   ζ λ = 2 + ξ λ 2 − ξ λ , (8a) A = ψ ζ λ − E λ , 0 j E λ , 0 dif , B = E λ , 0 dif ζ λ E λ , 0 j − ψ E λ , 0 dif , ψ = 2 + ν j 2 − ν j , (8b) E λ , 0 j = exp − ν j τ λ , 0 tr , E λ dif = ⁡ exp − ξ λ τ λ tr , E λ , 0 dif = ⁡ exp − ξ λ τ λ , 0 tr . (8c)

The relation between G ¯ λ and G ¯ λ dif also simplifies when r λ = 0 : G ¯ λ = G ¯ λ dif + E λ j .

The radiative transfer problem in water is as follows: μ ∂ I ¯ λ w ∂ τ λ w + I ¯ λ w = 0 ,   τ λ w = α λ w × z − d , z > d , (9a) I ¯ λ w d , μ = I ¯ λ d , μ ,   I ¯ λ w ∞ , μ = I ¯ λ w ∞ , − μ = 0 , μ > 0 . (9b) where α λ w = 4 π κ w / λ is the absorption coefficient of water. The solution to this problem is: I ¯ λ w z , μ = I ¯ λ d , μ exp − τ λ w / μ  when  μ > 0   0  when  μ < 0 . (10)

According to Eq. 10, the intensity of light in water at μ = 1 decreases most slowly. As a result, with increasing depth the light becomes less diffuse and closer to directed vertically downward.

The solar radiation power absorbed in ice and water can be calculated as follows: P z = ∫ λ uv ‐ vis λ * p λ , z d λ ,   p λ , z = α λ G λ z  when  z ≤ d α λ w G λ w z  when  z > d . (11)

The functions in the right-hand side of the second of these equations are defined as: G λ z = I λ , i inc × G ¯ λ dif z + ⁡ exp − ν j β λ tr z } , (12a) G λ w z = I λ , i inc × { G ¯ λ dif d exp − 2 τ λ w + ⁡ exp − ν j β λ tr d + τ λ w } . (12b)

Obviously, G λ w d = G λ d , whereas the function p λ , z is not continuous at z = d . The diffuse irradiation component in water (the first term in Eq. 12b) decreases with depth twice as fast as the directional one. It is clear that Bouguer’s law with a constant extinction coefficient in the exponent, which was applied for the light propagation in water under an ice layer in limnology studies, is not correct.

The spectral radiative flux at the surface of Ngoring Lake at the end of March (the time of the beginning of ice melting) at different values of the Sun’s zenith angle has been calculated in (Dombrovsky and Kokhanovsky, 2023). The results obtained are used below in calculating profiles of absorbed radiation power in ice and water. Typical profiles of this radiation power at the zenith angle 60^o for two values of the scattering parameter and different thicknesses of the ice layer are shown in Figure 1. A significant attenuation of light in the ice layer draws attention. Increasing the scattering parameter leads to an increase in the absorption of solar radiation in the upper part of the ice layer and a significant decrease in the absorption in the lake water. When ice melts and the thickness of the ice becomes smaller, the radiation absorption in ice decreases considerably and the absorption in water increases very strongly. The latter is the physical cause of an observation by Lazhu et al. (2021) for several lakes in the Tibetan Plateau where the water temperature at some distance from the ice-water interface increased rapidly during ice melting.

The temperature field in the lake water is largely determined by the non-monotonic dependence of the water density on temperature with a maximum density ρ w * = 1000  kg/m³ at temperature T * = 4 ℃ . At the beginning of ice formation, the water temperature at the ice-water interface becomes equal to T 0 = 0 ℃ , while the water at depth may remain warmer. However, even at the bottom of a deep lake, the water temperature cannot be higher than T * . In a severe winter, even a deep lake can freeze to the bottom (as in the case of Ngoring Lake), but this paper considers early spring, when the ice thickness does not exceed 1 m.

Field observations of Ngoring Lake and some other lakes in the Qinghai–Tibet Plateau have shown an interesting thermal regime, named “anomalous winter” by Kirillin et al. (2021), when at a depth of 1.5–3 m under the ice layer the water is heated by solar radiation to a temperature T max > T * . Solar radiation absorbed in the upper water layer leads to an increase in T max , and this effect was measured by Kirillin et al. (2021) in March and early April. An estimate of the heat flux from relatively warm water to the ice layer in (Dombrovsky and Kokhanovsky, 2023) gave the value q w ‐ ice ≈ 1.2  kW/m². The effect of this heat flux is much smaller than ice heating by solar radiation.

The heat transfer on the illuminated ice surface changes significantly during the day but this does not affect the temperature of ice at some distance from the surface because of the large heat capacity of the ice layer. A simple estimate based on the Fourier criterion confirms that a thermal relaxation time for the 0.5 m thick ice is about 10 days. This allows using a steady-state model with constant heat transfer parameters, varying according to weather changes from week to week.

The ice melting at Ngoring Lake in March is very slow, and only in April does ice melting accelerate, completed by April 16 (Kirillin et al., 2021). This result is clear from the observations by Zhou et al. (2022): the air temperature in March is almost constant and does not exceed − 10 ℃ , while it increases up to − 3 ℃ in the middle of April at the lake. Therefore, the analysis of ice melting on the lake surface should refer to the conditions of the first half of April, when the day-averaged solar illumination of the lake practically does not change and the air temperature is the only changing parameter of the problem. According to (Wang et al., 2022), the wind speed in the first 2 weeks of April is equal to 4 m/s and not changed during the day. The corresponding value of the heat transfer coefficient is about h = 20  W/(m²K) (Defraeye et al., 2011; Mirsadeghi et al., 2013).

The boundary-value problem for the quasi-steady temperature profile in the ice layer is as follows: k ice T ″ + P ¯ z = 0 ,   0 < z < d , (13a) − k ice T ′ 0 = q ¯ inf sol − q conv   − q rc ,   T d = T 0 , (13b) where P ¯ z = 1 t day ∫ 0 t day P z , t d t ,   q ¯ inf sol = 1 t day ∫ 0 t day q inf sol t d t ,   q conv   = h T 0 − T air (14a) q inf sol t = ∫ λ * ∞ q i , λ inc λ , t d λ ,   q rc = π ∫ λ w 1 λ w 2 I λ , b T 0 d λ . (14b) Here q i , λ inc is the incident radiative flux at the solar zenith angle θ i . To determine the function P ¯ z = P z × t dl / t day = P z / 2 one can use the average profiles of P z for daylight hours. The same relation is true for the solar infrared radiation flux, taking into account its contribution to the integral radiative flux. The value q ¯ inf sol = 37  W/m² is used in the calculations. It is also convenient to use the following approximation for the temperature dependence of q rc : q rc = q 0 + η × T 0 − T 0 ,   q 0 = 93   W / m 2 , η = 1.6   W / m 2   K . (15)

The radiative cooling of ice compensates the radiative flux to the ice surface in the opacity range and the infrared solar heating cannot lead to the surface ice melting.

The condition for the beginning of ice melting can be obtained using the analytical solution to the problem (13a)–(13b): T z = T 0 + f 2 d − f 2 z k ice − Q k ice d − z , (16a) f 1 z = 1 2 ∫ 0 z P z d z ,   f 2 z = ∫ 0 z f 1 z d z , Q = q conv   + q inf ice − q ¯ inf sol   . (16b)

The ice melting on its lower surface begins when T ′ d = 0 . This enables us to obtain the threshold temperature profile in the ice layer: T z = T 0 + f 2 d − f 2 z k ice − f 1 d k ice d − z . (17)

The calculated temperature profiles are shown in Figure 2. As one might expect, scattering plays a significant role in the case of a thick ice layer and the assumption of uniform distribution of gas bubbles in the ice is acceptable only for ice layers less than 0.5 m thick. One can also determine the temperature of the illuminated ice surface, at which melting begins at the ice-water interface: T surf * = T 0 + f 2 d − f 1 d d k ice . (18)

Typical dependences of T surf * d are plotted in Figure 3. As one might expect, the effect of light scattering by gas bubbles is more significant for thick ice layers.

In the absence of snow cover on the lake ice surface, solar radiation penetrates through the semitransparent ice and already in the beginning of spring significantly heats up the water under the ice. Interestingly, the lake ice, even with a thickness of about 1 m, begins to melt from the lower surface. Calculations have shown that this is not due to heating of the ice by the warmer water, but almost exclusively due to solar heating of the ice. This is so because the upper layer of ice is continuously cooled by the colder air, as well as by the thermal radiation in the middle-infrared window of transparency of the cloudless atmosphere. This cooling is not compensated by the daytime heating of the ice surface by the infrared radiation of the Sun but does not prevent heating of the lower part of the ice layer and the beginning of spring ice melting.

The discussed peculiarity of the thermal regime of ice on the lake surface and its melting under the action of spring solar heating changes radically in the presence of even a thin layer of snow on the ice surface. The effect of snow is due to two main factors: firstly, snow significantly reduces the solar radiative flux on the ice surface due to strong scattering of radiation and, secondly, the thermal conductivity of snow is so small that even a thin layer of snow protects the ice surface from convective and radiative cooling. The physical model of snow’s effect on ice melting should include not only the transfer of solar radiation in the snow and ice layer but also the heat conduction process. If we do not take into account changes in the structure and optical properties of snow during its heating (before the melting), the above physical problems can be solved sequentially: first to calculate radiative transfer, and then to solve the heat transfer problem taking into account the absorption of solar radiation both on the snow surface and in the volume of snow and ice. Of course, in the thermal part of the computational model, it is necessary to consider the convective heat transfer with the surrounding air, as well as the mid-infrared radiative cooling of the snow surface.

Let us first consider the propagation of direct solar radiation through the snow layer. In order not to complicate the solution, we will not take into account the scattered radiation from a cloudless sky. Nor will we take into account that there is ice under the snow, which also scatters sunlight. The ice with air bubbles scatters light much more weakly than snow, and only a small fraction of the light transmitted through the snow is scattered by the ice in the direction of the snow. In the case of a very thin snow layer, the scattering of radiation in the ice layer has some influence on the radiative transfer in the snow layer, but this does not affect the main results of the calculations presented below.

As usual, when calculating radiative transfer in a medium with multiple scattering, the transport approximation is used. In addition, the radiation intensity is represented in the form of two additive components: direct radiation and a diffuse component formed due to the scattering of radiation in the medium. The irradiation is also additive: G ¯ λ = E λ i + G ¯ λ dif ,   E λ i = ⁡ exp − ν i τ λ tr ,   ν i = 1 / μ i . (19)

Here τ λ tr is the optical coordinate measured from the illuminated surface and μ i is the cosine of the Sun’s zenith angle measured from the normal to the horizontal snow surface. The use of the two-flux method leads to the boundary-value problem for the normalized diffuse irradiation: − ( G ¯ λ dif ) ″ + ξ λ 2 G ¯ λ dif = 4 ω λ tr E λ j ,   ξ λ = 2 1 − ω λ tr (20a) τ λ tr = 0 ,   ( G ¯ λ dif ) ′ = 2 G ¯ λ dif ,   τ λ tr = τ λ , 0 tr ,   ( G ¯ λ dif ) ′ = − 2 G ¯ λ dif . (20b)

This problem is a little simpler than that for the ice layer since the sunlight is not reflected and refracted on the snow surface. For a snow layer with constant values of α λ and β λ tr , there is a following analytical solution, which looks similar to Eq. (8a–c), but with replacement of ν j by ν i : G ¯ λ dif = 4 ω λ tr ξ λ 2 − ν i 2 E λ i + 2 − ν i 2 − ξ A E λ dif + B / E λ dif E λ , 0 dif 2 − ζ λ 2 ,   ζ λ = 2 + ξ λ 2 − ξ λ , (21a) A = ψ ζ λ − E λ , 0 i E λ , 0 dif ,   B = E λ , 0 dif ζ λ E λ , 0 i − ψ E λ , 0 dif ,   ψ = 2 + ν i 2 − ν i , (21b) E λ , 0 i = exp − ν i τ λ , 0 tr ,   E λ dif = ⁡ exp − ξ λ τ λ tr ,   E λ , 0 dif = ⁡ exp − ξ λ τ λ , 0 tr . (21c)

In the limit of τ λ , 0 tr = β λ tr d snow → ∞ , this solution coincides with the solution derived by Dombrovsky et al. (2019).

The optical properties of snow were calculated in the same way as in (Dombrovsky et al., 2019). The analytical solution for single spherical ice grains obtained by Kokhanovsky and Zege (1995) in the geometrical optics approximation was used. As usual, it was assumed that the ice grains have a mean radius a = 100  µm and their volume fraction is equal to f v = 0.33 .

The real parameters of the problem may differ from the accepted values. However, the effect of this choice on the subsequent physical analysis is expected to be insignificant. The possible contamination of snow by particles suspended in the atmosphere is not considered. This effect was considered, for example, by He et al. (2018) and Dombrovsky and Kokhanovsky (2020a).

The effect of snow layer thickness on the directional-hemispherical transmittance of snow, T λ d − h , is shown in Figure 4. The value of T λ d − h was calculated at a solar zenith angle θ i = 60 o . One can see that even a 5 mm thick snow layer transmits only about 10%–14% of the incident radiation whereas a thicker snow layer practically does not transmit solar radiation to the ice surface. Note that this result agrees well with the measurements by Perovich (2007). In addition, a small reflection of light from the ice surface and also small scattering of light by gas bubbles in the ice allows a significant simplification of the radiative transfer problem: one can assume that the sunlight transmitted through the snow layer does not return back.

Consider now the heat transfer problem for solar heating and possible snow melting on the ice surface, taking into account the conductive heating of the ice layer. The temperature of the lower surface of the ice (at the ice-water interface) is 0°C. This boundary condition allows one not to consider the heat transfer under the ice, which is necessary when the ice layer is located, for example, on bare ground. At the same time, when analyzing the thermal state of the snow layer illuminated by the Sun, it is necessary to take into account the heat conducted away from the snow into the ice layer. In other words, the energy equation should be solved in the computational domain including both the snow cover and the ice layer: ρ c ∂ T ∂ t = ∂ ∂ z k ∂ T ∂ z + P ,   t > 0 ,   0 < z < d th = d snow + d ice , (22a) T 0 , z = T init z , (22b) z = 0 ,   − k ∂ T ∂ z = h   T air − T + q inf sol − ε   q rc , z = d th , T t , z = T 0 , (22c) q inf sol = ∫ λ * λ op q λ sol t d λ ,   q rc = π ∫ λ w 1 λ w 2 I λ , b T , λ   d λ . (22d) where T t , z is the medium temperature which depends on current time t and coordinate z measured from the upper surface of the snow layer, P t , z is the solar radiation power absorbed in the medium, T air t is the air temperature outside the thermal boundary layer, T 0 = 0 ℃ , q inf sol t is the integral infrared solar radiative flux in the opacity range of λ * < λ < λ op (the upper boundary of this range was taken equal to λ op = 2.4  µm). Specific volumetric heat capacity ρ c and thermal conductivity k of the medium are: ρ c = ρ c snow  when  0 < z < d snow ρ c ice  when  d snow < z < d th , k = k snow  when  0 < z < d snow k ice  when  d snow < z < d th . (23)

The initial temperature profile was assumed to be as follows: T init z = { T air  when  z ≤ d snow T air + T 0 − T air z − d snow / d ice  when  z > d snow . (24)

When modeling the effect of the snow layer on the opening of the lake from ice, it would be incorrect to consider the solar heating conditions typical of the high-mountain lakes of northeastern Tibet, when, due to low precipitation and constantly strong winds, the ice on these lakes is not covered by snow. In addition, at low altitudes, unlike in high mountains, not only direct solar radiation but also diffuse radiation from the light-scattering atmosphere should be taken into account. In this problem, we use the same data as in the papers by Dombrovsky et al. (2019) and Dombrovsky and Kokhanovsky (2022) for the summer solstice during the Arctic polar summer at latitude 70^o.

The motion of the Sun across the sky during the day was calculated in the same way as in (Dombrovsky et al., 2019). At the same time, taking into account the diffuse atmospheric radiation requires additional relations, which are obtained by solving the following boundary-value problem for spectral irradiation: G ¯ λ ″ − ξ 2 G ¯ λ = 0 , (25a) τ λ tr = 0 , G ¯ λ ′ = 2 G ¯ λ − 2 , τ λ tr = τ λ , 0 tr , G ¯ λ ′ = − 2 G ¯ λ . (25b)

The analytical solution to this problem is given by: G ¯ λ = 4 2 − ξ ζ E λ dif − E λ , 0 dif 2 / E λ dif ζ 2 − E λ , 0 dif 2 , (26)

In the limit of τ λ , 0 tr → ∞ , Equation 26 reduces to the following one : G ¯ λ = 4 E λ dif / 2 + ξ .

To exclude the influence of the initial temperature profile, we consider the calculated temperature profiles on day 14 from t = 0 . For certainty, the ice layer thickness is assumed to be equal to d ice = 40  cm, whereas the thickness of the snow layer varies widely. The computational results presented in Figure 5 for three values of d snow differ most markedly at midnight when the snow has time to cool by 0.3–0.4°C as a result of convective and radiative cooling. It is interesting that the maximum heating of snow occurs in the second half of the day and even at d snow = 200  mm only the surface layer of snow 2–3 mm thick melts. This result is qualitatively different from the case of thick snow cover considered in (Dombrovsky et al., 2019; Dombrovsky and Kokhanovsky, 2022) due to heat sinking into the ice layer. Of course, this heat sink is more noticeable in case of a thin snow layer, but remains significant even at d snow = 200  mm. Calculations have shown that snow less than d snow * ≈ 15  mm thick does not melt at all due to cooling by the ice layer. The found threshold value d snow * of the non-melting snow layer is another characteristic parameter of the considered problem. It is important that in all the variants, the ice layer is far from melting.

The above computational study gives the following general picture of solar heating of the ice layer covered by snow. When the upper surface of ice on a lake is illuminated by the spring Sun, the thick ice layer begins to melt even in very cold air and melting occurs from the lower surface of the ice layer. But if there is a thin layer of snow on the surface of the ice (even less than 1 cm thick) the snow does not allow much of the sunlight to penetrate into ice and makes melting of the ice surface underneath impossible. At the same time, the snow itself on the ice surface does not heat up because it scatters a significant part of the visible solar radiation, and the absorbed solar heat is almost immediately transferred to the relatively thick layer of ice. Snow melting can only begin when the snow layer is thicker than d snow * ≈ 15  mm and only when the snow layer is d snow ≈ 200  mm thick does melting become significant. Of course, the produced water flows down through the pores in the snow and freezes again at some distance from the ice surface. If the initial thickness of the snow layer is more than 300 mm, melting of the snow near the sunlit surface continues and leads to the formation of a kind of mushy zone on the ice with ice particles suspended in water, which tend to float to the water surface. This stage of the process can be quite long, but as a result the ice surface turns out to be covered with a layer of water.

Meltwater, formed initially by snowmelt, does not cover the ice surface uniformly. This is evidenced by numerous observations both for ice-covered lakes and for the much better-studied Arctic Sea ice. Note that the abnormally strong melting of sea ice in the Arctic during the polar summer is an important process that has been intensively analyzed over the last decade. One should name a number of studies that are specifically devoted to the formation and evolution of melt ponds on the Arctic Sea ice. In chronological order, these are experimental and analytical works by Polashenski et al. (2012), Hudson et al. (2013), Schröder et al. (2014), Popović et al. (2018), Malinka et al. (2018), Ma et al. (2019), König and Oppelt (2020), Perovich et al. (2021), Sterlin et al. (2021), Lei et al. (2022), and Rosenburg et al. (2023). Photographs in the literature show numerous melt ponds and the change in the ice surface area occupied by these ponds in the polar summer. Note that similar melt ponds are also observed on glaciers (Rockström and Gaffney, 2021). In most of the above-mentioned works on sea ice, the authors are focused on positive feedback that leads to an enormous increase in the total area of the melt ponds. The point is that the albedo of snow in the visible range of the spectrum is very high, and when snow melts in polar summer, forming melt ponds, the reflection of solar radiation is radically reduced. As a result, sea ice receives much more solar heat, causing it to melt and increasing the area occupied by melt ponds. This positive feedback over a large surface area of the Arctic Sea is accompanied by an increase in the water vapor content of the atmosphere and is important for the continuation of global warming.

Even visual observations allow us to distinguish between two types of melt ponds: the so-called bright ponds and dark ponds. The ice on the bottom of the bright pond is mostly smooth and dense but with small cracks and highly light-scattering porous areas with fine pores. The bottom of the dark pond is more heterogeneous and has relatively large cracks and voids (König and Oppelt, 2020). This is related to the different structure of the first-year ice and multiyear ice (Li et al., 2020). The experimental work by König and Oppelt (2020) confirmed the qualitative results of the analytical study of melt pond reflectivity by Malinka et al. (2018). Note that the work by Malinka et al. (2018) was based on the same methodological framework as the present paper, including the transport approximation and the two-flux model for radiative transfer.

For simplicity, we will focus on a more simple melting problem for the first-year ice. The formation of a layer of water on the ice surface and the subsequent melting of ice are complex processes that deserve special modeling. Nevertheless, it is possible to suggest a simple physical model for ice melting under a meltwater layer. It is obvious that on both surfaces of the ice layer, the temperature is constant and equal to 0°C. Of course, the total radiation power absorbed in the ice volume leads not to increase the ice temperature, but to its melting. The heat flux from water under the ice, as was shown above, is insignificant. Indeed, there is a stable temperature stratification of water directly under the ice, and the thermal conductivity of water is rather small. On the contrary, the natural convection of water in a melt pond can give a noticeably larger contribution to ice melting, and this should be taken into consideration. This effect is more significant in the case of a small ice layer thickness.

The temperature profile along the depth of the melt pond is non-monotonic and has a maximum at some depth below the surface. Indeed, the heating of the water surface by infrared solar radiation is compensated by convective cooling and radiative cooling in the middle-infrared transparency window of the cloudless atmosphere. Even at temporary heating of the surface layer of water to the temperature T * = 4 ℃ , this water, as denser, goes down to some depth d * < d p , where d p is the depth of the melt pond. At the bottom of the melt pond, at a depth of d * < d < d p , the water temperature decreases from T * at d = d * to T 0 = 0 ℃ at d = d p (at the ice surface). Intensive natural convection takes place in this lower layer. Obviously, most of the solar radiative flux absorbed in the water volume is transferred to the ice in this lower layer of the melt pond. This is so because in the upper layer of the pond water (at a depth of d < d * ) due to stable temperature stratification, the only heat transfer mode–heat conduction–operates, and the thermal conductivity of water is small: k w = 0.57  W/(m K). As a result, for example, at d − d * = 10  cm and temperature difference T * − T 0 = 4 ℃ the upward heat flux is very small: q cond ≈ 23  W/(m² K). Therefore, in the suggested approximate model, it is assumed that the main part of solar radiation in the spectral range of water semitransparency absorbed in the volume of the melt pond is transferred to the ice layer. We will consider only the day-averaged radiative flux at the melt pond surface. The change of heat transfer conditions on the surface of the melt pond during the day does not matter for the considered ice melting because of the large heat capacity of the ice layer.

Small reflection of sunlight from water can be neglected. In subsequent calculations, the cosine of the refraction angle is assumed constant and equal to μ j = 0.75 , which approximately corresponds to the zenith angle of the Sun θ i = 60 ° . As before, the insignificant difference between the refractive indices of water and ice makes it possible to neglect the reflection and refraction of light both at the bottom of the melt pond and at the lower surface of the ice layer. To simplify the calculations, the propagation of only the direct solar radiation is considered. Of course, a 1D approach that does not take into account the size of the melt pond and its unequal depth in the central part and at the periphery is acceptable for the model problem. Perovich et al. (2021) found that the depth of melt ponds during ice melt varies little with time because it is regulated by the flow of horizontal and vertical water drainage to the ocean. As a result, some of the solar heat absorbed by the water in the pond is not transferred to the ice layer. This should be taken into account, especially for small ponds.

To calculate the transfer of solar radiation in the ice layer under the melt pond, the above suggested method with the spectral irradiation G ¯ λ = G ¯ λ dif + E λ j and the diffuse component calculated by analytical solution (Eq. 10a–10c) is used. Since we are not interested in the distribution of absorbed radiation over the thickness of the ice layer, it is sufficient to consider the spectral irradiation integrated over the optical thickness of the ice. However, the profile of the ice scattering parameter is not known. In addition, the solution depends on physical parameters related to solar illumination and the depth of the melt pond. The natural uncertainty of these data does not allow one to expect accurate results of detailed calculations. Under these conditions, relying on the field measurement data available in the literature is reasonable. The experimental data on geometric parameters and heat fluxes for first-year Arctic Sea ice in advanced stages of melt are given in Hudson et al. (2013).

The following approximation of the radiation power absorbed in the ice layer is used: P = P max ⁡ exp − z / d ice * . (27)

According to Figure 1, at the scattering parameter S = 2  m⁻¹ this dependence with the value d ice * = 0.8  m is universal for ice layers of different thicknesses in the range 0.2 ≤ d ice ≤ 0.8  m. The integration of Eq. 30 gives the following formula for the radiative flux absorbed in an ice layer: q ice = q ice inc 1 − ⁡ exp − d ice / d ice * . (28)

Of course, the dependences of q ice d ice at the real scattering parameter, which, in general, varies with the ice thickness, may differ significantly from function (Equation 28), but retrieval of profiles S z from the available limited experimental data is impossible. Therefore, Equation 28 should be considered an essential assumption adopted in the approximate model.

The integral radiative flux absorbed in a melt pond is given by the following equation: q w = ∫ λ uv ‐ vis λ * q λ w d λ ,   q λ w = q λ , w inc 1 − ⁡ exp − α λ w d p / μ j . (29)

It is assumed that q λ , w inc is constant during a large part of the polar summer and its spectral dependence roughly coincides with the Planck function at temperature T sol = 5800 K. A day-averaged integral incident flux of solar radiation in the range of water and ice semitransparency q inc = q w + q ice inc can be taken from Hudson et al. (2013).

Using the approximate value of q ice , one can estimate the linear rate of ice melting using the heat balance equation: ρ ice L d ˙ ice = − 1 − K loss q w + q ice ,   d ice 0 = d ice , 0 . (30)

The value K loss < 1 is introduced to account the heat losses with horizontal drainage of water. The results of calculations at K loss = 0.2 presented in Figure 6 show two periods of the process: first, the ice melts faster due to the absorption of solar radiation passing through it, and then the heat is transferred to ice mainly from the melt pond. The characteristic time of melting of thick ice is comparable to the duration of the warmest period of the polar summer.

The suggested model can be considered only as a simple physical assessment. The real process needs a more sophisticated analysis which is beyond the scope of the present paper. However, it is clear that a relatively simple approach to solving the radiative transfer problem may be useful in the physical modeling of one of the stages of ice melting and the complex evolution of melt ponds.

As a result of the strong melting of sea ice during the polar summer, significant areas of the Arctic Sea may become ice-free. Most likely, this will not prevent the restoration of the ice cover during the polar winter, which would be an extremely undesirable effect of global warming. The fact is that the surface layer of water in the Arctic Sea cannot heat up above the melting point of ice due to the natural convection of water in its unstable temperature stratification. Thus, the unique physical properties of water allow us to count on serious negative feedback and retain some optimism about the rate of global climate change accompanied by intensive seasonal melting of polar ice.

Physical and computational models of solar radiative heating and melting of ice on water surface were developed. In the semitransparency range of ice and snow, the multiple scattering of radiation in the ice which contains gas bubbles and in the snow cover that may be present on the ice surface is described using the transport approximation for the scattering phase function and two-flux model for the diffuse component of the spectral radiation intensity. Analytical solution derived for the radiative transfer are coupled with a numerical solution for the transient energy equation. The heat transfer model takes into account convective heat transfer with the atmospheric air, infrared solar radiation absorbed at the illuminated surface, and radiative cooling in the mid-infrared transparency window of the cloudless atmosphere.

It was shown that a thick layer of ice not covered with snow begins to melt at the ice-water interface due to solar heating of ice. The computational data are in good agreement with the field observations for Ngoring Lake in the Qinghai–Tibet Plateau. The theoretical analysis showed a dramatic change in the process if there is a layer of snow on the ice. Even in the case of a thin snow layer less than 1 cm thick, the snow does not transmit most of the sunlight and makes ice melting impossible. At the same time, the snow itself on the ice surface does not heat up because it scatters a significant part of visible solar radiation, and the absorbed solar heat is transferred almost instantly to the relatively thick layer of ice.

Snow melting can only begin when the snow layer is about 15 mm thick, and only when the snow layer is about 200 mm thick does melting become significant. If the initial thickness of the snow layer exceeds 300 mm, snow melting near the sunlit surface leads to the formation of a melt pond on the ice surface. Such melt pools are regularly observed in polar summer on the ice of the Arctic Sea. Abnormally strong melting of sea ice in the Arctic during polar summer is an extremely important process that has been intensively studied in recent years. The last part of the paper provided an estimate of ice melting under a melt pond. The results obtained are in qualitative agreement with in-situ observations.

The general approach and particular solutions suggested for approximate calculations of solar heating and melting of lake and sea ice with possible snow cover on the ice surface can be used as the basis of more specific computational models for a variety of combined heat transfer problems in geophysics and other research areas.