## SOLAR HEATING OF ICE-COVERED LAKES AND ICE MELTING

**Following from: **
Solar heating of ice containing bubbles

The study of solar heating of large ice-covered lakes is of interest for two reasons. First, freeing such lakes from ice significantly affects the atmosphere in local regions as well as much larger areas (Su et al., 2020). Second, the melting of ice in lakes in recent years is an indicator of global climate change (Su et al., 2019; Zhang and Duan, 2021). It should also be noted that sunlight favorably affects microalgae life present under ice (Song et al., 2019).

In an earlier work, Matthews and Heaney (1987) studied the solar heating of an ice-covered lake and proposed a model that takes into account the natural convection of water. Unfortunately, in that study, Bouguer's law was used to describe the transfer of solar radiation in water and a layer of ice. This drawback remained in subsequent studies (Mironov et al., 2002; Kirillin et al., 2012, 2021; Aslamov et al., 2014; Leppäranta, 2015). The simplified approach was partially compensated for by the selection of an extinction coefficient, which provided satisfactory agreement with field measurements. The main efforts in these previous works were focused on the convection of water (Mironov et al., 2002). The traditional approach to radiative transfer used in limnology studies is insufficient. The latter approach is especially important in ice layers containing microcracks and small gas bubbles, which intensely scatter radiation. Obviously, multiple scattering leads to stronger solar heating of ice layers. A recent study by Dombrovsky and Kokhanovsky (2023) has shown that taking into account the effects of multiple scattering of light changes the common belief that ice melts due to heating from warmer water. It turns out that ice melting at the ice–water interface is mainly due to solar heating of the ice layer.

In this study, we considered the effects of solar radiation in the visible and near-infrared regions that contribute significantly to heating an ice layer and the water beneath the ice. This radiation is partially reflected from the ice surface; however, the main part of the radiation is scattered and absorbed in both the ice layer and the lake water.

The radiative transfer problem under consideration is more complicated than the one solved by Dombrovsky and Kokhanovsky (2020) for an optically thick layer of ice. However, the problem can be simplified for lakes similar to those found in the Qinghai–Tibet Plateau, which is located about 4 km above sea level. The sky in this region is usually clear and the light scattered by the atmosphere can be neglected. Insignificant atmospheric precipitation and strong winds result in the absence of snow cover on the ice surface. Another simplification of the problem is due to the small differences between the refractive indices of ice and water in the range of \(\lambda_{{\rm uv-vis}} \lt\lambda \lt\lambda_{*}\), where \(\lambda_{{\rm uv-vis}} =\) 0.4 μm and \(\lambda_{*} =\) 1.2 μm (Hale and Querry, 1973; Warren and Brandt, 2008); therefore, the reflection and refraction of light at the ice–water interface can be neglected. According to Kirillin et al. (2012), the ice layer on a lake's surface usually contains gas bubbles and microalgae. Gas bubbles scatter radiation much more strongly than particles of the same size (Dombrovsky and Baillis, 2010) and scattering by small microalgae can be neglected. The bubbles are usually unevenly distributed in the ice layer. However, in this study we limited ourselves to two bubble parameters: the average volume concentration and size. Moreover, following the study by Dombrovsky and Kokhanovsky (2020), we used only one parameter to describe light scattering by polydisperse gas bubbles: \(S=f_{v} /a_{32}\), where \(f_{v}\) is the volume fraction of the bubbles in ice and \(a_{32}\) is the Sauter average radius of the bubbles.

The shortwave solar radiative transfer model is based on the radiation transfer theory. We assumed that the optical properties of ice do not change along the ice surface. Therefore, it is sufficient to calculate only the vertical profiles of the absorbed radiation power at each time point, i.e., under different illumination conditions. The method for calculating the propagation of obliquely incident radiation by solving an equivalent one-dimensional problem is described in detail in Dombrovsky and Kokhanovsky (2022). Since ice weakly absorbs shortwave radiation, this results in the multiple scattering of radiation by the gas bubbles and the transport approximation can be used (Dombrovsky, 2012). According to this approach, the scattering phase function is replaced by the sum of the isotropic component and the term describing the peak of forward scattering. The sufficient accuracy of this approach has been demonstrated by comparison with calculations using the Heyney–Greenstein model (Dombrovsky and Baillis, 2010; Dombrovsky et al., 2011).

