SNOWTRAN: A Fast UV-VNIR-SWIR Radiative Transfer Model for Polar Hyperspectral Remote Sensing Applications

1. Introduction

Radiative transfer modelling is an important tool for simulating and studying the interaction of radiation with different underlying surfaces and the atmosphere in order to make better use of remote sensing data. Multiple radiative transfer codes exist for the simulation of satellite signals at the top-of-atmosphere (TOA) (see, e.g., [1,2,3,4,5,6,7,8]). Papers [2,3,4,5,6,7,8]) are aimed at the solution of the general radiative transfer in various turbid media. In contrast, [1] has a special target - the spectral TOA reflectance over snow. The authors of [1] rely on the numerical solution of the radiative transfer equation (RTE), which slows down the calculation procedure considerably.

A new generation of spaceborne imaging spectroscopy sensors have been appearing in the recent years (e.g., EnMAP, PRISMA[9,10]) that have observational capacities to observe on demand polar areas with high spatial and spectral resolution covering the broad spectral range 400-2500nm. Furthermore, the next generation is in preparation for the operational, regularly repeating satellite measurements including polar areas such as CHIME[11] and SBG[12], which will significantly increase the amount of observation data. In this frame, the task of this paper is to develop a fast and accurate radiative transfer code for the simulation of hyperspectral satellite measurements over Antarctica and other polar regions with underlying snow surfaces based on the approximate analytical solutions of the RTE, that will be applicable to current and future nadir looking imaging spectroscopy missions. The Sobolev approximation [13] is used to simulate the radiative transfer in the atmosphere. This approximation is valid at small values of the average cosine of scattering angle g, which is a valid approximation, e.g., over Antarctica, where the aerosol load is extremely low and, therefore, the light scattering is mostly due to molecules. It is known that g=0 for molecular scattering. The spherical albedo of semi–infinite homogeneous snow surfaces is considered in the framework of the van de Hulst approximation valid at any level of light absorption in snow [14,15] We have extended the van Hulst approximation to the case of the nadir reflectance as observed by hyperspectral instruments orbiting our planet at present (see, e.g., [9,10]. Furthermore, we have used the analytical representation of both aerosol phase function and aerosol spectral extinction coefficient, which makes it possible to avoid time-consuming Mie calculations or the usage of respective look-up-tables. We have found that the shortwave infrared (SWIR) TOA reflectance spectra closely coincide with the spectra of the underlying surface modulated by the atmospheric gaseous absorption processes in this case. This is due to the weak aerosol load in polar regions.

The first version of a software suite for studies of atmosphere – underlying SNOW- spaceborne optical receiver light TRANsmission calculations (SNOWTRAN) has been proposed in [16] $\mathrm{f}\mathrm{o}\mathrm{r} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{l} \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e} 400-1000 \mathrm{n}\mathrm{m}.$ In this paper, we extend the range of the applicability of the SNOWTRAN to the spectral region 320-2500 nm using the extended databases of absorption cross sections of atmospheric molecules such as $O_2$ $,O_3\mathrm{a}\mathrm{n}\mathrm{d}{H}_{2}O$. The local (single scattering albedo, asymmetry parameter) and global (reflectance, albedo) snow optical properties are also calculated at any level of light absorption, making it possible to extend the range of the applicability of the respective equations and software to the shortwave infrared (SWIR) region. The main advantage of SNOWTRAN is the possibility to perform fast and accurate radiative transfer calculations for atmosphere – underlying snow system both for clean and polluted snow cases.

2. The Analytical Representation of the Top-of-Atmosphere Reflectance over Snow

2.1. The Radiative Transfer in a Snow Layer

The top-of-atmosphere (TOA) reflectance over snow $R_{TOA}$as used in version 1.1 of the SNOWTRAN software, is given by the following analytical expression [16]:

$$R_{TOA}=R{T}_g,$$

(1)

where $T_g$ is the gaseous transmittance, R is the TOA reflectance of the gaseous absorption free atmosphere - underlying snow surface derived as

$$R{=R}_a+\gamma T_a{R}_s,$$

(2)

$R_a$ is atmospheric reflectance for the case of a black underlying surface, $T_a$is the two-way atmospheric transmittance, ${\mathrm{a}\mathrm{n}\mathrm{d}R}_s$ is the snow reflectance (bottom of atmosphere (BOA) reflectance). The parameter $\gamma$ accounts for the atmosphere–snow retro-reflections. The value of $\gamma$ can be calculated using the following approximation strictly valid for the Lambertian underlying surfaces [16]:

$$\gamma ={(1-r_a{r}_s)}^{-1},$$

(3)

where $r_a$ is the atmospheric albedo, $r_s$ is the snow spherical albedo.

It is assumed that the snow can be represented as a semi-infinite layer with irregularly shaped ice crystals. The spherical albedo of such a layer with the single scattering albedo ${\omega }_0$ and the asymmetry parameter $g$ can be estimated using the following approximation [14]:

