In recent years, there has been a marked melting of snow and ice in the polar and circumpolar regions, which is a significant factor in overall climate change (Barry and Hall-McKim 2018; Kokhanovsky and Tomasi 2020). One of the problems deserving special attention is the solar heating and even melting of snow at a depth below the surface of the snow cover. This effect was probably first discussed by Brandt and Warren (1993) and the term “solid-state greenhouse” was suggested for this phenomenon which was discussed also by Liston and Winther (2005). The physical explanation for this effect, typical of a long polar summer, is based on the interaction of two factors. First of all, it is the very weak absorption of visible radiation by the ice particles that form the snow. Therefore, a part of the visible solar radiation is absorbed several centimeters from the surface of the snow cover (Kokhanovsky, 2021a; 2021b). The absorption maximum turns out to be at some distance from the surface, and a special attention is paid to this issue in the article. Second, due to the low thermal conductivity of snow, much of this heat does not escape to the surface. This heat is continuously transferred by heat conduction to a considerable depth, where it accumulates and can even lead to snow melting away from the surface, which is usually cooled by convective heat transfer to the cold air above snow, as well as by radiative cooling in the infrared transparency window of the cloudless atmosphere.

Let us first consider the conditions under which the absorption maximum of the shortwave external radiation is located at some distance from the surface of the scattering medium illuminated by the directional radiation. To solve this problem, we can use the approach for arbitrary oblique illumination proposed in (Dombrovsky and Randrianalisoa 2018) and used in (Dombrovsky et al., 2019) to calculate the radiative transfer in a layer of snow. The computational results obtained by Dombrovsky et al. (2019) make clear the role of convective cooling or heating of snow, mid-infrared radiative cooling of snow surface, and continuous heating of deep snow layers due to heat conduction. Some numerical results for the effect of a sloping snow surface on solar heating of snow were also obtained.

In the layer of isotropic snow, the general radiative transfer equation (RTE) is as follows (Dombrovsky and Baillis 2010; Howell et al., 2021; Modest and Mazumder 2021): Ω ∇ I λ ( r , Ω ) + β λ ( r ) I λ ( r , Ω ) = σ λ ( r ) 4 π ∫ ( 4 π ) I λ ( r , Ω ′ )   Φ λ ( r , Ω ′ × Ω )   d Ω ′ (1) where I λ is the spectral radiation intensity at point r in direction Ω , β λ = α λ + σ λ is the extinction coefficient, α λ and σ λ are the coefficients of absorption and scattering, respectively, and Φ λ is the scattering phase function. The thermal radiation of snow in the middle-infrared is responsible for the so-called radiative cooling. Obviously, this radiation is not contained in Eq. 1. It should be taken into account only in the boundary condition of the energy equation.

The weak absorption of short-wave radiation in an optically thick layer of snow leads to the strong multiple light scattering. As a result, it is sufficient to use the simplest transport model for the single scattering (Dombrovsky 2012): Φ λ ( r , μ 0 ) = ( 1 − μ ¯ λ ( r ) ) + 2 μ ¯ λ ( r ) δ ( 1 − μ 0 ) (2) where μ 0 = Ω ′ × Ω is the cosine of the angle of scattering, δ is the Dirac delta-function, and μ ¯ λ is the asymmetry factor of scattering: μ ¯ λ ( r ) = 1 2 ∫ − 1 1 μ 0   Φ λ ( r , μ 0 ) d μ 0 (3)

Of course, the normalization condition 1 2 ∫ − 1 1 Φ λ ( μ 0 )   d μ 0 = 1 is satisfied by function (2). The transport model proved to be sufficiently accurate for many problems of radiative transfer in weakly absorbing and strongly scattering media (Dombrovsky and Baillis 2010; Dombrovsky and Lipiński 2010). Interestingly, the computational results obtained using the transport model are normally very close to those obtained with the more complex Henyey–Greenstein analytical model (Dombrovsky et al., 2011a; Dombrovsky et al., 2013). This statement is also correct in the case under study.