Insignificant scattering of radiation by bubbles and plankton in water can be neglected. Therefore, the radiative transfer in the ice layer can be considered independently. The equations for the radiative transfer and oblique illumination boundary conditions of the ice surface are the following:

\(\mu \dfrac{\partial \bar{I}_{\lambda ,i} }{\partial \tau_{\lambda }^{\rm tr} } +\bar{I}_{\lambda ,i} =\dfrac{\omega_{\lambda }^{\rm tr} }{2} \bar{G}_{\lambda ,i} ,\quad \bar{G}_{\lambda ,i} \left(\tau_{\lambda }^{\rm tr} \right)=\int\limits_{-1}^{1} \bar{I}_{\lambda ,i} \left(\tau_{\lambda }^{\rm tr} ,\mu \right)d\mu\) | (1a) |

\(\bar{I}_{\lambda ,i} (0,\mu)=r_{\lambda ,i} \bar{I}_{\lambda ,i} (0,-\mu)+(1-r_{\lambda ,i})\; \delta (\mu_{j} -\mu),\quad \bar{I}_{\lambda ,i} \left(\tau_{\lambda 0}^{\rm tr} ,-\mu \right)=0,\quad \; \mu ,\; \mu_{j} \gt 0\) | (1b) |

where \(\bar{I}_{\lambda ,i} =I_{\lambda ,i} / I_{\lambda ,i}^{\text{inc}}\); \(\bar{G}_{\lambda ,i} =G_{\lambda ,i} / I_{\lambda ,i}^{\text{inc}}\); \(I_{\lambda ,i}^{\text{inc}}\) is the spectral intensity of the incident radiation in direction \(\mu_{i} =\cos \theta_{i}\) (where \(\theta_{i}\) is measured from the external normal); \(r_{\lambda ,i} =r_{\lambda } (n,\mu_{i})\) is the Fresnel reflection coefficient; \(\mu_{j} =\sqrt{1-\left(1-\mu_{i}^{2} \right)/ n^{2}}\) is the cosine of the refraction angle; \(\tau_{\lambda }^{\rm tr} (z)=\int_{0}^{z}\beta_{\lambda }^{\rm tr} (z)dz\) is the optical depth measured from the upper ice surface (where \(z=\) 0); \(\tau_{\lambda ,0}^{\rm tr} =\tau_{\lambda }^{\rm tr} (d)\); \(d\) is the ice layer thickness; \(\omega_{\lambda }^{\rm tr} (z)=\sigma_{\lambda }^{\rm tr} (z)/ \beta_{\lambda }^{\rm tr} (z)\) is the transport albedo of the medium; \(\sigma_{\lambda }^{\rm tr} =\sigma_{\lambda } \times (1-\bar{\mu })\) and \(\beta_{\lambda }^{\rm tr} =\alpha_{\lambda } +\sigma_{\lambda }^{\rm tr}\) are the scattering and extinction transport coefficients, respectively; \(\bar{\mu }\) is the asymmetry factor of scattering; and \(\alpha_{\lambda}\) is the spectral absorption coefficient.

The linearity of the radiative transfer problem allows us to represent the radiation intensity as the sum of the diffuse component (which is due to light scattering in the medium) and the directional sunlight:

\(\bar{I}_{\lambda ,i} =\bar{J}_{\lambda ,i} +(1-r_{\lambda ,i})E_{\lambda ,j} \delta (\mu_{j} -\mu ),\quad E_{\lambda ,j} =\exp \left({-\tau_{\lambda }^{\rm tr} / \mu_{j} } \right)\) | (2a) |

\(\bar{G}_{\lambda ,i} =\bar{G}_{\lambda ,i}^{\text{diff}} +(1-r_{\lambda ,i})E_{\lambda ,j} ,\quad \bar{G}_{\lambda ,i}^{\text{diff}} =\int\limits_{-1}^{1}\bar{J}_{\lambda } d\mu\) | (2b) |