$$r_s=\frac{(1-as)(1-s)}{(1+bs)},$$

(4)

where a=0.139, b=1.17 and

$$s=\sqrt{\frac{(1-{\omega }_0)}{(1-g{\omega }_0)}}$$

(5)

is the similarity parameter. It follows from Eq. (4) that various turbid media with close values of the similarity parameter s have similar values of spherical albedo. The approximation given by Eq. (4) is very accurate and can be used at any level of light absorption in snow, which is relevant for the modelling of visible-near infrared (VNIR) and SWIR snow spectra.

In this paper we derive the relevant approximation for the value of the nadir snow reflectance $R_s,\mathrm{w}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{V}\mathrm{N}\mathrm{I}\mathrm{R}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}\mathrm{r}\mathrm{e}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}.$In principle, $R_s$ can be found solving the respective integro-differential radiative transfer equation [13,15]. It is assumed in the framework of SNOWTRAN that the nadir snow reflectance can be calculated using the following approximation derived from the parameterization of the simultaneous radiative transfer calculations of the snow reflectance $R_s$ and spherical albedo $r_s$as presented in [15]:

$$R_s=a_0+a_1{r}_s+a_2{r}_s^2,$$

(6)

where the parameters $a_0,a_1,a_2$ depend on the cosine of the solar zenith angle${\mu }_0$ (see Appendix A). The accuracy of Eqs. (4), (6) is demonstrated in Figure 1 for the case of the Henyey - Greenstein phase function [13,15] with the asymmetry parameters g equal to 0.75 and 0.875 (typical for snow [17]) and several values of the single scattering albedo ${\omega }_0$.

Eqs. (6) and (4) make it possible to calculate the global radiative transfer characteristics of snow layers for a given value of the similarity parameter s. It follows from Eq. (5) that the similarity parameter s depends on just two local optical characteristics of snow (${\omega }_0$ and $g$). Therefore, an important point is to find analytical relationships of these parameters with the size of ice grains. This also makes it possible to reduce the number of input parameters for the model proposed in this paper. Under assumption of clean snow, we shall use the following approximate equations for the relevant parameters [17,18]:

$${\omega }_0=1-{\beta }_i,{\beta }_i=\frac{1}{2}(1-\rho )(1-e^{-\sigma z}),g=g_{\infty }-(g_{\infty }-g_0)e^{-\epsilon z}$$

(7)

derived using the geometrical optics approximation for the fractal ice grains [17,19]. Here, ${\beta }_i$ is the probability of photon absorption (PPA) by pure snow at the wavelength $\lambda$, $\alpha =4\pi \chi /\lambda ,$ $z=\alpha d_{ef}$ is the bulk ice absorption coefficient, $\chi$ is the imaginary part of the ice refractive index, $d_{ef}$ is the effective diameter of ice grains and parameters $\rho ,\sigma ,\epsilon ,g_0,g_{\infty }$depend on the assumed type/shape of ice crystals in snow. The effective grain diameter is defined as[17,20,21]:

$$d_{ef}=\frac{3V}{2\Sigma },$$

(8)

where V is the average volume of the ice grains and $\Sigma$ is their average cross section perpendicular to the incident beam. It is assumed that particles are randomly oriented.

We have found that the volume of fractal grains used in simulations (Koch fractals of second generation [19]) is related to the side length $\mathcal{l}$ by the following formula:

$$V=b{\mathcal{l}}^3,$$

(9)

where $b=27/52\sqrt{2}\approx \mathrm{0.367.}$ Then it follows: $d_{ef}=\varsigma {\mathcal{l}}^3/\Sigma$, where $\varsigma =1.5b$. The numerical calculations show that the value of $d_{ef}$ is close to $\mathcal{l}$/2.

We use the following equations for the parameters$\rho ,g_0,g_{\infty }$ derived for the real part of the refractive index n in the range 1.25-1.35, which is the case for ice in the spectral region 0.3-2.5 μm [17,22]:

$$\mathsf{\rho }=0.0123+0.1622(n-1),g_0=0.9919-0.769(n-1),g_{\infty }=1.008-0.11(n-1).$$

(10)

In addition, when fitting the numerical Monte–Carlo geometrical optical calculations with the use of Eqs. (7), we found that σ=0.9045, ε=0.8571. These parameters depend on the assumed shapes of ice crystals in snow.

The extinction cross section of large nonspherical randomly oriented particles can be calculated as follows [20,21]:

$$C_{ext}=2\Sigma .$$

(11)

Therefore, it follows for the extinction coefficient [17]:

$$k_{ext}=\frac{3c}{d_{ef}},$$

(12)

where c is the volumetric concentration of particles. Taking into account that c is close to 1/3 for natural snow$\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}$ [17], it is concluded that the extinction length in snow $L_{ext}=k_{ext}^{-1}$ is close to the effective diameter of ice grains.

