COMPUTATIONAL MODELS FOR COMBINED HEAT TRANSFER IN SNOW AND ICE
Following from: Spectral optical properties of pure ice and snow
The transport radiative transfer equation (RTE) for a flat layer of isotropic snow can be written as follows (Dombrovsky and Baillis, 2010; Dombrovsky, 2019):
where Iλ is the spectral radiation intensity at point r⃗ in direction Ω⃗, βλ = αλ + σλ is the extinction coefficient, and αλ and σλ are the coefficients of absorption and scattering, respectively. The physical meaning of Eq. (1) is evident: variation of the radiation intensity takes place due to extinction and the scattering from other directions.
The thermal radiation of snow and ice should also be taken into account, but only in the infrared part of the spectrum, outside the near-infrared range. Both snow and ice are opaque in the middle-infrared range, and thermal radiation is emitted only by the surface of the snow cover or ice sheet. Therefore, the intrinsic thermal radiation of the medium does not appear in the radiative transfer equation.
The multiple scattering of solar radiation takes place in snow and scattering ice, which contains numerous small gas bubbles. This makes it possible to use the so-called transport approximation for the single-scattering phase function (Dombrovsky, 2012). In this case, it is sufficient to know the asymmetry factor of scattering (the average cosine of the scattering angle), which is independent of polarization:
The transport RTE can be written as follows:
where Gλ is the irradiation function and the transport scattering and extinction coefficients are introduced as follows:
An accurate numerical solution to the RTE in scattering media is a very complicated task, even in the case of transport approximation for the single-scattering phase function. One can find a number of studies in the literature on specific numerical methods developed to obtain more accurate spatial and angular characteristics of the radiation intensity field. Several modifications of the discrete ordinates method, the finite-volume method, and statistical Monte Carlo methods are the most popular tools employed by many authors. However, simple and physically clear differential approximations are widely used for solving the radiative transfer problems in scattering media, particularly in combined heat transfer problems (Dombrovsky, 2019).
The angular dependence of the radiation intensity is the main difficulty in solving the RTE. Fortunately, this angular dependence is rather simple in many applied problems. This enables one to derive the so-called differential approximations. These approximations have a long history. The latter is reflected in their well-known names, such as the Eddington method, Schwarzschild–Schuster method, etc. It should be emphasized that all the differential approximations are based on one or another simple presentation of angular dependence of the radiation intensity. As a result, it is sufficient to consider a limited number of coordinate functions and turn to the coupled ordinary differential equations by the use of integration of RTE.
The simplest differential approximations, brought together by the general term “diffusion approximation,” give the following representation of the spectral radiative flux:
and differ only by the expression for the radiation diffusion coefficient Dλ. It is known that Eq. (5) can be also derived based on some assumptions concerning the angular dependence of radiation intensity. Substituting Eq. (5) into the radiation energy balance
we obtain the modified Helmholtz equation for the spectral irradiation
Note that the radiation power absorbed in a medium is expressed as follows:
To clarify the physical sense of the diffusion approximation, consider the linear angular dependence of radiation intensity in the popular P1 approximation
By multiplying the transfer equation (3) by Ω⃗ and integrating it over a solid angle, one can find that the spectral radiative flux is simply related to the spectral irradiation
Note that Eq. (10) is obtained for an arbitrary scattering function, but it is the same as that for transport approximation.
The two-flux approximation is also often used to solve 1D radiative transfer problems, including those considered in the present article. Consider the transport RTE
where 0 < z < d is the coordinate across the medium layer, d is the layer thickness, and θ is the angle measured from the z-axis direction. The spectral radiative flux is obtained as follows:
Note that Gλ = Gλ– + Gλ+. With the use of two-flux approximation, the radiation intensity is assumed to be constant over the backward and forward hemispheres
The integration of RTE over the hemispheres leads to the following coupled equations:
These equations and the expression for the spectral radiative flux can be written as follows:
where Dλ = 1/(4βλtr). One can see that two-flux approximation leads to the typical relations of the diffusion approximation. We will not discuss the boundary conditions for the P1 and two-flux equations. The details of these approaches, including the accuracy analysis based on a comparison with reference solutions for the model problems, can be found in Dombrovsky and Baillis (2010). However, it is important to recall that P1 is insufficiently accurate at the medium boundaries in the case of discontinuous angular dependence of radiation intensity. Therefore, the two-flux approximation is preferable when the solar light irradiates a snowpack or ice sheet.
Consider first the radiative transfer in snow. This problem is relatively simple because there is no refracting host medium. The solution considered in this section is applicable to the snowpack irradiated by collimated solar radiation and also by diffuse radiation from the sky. In the case of oblique solar irradiation, the radiation field in a snow layer is three-dimensional. However, the distribution of the absorbed power over the thickness of a medium layer depends only on the solar zenith angle and it does not matter which side of the normal the Sun is located. Let us imagine that the radiative flux is uniformly distributed over the surface of a cone with the same zenith angle at the vertex. Of course, in this case the power profile of the absorbed radiation will be the same. In other words, the original problem is equivalent to the 1D axisymmetric problem of the oblique irradiation along the cone surface. This idea has been used by Dombrovsky et al. (2019). The axial symmetry of the equivalent problem enables one to integrate the original transport RTE over the azimuth angle. The resulting equations can be written as follows:
where μi = cos θi corresponds to the direction of incident collimated light, Iλsol is the intensity of solar radiation transmitted through the atmosphere, Iλsky is the intensity of diffuse radiation from the sky, and d is the thickness of a medium layer. The conditions of solar irradiation depend on current time, t, but the value of t is just a parameter of the radiative transfer problem.
The linearity of the RTE makes possible a separate consideration of additive contribution of the direct solar radiation and radiation from the sky. Consider first the transfer of direct solar radiation using the dimensionless form of Eqs. (16a) and (16b):
There are some simplifications in these equations.The boundary condition at z = d is replaced by the same condition at τλtr → ∞ (this can be done for optically thick layers), and the following dimensionless variables are introduced: the ratio of I¯λ,i = Iλ /Iλsol(μi), the current optical thickness τλtr = ∫0z βλtr(z) dz, and transport albedo of single scattering ωλtr(z) = σλtr(z)/βλtr(z).
Following the usual technique, the radiation intensity is represented as a sum of the diffuse component J¯λ,i and the term corresponding to the direct solar radiation:
The spectral irradiance can be also written as a sum of two components
The mathematical problem statement for the diffuse component of radiation intensity is as follows:
The source term on the right-hand side of Eq. (20a) does not depend on angular variable μ. This enables further simplification of the problem with the use of two-flux approximation.
For the diffuse radiation from the sky, one should use another normalization, I¯λ = Iλ /qλsky and replace the term δ(μi – μ) by 1 in Eq. (17b). As a result, the following equations are true:
The spectral component of the absorbed radiation power taking into account both the direct solar radiation and the diffuse radiation of the sky, can be written as follows:
According to the two-flux method, the following approximation of the diffuse component of radiation intensity is considered:
Integrating Eq. (20a) separately over the intervals –1 < μ < 0 and 0 < μ < 1, one can obtain the boundary-value problem for the diffuse irradiance G¯λ,idiff = J¯λ,i– + J¯λ,i+ :
Equations (24a) and (24b) are true at arbitrary variations of ξλ and ωλtr with the current optical thickness τλtr, and they can be easily solved numerically. However, it is interesting to obtain an analytical solution of the problem in the case of uniform optical properties of snow. The analytical solution to Eqs. (24a) and (24b) at ξλ ≠ 1/μi can be written as follows (Dombrovsky and Randrianalisoa, 2018; Dombrovsky et al., 2019):
where μi should be considered as a parameter. There are two different exponential functions in Eq. (25). The first one, Eλ,i, is related with a contribution of the collimated solar radiation, whereas the second one, Eλdiff, corresponds to the diffuse component of the radiation field.
The following simple expression can be easily obtained for the diffuse radiation of the sky:
A comparison of functions Eλ,i(τλtr) and Eλdiff(τλtr) indicates that the propagation depth of collimated radiation (defined as a distance from the snow surface at which the irradiance decreases e times) is less than the propagation depth of diffuse radiation component when ωλtr > 0.75.
In the case of a semi-transparent scattering ice sheet, the abovementioned radiative transfer model should be modified to take into account the refraction of solar light in ice. The generalized boundary condition leads to the following equations for the direct solar radiation instead of Eqs. (16a) and (16b):
where I¯λ,i = Iλ,i / Iλ,iinc, G¯λ,i = Gλ,i / Iλ,iinc, and Iλ,iinc is the incident radiation intensity in direction μi, rλ,i = rλ(n,μi) is the Fresnel reflection coefficient (Born and Wolf, 1999), and μj is the directional cosine of the angle of refraction. The refraction leads to the change of direction of the incident solar radiation at the medium boundary. The new direction of the refracted solar light under the horizontal surface of a refracting medium is . In particular, sunlight at sunrise and sunset propagates in a horizontal layer of a refracting and nonscattering medium along the surface of a cone with a half-angle at the apex θ = arccos (1 – 1/n2). Note that simple equation of rλ,i = rλ,n = (n – 1)2/(n + 1)2 can be used instead of the general Fresnel equations in the case of normal irradiation (μi = 1) of a weakly absorbing host medium, when κ ≪ n.
As before, the radiation intensity is presented as a sum of the diffuse and collimated components:
and the problem statement for the diffuse component is as follows:
According to the modified two-flux approximation (Dombrovsky et al. 2006), the angular dependence of the radiation intensity takes into account the total internal reflection at τλtr = 0:
The term “two-flux” is applicable to this approximation because the interval of –μc < μ < μc does not contribute to the radiative flux. Integrating Eq. (27a) over the angles, one can obtain the boundary-value problem for the normalized diffuse irradiance:
where r¯λ is the angle-averaged reflectance for the diffuse radiation at the interface τλtr = 0.
The analytical solution to the Eqs. (31a)–(31c) at ξλ ≠ 1/μj is:
In the case of a nonrefracting medium (rλ,i = 0, μj = μi, and γ = 1), Eqs. (31a)–(31c) and (32) were used by Dombrovsky et al. (2019) and Dombrovsky and Kokhanovsky (2019) for solar radiation transfer in a snowpack.
It should be recalled that a comparison to the exact numerical results obtained using the high-order composite discrete ordinates method performed by Dombrovsky et al. (2006) and also a comparison to the Monte Carlo simulation showed that the modified two-flux approximation is sufficiently accurate for the use in diverse applications (Dombrovsky, 2019).
The ordinary or modified two-flux approximation is usually a good approach for estimating the irradiation field. With the use of transport approximation, the irradiation field is sufficient to obtain the source function (the right-hand side) of RTE. As a result, there is a possibility to improve the solution using the second iterative step of a computational procedure. It is sufficient to solve the RTE with the known source function using one of the ray-tracing procedures. A study of applicability and accuracy of the combined two-step procedure with the use of the Monte Carlo method at the second step of the iterative solution has been reported by Dombrovsky and Baillis (2010).
Consider now the 1D heat transfer problem. Generally speaking, there are various thermal or related processes in a snowpack or scattering ice sheet and many of these processes should be involved in a complete transient computational model for snow heating. As an example, one should recall the processes of ice sublimation and diffusion of water vapor through a snow layer. These processes appear to be important for the snow microstructure, which determines macroscopic snow properties. However, this is beyond the scope of this paper.
The 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:
where z is the coordinate measured from the irradiated surface of the layer; ρ, c, and k are the density, the specific (per unit mass) heat capacity, and the thermal conductivity of snow or ice; Tair(t) is the temperature of ambient air; and h(t) is the convective heat transfer coefficient. The adiabatic condition at the boundary z = dth ≫ d means that we neglect heat transfer at z > dth. Of course, the value of dth increases with time and should be estimated using additional calculations. The time dependence of heat transfer coefficient is determined by wind speed.
The last term on 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 surface in the atmospheric window of λw1 < λ < λw2 (λw1 = 8 μm, λw2 = 13 μm) (Hossain and Gu, 2016). The radiative cooling is limited during the nighttime when it is not compensated by solar mid-infrared radiation. Therefore, the coefficient ϵ varies from zero in the daytime to the unity at night. The absorbed radiation power, P(z), at arbitrary conditions of solar irradiation should be recalculated time by time during the combined problem solution.
Of course, the great value of the latent heat of ice melting, L = 0.34 MJ/kg, should be taken into account. This can be done using an equivalent additional heat capacity, Δc, in narrow temperature range of Tm – ΔT < T < Tm + ΔT (where Tm = 273 K is the ice melting temperature, ΔT ≪ Tm – T0). The following simple variant of this approach was used by Dombrovsky et al. (2019):
A correct choice of ΔT value is determined by the interval of the computational grid and the time step of numerical solution to Eqs. (33a)–(33c). Obviously, it is difficult to choose a realistic initial profile of temperature, T0(z), for the heat transfer calculations. Fortunately, the effect of this temperature profile decreases with time. It was shown by the authors that the choice of T0(z) makes no difference for the snow temperature after about 4 hr from the conventional initial time moment. An implicit finite-difference scheme of the second order of spatial approximation can be recommended for numerical calculations. The specific combination of a strong decrease of the heat-generation rate with the distance from the irradiated surface and rather deep propagation of heat in the snow or ice layer due to longtime conduction makes it necessary to use a detailed nonuniform computational grid and double precision calculations.
Note that the abovementioned model of heat transfer is incomplete in the case of a considerable melting of snow, because water penetrates downward through porous snow and solidifies there. This process should be taken into account in a more general, physical and computational model. The accompanying variation of optical properties of a nonuniform snowpack should be also analyzed on the basis of more sophisticated models.
It is important that the heat transfer problem considered is not a classical problem of radiative-conductive heat transfer (Dombrovsky and Baillis, 2010). In the case of solar heating of a snow cover, the snow temperature does not directly affect the radiative transfer. This simplifies significantly the solution of the problem, since a single spectral calculation of the radiative transfer is sufficient.
Born, M. and Wolf, E. (1999) Principles of Optics, 7th (expanded) ed., New York: Cambridge University Press.
Dombrovsky, L.A. (2012) The Use of Transport Approximation and Diffusion-Based Models in Radiative Transfer Calculations, Comput. Therm. Sci., 4(4): 297–315.
Dombrovsky, L.A. (2019) Scattering of Radiation and Simple Approaches to Radiative Transfer in Thermal Engineering and Bio-Medical Applications, In A. Kokhanovsky (Ed.), Springer Series in Light Scattering, vol. 4, Cham, Switzerland: Springer, pp. 71–127.
Dombrovsky, L.A and Baillis, D. (2010) Thermal Radiation in Disperse Systems: An Engineering Approach, Danbury, CT: Begell House.
Dombrovsky, L.A. and Randrianalisoa, J.H. (2018) Directional Reflectance of Optically Dense Planetary Atmosphere Illuminated by Solar Light: An Approximate Solution and Its Verification, J. Quant. Spectrosc. Radiat. Transf., vol. 208: 78–85.
Dombrovsky, L.A., 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., Kokhanovsky A.A., Randrianalisoa J.H. (2019) On Snowpack Heating by Solar Radiation: A Computational Model, J. Quant. Spectrosc. Radiat. Transf., 227: 72–85.
Hossain, M.M. and Gu, M. (2016) Radiative Cooling: Principles, Progress, and Potentials, Adv. Sci., 3(7): 1500360.
- Born, M. and Wolf, E. (1999) Principles of Optics, 7th (expanded) ed., New York: Cambridge University Press.
- Dombrovsky, L.A. (2012) The Use of Transport Approximation and Diffusion-Based Models in Radiative Transfer Calculations, Comput. Therm. Sci., 4(4): 297–315.
- Dombrovsky, L.A. (2019) Scattering of Radiation and Simple Approaches to Radiative Transfer in Thermal Engineering and Bio-Medical Applications, In A. Kokhanovsky (Ed.), Springer Series in Light Scattering, vol. 4, Cham, Switzerland: Springer, pp. 71–127.
- Dombrovsky, L.A and Baillis, D. (2010) Thermal Radiation in Disperse Systems: An Engineering Approach, Danbury, CT: Begell House.
- Dombrovsky, L.A. and Randrianalisoa, J.H. (2018) Directional Reflectance of Optically Dense Planetary Atmosphere Illuminated by Solar Light: An Approximate Solution and Its Verification, J. Quant. Spectrosc. Radiat. Transf., vol. 208: 78–85.
- Dombrovsky, L.A., 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., Kokhanovsky A.A., Randrianalisoa J.H. (2019) On Snowpack Heating by Solar Radiation: A Computational Model, J. Quant. Spectrosc. Radiat. Transf., 227: 72–85.
- Hossain, M.M. and Gu, M. (2016) Radiative Cooling: Principles, Progress, and Potentials, Adv. Sci., 3(7): 1500360.