The mathematical formulation of the problem for the diffuse component is as follows:

\(\mu \dfrac{\partial \bar{J}_{\lambda ,i} }{\partial \tau_{\lambda }^{\rm tr} } +\bar{J}_{\lambda ,i} =\dfrac{\omega_{\lambda }^{\rm tr} }{2} \bar{G}_{\lambda ,i}\) | (3a) |

\(\bar{J}_{\lambda ,i} (0,\mu )=r_{\lambda ,i} \bar{J}_{\lambda ,i} (0,-\mu),\quad \bar{J}_{\lambda ,i} (d,-\mu)=0,\quad\mu \gt 0\) | (3b) |

The asymmetry of the boundary conditions on the upper and lower surfaces of an ice layer does not allow us to use the modified differential approximation (Dombrovsky et al. 2006), as was done in Dombrovsky and Kokhanovsky (2020). Therefore, we use the ordinary two-flux method. We obtain the following boundary-value problem by integrating Eqs. (3a) and (3b) over the angle:

\(-\left(\bar{G}_{\lambda ,i}^{\text{diff}} \right)'' +\xi_{\lambda }^{2} \bar{G}_{\lambda ,i}^{\text{diff}} =4\omega_{\lambda }^{\rm tr} (1-r_{\lambda ,i})E_{\lambda ,j}\) | (4a) |

\(\tau_{\lambda }^{\rm tr} =0,\quad\left(\bar{G}_{\lambda ,i}^{\text{diff}} \right)' =2\gamma \bar{G}_{\lambda ,i}^{\text{diff}},\quad \tau_{\lambda }^{\rm tr} =\tau_{\lambda ,0}^{\rm tr},\quad \left(\bar{G}_{\lambda ,i}^{\text{diff}} \right)' =-2\bar{G}_{\lambda ,i}^{\text{diff}}\) | (4b) |

where \(\bar{r}_{\lambda }\) is the angle-averaged reflection coefficient; \(\xi_{\lambda } =2\sqrt{1-\omega_{\lambda }^{\rm tr} }\); and \(\gamma ={\left(1-\bar{r}_{\lambda } \right)/\left(1+\bar{r}_{\lambda } \right)}\). For simplicity, we assume that the optical properties of ice are uniform. This makes it possible to obtain an analytical solution to the problem. It will be shown below that varying the optical properties of ice is sufficient to estimate the small effect of changing these properties along the ice layer thickness. The analytical solution to this problem at \(\xi_{\lambda } \ne \nu_{j} ={1/\mu_{j}}\) is as follows:

\(\bar{G}_{\lambda ,i}^{\text{diff}} =C\left(E_{\lambda ,j} -AE_{\lambda ,i}^{\text{diff}} +\dfrac{B}{E_{\lambda ,i}^{\text{diff}} } \right),\quad C=\dfrac{4\omega_{\lambda }^{\rm tr} \left(1-r_{\lambda ,i} \right)}{\xi_{\lambda }^{2} -\nu_{j}^{2} },\quad E_{\lambda ,i}^{\text{diff}} =\exp \left(-\xi_{\lambda } \tau_{\lambda }^{\rm tr} \right)\) | (5a) |

\(\begin{array}{c} {A} = \dfrac{(2+\xi_{\lambda })(2\gamma +\nu_{j} )\Big/(2\gamma -\xi_{\lambda })-(2\gamma -\nu_{j})\exp\!\left[-(\xi_{\lambda } +\nu_{j})\tau_{\lambda ,0}^{\rm tr} \right]}{(2+\xi_{\lambda })\zeta_{\lambda } -(2-\xi_{\lambda } )\exp \!\left(-2\xi_{\lambda } \tau_{\lambda ,0}^{\rm tr} \right)}, \\[5pt] B = (A-1)\zeta_{\lambda },\quad \zeta_{\lambda } =\dfrac{2\gamma +\xi_{\lambda } }{2\gamma -\xi_{\lambda } } \end{array}\) | (5b) |