The transport form of the RTE, which follows from Eqs 1–3, is similar to that for a hypothetical isotropic scattering (Dombrovsky 2012; Dombrovsky 2019): Ω ∇ I λ ( r , Ω ) + β λ tr I λ ( r , Ω ) = σ λ t r 4 π G λ ( r ) , G λ ( r ) = ∫ ( 4 π ) I λ ( r , Ω )   d Ω (4) where G λ is the irradiance and the modified coefficients are: σ λ t r = σ λ ( 1 − μ ¯ λ ) , β λ t r = α λ + σ λ t r = β λ − σ λ μ ¯ λ (5)

Consider a semi-infinite plane-parallel layer of a scattering medium, such as snow, illuminated directly by solar radiation. In the case of oblique incidence, the radiation field in the medium is three-dimensional. Fortunately, the 3D problem can be radically simplified. Let us note that the irradiance field in snow at illumination at the same angle relative to the normal to the surface, but at different azimuthal angle differs only in the rotation relative to the normal. As a result, the irradiance field is the same as in the case of uniform radiation along a conical surface. This cone with a vertical axis is formed by solar rays with the same zenith angle. Hence, we can solve the equivalent axisymmetric problem as it was done by Dombrovsky et al. (2011b), Dombrovsky and Randrianalisoa (2018) and Dombrovsky et al. (2019). The axial symmetry allows us to integrate the RTE over the azimuthal angle: μ ∂ I λ , i ∂ z + β λ t r I λ , i = σ λ t r 2 G λ , i , G λ , i ( z ) = ∫ − 1 1 I λ , i ( z , μ )   d μ , μ = cos ⁡ θ , z > 0 (6a) I λ , i ( 0 , μ ) = I λ sol ( μ i )   δ ( μ i − μ ) , I λ , i ( ∞ , − μ ) = 0 , μ , μ i > 0 (6b) where μ i = cos ⁡ θ i , θ i is the incidence angle, I λ s o l is the solar radiation intensity. The diffuse radiation from the sky is relatively small, and it is neglected in boundary conditions (6b).

It is important that the radiative transfer equation is linear. This allows us to separately consider the contribution of direct radiation from the Sun and radiation scattered by the atmosphere. For the direct irradiation we have: μ ∂ I λ , i ∂ τ λ t r + I λ , i = ω λ t r 2 G λ , i , G λ , i ( τ λ t r ) = ∫ − 1 1 I λ , i ( τ λ t r , μ )   d μ (7a) I λ , i ( 0 , μ ) = I λ sol ( μ i ) δ ( μ i − μ ) , I λ , i ( ∞ , − μ ) = 0 , μ , μ i > 0 (7b) where τ λ t r ( z ) = ∫ 0 z β λ t r ( z )   d z is the dimensionless optical coordinate, ω λ t r ( z ) = σ λ t r ( z ) / β λ t r ( z ) is the transport albedo. As usual, the total values of radiation intensity and irradiance are: I λ , i ( τ λ t r , μ ) = J λ , i ( τ λ t r , μ ) + I λ s o l ( μ i ) E λ , i δ ( μ − μ i ) , E λ , i = exp ( − τ λ t r / μ i ) (8) G λ , i ( τ λ t r ) = G λ , i d i f ( τ λ t r ) + I λ s o l ( μ i ) E λ , i , G λ , i d i f = ∫ − 1 1 J λ , i   d μ (9)

The radiative transfer problem for the diffuse radiation component can be written as: μ ∂ J λ , i ∂ τ λ t r + J λ , i = ω λ t r 2 ( G λ , i d i f + I λ s o l ( μ i ) E λ , i ) (10a) J λ , i ( 0 , μ ) = J λ , i ( ∞ , − μ ) = 0 , μ > 0 (10b)