The intercomparison of results derived using simple Eqs. (7) and time consuming geometrical optical Monte - Carlo calculations [18] is shown in Figure 2 a,b for the case of large nonspherical particles (Koch fractals of the second generation [18] with the side length of the initial tetrahedron $\mathcal{l}$ equal to 50, 100 and 300 μm ( or $d_{ef}$=0.027, 0.054, and 0.162 mm, respectively). The ice spectral refractive index proposed in [22] has been used. The high accuracy of Eqs. (7) is clearly demonstrated. Some deviations could be actually due to the statistical error of Monte – Carlo calculations and not due to the errors of Eqs. (7). Figure 2c shows the results of calculation of the similarity parameter s using Eqs. (5), (7) and Monte – Carlo geometrical optics calculations. Both results almost coincide, which makes it possible to derive analytical relationship between the value of the nadir reflectance (see Eqs. (6), (4)) and the effective diameter of ice grains, which is of importance for the solution of the inverse snow optics problems.

Figure 2. a. The spectral dependence of the probability of photon absorption derived using the geometrical optics Monte Carlo code (lines) and simple Eq. (7) (symbols) for the length $\mathcal{l}$ of the side of the initial ice tetrahedron for the Koch fractal of the second generation [17] equal to 50, 100 and 300µm.

Figure 2. a. The spectral dependence of the probability of photon absorption derived using the geometrical optics Monte Carlo code (lines) and simple Eq. (7) (symbols) for the length $\mathcal{l}$ of the side of the initial ice tetrahedron for the Koch fractal of the second generation [17] equal to 50, 100 and 300µm.

Figure 2. b. The same as in Figure 2a except the asymmetry parameter is considered. The lower and upper lines give the spectral dependencies of the parameters $g_0\left(\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\right)$and $g_{\infty }$ (upper line).

Figure 2. b. The same as in Figure 2a except the asymmetry parameter is considered. The lower and upper lines give the spectral dependencies of the parameters $g_0\left(\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\right)$and $g_{\infty }$ (upper line).

Figure 2. c. The same as in Figure 2a except for the similarity parameter.

Figure 2. c. The same as in Figure 2a except for the similarity parameter.

We have studied the case of type 1 snow (clean snow) in the discussion given above. Let us consider now the case of type 2 snow (snow with impurities). The level of snow pollution is usually low and, therefore, in a good approximation, we can assume that the light scattering and extinction processes in snow are dominated by ice grains. The absorption coefficient of snow with impurities can be presented as:

$$k_{abs}=c_i{\kappa }_i+c_p{\kappa }_p,$$

(13)

where $c_i$ is the volumetric concentration of ice grains and $c_p$is the volumetric concentration of impurities. We assume that just one type of impurity is present in the snow, although the technique can easily take into account multiple types of impurities by adding additional terms in Eq. (13). The values of ${\kappa }_i\mathrm{and}{\kappa }_p$ are volumetric absorption coefficients of ice grains and impurities, respectively. They can be calculated as follows:

$${\kappa }_i=\frac{〈C_{abs,i}〉}{〈V_i〉},{\kappa }_p=\frac{〈C_{abs,p}〉}{〈V_p〉},$$

(14)

where $〈C_{abs,i}〉$ is the average absorption cross section of ice grains, $〈V_i〉$is the average volume of ice grains, $〈C_{abs,p}〉$ is the average absorption cross section of impurities, $〈V_p〉$is the average volume of impurities. It follows for the probability of photon absorption:

$$\beta \equiv \frac{k_{abs}}{k_{ext}}={\beta }_i+c_p\frac{{\kappa }_p}{k_{ext}},$$

(15)

where ${\beta }_i$is the PPA for snow without impurities (see Eq. (7)), and thus (see Eqs. (15), (12)):

$$\beta ={\beta }_i+\frac{c{d}_{ef}}{3}{\kappa }_p,$$

(16)

where $c$=$c_p/c_i$. One can see that the probability of photon absorption for the snow with impurities depend not only on the value of ${\beta }_i$ but also on the relative volumetric concentration of impurities c and their spectral volumetric absorption coefficient ${\kappa }_p$. We shall parameterize the spectrum ${\kappa }_p$ in the following way:

$${\kappa }_p(\lambda )={\kappa }_p({\lambda }_0){(\lambda /{\lambda }_0)}^{-m},$$

(17)

where ${\kappa }_0={\kappa }_p({\lambda }_0)$ and m are the input parameters of the discussed analytical model. We shall assume that ${\lambda }_0$= 550 nm.

2.2. The Atmospheric Radiative Transfer

Let us consider the atmospheric radiative transfer processes now. It is assumed that the polar atmosphere is clean and the main contribution to light scattering and extinction in atmosphere is due to molecular scattering. Therefore, we use the same approach for the calculation of the atmospheric characteristics (${{r_a,R}_a,T}_a$) as described in [16] The input for the atmospheric model is the spectral aerosol optical thickness ${\tau }_{aer}$ presented as ${\tau }_{aer}={(\lambda /{\lambda }_0)}^{-B}{\tau }_{aer}({\lambda }_0)$, where ${\tau }_{aer}({\lambda }_0)$is the aerosol optical thickness at the wavelength ${\lambda }_0$= 550 nm[16] and B is the Angström parameter. The molecular optical thickness is calculated as suggested in [16]. The gaseous transmittance $T_g$ in Eq. (1) is found using the following approximation[16]:

$$T_g=T_{H_{2}O}{T}_{O_2}{T}_{O_3},$$

(18)

where gaseous transmittances for water vapor, oxygen and ozone are given by the following expressions:

$$T_{H_{2}O}=\exp (-{(M{f}_{H_{2}O}{\tau }_{H_{2}O})}^{\varsigma }),T_{O_2}=\exp (-{(M{f}_{O_2}{\tau }_{O_2})}^{\upsilon }),T_{O_3}=\exp (-M{\tau }_{O_3}).$$

(19)

Here M is the airmass factor, $f_{H_{2}O}={(P/P_0)}^{m_{H_{2}O}}{(T_0/T)}^{p_{H_{2}O}}$, $f_{O_2}={(P/P_0)}^{m_{O_2}}{(T_0/T)}^{p_{O_2}}$. The pair $(T_0,P_0)$represents normal temperature (273.16K) and pressure (1013.25hPa) and the values $(T,P)$ give the average values of temperature and pressure over vertical for a given location. The parameters ($\varsigma ,\upsilon ,m_{H_{2}O},p_{H_{2}O},m_{O_2},p_{O_2}$ ) depend on the absorption band. The optical thicknesses for various gaseous absorbers are calculated as follows:

$${\tau }_{H_{2}O}=k_{abs,H_{2}O}{\mathcal{l}}_{H_{2}O},{\tau }_{O_2}=k_{abs,O_2}{\mathcal{l}}_{O_2},{\tau }_{O_3}=C_{abs,O_3}{N}_{O_3}$$

(20)

where $k_{abs,H_{2}O}$ is the absorption coefficient of water vapor$(c{m}^{-1})$, ${\mathcal{l}}_{H_{2}O}$is the precipitable water vapor (PWV) measured in cm, $k_{abs,O_2}$is the molecular oxygen absorption coefficient $(c{m}^{-1})$,${\mathcal{l}}_{O_2}$ is the total molecular oxygen column (TMOC) (cm),$N_{O_3}$ is the total ozone column (in Dobson Units (DU)), $C_{abs,O_3}$ is the absorption cross section of ozone molecule in the inverse DU (the multiplier 2.69$\times$ ${10}^{16}$ should be applied to transfer $C_{abs,O_3}$ (usually provided in the units ${cm}^2/molec)$ to the Dobson Units).

The vertical profile of molecular oxygen concentration could be considered as a constant one as compared to highly variable profiles of ozone and water vapor concentrations. Let us introduce the effective oxygen vertical column: $L_{O_2}$=$f_{O_2}$${\mathcal{l}}_{O_2}$. Then it follows: $T_{O_2}=\exp (-{(M{L}_{O_2}{k}_{abs,O_2})}^{\upsilon })$. Therefore, for the calculation of the oxygen transmittance we use the value of $L_{O_2}$, which reduces the number of parameters needed for the calculations. More precisely, the Normalized effective OXygen amount (NOX) $K_{O_2}$is used as input in SNOWTRAN. It is defined as follows: $K_{O_2}$=$L_{O_2}$/${\mathcal{l}}_{O_2}^{st}$, where ${\mathcal{l}}_{O_2}^{st}=87068.53cm$ is the molecular oxygen vertical column fo the standard atmosphere$\left[8\right]$. The value of $K_{O_2}$ is equal to the ratio ${\mathcal{l}}_{O_2}$/${\mathcal{l}}_{O_2}^{st}$at $f_{O_2}=1.$

We use the spectra $k_{abs,H_{2}O}(\lambda )$, $k_{abs,O_2}$, $C_{abs,O_3}$and also parameters $\varsigma ,\upsilon ,m_{H_{2}O},p_{H_{2}O}$ provided in $\left[8\right]$ at the spectral resolution of 7 nm (the triangular slit function has been used) in steps of 1 nm. In addition, we have calculated the oxygen absorption coefficient using the approximation provided in [23] for the oxygen absorption band located at the wavelength 1.27 μm, which is not provided in $\left[8\right]$.

The geometrical approximation for the airmass factor has been used:

$$M={\xi }^{-1}+{\eta }^{-1}.$$

(21)

Here (ξ, η) are the cosines of the solar and observation zenith angles, respectively.

The example of calculations using SNOWTRAN for the case of polluted snow is given in Figure 3. The values of parameters used in the calculations are given in Table 1.

3. The Comparison of Theoretical Calculations with Satellite Measurements

Let us compare SNOWTRAN calculations with spaceborne hyperspectral measurements over Antarctica. We will use data of the German spaceborne imaging spectrometer mission EnMAP[10,24] and also the Italian Space Agency imaging spectrometer mission PRISMA [29].

The EnMAP performs the spectral radiance measurements with the spatial resolution 30m in the spectral range 418-2450 nm with the average spectral sampling distance 6.5nm in VNIR and 10 nm in SWIR (https://www.enmap.org/). The measurements are performed using two separate spectrometers (VNIR and SWIR). PRISMA (https://www.eoportal.org/satellite-missions/prisma-hyperspectral) has the same spatial resolution and the spectral range 407 – 2469 nm is covered with a spectral sampling interval ≤ 12 nm in both the VNIR and SWIR (https://www.prisma.asi/). The conversion of the spectral radiance I provided in L1C EnMAP data and also PRISMA top-of-atmosphere radiance data to the spectral reflectance R has been performed using the following equation:

$$R=\frac{\pi I}{\xi d_{sol}{F}_{sol}},$$

(22)

where the parameter $d_{sol}$ accounts for the Sun-Earth distance and $F_{sol}$ is the solar irradiance at the top-of-atmosphere, which was also considered for the radiometric calibration of the sensor. The latter is important because TOA reflectance should be independent of the solar model, which is not the case if two different models were used. Therefore, we have applied the solar irradiance model given in [25] in case of EnMAP and the model provided in [26] in case of PRISMA.

Figures 4 and 5 show SNOWTRAN results compared with respectively EnMAP and PRISMA at different spectral resolutions, and Table 2 resumes the simulation parameters for both cases. The spectral reflectance derived from the EnMAP measurements for a single pixel located at (75.119S, 123.902E) is given in Figure 4a together with the calculations of the top-of atmosphere and bottom-of-atmosphere (BOA) snow reflectance. The EnMAP measurements have been performed on February 15, 2023, 00:36 UTC. We have used the values of surface pressure 651hPa, surface height of 3233m and average values of pressure and temperature along vertical equal to 491 hPA and 229 K as in our previous publication[27]. The values of TOC (250DU) and PWV (0.175cm) were taken coincident with those measured at the closely located Concordia research station ( Dome C, 75S, 123.3E). Other parameters used in SNOWTRAN simulations are given in Table 2. They have been selected by fitting SNOWTRAN and satellite measurements at several spectral points. It follows from Figure 4a that SNOWTRAN is indeed capable to describe efficiently the EnMAP hyperspectral satellite measurements.

Figure 4. a. The intercomparison of SNOWTRAN spectra with the EnMAP measured spectra (75.119S,123.902E, February 15, 2023, 00:36UTC). The EnMAP data in the spectral range with strong water vapor absorption (e.g., around 1800nm) are not provided by the satellite ground system.

Figure 4. a. The intercomparison of SNOWTRAN spectra with the EnMAP measured spectra (75.119S,123.902E, February 15, 2023, 00:36UTC). The EnMAP data in the spectral range with strong water vapor absorption (e.g., around 1800nm) are not provided by the satellite ground system.

Figure 4. b. The relative difference (in percent) between theoretical calculations and experimental data in the VNIR region.

Figure 4. b. The relative difference (in percent) between theoretical calculations and experimental data in the VNIR region.

Figure 4. c. The absolute difference between theoretical calculations and experimental data in the VNIR (upper panel) and SWIR (lower panel).

Figure 4. c. The absolute difference between theoretical calculations and experimental data in the VNIR (upper panel) and SWIR (lower panel).

The reason for the deviations of the theoretical calculations and EnMAP measurements in the range 418-500 nm for bright snow surfaces must be further investigated. It is unlikely that they are related to the inaccuracy of the SNOWTRAN model. Some deviations of theoretical calculations in gaseous absorption bands are due to different spectral sampling of SNOWTRAN (1 nm) and EnMAP (6.4- 10 nm) measurements.

We have re-sampled the SNOWTRAN data derived with 1nm spectral resolution to match the spectral sampling and spectral instrument response function of EnMAP. The difference between the re-sampled SNOWTRAN data and EnMAP spectra in the spectral range 418-920nm and 920-2500nm is given in Figure 4b,c. One can see that the difference is smaller than 2% for most of points not contaminated by gaseous absorption, which is an excellent result both with respect to the SNOWTRAN and EnMAP performance. The differences are somewhat larger in SWIR (see Figure 4c), which is explained by the influence gaseous absorption bands, especially water vapor. The difference in gaseous absorption bands needs more thorough analysis with possible account for other gaseous absorbers and change in the assumed abundance of atmospheric gaseous components in the SNOWTRAN run. The errors related to the integration of the oscillating spectral radiance in the gaseous absorption bands with account for the EnMAP spectral response function could be also of importance.

The fit of PRISMA measurements at the location (76.14S, 129.93E) situated at somewhat larger distance from the ground validation site located at Dome C is given in Figure 5a. The measurements have been taken on 21.12.2022, 00:14:07UTC. We also see that the SNOWTRAN can be successfully used to simulate the TOA PRISMA spectral reflectance. The ground observations at DOME C give TOC=269 DU and PWV=0.045 cm at the same date and time. These parameters are somewhat smaller than the parameters used in the simulations (289 DU and 0.055 cm, respectively). The parameters used in simulations have been selected to provide the best fit of experimental and theoretical curves in the ozone and water vapor absorption bands. Therefore, they can be considered as retrieved characteristics. The possible explanation for the differences between ground and satellite - derived parameters could be due to not exact collocation of ground and satellite measurements and also due to corresponding experimental errors of related retrieval techniques. The spectral relative and absolute difference of PRISMA data and re-sampled SNOWTRAN calculations are given in Figure 5b,c. The resampling has been performed to match the spectral sampling and spectra instrument response function for PRISMA. One can see that the relative difference is not larger than 5% in the deep blue region in contrast to the difference with EnMAP data shown in Figure 4b.

Figure 5. a. The intercomparison of SNOWTRAN spectra with the PRISMA measured spectra (76.14S,129.93E, December 21, 2022, 00:14UTC).

Figure 5. a. The intercomparison of SNOWTRAN spectra with the PRISMA measured spectra (76.14S,129.93E, December 21, 2022, 00:14UTC).

Figure 5. b. The spectral relative difference (in percent) between theoretical calculations and PRISMA data in the VNIR region.

Figure 5. b. The spectral relative difference (in percent) between theoretical calculations and PRISMA data in the VNIR region.

Figure 5. c. The spectral absolute difference between theoretical calculations and PRISMA data.

Figure 5. c. The spectral absolute difference between theoretical calculations and PRISMA data.

It follows from Figure 4a that EnMAP measurements can be well fitted using SNOWTRAN both at 1030 nm and 2200 nm bands almost free of gaseous absorption and atmospheric scattering effects. This is not the case for the PRISMA measurements (see Figure 5a). This points to the fact that at the time and location of PRISMA measurements, the snow vertical inhomogeneity cannot be ignored. This brings us to the conclusion that snow vertical inhomogeneity could be potentially assessed using hyperspectral top-of-atmosphere reflectance measurements [27,28,29]. Also we can see that SNOWTRAN calculations and PRISMA measurements at 1500-2500 nm (and also for EnMAP, see Figure 4a) matches very well. This points out that this spectral range can be used to derive the snow grain size in upper snow layer due to the small penetration depth of the SWIR radiation in snow. This subject is considered in the next section.

One concludes from Figures 4a and 5a that the TOA reflectance is smaller than BOA reflectance in the gaseous absorption bands. At some channels the BOA reflectances are a little bit smaller than TOA reflectances. This is related to the uncertainties of the SNOWTRAN model. However, it should be underlined that the respective errors are smaller as compared to the calibration errors of the hyperspectral measurements.

Currently investigations are on-going for the differences in the spectral range 418-464nm (see Figure 4a) between SNOWTRAN and EnMAP data. Such differences and small – scale oscillations as seen in Figure 4b have been found for water surfaces as well. They might be related to the calibration of EnMAP data, which is indirectly supported by data shown in Figure 5a for PRISMA, where respective differences in the deep blue region of the electromagnetic spectrum are much smaller.

Table 2. The parameters of the model used in simulations for the case of nadir observations shown in in Figs.4a, 5a. The ground-observed TOC and PWV at Dome C (75S, 123.3E) are given in the brackets. It has been assumed that snow is free of impurities. The small difference of satellite measurements ( up to 5 degrees) from the exact nadir direction is ignored.

Parameter	EnMAP 15.02.2023	PRISMA 21.12.2022
$d_{ef}$, mm	0.11	0.2
TOC, DU	250	289 (269)
PWV, cm	0.175	0.055 (0.045)
TOX	0.45	0.13
${\tau }_{aer}(550nm)$	0.14	0.05
B	1.0	1.35
${\vartheta }_o,degrees$	68.0	57.6
Latitude, degrees	75.119S	76.140S
Longitude, degrees	123.902E	129.930E
Date/time (UTC)	15.02.2023/ 00:36UTC	21.12.2022/ 00:14UTC

Table 2. The parameters of the model used in simulations for the case of nadir observations shown in in Figs.4a, 5a. The ground-observed TOC and PWV at Dome C (75S, 123.3E) are given in the brackets. It has been assumed that snow is free of impurities. The small difference of satellite measurements ( up to 5 degrees) from the exact nadir direction is ignored.

Parameter	EnMAP 15.02.2023	PRISMA 21.12.2022
$d_{ef}$, mm	0.11	0.2
TOC, DU	250	289 (269)
PWV, cm	0.175	0.055 (0.045)
TOX	0.45	0.13
${\tau }_{aer}(550nm)$	0.14	0.05
B	1.0	1.35
${\vartheta }_o,degrees$	68.0	57.6
Latitude, degrees	75.119S	76.140S
Longitude, degrees	123.902E	129.930E
Date/time (UTC)	15.02.2023/ 00:36UTC	21.12.2022/ 00:14UTC

4. The Determination of the Snow Grain Size

In this section we consider the application of the proposed theoretical framework (forward model) for the determination of the snow grain size from EnMAP measurements. It is an important fact that both albedo and nadir reflectance primarily depend just on one parameter characterizing the snow microstructure (the similarity parameter s). Therefore, the inverse problem is reduced to the determination of snow microstructure from the value of the similarity parameter s. This parameter can be directly retrieved from the measurements of $r_s\mathrm{or}{R}_s$, i.e., it follows from Eq. (4):

$$s=\frac{\psi -\sqrt{{\psi }^2-4a(1-r_s)}}{2a},$$

(23)

where $\psi =1+a+b{r}_s$. One can see that s=0 at $r_s$=1 and s=1 at $r_s$=0 as it should be (see Figure 1 and Eq. (4)). The value of $r_s$can be measured directly or derived from the nadir reflectance (see Eq. (6)):

$$r_s=\frac{\sqrt{a_1^2-4a_2(a_0-R_s)}-a_1}{2a_2}.$$

(24)

Eqs. (5), (7) can be used to determine the value of the effective grain size from the value of the similarity parameter s analytically. Namely, it follows:

$$d_{ef}=\frac{1}{\alpha \sigma }\ln \left[\frac{\varphi /s^2+{\gamma }_1}{\varphi /s^2-{\gamma }_2}\right],$$

(25)

where $\alpha$ is the bulk ice absorption coefficient at the wavelength λ, $\varphi =\frac{1}{2}(1-\rho ),{\gamma }_1=g_{\infty }-g_0,{\gamma }_2=1-g_{\infty }$. The derivation of Eq. (25) has been performed at the assumption $\epsilon =\sigma$ (see Eq. (7)), which is a valid approximation for the case at the study.

One can see that the value of $d_{ef}$ can be derived analytically from the values of $r_s$or $R_s$ (see Eqs. (23) - (25)). It is of importance that the values of s, $r_s$and $R_s$depend on the wavelength and $d_{ef}$ is not (for vertically homogeneous snow layers). Therefore, the spectral dependence of $d_{ef}$ is an indication of the snow vertical inhomogeneity. We select the measured reflectance at the wavelengths of 1030 nm, 1235 nm and 2200 nm, which are least affected by atmospheric scattering and absorption processes (see Figure 4a) and also sensitive to the snow grain size. Then it follows: $R_{TOA}\approx R_s$ ( see Eqs. (1), (2) and Figure 3).

Due to the different penetration of radiation in vertically inhomogeneous snow at different wavelengths$\left[17\right]$, the values of $d_{ef}$ derived at different wavelengths may differ. Let us introduce the ratio of effective diametersretrieved at 1235nm ($d_{ef,1235nm}$) and 2200nm ($d_{ef,2200nm}$) to that at 1030nm ($d_{ef,1030nm}$):

$$K_1=\frac{d_{ef,2200nm}}{d_{ef,1030nm}},K_2=\frac{d_{ef,1235}}{d_{ef,1030nm}}.$$

(26)

It should be noted that it can generally be assumed that $K_1<K_2$because the SWIR radiation with longer wavelengths generally penetrates to the smaller distances from the snow top, where new fresh fallen snow with smaller crystals are located. The radiation with the wavelength 1235 nm can penetrate into somewhat deeper snow layers as compared to that at the wavelength 2200nm.

We have derived the values of $d_{ef,1030nm}$, $d_{ef,1235nm}$and also $d_{ef,2200nm}$ for the EnMAP 30 km x 30 km scene located at the Aviator Glacier in Antarctica (16.02.2023, 22:11UTC). A glacier is a mass of ice with sufficient thickness to flow away from the source area in lobes, tongues or masses. The results are given in Figure 6a-c.

The valid results are given in the central part of the picture (Aviator Glacier Tonque) between two mountainous areas, where the ice flow to the ocean occurs. One can see that the value of $d_{ef,1235nm}\approx d_{ef,1030nm}$ and $d_{ef,2200nm}<d_{ef,1030nm}$ as one may expect. In particular, we have found that $d_{ef,1030nm}\approx$ 0.52mm, $d_{ef,1235nm}\approx$0.58mm, $d_{ef,2200nm}\approx$0.21mm, $K_1=0.4$ and $K_2=1.1$ for the area A (73.75S-73.85S, 165-165.1E). Hence, the sizes of ice crystals at the top ($d_{ef,2200nm}$) are more than two times smaller compared to the sizes inside the snow layer ($d_{ef,1030nm}$) as one might expect due to the snow metamorphism processes. The spatial distribution of the parameter $K_1$ is given in Figure 7. It should be pointed out that although the parameters $d_{ef}$are affected by the assumed shape of light scattering particles$\left[17\right]$, this is less the case for the ratios $K_1\mathrm{and}{K}_2$ due to the cancellation of a common multiplier.

The retrievals for the case of vertically inhomogeneous snow requires the change of the forward model, which is out of scope of this paper. Because the value of $K_1$ is close to one, one may expect that the snow observed using the couple of wavelengths 1030 and 1235 nm can be considered as a homogeneous one. Therefore, the algorithm proposed in [27] is valid for the case considered here. Indeed, we have derived the value of $d_{ef,0}=d_{ef,1}=$0.53mm from the algorithm described in [27] for the area A (see above), which is close to the value of $d_{ef,1030nm}$ derived in the framework of the algorithm discussed in this paper. It should be pointed out that the presented algorithm and the one described in [27] are based on different system of equations. Therefore, our conclusions confirm the consistency of both algorithms even though the current algorithm is able to determine the grain sizes at the very top of the snow surface, which is not the case for the algorithm described in [27]. Although we do not have grain sizes measured in situ and collocated with satellite measurements, our results are consistent with those presented in [31], where the most frequent grain size diameters in the range 0.2-0.4 mm at another coastal Antarctic site (67.7S, 139.72E, Adelie Land). More validation work is needed to relate the satellite - measured snow grain size with in situ derived values of the same parameter.

The snow grain size derived from EnMAP measurements at another location (closer to DOME C) using band 1030 nm is shown in Figure 8. It is obvious that the spatial distribution of the retrieved snow grain size is not so pronounced as in the case of the Aviator Glacier due to the reduced number of factors influencing the snow metamorphism processes (e.g., precipitation, surface inclination, proximity to the ocean and mountains). The average grain diameter derived from 1030 nm EnMAP measurements is 0.13 mm, which is close to the range of diameters (0.15-0.3 mm) measured at the site close to Dome C [31]. The spatial distribution of the ice grain diameter derived from the EnMAP measurements at 2200 nm is somewhat smaller (0.1 mm) with the spatial distribution similar to that shown in Figure 8. Basically, we can conclude that the vertical inhomogeneity of the snow is less pronounced as compared to the Aviator Glacier case shown in Figure 6, where the sizes at the top and inside the snow differ by more than two times, being larger at larger depths inside the snow. The small vertical variability of snow grain size around Dome C has been found using ground measurements as well [31].

It should be pointed out that the application of the channel 2200nm to the snow grain size retrieval must be performed only for the case of small ice crystals. With the increase of the crystal size, the reflectance at this channel saturates and the information on the snow grain size is lost (see Figure 9). The same applies to the channel 1640nm.

The EnMAP data can be used not only for the determination of snow parameters but also for the characterization of atmospheric state over the snow field as discussed in [27]. In particular, we have found a significant difference in the value of the Precipitable Water Vapor at Dome C on February 15, 2023 as compared to December 21, 2022 (see Table 2). This peculiarity is explained by the advective warming event at Dome C on February 15, 2023 reported in [32].

5. Conclusions

We have proposed a fast and simple fully analytical model for the calculation of the TOA reflectance over snow surfaces. The derived equations have been also used to derive the snow grain size, which can be used to estimate other snow parameters such as snow spectral and broadband albedo, bottom of atmosphere reflectance and snow specific surface area. It has been assumed that the snow is an ideally horizontal semi-infinite homogeneous surface without sastrugi and any inclination. The model can be easily extended considering Mie scattering by atmospheric particles and inclined surfaces. The case of nadir observations common for current and future hyperspectral missions is considered. The technique can be easily extended to the case off-nadir observations, if needed. The simplicity of the model makes it possible to apply the developed software for the solution of inverse problem, i.e. the determination of atmosphere and underlying snow properties from hyperspectral satellite observations as demonstrated here for EnMAP measurements. In particular, we have found a significant difference in the value of the Precipitable Water Vapor (PWV) at Dome C on February 15, 2023 as compared to December 21, 2022 (see Table 2).

Although the derived expression for the nadir snow reflectance as the function of the effective grain size is only approximation to the exact radiative transfer solution for the case of clean Antarctic atmosphere [33], the radiative transfer in natural snow is such a complicated process that exact solution can not be derived in principle taking into account irregular shape of particles and layered nature of snow layers with possible areas covered by frost and sastrugi. Therefore, we believe that our approach for modelling of radiative transfer in snow – underlying surface system will be useful for the various applied problems including the snow grain size and new fallen snow detection in polar regions using hyperspectral SWIR measurements. The technique can be also applied to other regions and surfaces (e.g., deserts). In this case the atmospheric correction procedure and respective expression for the similarity parameter as the function of the particle size must be modified taking into account spectral refractive index of light scattering particles and updates of several constants used in Eq. (7).

We have found that the snow grain sizes around Aviator Glacier are about 4 times larger as compared to those at Dome C. This leads to the increased light absorption by snow at this coastal site, which contributes to the increased rate of snow melting and formation of glacier rivers and lakes at favourable temperatures of surrounding air [34].

The SNOWTRAN can be coupled with the fast infrared radiative transfer codes (see, e.g., [35]) to provide users with a tool, which can be used to study various applied problems related to hyperspectral remote sensing observations in a broad spectral range from ultraviolet to the far infrared. The software is distributed upon request.