Radiative transfer in a layer of water is a much simpler problem:

\(\mu \dfrac{\partial \bar{I}_{\lambda ,i}^{w} }{\partial \tau_{\lambda }^{(w)} } +\bar{I}_{\lambda ,i}^{w} =0,\quad \tau_{\lambda }^{w} =\alpha_{\lambda }^{w} \times (z-d),\quad z\gt d\) | (6a) |

\(\bar{I}_{\lambda ,i}^{w} (d,\mu)=\bar{I}_{\lambda ,i} (d,\mu),\quad \bar{I}_{\lambda ,i}^{w} (\infty ,\mu )=\bar{I}_{\lambda ,i}^{w} (\infty ,-\mu )=0,\quad \mu \gt 0\) | (6b) |

The solution to this problem is obvious:

\(\bar{I}_{\lambda ,i}^{w} \left(z,\mu \right)=\left\{\begin{array}{ccc} {\bar{I}_{\lambda ,i} \left(d,\mu \right)\exp \left({-\tau_{\lambda }^{w} / \mu } \right)\; } & \text{when} & {\mu \gt 0} \\[4pt] {0} & \text{when} & {\mu \lt 0} \end{array}\right.\) | (7) |

\(P(z)=\int\limits_{\lambda_{{\rm uv-vis}} }^{\lambda_{*} } p\left(\lambda ,z\right)d\lambda,\quad p\left(\lambda ,z\right)=\left\{\begin{array}{ccc} \alpha_{\lambda } G_{\lambda ,i} (z)\; & \text{when} & z\le d \\[4pt] \alpha_{\lambda }^{w} G_{\lambda ,i}^{w} (z) & \text{when} & z\gt d \end{array}\right.\) | (8a) |

\(G_{\lambda ,i} (z)=I_{\lambda ,i}^{\text{inc}} \times \left\{\bar{G}_{\lambda ,i}^{\text{diff}} (z)+\left(1-r_{\lambda ,i} \right)\exp \left(-\beta_{\lambda }^{\rm tr} z\right)\right\}\) | (8b) |

\(G_{\lambda ,i}^{w} (z)=I_{\lambda ,i}^{\text{inc}} \times \left\{\bar{G}_{\lambda ,i}^{\text{diff}} (d)\exp \left(-2\tau_{\lambda }^{w} \right)+\left(1-r_{\lambda ,i} \right)\exp \left(-\beta_{\lambda }^{\rm tr} d-\tau_{\lambda }^{w} \right)\right\}\) | (8c) |

The first term in Eq. (8c) is written as a two-flux approximation for the incoming diffuse radiation, and the second term corresponds to Bouguer's law for directional radiation.

The experimental data for the refractive index, \(n(\lambda)\), and absorption index, \(\kappa(\lambda)\), of water and ice are presented in Fig. 1. In general, the spectral dependences of both indices for water and ice are close to each other; however, there are significant discrepancies between the data obtained by different authors for the index of absorption of ice in the wavelength range of 0.4 \(\lt\lambda\lt\) 0.6 μm [Fig. 1, right, curves 1 and 2]. Further calculations for ice are performed both for data of Warren and Brandt (2008) and Picard et al. (2016).

**Figure 1. **Indices of refraction and absorption of freshwater and ice in the semi-transparency range [for water, data from Hale and Querry (1973) are shown; for ice, data from Warren and Brandt (2008) (1) and Picard et al. (2016) (2) are shown; reprinted from Dombrovsky and Kokhanovsky (2023) with permission from Elsevier, copyright 2023]

Water and ice absorption coefficients are calculated using simple formulas linking the absorption coefficient of a homogeneous medium to the absorption index, since gas bubbles in ice have almost no effect on the absorption coefficient:

\(\alpha_{\lambda }^{(w)} =\dfrac{4\pi \kappa_{w} (\lambda)}{\lambda } ,\quad \; \alpha_{\lambda } =\dfrac{4\pi \kappa_{\text{ice}} (\lambda)}{\lambda }\) | (9) |

At wavelength \(\lambda_{*}\), ice absorption coefficient \(\alpha_{\lambda } =\) 70 m^{–1}. It can be shown that at \(\lambda =\) 1.4 μm this value increases to 178 m^{–1}. This means that when an ice layer is illuminated along the normal to the surface, the radiation intensity at a depth of 1 cm at wavelength \(\lambda_{*}\) decreases by two times, and at \(\lambda =\) 1.4 μm it decreases by six times. The ice layer on the lake surface is several times thicker and the radiation at wavelength \(\lambda \gt \lambda_{*}\) is absorbed in a relatively thin surface layer. The transport scattering coefficient of ice with gas bubbles is determined as follows (Dombrovsky and Kokhanovsky, 2020):

\(\sigma_{\lambda }^{\rm tr} =0.675\left[n(\lambda)-1\right]S\) | (10) |

As a rule, the value of \(\sigma_{\lambda }^{\rm tr}\) is significantly greater than the absorption coefficient. Therefore, the extinction of sunlight in an ice layer is determined by the scattering of light by gas bubbles.

The typical spectral dependencies of the solar radiative flux in the Ngoring Lake area at the end of March are shown in Fig. 2. The time dependencies of the solar zenith angle and integral radiative flux are well approximated by the following functions (Dombrovsky and Kokhanovsky, 2023):

\(\theta_{i} (t)=\theta_{\min } +\left(\dfrac{\pi }{2} -\theta_{\min } \right)\left(1-\cos \psi \right),\quad q_{\text{int}}^{\text{sol}} (t)=q_{\text{int},\max }^{\text{int}} \cos \psi ,\quad \psi =\dfrac{t-t_{*} }{t_{dl} } \dfrac{\pi }{4}\) | (11) |

where \(t_{dl} =\) 12 h is the duration of daylight; \(t_{*}\) is the local time at which the zenith angle of the sun is minimal; \(\theta_{\min} =\pi / 6\); and \(q_{\text{int},\max }^{\text{int}} =\) 940 W/m^{2}. One can choose a value of \(\theta_{i}\) corresponding to the average irradiation intensity. Taking \(\cos \psi =\) 0.5, we obtain \(\theta_{i} =\pi / 3\). The calculated profiles of the absorbed radiation power for two values of the scattering parameter are shown in Fig. 3.

**Figure 2. **The spectral solar radiative flux at the lake surface [reprinted from Dombrovsky and Kokhanovsky (2023) with permission from Elsevier, copyright 2023]

**Figure 3. **Profiles of absorbed radiation power calculated using various data for the ice absorption index: (a) \(S =\) 2 m^{–1}; (b) \(S =\) 10 m^{–1} [solid lines, data from Warren and Brandt (2008); dashed lines, data from Picard et al. (2016); reprinted from Dombrovsky and Kokhanovsky (2023) with permission from Elsevier, copyright 2023]

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 ice layer and a decrease in the absorption in water under the ice [compare Figs. 3(a) and 3(b)]. When ice melts, the absorption of radiation in the ice layer decreases significantly; at the same time, the absorption of radiation in water increases very strongly. The latter is the physical cause of a recently reported observation (Lazhu et al., 2021) at several lakes in the Tibetan Plateau, where the water temperature at some distance from the ice–water interface increased rapidly during ice melting. The use of different data for the ice absorption index does not lead to a significant change in the calculated absorption of radiation in the ice layer.

The lake water temperature field is largely determined by the known non-monotonic dependence of the freshwater density on the temperature with a maximum density at \(T_{*} =\) 4°C. Water weakly absorbs sunlight and the radiation is absorbed in a layer of water several meters thick. Natural convection moves water heated up to a temperature of \(T_{0} \lt T\lt T_{*}\) closer to the bottom of the lake. Gradually, the region at temperature \(T_{*}\) covers almost the entire volume, and only in a relatively thin layer of water under the ice does the temperature remain lower. The period during which almost all of the water in a lake is heated by solar radiation to temperature \(T_{*}\) is called normal winter. After that, an interesting thermal regime, called anomalous winter (see Kirillin et al., 2021), occurs when the water under the ice layer is heated by solar radiation to temperature \(T_{\max } \gt T_{*}\) at a depth of 1.5–3 m. Water at temperatures in the range of \(T_{*} \lt T\lt T_{\max }\) is less dense than at \(T_{*}\) and there is stable temperature stratification below the maximum temperature position. At the same time, the upper layer of water at temperatures less than the maximum can be divided into two zones: a thin lower part (where \(T\gt T_{*}\)) and an upper part (where there should be no convection). However, solar radiation is also absorbed in the upper part and this heat is transferred downward by natural convection. As a result, the maximum water temperature in March and early April increases rapidly.

The heat flux from water to the ice layer can be easily estimated. Dombrovsky and Kokhanovsky (2023) showed that this heat flux can be neglected when calculating the temperature profile in the ice layer. The complete computational model of the transient thermal state of the ice layer on the lake surface can be based on the transient energy equation by taking into account the volumetric absorption of solar radiation, heat conduction, and latent heat of melting. The boundary conditions at the upper surface should consider infrared solar radiation, convective heat transfer to the air, and radiative cooling of ice in the infrared transparency window of a cloudless atmosphere.

The conditions for heat transfer on the upper surface of an ice layer change significantly during the day due to variations in solar illumination and changes in the wind speed. This affects the surface temperature; however, the thermal state of ice at some distance from the surface is not sensitive to these conditions. The volumetric absorption of shortwave radiation also mainly affects the heat generation in the upper part of the ice layer. Since the longest initial period of heating and melting of a thick ice layer is of the greatest interest, the changes in the heat transfer parameters during the day are not important. Therefore, an approximate computational method based on the averaged heat flux was suggested in Dombrovsky and Kokhanovsky (2023). The applicability of this approach can be evaluated using the Fourier criterion for an ice layer. It turns out that even at ice layer thickness \(d =\) 0.5 m the thermal relaxation time is about 10 days. Therefore, a steady-state model can be used with heat transfer parameters that vary from week to week.

The initial period of ice melting at Ngoring Lake, when the ice layer thickness decreases by 10 cm from 1.1 to 1 m, takes all of March, and only with the onset of April does ice melting accelerate, with all of the ice having completely melted by April 16. This result is clear since the air temperature in March is approximately constant and does not exceed –10°C, while in April it increases significantly and reaches –3°C by the middle of the month. Analyses of ice melting conditions should refer to the conditions in the first half of April, when solar illumination does not change, the wind speed is 4 m/s [the convective heat transfer coefficient is \(h =\) 20 W/(m^{2} K)], and the air temperature is the only changing parameter. The boundary-value problem for the temperature profile in the ice layer is the following:

\(k_{\text{ice}} T''+\bar{P}(z)=0,\quad 0\lt z\lt d\) | (12a) |

\(-k_{\text{ice}} T'(0)=\bar{q}_{\text{inf}}^{\text{sol}} -q_{\text{conv}} -q_{\text{inf}}^{\text{ice}} ,\quad T(d)=T_{0}\) | (12b) |

\(\bar{P}(z)=\dfrac{1}{t_{\text{day}} } \int\limits_{0}^{t_{\text{day}} } P(z,t)dt,\quad \; \bar{q}_{\text{inf}}^{\text{sol}} =\dfrac{1}{t_{\text{day}} } \int\limits_{0}^{t_{\text{day}} } q_{\text{inf}}^{\text{sol}} (t)dt,\quad q_{\text{conv}} =h\left[T(0)-T_{\text{air}} \right]\) | (12c) |

\(q_{\text{inf}}^{\text{sol}} (t)=\int\limits_{\lambda_{*} }^{\infty } q_{i,\lambda }^{\text{inc}} (\lambda ,t)d\lambda \; q_{\text{inf}}^{\text{ice}} =\int\limits_{\lambda_{c,1} }^{\lambda_{c,2} } I_{b} [T(0),\lambda]d\lambda ,\quad \lambda_{c,1} =8\; \text{μm},\quad \lambda_{c,2} =13\; \text{μm}\) | (12d) |

where \(q_{i,\lambda }^{\text{inc}}\) is the incident radiative flux at solar zenith angle \(\theta_{i}\); and \(q_{\text{inf}}^{\text{ice}}\) is the integral flux of infrared radiation from the ice surface in the transparency window of the atmosphere (the so-called radiative cooling). Functions \(\bar{P}(z)=P(z) t_{dl} / t_{\text{day}} =P(z)/2\) can be calculated using the average profiles of absorbed radiation, \(P(z)\). The same relation is true for the solar infrared radiative flux. The resulting value \(\bar{q}_{\inf}^{\text{sol}} =\) 37 W/m^{2} is used in the calculations. It should be noted that radiative cooling completely compensates for the solar radiative flux in the opacity range and infrared solar radiation cannot lead to surface ice melting. The analytical solution to the problem is as follows:

\(T(z)=T_{0} +\dfrac{f_{2} (d)-f_{2} (z)}{k_{\text{ice}} } -\dfrac{Q}{k_{\text{ice}} } (d-z)\) | (13a) |

\(f_{1} (z)=\dfrac{1}{2} \int\limits_{0}^{z} P(z)dz,\quad f_{2} (z)=\int\limits_{0}^{z} f_{1} (z)dz,\quad Q=q_{\text{conv}} +q_{\inf }^{\text{ice}} -\bar{q}_{\inf }^{\text{sol}}\) | (13b) |

The condition for the beginning of ice melting, \(T'(d) =\) 0, enables us to obtain the threshold profile as follows:

\(T(z)=T_{0} +\dfrac{f_{2} (d)-f_{2} (z)}{k_{\text{ice}} } -\dfrac{f_{1} (d)}{k_{\text{ice}} } (d-z)\) | (14) |

The calculations showed that temperature profiles at the beginning of melting are weakly dependent on the scattering parameter. This result indicates that the model for uniform ice structure is acceptable. Of interest is the low temperature of the upper ice surface (\(T_{{\rm surf}}^{*}\)) at which melting begins. This temperature can be determined by remote sensing to estimate the beginning of ice melting. The \(T_{{\rm surf}}^{*} (d)\) curves shown in Fig. 4 illustrate that the effect considered is not sensitive to the scattering.

**Figure 4. **Dependencies of the surface temperature threshold value on the ice thickness [reprinted from Dombrovsky and Kokhanovsky (2023) with permission from Elsevier, copyright 2023]

The presented methodological developments are general and can be useful in computational modeling of solar heating of ice-covered lakes under different climatic conditions. Similar models can be developed for large ice-covered rivers as well as for sea ice melt analysis.

#### REFERENCES

Aslamov, I.A., Kozlov, V.V., Kirillin, G.B., Mizandrontsev, I.B., Kucher, K.M., Makarov, M.M., Gornov, A.Yu., and Granin, N.G. (2014) Ice-Water Heat Exchange during Ice Growth in Lake Baikal, *J. Great Lakes Res.*, 40(3): 599–607.

Dombrovsky, L.S. (2012) The Use of Transport Approximation and Diffusion-Based Models in Radiative Transfer Calculations, *Comput. Thermal Sci.: Int. J.*, 4(4): 297–315.

Dombrovsky, L.A. and Baillis, D. (2010) *Thermal Radiation in Disperse Systems: An Engineering Approach*, New York: Begell House.

Dombrovsky, L.A. and Kokhanovsky, A.A. (2020) Solar Heating of Ice Sheets Containing Gas Bubbles, *J. Quant. Spectrosc. Radiat. Transf.*, 250: 106991.

Dombrovsky, L.A. and Kokhanovsky, A.A. (2022) Deep Heating of a Snowpack by Solar Radiation, *Front. Therm. Eng.*, 2: 882941.

Dombrovsky, L.A. and Kokhanovsky, A.A. (2023) Solar Heating of Ice-Covered Lake and Ice Melting, *J. Quant. Spectrosc. Radiat. Transf.*, 294C: 108391.

Dombrovsky, L., Randrianalisoa, J., and Baillis, D. (2006) Modified Two-Flux Approximation for Identification of Radiative Properties of Absorbing and Scattering Media from Directional-Hemispherical Measurements, *J. Opt. Soc. Am. A*, 23(1): 91–98.

Dombrovsky, L.A., Randrianalisoa, J.H., Lipiñski, W., and Baillis, D. (2011) Approximate Analytical Solution to Normal Emittance of Semi-Transparent Layer of an Absorbing, Scattering, and Refracting Medium, *J. Quant. Spectrosc. Radiat. Transf.*, 112(12): 1987–1994.

Hale, G.M. and Querry, M.P. (1973) Optical Constants of Water in the 200-nm to 200-μm Wavelength Region, *Appl. Opt*., 12(3): 555–563.

Kirillin, G., Leppäranta, M., Terzhevik, A., Granin, N., Bernhard, J., Engelhardt, C., Efremova, T., Golosov, S., Palshin, N., Sherstyankin, P., Zdorovennova, G., and Zdorovennov, R. (2012) Physics of Seasonally Ice-Covered Lakes: A Review, *Aquat. Sci.*, 74(4): 659–682.

Kirillin, G.V., Shatwell, T., and Wen, L. (2021) Ice-Covered Lakes of Tibetan Plateau as Solar Heat Collectors, *Geophys. Res. Lett.*, 48(14): e2021GL093429.

Lazhu, Y., Yang, K., Hou, J., Wang, J., Lei, Y., Zhu, L., Chen, Y., Wang, M., and He, X. (2021) A New Finding on the Prevalence of Rapid Water Warming during Lake Ice Melting on the Tibetan Plateau, *Sci. Bull.*, 66: 2358–2361.

Leppäranta, M. (2015) *Freezing of Lakes and the Evolution of Their Ice Cover*, Berlin: Springer, 2015.

Matthews, P.C. and Heaney, S.I. (1987) Solar Heating and Its Influence on Mixing in Ice-Covered Lakes, *Freshwater Biol.*, 18: 135–149.

Mironov, D., Terzhevik, A., Kirillin, G., Jonas, T., Malm, J., and Farmer, D. (2002) Radiatively Driven Convection in Ice-Covered Lakes: Observations, Scaling, and a Mixed Layer Model, *J. Geophys. Res.*, 107(C4): 3032.

Picard, G., Libois, Q., and Arnaud, L. (2016) Refinement of the Ice Absorption Spectrum in the Visible Using Radiance Profile Measurements in Antarctic Snow, *Cryosphere*, 10(6): 2655–2672.

Song, S., Li, C., Shi, X., Zhao, S., Tian, Q., Li, Z., Bai, Y., Cao, X., Wang, Q., Huotari, J., Tulonen, T., Uusheimo, S., Leppäranta, M., Loehr, J., and Arvola, L. (2019) Under-Ice Metabolism in a Shallow Lake in a Cold and Arid Climate, *Freshwater Biol.*, 64(10): 1710–1720.

Su, D., Hu, X., Wen, L., Lyu, S., Gao, X., Zhao, I., Li, Z., Du, J., and Kirillin, G. (2019) Numerical Study on the Response of the Largest Lake in China to Climate Change, *Hydrol. Earth Syst. Sci.*, 23: 2093–2109.

Su, D., Wen, L., Gao, X., Leppäranta, M., Song, X., Shi, Q., and Kirillin, G. (2020) Effects of the Largest Lake of the Tibetan Plateau on the Regional Climate, *J. Geophys. Res.: Atmos.*, 125(22): e2020JD033.

Warren, S.G. and Brandt, R.E. (2008) Optical Constants of Ice from the Ultraviolet to the Microwave: A Revised Compilation, *J. Geophys. Res.*, 113(D14): D14220.

Zhang, G. and Duan, S. (2021) Lakes as Sentinels of Climate Change on the Tibetan Plateau, *All Earth*, 33(1): 161–165.