This problem can be further simplified because the source function in Eq. 10a is independent of μ and the two-flux (Schwarzschild–Schuster) model is applicable. The detailed analysis of this model in (Dombrovsky and Baillis 2010) showed its sufficiently high accuracy. For the considered problem of radiative transfer in a snowpack this approach was proved to be also applicable (Dombrovsky et al., 2019). The two-flux method is based on the following approximation of the diffuse radiation intensity: J λ , i ( τ λ t r , μ ) = { J λ , i − ( τ λ t r )   when   − 1 < μ < 0 J λ , i + ( τ λ t r )   when   0 < μ < 1 (11)

Integrating Eq. 10a, one can obtain the boundary-value problem for the irradiance G λ , i d i f = J λ , i − + J λ , i + (Dombrovsky and Baillis 2010). It the case of vertically homogeneous snow, one can obtain the analytical solution: G λ , i d i f = 4 ω λ t r ξ λ 2 − 1 / μ i 2 ( E λ , i − 2 + 1 / μ i 2 + ξ λ E λ d i f ) I λ s o l ( μ i ) (12a) E λ d i f = e x p ( − ξ λ τ λ t r ) , ξ λ = 2 1 − ω λ t r ≠ 1 / μ i (12b) where μ i is a parameter. There are two different exponential functions in Eq. 12a: E λ , i is related to the not-scattered collimated radiation, whereas E λ d i f corresponds to the diffuse radiation. It is not difficult to see that the depth of propagation of collimated radiation is smaller than that of the diffuse one when ω λ t r > 0.75 . This makes clear the evolution of the radiation field with the depth. After entering a strongly scattering medium, the directional solar radiation in the vicinity of snow surface is converted into diffuse (but not isotropic) radiation. Due to multiple scattering, absorption of diffuse radiation is enhanced in a medium with a very small absorption coefficient. As a result of this absorption, the radiation flux decreases, and with it the absorbed power of radiation also decreases. The absorbed radiation power is determined using the calculated irradiance profiles for direct solar radiation: P ( z ) = ∫ λ m i n λ m a x p ( λ , z ) d λ , p ( λ , z ) = α λ G λ , i ( z ) (13) where λ min and λ max are the boundaries of the main wavelength range that contributes to the radiative transfer in the layer.

In the case of an optically thick homogeneous medium layer, the position of the maximum (along the optical coordinate τ λ t r ) of the spectral component of the absorbed radiation energy, p ( λ , τ λ t r ) , inside the layer is determined from the condition d G ¯ λ , i d τ λ t r = 0 , where G ¯ λ , i ( τ λ t r ) = E λ , i + 4 ω λ t r ξ λ 2 − 1 / μ i 2 ( E λ , i − 2 + 1 / μ i 2 + ξ λ E λ d i f ) (14)

It the case of a strongly scattering medium, it follows that the maximum of the absorbed radiation power occurs at: τ λ , m a x t r = μ i ln ( 2 μ i − 1 ξ λ ( 2 − ξ λ ) ) when ξ λ ≪ 1 / μ i (15)

The value ω λ t r was excluded using the formula 4 ω λ t r = 4 − ξ λ 2 . Equation 15 makes sense only when μ i > 0.5 ( θ i < 60 ° ) , but this restriction suits us fine. The formal use of this equation may give zero or negative values of τ λ , m a x t r . The latter means that there is no maximum of absorbed power inside the medium layer at the angle of incidence less than 60 ° (e.g., at the Sun in the zenith). The condition of τ λ , m a x t r = 0 gives the boundary of the region of ( ω λ t r , μ i ) where the maximum of the absorbed radiation power exists at some distance under the irradiated surface of the medium layer. This boundary can be obtained analytically as a solution for the second-order algebraic equation: ξ λ 2 − 2 ξ λ − ( 2 μ i − 1 ) = 0 (16)

The resulting solution is as follows: ξ λ = 1 − 2 ( 1 − μ i ) (17)

Obviously, the deepest penetration of the collimated radiation into the medium takes place at the normal incidence ( μ i = 1 ) when Eq. 17 gives ξ λ = 1 and ω λ t r = 0.75 . The large boundary value of transport albedo shows that the internal maximum of absorption takes place for highly-scattering media only. The results of calculations in the most interesting range of 0 ≤ θ i ≤ 45 ° are presented in Figure 1, where the curve ω λ t r ( θ i ) separates the regions corresponding to the absorption maximum directly at the surface of the illuminated medium and below this surface. Interestingly, the absorption maximum can occur only at a distance from the irradiated surface when the angle of incidence is less than 60^o even the case of a weakly scattering medium.

It is also interesting to calculate the ratio of the maximum local absorption of radiation to that on the surface of the medium illuminated along the normal: p ¯ m a x ( ω λ t r ) = ( 1 + ζ ) ( 2 + ξ λ ) e x p ( − τ λ , m a x t r ) − 3 ζ e x p ( − ξ λ τ λ , m a x t r ) 2 + ξ λ − ζ ( 1 − ξ λ ) , ζ = 4 ω λ t r / ( ξ λ 2 − 1 ) (18) where τ λ , m a x t r is determined by Eq. 15. The increase at the optical depth τ λ , m a x t r where the absorption maximum is located, as well as the relative value of p ¯ m a x of this maximum for media with high values of the transport albedo of single scattering ω λ t r is shown in Figure 2. It can be seen that the effect of deep absorption of radiation is most pronounced for media with ω λ t r > 0.99 . For example, when ω λ t r = 0.999 9, the absorption maximum is at the optical depth τ λ , m a x t r = 3.23 and its relative value p ¯ m a x = 1.85 . This seemingly exotic situation takes place for pure snow in the visible part of the solar spectrum. Against this background, the application of Bouguer’s law in some publications on radiation propagation in scattering and weakly absorbing media looks especially inappropriate.

Figure 3 shows several profiles of the absorption of external radiation in highly scattering media illuminated along the normal. The absorption of radiation is large at optical depths much greater than the depth of the absorption maximum. For example, at ω λ t r = 0,9999 even at τ λ t r = 30 the absorbed power of radiation is higher than near the illuminated surface.

To calculate the absorption of solar radiation in snowpack, one needs the spectral local optical properties of snow. These properties, in turn, can be determined only with the known spectral optical constants of ice as the only substance for pure snow. The optical constants are considered as real and imaginary parts of the complex index of refraction, m ( λ ) = n ( λ ) − i κ ( λ ) , where n is the index of refraction and κ is the index of absorption. The optical constants of ice obtained by Warren and Brandt (2008) are plotted in Figure 4. The presented wavelength range is wider than that required to solve the problem of radiative heating and includes the near-ultraviolet spectral range considered in (Dombrovsky et al., 2022). Note that the measurements by Picard et al. (2016) showed larger values of κ in the wavelength range of 0.39 < λ < 0.6 µm than those obtained by Warren and Brandt (2008) (see Figure 4). In this paper, for the first time, we compare the calculated data on the transfer of shortwave radiation in the snow layer for the different dependences of κ ( λ ) in the visible.

The transition from very small values of the absorption index in the visible range to its large values in the near-infrared determines quite different optical properties of ice and snow in these spectral ranges. In particular, the high value of snow albedo in the visible range is a result of high transparency of pure ice. Note that even very small impurities such as aerosol particles may strongly affect the albedo of snow (Kokhanovsky, 2021c). In this regard, one can recommend the following recent papers: Kokhanovsky et al. (2018), He and Flanner (2020), Kokhanovsky (2021a), Pu et al. (2021). Note that the optical properties of snow polluted with atmospheric dust or soot particles (He et al., 2018; Dombrovsky and Kokhanovsky 2020a; Shi et al., 2021) do not depend on the choice between the absorption data for pure snow in the short-wave range. However, in this paper, we consider only pure snow, for which this choice may matter.

Modeling of the optical properties of various scattering media is often based on classical Mie theory (Bohren and Huffman 1998). However, in many cases one can use the geometric optics (GO), which works well for large particles, as well as some combined methods (Bi et al., 2011; Lindqvist et al., 2018; Kokhanovsky, 2021c). Such a possibility is important for ice grains of complex shape (Libois et al., 2013; Liou and Yang 2016; Ishimoto et al., 2018). Interestingly, in some cases reliable results can be obtained from physically sound approximate solutions (Kokhanovsky 2021b).

Dombrovsky et al. (2019) have considered the spherical ice grains of different size instead of nonspherical ice particles in random orientation. It was assumed that the complex shape of the ice grains has little effect on solar heating of the snowpack. The different sizes of ice grains were treated as those corresponding to different snow morphology. It was shown that the monodisperse model can be used instead of the detailed calculations for polydisperse grains. Therefore, we restrict ourselves to the calculation results for the case of monodispersed spheres with a = 100 µm. A comparison with the Mie theory in (Dombrovsky et al., 2019) showed that the GO approximation is very accurate and one can use a relatively simple analytical solution derived in the early paper by Kokhanovsky and Zege (1995). In some cases, it is convenient to use the simplest analytical approximations for the main optical properties of weakly absorbing spherical particles suggested in (Dombrovsky 2002) and in the recent paper by Kokhanovsky (2021b). In particular, following (Dombrovsky et al., 2019) we use simpler relations in the region of weakest absorption, where a formal application of a more general solution yields erroneous results.

For calculations of radiation transport in the transport approximation we need only two dimensionless characteristics of absorption and scattering of light by individual particles: the absorption efficiency factor, Q a , and the transport efficiency factor of scattering, Q s t r = Q s × ( 1 − μ ¯ ) , where μ ¯ is the asymmetry factor of scattering. Note that the values of Q a and Q s are introduced as absorption or scattering cross sections divided by π a 2 , where a is the particle radius (Bohren and Huffman 1998). We consider also the transport efficiency factor of extinction, Q t r = Q a + Q s t r , and transport albedo, ω t r = Q s t r / Q t r . It is known that these characteristics depend not only on optical constants but also on the dimensionless diffraction parameter x = 2 π a / λ .

Calculations show that light scattering significantly dominates over absorption in the visible range (Figure 5). This leads to the reflection of most of the shortwave solar radiation from the snow. At the same time, the remaining part of the solar radiation is scattered in a snowpack (Kokhanovsky 2021b; 2021c). As a result, visible light contributes significantly to heating the deep layers of snow, while infrared radiation from the Sun is absorbed in a relatively thin surface layer.

Note that the hypothesis of independent scattering by single ice grains in snow is true (Mishchenko 2018). It is really the case because of a high porosity of snow and large sizes of ice grains. Therefore, the optical properties of snow can be calculated as follows (Dombrovsky and Baillis 2010): { α λ , σ λ t r , β λ t r } = 0.75 f v { Q a , Q s t r , Q t r } / a (19)

The transport albedo of scattering for snow coincides with the same value for a single ice grain: ω λ t r = ω t r . Returning to Figure 5, we note that the calculated value of the snow albedo in the short-wave region is very close to unity for both the absorption data (Warren and Brandt 2008) and the data (Picard et al., 2016). For example, at wavelength λ = 0.4 µm in the first case we have ω λ t r = 1.4 × 10 − 5 and in the second case – 1 − ω λ t r = 2.9 × 10 − 5 . This means that regardless of the initial absorption data, we are dealing with a deep and strong absorption maximum for the visible radiation.

In this section, we limit ourselves to calculating the absorption of direct solar radiation, neglecting the relatively small contribution of radiation scattered by the atmosphere. The vertical profile of the normalized absorbed power is illustrated in Figure 6 for the normal incidence separately for the visible and near-infrared spectral ranges. It is not necessary to consider the real spectrum of solar radiation that has passed through the atmosphere in order to make a qualitative estimate. Instead, we will use the Planck spectrum at temperature T s o l = 6000 K . It is interesting to calculate separately the contributions of the visible and the near-infrared ranges to the absorbed power of radiation: P ¯ ( z ) = P ¯ v i s ( z ) + P ¯ n i r ( z ) , P ¯ v i s , n i r = P v i s , n i r / ∫ λ m i n λ m a x I b ( T sol , λ )   d λ (20) where λ m i n = 0.4 µm and λ m a x = 2.8 µm and   P v i s , n i r is calculated using Eq. 13 for the spectral ranges 0.4–0.78 and 0.78–2.8 µm, respectively. Of course, one can use more accurate values for the conventional solar temperature and the boundaries of the visible and near-infrared ranges as it was recommended in the monograph by Kokhanovsky (2021c). However, it is not important for the presented estimations. According to (Dombrovsky et al., 2019), we introduce the normalized coordinates f v z and P ¯ / f v in Figure 6A, where f v is the volume fraction of ice is a snow layer.

As one might expect, infrared radiation is absorbed mostly in a thin surface layer, while visible light is absorbed almost uniformly in a relatively thick layer. The result is in good agreement with the work of Munneke et al. (2009), which discusses the deep penetration of shortwave solar radiation and its role in the thermal balance of the snow cover.

For better clarity, Figure 6B shows similar dependences for the two parts of the visible spectrum: 0.4 < λ < 0.6 µm (A) and 0.6 < λ < 0.78 µm (B). In the calculations for the sub-range A, different data on the ice absorption index are used, which give similar dependences of the absorbed power on depth, differing slightly only in the absolute value. As we have already discussed, the position and relative magnitude of the absorption maximum are the same for the two variants of the ice absorption index in the visible range. The calculations using the data by Picard et al. (2016) give a slightly greater absorption of the shortwave radiation.

Of course, knowing the local value of the absorbed radiation power is not enough to calculate the snow temperature. The next step in solving the problem requires the transient energy equation with appropriate initial and boundary conditions. The wide-range spectral dependence of κ ( λ ) as given by Warren and Brandt (2008) will be used in these calculations.

Consider now the 1D problem of heat transfer. As in (Dombrovsky and Kokhanovsky 2021), we will not consider such processes as ice sublimation and water vapor diffusion in the snow layer, which may be of interest for studying the changing microstructure of snow. The problem statement for the transient temperature field, T ( t , z ) , in a layer of snow is: ρ c ∂ T ∂ t = ∂ ∂ z ( k ∂ T ∂ z ) + P , t > 0,0 < z < d t h (21a) T ( 0 , z ) = T 0 ( z ) (21b) z = 0 , − k ∂ T ∂ z = h   ( T a i r − T ) − ε   π ∫ λ w 1 λ w 2 I b ( T , λ )   d λ z = d t h , ∂ T ∂ z = 0 (21c)

where z is the coordinate measured from the irradiated surface, ρ , c , and k are the density, the specific heat capacity, and the thermal conductivity of snow or ice, T a i r ( t ) is the temperature of ambient air (outside the thermal boundary layer), and h ( t ) is the convective heat transfer coefficient. The adiabatic condition at the boundary z = d t h > > d means that we neglect heat transfer at z > d t h . Of course, the value of d t h increases with time and should be estimated using additional calculations. The heat transfer coefficient is determined by a wind speed. The last term in the right-hand side of the boundary condition at z = 0 is the mid-infrared radiative cooling due to thermal radiation of snowpack or ice sheet in the atmospheric transparency window of λ w 1 < λ < λ w 2 ( λ w 1 = 8 µm, λ w 2 = 13 µm) (Hossain and Gu 2016). The radiative cooling is limited during the night time when it is not compensated by solar mid-infrared radiation. Therefore, the coefficient ε varies from zero in the day time to the unity at night. The absorbed radiation power, P ( z ) , at arbitrary conditions of solar illumination should be recalculated during the combined problem solution.

The latent heat of ice melting was accounted for by an equivalent increase in the heat capacity in a narrow temperature range near the melting temperature. This technique was successfully used in (Dombrovsky et al., 2009; Dombrovsky et al. 2015; Dombrovsky et al. 2019). It is interesting that the three cited papers considered completely different problems: the first of them dealt with the solidification of core melt droplets in a severe nuclear reactor accident, the second one with the freezing and thawing of biological tissues in cryosurgery, and the third one with the melting of freezing snow under periodic solar heating.

In numerical calculations, an implicit finite-difference scheme of the second order of approximation was used. Note that the model should be generalized in the case of significant snow melt or rain over snow, as water penetrates through the snow and solidifies at a greater depth.

We used the algorithm and input data from (Dombrovsky et al., 2019) to calculate the solar heating of snow. Because of the considerable uncertainty in the snow thermal conductivity (Arenson et al., 2015), there may be a systematic error. At the same time, changes in solar illumination conditions over time in the summer solstice and the local altitude of 70^o were calculated accurately and the light scattered by cloudless atmosphere was taken into account as it was done by Dombrovsky et al. (2019).

As an example, we consider the case of windless weather, when the heat transfer coefficient can be assumed to be h = 6 W/(m² K) (Defraeye et al., 2011; Mirsadeghi et al., 2013). The temperature profiles in a snowpack calculated for the conventional first day of solar irradiation at T 0 = T a i r ( 0 ) = − 10 °C are presented in Figure 7.

The calculations showed that the melting temperature is reached at noon near the surface (it is not shown in the figure) and temperature maximum is shifted to the depth about 5 mm at 6 p.m. due to convective cooling (see Figure 7). Further moving of the temperature maximum to the larger depth is accompanied by cooling of the surface layer of the snow cover. Convective cooling predominates in the evening hours, and radiative cooling is strongly manifested closer to midnight. Obviously, radiative cooling under cloudless skies should be taken into account during the period of weak solar irradiation.

Following (Dombrovsky et al., 2019), similar calculations can be made for solar heating of snow, for example, over 3 days. However, the calculation of snow cover heating over many days and, even more so, over a long polar summer is associated with serious difficulties and requires the development of a special physical model. This task is beyond the scope of the present study. A simple physical estimate using the Fourier number for snow shows that it takes about 10 years for snow to warm to a depth of 100 m. Note that the heating time is directly proportional to the square of the depth and the solar heat “stored” at depth during the summer continues to spread deep into the depths even during the polar winter.

Of course, the real process is more complicated and is not reduced to simple heat conduction. It is interesting to discuss briefly the main stages of the process for the case of a very thick snow layer. Most likely, at a not very deep depth, the snow will begin to melt. The resulting water will flow down between the snow particles. This flow is very slow because of the high hydraulic resistance of the dispersed medium with small pores and water will freeze at a slightly deeper depth, where the snow temperature is noticeably lower. This is how the ice layer will start to form, probably with gas bubbles scattering the radiation (Dombrovsky and Kokhanovsky 2020b). Solar heating in several years increases the thickness of the ice layer. Over many years, conditions may develop in which ice melting will lead to the formation of a gas cavity and a lake under the ice sheet. It is possible that deep solar heating is one of the reasons (along with ice melting under pressure and geothermal heat) for the formation of some of the many discovered lakes under the Greenland and Antarctic ice sheets.

A theoretical analysis of the conditions under which the absorption of collimated radiation penetrating into a weakly absorbing scattering medium has a maximum at a considerable distance from the illuminated surface was presented for the first time. It is shown that the maximum absorption of radiation inside an optically thick medium exists only at illumination angles less than 60° from the normal. Both the depth at which the local absorption of solar radiation is maximum and the magnitude of this maximum are greatest at illumination along the normal. The greatest optical depth of the maximum location (by transport attenuation coefficient) is 3.23. In this case, the maximum local power of absorbed radiation energy is 85% higher than that at the illuminated surface of the medium. The obtained approximate analytical solution is general and applicable to various scattering media.

The possibility of deep heating of the pure snow cover by solar radiation during the polar summer is confirmed by calculations based on a general physical model, which is not limited to determining the solution for the absorbed radiation power. The transient problem of combined heat transfer was solved taking into account the change in the zenith angle of the Sun with time as well as the heat conduction in snow and both convective and radiative cooling of the snow surface.

The obtained numerical solution can be used for studies of the propagation of heat into the depth of the snowpack. It was shown that the computational results are not very sensitive to the choice between the published data for the absorption index of pure ice in the visible range. It is interesting that by the evening of the first sunny day, the snow temperature at a depth of 4 mm becomes about 2° higher than on the surface, and by midnight this temperature difference can reach 6°, while the zone of maximum heating shifts to a depth of about 12 cm. This process continues during the long polar summer, leading to heating and possible melting of snow at a considerable depth.

The calculations have shown that the radiative cooling of snow in the transparency window of the cloudless atmosphere should be taken into account during periods of relatively low solar illumination, even during the polar summer.